作為表格文本處理命令awk的延伸消请，Heng Li開發(fā)了專門處理生物信息學(xué)的文件（例如bed, sam, vcf, gtf等）的命令bioawk虐骑，我們了解一下叽奥。
簡介
bioawk發(fā)布在github上，可以通過如下命令進行下載饱狂，編譯與安裝：
git clone git://github.com/lh3/bioawk.git && cd bioawk && make && sudo cp bioawk /usr/local/bin/
bioawk與awk命令的使用方式類似（若不了解awk可以見上節(jié)內(nèi)容）曹步，命令最大的特點在于-c參數(shù)，通過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:strand 8:frame 9:attribute
fastx:
1:name 2:seq 3:qual 4:comment
通過-c輸入文件類型休讳，bioawk會智能地將各列賦給指定的變量（與awk類似讲婚，使用$1，$2也可以）俊柔。
實例
- 以
Mus_musculus.GRCm38.75_chr1.gtf為例筹麸，假如我們想要統(tǒng)計所有蛋白編碼基因的長度的話，可以采用如下命令：
$ bioawk -c gff '$source~/protein_coding/ && $feature~/gene/ {print $seqname"\t"($end - $start)}' Mus_musculus.GRCm38.75_chr1.gtf | head -n4
1 465597
1 16807
1 5485
1 12533
-
-c設(shè)置為fastx的話雏婶，bioawk會自動判斷輸入文件類型為fastq還是fasta物赶。例如我們統(tǒng)計fastq文件的行數(shù)：
$ bioawk -c fastx 'END{print NR}' contam.fastq
8
- 可以通過如下命令將
fastq文件轉(zhuǎn)換為fasta文件：
$ bioawk -c fastx '{print ">"$name"\n"$seq}' contam.fastq | head -n4
>DJB775P1:248:D0MDGACXX:7:1202:12362:49613
TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
>DJB775P1:248:D0MDGACXX:7:1202:12782:49716
CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG
還可以使用revcomp函數(shù)來確定序列的逆向互補序列：
$ bioawk -c fastx '{print ">"$name"\n"revcomp($seq)}' contam.fastq | head -n4
>DJB775P1:248:D0MDGACXX:7:1202:12362:49613
TTCAGACGTGTGCTCTTCCGATCTAAGCAGTGGTATCAACGCAGAGTAAGCA
>DJB775P1:248:D0MDGACXX:7:1202:12782:49716
CCGATCTAAGCAGTGGTATCAACGCAGAGTAAGCAGTGGTATCAACGCAGAG
- 最后我們統(tǒng)計一下之前下載的小鼠參考基因組fa壓縮文件的序列長度：
$ bioawk -c fastx '{print $name,length($seq)}' Mus_musculus.GRCm38.74.dna.toplevel.fa.gz | head
1 195471971
10 130694993
11 122082543
12 120129022
13 120421639
14 124902244
15 104043685
16 98207768
17 94987271
18 90702639
- 其實，
bioawk同樣可以處理tsv文件留晚，使用-t或者-c hdr作為參數(shù)酵紫。前者代表輸入為純tsv文件，而后者則代表此tsv的第一行為表頭倔丈，bioawk為自動將每一列的信息賦給以表頭命名的變量憨闰。以文件genotypes.txt為例進行說明：
id ind_A ind_B ind_C ind_D
S_000 G/G A/G A/A A/G
S_001 A/T A/T T/T T/T
S_002 C/T T/T C/C C/T
S_003 C/T C/T C/T C/C
S_004 C/G G/G C/C C/G
S_005 A/T A/T A/T T/T
S_006 C/G C/C C/G C/G
S_007 A/G G/G A/A G/G
S_008 G/T G/T T/T G/T
S_009 C/C C/C A/A A/C
這里通過bioawk尋找ind_A與ind_B列值相同的id：
$ bioawk -c hdr '$ind_A == $ind_B{print $id}' genotypes.txt
S_001
S_003
S_005
S_008
S_009
