李恒大神的又一杰作:bioawk让你轻松搞定生物序列分析






李恒大神的又一杰作:bioawk让你轻松搞定生物序列分析

小果  生信果  2023-12-22 19:00:22

大家好,小果今天想和大家分享一个很有用的工具,叫做bioawk。
Bioawk是生物信息学大神李恒基于awk编辑的。在此必须膜拜以下!做过基因组二代测序分析的生信人几乎都用过BWA,他就是BWA的第一作者,他还是SAMtools、MAQ、TreeSoft和TreeFam等著名生物信息软件的核心作者!
李恒大神不仅在软件开发方面有着卓越的能力,还在生物信息学理论和应用方面有着深厚的造诣。他发表了超过200篇高水平的论文,其中包括Nature、Science、Cell等顶级期刊。他的论文被引用了超过2万次,H指数达到了70。他还获得了很多荣誉和奖励,例如:中国青年科技奖、国家杰出青年科学基金、长江学者特聘教授、中国生物信息学会理事长等。    
bioawk是一个基于awk的文本处理程序,专门用于处理生物信息学相关的数据格式,比如FASTA,FASTQ,GFF,SAM、BED、VCF等和带有列名的TAB分隔格式。bioawk可以让你快速地对这些数据进行过滤,转换,统计和分析,而不需要写复杂的脚本或者使用其他的软件。bioawk的语法和awk基本一致,但是增加了一些针对生物信息学数据的特性,比如可以直接读取序列的长度,质量,注释等信息,也可以方便地对序列进行反向互补,翻译等操作。bioawk还支持多种输出格式,比如CSV,TSV,JSON等,方便你将结果导入其他的程序或者平台进行进一步的分析。
接下来,小果将给大家介绍一些bioawk的基本用法和常见的应用场景。如果你对bioawk感兴趣,或者想要提高你的文本处理能力,那么请继续阅读吧。    

bioawk的安装

bioawk是一个开源的软件,你可以从GitHub上克隆它的源代码,并且按照说明进行编译和安装。你也可以使用conda或者brew等包管理器来安装bioawk。
源代码下载网址:https://github.com/lh3/bioawk
安装命令:
$ git clone git://github.com/lh3/bioawk.git $ cd bioaw$ make或者connda安装:$:connda install -p ~/YOUR/conda path/ bioawk -y

bioawk的基本语法
 

Bioawk是一个命令行工具,它的用法是:./bioawk [-F fs] [-v var=value] [-c fmt] [-tH] [-f progfile | ‘prog’] [file …]。下面是各个选项的解释:
·-F fs:指定输入字段分隔符。默认为“空格”。    
·-v var=value:在程序开始运行之前,将变量var设置为值value。
·-c fmt:指定输入格式。当使用此选项时,bioawk将根据格式或输入的第一行无缝添加命名字段的变量,具体取决于fmt。此选项还使bioawk能够读取gzip压缩的文件。
·-t:使用此选项等同于bioawk -F’t’ -v OFS=”t”。
·-H:打印输出标题行。
·-f progfile:从文件progfile中读取程序。
·’prog’:直接在命令行中指定程序。
·[file …]:指定要处理的文件。如果未指定,则从标准输入读取。
我们可以用bioawk -c -h命令列出支持的格式:
bed:     1:chrom 2:start 3:end 4:name 5:score 6:strand 7:thickstart 8:thickend 9:rgb 10:blockcount 11:blocksizes 12:blockstartssam:    1:qname 2:flag 3:rname 4:pos 5:mapq 6:cigar 7:rnext 8:pnext 9:tlen 10:seq 11:qualvcf:    1:chrom 2:pos 3:id 4:ref 5:alt 6:qual 7:filter 8:infogff:    1:seqname 2:source 3:feature 4:start 5:end 6:score 7:filter 8:strand 9:group 10:attributefastx:       1:name 2:seq 3:qual 4:comment
bioawk的基本语法和awk一样,都是`bioawk ‘pattern {action}’ file`。pattern是一个匹配条件,用于选择符合条件的记录(record)。action是一个执行动作,用于对选择的记录进行处理。file是一个输入文件,通常是一个生物信息学相关的数据文件。如果没有指定file,那么bioawk会从标准输入中读取数据。    

bioawk的常见应用
 

