方法一：linux下用awk计算fasta序列的长度
前面发表一篇文章《用awk和sed快速将fasta格式的序列改成一行显示》，其实我的这种方法就是在这基础上进行的。加入已经有一个fasta文件为contig.fa，文件中的序列如下：
04 |
GAATAGTGATACAAATTATATAGTTTCAAGTATGTGACTTGAACATGAGATTAT |
06 |
TAATCTAGGCTTGAAACTATATAATTTGTATCACTATTCTAAGGATTTTTTT |
08 |
TATTCATCTTTGCACTACGTTCATCTCAA |
10 |
TCCGTTGTGGGGTCCACCAATGATTAAAACGAATATTCCC |
首先通过上面的命令将fasta序列转换成一行显示，命令如下：
1 |
awk '/^>/&&NR>1{print
"";}{ printf "%s",/^>/ ? $0" ":$0
}' contig.fa |
得到如下结果：
http://upload.plob.ybzhao.com/wp-content/uploads/2012/10/f1-600x115.jpg
如果想直接显示每条序列的长度，可以运行如下命令：
1 |
awk '/^>/&&NR>1{print
"";}{ printf "%s",/^>/ ? $0" ":$0
}' contig.fa
|awk '{print
$1"\t"length($3)}' |
得到结果如下：
方法二：利用bioperl计算fasta序列长度
上面的方法是基于linux计算的，直接输出结果。但是有是有计算fasta序列的长度只是程序某一个小的操作步骤，那我们可以采用下面的方法.
首先，确定bioperl正确安装了。
然后再perl中利用如下的代码：
05 |
my $in=Bio::SeqIO->new(-file=>"$file",-format=>"fasta"); |
06 |
while ($seq=$in->next_seq()) |
08 |
$hash{$seq->id}=length($seq->seq()); #
length($seq->seq()) 计算的是序列长度，序列的长度被存入hash表中 |
09 |
print $seq->id."\t".$seq->seq()."\n";#
直接输入，输出的结果与上面awk的方法是一致的 |
这样每一条序列的长度就被存入以其序列名字为key的hash表中
加载中，请稍候......
加载中…
(2016-11-16 22:42:18) 
