bioawk是由生信大牛李恒用C开发,不仅涵盖了awk的功能,在此基础上添加了很多额外的生物大数据文件的支持,如:BED、 GFF、 SAM、VCF、FASTA/Q等等。

安装:conda install bioawk

下面简单介绍其使用方法:

基础应用


# 查看帮助
bioawk [command] help
# 如:bioawk -c help

# 输出:
bed:
	1:chrom 2:start 3:end 4:name 5:score 6:strand 7:thickstart 8:thickend 9:rgb 10:blockcount 11:blocksizes 12:blockstarts 
sam:
	1:qname 2:flag 3:rname 4:pos 5:mapq 6:cigar 7:rnext 8:pnext 9:tlen 10:seq 11:qual 
vcf:
	1:chrom 2:pos 3:id 4:ref 5:alt 6:qual 7:filter 8:info 
gff:
	1:seqname 2:source 3:feature 4:start 5:end 6:score 7:filter 8:strand 9:group 10:attribute 
fastx:
	1:name 2:seq 3:qual 4:comment 

1)提取未mapping上的Reads(不含文件头):bioawk -c sam 'and($flag,4)' aln.sam.gz

2)提取mapping上的Reads:bioawk -Hc sam '!and($flag,4)'

3)得到反向互补序列:bioawk -c fastx '{print ">"$name;print revcomp($seq)}' seq.fa.gz

4)从SAM创建FASTA序列


samtools view aln.bam | \
     bioawk -c sam '{s=$seq; if(and($flag, 16)) {s=revcomp($seq)} print ">"$qname"\n"s}' > aln.fasta

5)输出样本foo以及bar的基因型:grep -v ^## in.vcf | bioawk -tc hdr '{print $foo,$bar}'

更多应用

1)统计Reads:bioawk -t -c fastx 'END {print NR}' test.fastq

2)统计GC含量:bioawk -c fastx '{print $name, gc($seq)}' test.fa

3)根据ID提取Fasta文件中相应的序列信息:


$ cat test.fa
>blah_C1
ACTGTCTGTC
ACTGTGTTGTG
ATGTTGTGTGTG
>blah_C2
ACTTTATATATT
>blah_C3
ACTTATATATATATA
>blah_C4
ACTTATATATATATA
>blah_C5
ACTTTATATATT

$ cat IDs.txt
blah_C4
blah_C5

命令:bioawk -cfastx 'BEGIN{while((getline k <"IDs.txt")>0)i[k]=1}{if(i[$name])print ">"$name"\n"$seq}' test.fa

结果如下:


>blah_C4
ACTTATATATATATA
>blah_C5
ACTTTATATATT

最后补充一句,awk的功能和语法biowak也是支持的。

参考资料:

1.https://github.com/lh3/bioawk

2.https://www.biostars.org/p/127141/