35.《Bioinformatics Data Skills》之

作为表格文本处理命令awk的延伸，Heng Li开发了专门处理生物信息学的文件（例如bed, sam, vcf, gtf等）的命令bioawk，我们了解一下。

---

简介

bioawk发布在github上，可以通过如下命令进行下载，编译与安装：

git clone git://github.com/lh3/bioawk.git && cd bioawk && make && sudo cp bioawk /usr/local/bin/

bioawk与awk命令的使用方式类似（若不了解awk可以见上节内容），命令最大的特点在于-c参数，通过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也可以）。

---

实例

本节例子数据地址

1. 以Mus_musculus.GRCm38.75_chr1.gtf为例，假如我们想要统计所有蛋白编码基因的长度的话，可以采用如下命令：

$ 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

2. -c设置为fastx的话，bioawk会自动判断输入文件类型为fastq还是fasta。例如我们统计fastq文件的行数：

$ bioawk -c fastx 'END{print NR}' contam.fastq
8

3. 可以通过如下命令将fastq文件转换为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函数来确定序列的逆向互补序列：

$ 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

4. 最后我们统计一下之前下载的小鼠参考基因组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

---

5. 其实，bioawk同样可以处理tsv文件，使用-t或者-c hdr作为参数。前者代表输入为纯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