linux系统fasta程序,快速计算fasta序列长度的方法

最近看了一下进入PLoB的网页来路分析，看到有同学搜索计算fasta序列长度。其实自己在之前的数据分析中也遇到过相关的问题，这里给大家分享两种我常用的方法。

方法一：linux下用awk计算fasta序列的长度

前面发表一篇文章《用awk和sed快速将fasta格式的序列改成一行显示》，其实我的这种方法就是在这基础上进行的。加入已经有一个fasta文件为contig.fa，文件中的序列如下：>1 cvg_0.0_tip_0

ATTTTGGCTTTGGAAGGGC

>3 cvg_0.0_tip_0

GAATAGTGATACAAATTATATAGTTTCAAGTATGTGACTTGAACATGAGATTAT

>5 cvg_0.0_tip_0

TAATCTAGGCTTGAAACTATATAATTTGTATCACTATTCTAAGGATTTTTTT

>7 cvg_0.0_tip_0

TATTCATCTTTGCACTACGTTCATCTCAA

>9 cvg_0.0_tip_0

TCCGTTGTGGGGTCCACCAATGATTAAAACGAATATTCCC

>11 cvg_0.0_tip_0

GGAATATTCGTTTTAACAGGGAATATTCGTAGATGGCACAA

>13 cvg_0.0_tip_0

AGAAATAAATAAATTAAATAAAGTGATGTTTCTAATTTATTAAGGAAATTAA

>15 cvg_0.0_tip_0

GAAAGGACCAGACATCAATTATTATTGAAATAAATGTCAATTTT

>17 cvg_0.0_tip_0

GTTAATTACCCGATTGGTCAATATAACCTCCAGACATCAATTATTATTG

>19 cvg_0.0_tip_0

GATTATTTTTTATAACCTCCAGACA

首先通过上面的命令将fasta序列转换成一行显示，命令如下：awk '/^>/&&NR>1{print "";}{ printf "%s",/^>/ ? $0" ":$0 }' contig.fa

得到如下结果：

如果想直接显示每条序列的长度，可以运行如下命令：awk '/^>/&&NR>1{print "";}{ printf "%s",/^>/ ? $0" ":$0 }' contig.fa |awk '{print $1"\t"length($3)}'

得到结果如下：>1 19

>3 54

>5 52

>7 29

>9 40

>11 41

>13 52

>15 44

>17 49

>19 25

方法二：利用bioperl计算fasta序列长度

上面的方法是基于linux计算的，直接输出结果。但是有是有计算fasta序列的长度只是程序某一个小的操作步骤，那我们可以采用下面的方法.

首先，确定bioperl正确安装了。

然后再perl中利用如下的代码：use Bio::SeqIO;

my $file;

my $seq;

my %hash

my $in=Bio::SeqIO->new(-file=>"$file",-format=>"fasta");

while ($seq=$in->next_seq())

{

$hash{$seq->id}=length($seq->seq()); # length($seq->seq()) 计算的是序列长度，序列的长度被存入hash表中

print $seq->id."\t".$seq->seq()."\n";# 直接输入，输出的结果与上面awk的方法是一致的

}

这样每一条序列的长度就被存入以其序列名字为key的hash表中