快速计算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序列转换成一行显示,命令如下:
1 | awk '/^>/&&NR>1{print "";}{ printf "%s",/^>/ ? $0" ":$0 }'contig.fa |
得到如下结果:
如果想直接显示每条序列的长度,可以运行如下命令:
1 | 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表中
相关推荐:
- NCBI在线BLAST使用方法与结果详解 2938
- 神经网络术语:Epoch、Batch Size和迭代 527
- Consed的安装与使用教程 465
- 陈连福的NGS生物信息学培训教材V2.1 277
- WGCNA分析使用教程 272
最新创建圈子
-
原料药研发及国内外注册申报
2019-01-25 10:41圈主:caolianhui 帖子:33 -
制药工程交流
2019-01-25 10:40圈主:polysciences 帖子:30 -
健康管理
2019-01-25 10:40圈主:neuromics 帖子:20 -
发酵技术
2019-01-25 10:39圈主:fitzgerald 帖子:17 -
医学肿瘤学临床试验
2019-01-25 10:39圈主:bma 帖子:58
awk是linux自己带的命令应该是用C或者C++写的。我觉得要快一些,最主要的是方便。当然尽量少用管道符。多次套用的话,肯定会有点慢。