每一生信perl练习

Posted 2021-04-17 沈梦圆

1 数据的准备

• hg19基因组序列下载

• 测试数据 test.fa

存放在我的工作目录：/d/生信技能树-视频直播/第二讲

2 题目解读

 
求：hg19 每条染色体长度，每条染色体N的含量，GC含量。 
解题要点：G(g)C(c)N的计数；序列长度获取；哈希。

3 解题过程

这道题目我八月份的时候做过，但是是对一条序列的统计，分享在论坛上， 
我尝试了不同的碱基计数方法，结论是方法二更快点，但是好像忽略了小写碱基了。因为我以前做过这道题下面的方法肯定是对的，所以跳过测试数据，直接采用chr1.fa来重新进行用时测试。（写脚本来解决数据处理问题，一定要用已知答案的测试来测试调试代码，要不然直接上实验数据，报错很难找到原因的，切身感悟。）

• 方法1：1.pl

• 方法2：2.pl

• 方法三：3.pl

运行上面3种计算碱基的方法，还是方法2最快啦~ 
如果只需要看一条序列的长度、N及GC含量，写到这里就够了稍微修改下，如下：

• 统计单条序列：single.pl

• 然后改成perl单行命令行

perl -ne  'if(/>/){print $_;}else{chomp;$count_G=$count_G+($_=~tr/Gg//);$count_C=$count_C+($_=~tr/Cc//);$count_N=$count_N+($_=~tr/Nn//);$total_length+=length($_);}END{print qq{total count is $total_length bp
N%:},$count_N/$total_length,qq{
GC%:},($count_G+$count_C)/$total_length,qq{
} }' chr1.fa

显然题目是算多条序列的，所以要使用上哈希，如下：

运行时间大概是1分半吧，用普通笔记本电脑。

结果如下：不知道是不是完全正确的~_~这里的GC含量的计算，不知道分母要不要加上N，我一直都是算的总长度的。

 4 学习别人的代码

###############
#  直接统计一个序列的GC含量
###############
perl -e '$seq=shift;chomp $seq;$GC=$seq=~tr/GC/GC/;print $GC/length($seq),qq{
}' CTTTCATCTTCCCTGCAAAGAGTGTCAATAGATCTCCATGGCCGCGATTTCT

##########
#  输出fasta文件每个序列对应的长度
############
perl -0076 -ane '@F=map{s/[>
]//gr}@F;$id=shift @F;print $id,qq{	},length (join q{},@F),qq{
} if $id' test.fa

##########
#  输出fasta文件每个记录的A T G C 字数统计
###########
perl -0076 -ane '@F=map{s/[>
]//gr}@F;$id=shift @F;map {$aC+=($_=~tr/Aa//);$cC+=($_=~tr/Cc//);$gC+=($_=~tr/Gg//);$tC+=($_=~tr/Tt//)} @F;print qq{$id	$aC	$tC	$gC	$cC
} if($id);' test.fa

后面两条代码的做法说是比较危险，比较耗内存？我用第三条代码试了一下发现时间比上面的方法多了一分钟左右，稍微那么慢了点。不采用也可学习解题思想嘛~map函数~

5 我录了个视频

请忽略

~ 我自己录完也没看过 ~

有问题请联系我

个人微信ID：
Shenmengyuan1993