下面我将给大家展示一些bioawk在生物信息学中的常见应用。这些只是一些简单的例子,你可以根据你自己的需求和想象力来发挥bioawk的更多功能。
1. 统计文件信息。
当生信小白辛辛苦苦终于拿到我们的测序结果时,会发现得到了一堆看不懂的东西/(ㄒoㄒ)/~~,如何将这一堆让人脑壳痛的数据转换成我们可以看懂的数据呢?bioawk可以实现fastx数据的统计。如果我们有一个名为test.fastq的FASTQ文件,内容如下:          
如果我们想统计测序数据的基本信息,我们可以使用下面的命令来打印出每条序列的名称、GC比和长度:          
$ bioawk -c fastx '{ print $name, gc($seq),length($seq) }' test.fastq
输出结果如下:          
这里我们使用了`-c fastx`选项来告诉bioawk我们的输入文件是一个FASTA或者FASTQ格式的文件。这样bioawk就会自动将每条序列作为一个记录,并且提供一些内置的变量来访问序列的信息,比如$name, $seq, $qual等。我们还使用了length()函数来计算序列的长度,并且用print语句来打印出结果。    
我们还可以对gff文件进行信息统计。比如说,如果我们想要计算一个GFF文件中每种特征(feature)的个数和总长度。我们可以使用下面的命令来实现:          
bioawk -c gff '{count[$3]++; len[$3] += $5 - $4 + 1} END {for (f in count) print f, count[f], len[f]}' test.gff
输出结果:
这里我们使用了-c gff选项来告诉bioawk我们的输入文件是一个GFF格式的文件(gff3或者gtf为文件)。这样bioawk就会自动将每条注释作为一个记录,并且提供一些内置的变量来访问注释的信息,比如$1, $2, $3等。我们还使用了数组(array)来存储每种特征的个数和总长度,并且用END块来在处理完所有记录后打印出结果。          
2. 过滤序列
我们得到的测序数据是存在参差的,为例得到更好的结果,根据一些条件来过滤掉一些不需要的序列是很必要的。比如说,我们想要过滤掉长度小于100bp或者质量低于20的序列。我们可以使用下面的命令来实现:          
$ bioawk -c fastx 'length($seq) >= 100 && meanqual($qual) >= 35' test.fastq > filtered_test.fastq
这里我们使用了meanqual()函数来计算序列的平均质量,仅保留平均质量不低于35的序列,并且用&&符号来表示逻辑与(and)的关系。我们还使用了重定向符号(>)来将结果输出到一个新的文件中。          
3. 转换格式    
在生信分析的过程中,,我们很多时候需要将一种格式的文件转换成另一种格式的文件。比如说,我们想要将一个FASTQ文件转换成一个CSV文件,每行包含序列的名称,长度和GC含量。我们可以使用下面的命令来实现:          
bioawk -c fastx -v OFS=',' '{print $name, length($seq), gc($seq)}' test.fastq > test.csv
结果如下:          
这里我们使用了-v选项来设置输出字段分隔符(OFS)为逗号(,)。我们还使用了print语句来打印出结果。
Bioawk新增加了一个功能,可以直接从fasta/q文件中获得序列的反向互补序列,命令如下:
bioawk -c fastx '{ print ">"$name ORS revcomp($seq) }' test.fastq > rtest.fasta
结果如下:
在此命令中,'{ print “>”$name ORS revcomp($seq) }’ 是一个 awk 脚本,它用于提取每个记录的名称和序列,并计算序列的互补序列。$name 变量表示记录的名称,$seq 变量表示序列。revcomp 函数用于计算序列的互补序列。ORS 是输出记录分隔符,它的默认值是换行符。
Bioawk也可以从BAM文件中提取FASTA序列。还可以添加一些筛选条件,获得我们想要的高质量序列。比如如果想要得到FLAG为16的序列的方向互补序列,则使用revcomp。命令如下:    
samtools view test.bam | bioawk -c sam '{s=$seq; if(and($flag, 16)) {s=revcomp($seq)} print ">"$qname"n"s}'
结果文件:
      总之,bioawk是一个非常强大的文本处理工具,它可以帮助我们快速地处理各种生物信息学相关的数据格式。它的语法简单易学,功能强大,可以让我们轻松地完成各种复杂的数据分析任务。如果您是一名生物信息学工作者,那么bioawk一定是您不可或缺的工具之一。希望这篇文章能够帮助你提高你的文本处理能力,让你在生物信息学中更加游刃有余。如果你想要了解更多关于bioawk的信息,你可以访问它的GitHub页面,或者查看它的手册。如果你有任何问题或者建议,欢迎在评论区留言或者联系我。谢谢你的阅读,下次再见。
参考文档:https://github.com/lh3/bioawk/blob/master/README.md


往期推荐

1.搭建生信分析流水线,如工厂一样24小时运转Snakemake——进阶命令
2.比blast还优秀的序列比对工具?HMMER来了
3.对单细胞分析毫无头绪?让popsicleR领你入门
4.小果带你绘制ROC曲线评估生存预测能力
5.软件包安装、打怪快又好,1024G存储的生信服务器;还有比这更省钱的嘛!!!