搜索
查看: 7377|回复: 11

【Panda姐-perl练习题9】统计碱基个数及GC%

[复制链接]

58

主题

103

帖子

756

积分

版主

Rank: 7Rank: 7Rank: 7

积分
756
QQ
发表于 2016-8-30 19:30:15 | 显示全部楼层 |阅读模式
本帖最后由 Panda姐 于 2016-8-30 19:32 编辑

题目:
编写perl脚本,用于fastq、fasta格式文件的碱基数目和GC含量的统计。

代码记录:
## 方法一、二、三、四在统计的时候会把文件中的id行出现的A、T、C、G也统计进去,所以不怎么对,但是是我的一个学习过程还是放在这里。
方法一:
[Perl] 纯文本查看 复制代码
#! /usr/bin/perl
use strict;
open INPUT,"../4_Removehost/clean.fa";
my($GC,$count_A,$count_C,$count_G,$count_T,$total);
while(<INPUT>){
while ($_=~/A/ig){$count_A++} 
while ($_=~/C/ig){$count_C++} 
while ($_=~/G/ig){$count_G++} 
while ($_=~/T/ig){$count_T++}
}
$total=$count_A+$count_T+$count_G+$count_C;
$GC=($count_G+$count_C)/$total;
print "total count is $total \nGC is $GC";

#total count is 2875363107
#GC is 0.430857696888437
#real   13m52.721s
#user   13m51.099s
#sys    0m1.520s


方法二:

[Perl] 纯文本查看 复制代码
#! /usr/bin/perl
use strict;
open INPUT,"../4_Removehost/clean.fa";
my($GC,$count_A,$count_C,$count_G,$count_T,$total);
while(<INPUT>){
$count_A=$count_A+($_=~tr/A//); 
$count_T=$count_T+($_=~tr/T//); 
$count_G=$count_G+($_=~tr/G//); 
$count_C=$count_C+($_=~tr/C//);
}
$total=$count_A+$count_T+$count_G+$count_C;
$GC=($count_G+$count_C)/$total;
print "total count is $total \nGC is $GC";
 
#这个结果比4_2_Remocehost_count.pl快多了。
#total count is 2875363107
#GC is 0.430857696888437
#real   1m3.280s
#user   1m1.762s
#sys    0m1.509s


方法三:
[Perl] 纯文本查看 复制代码
#! /usr/bin/perl
use strict;
open INPUT,"../4_Removehost/clean.fa";
my($GC,$count_A,$count_C,$count_G,$count_T,$total);
while(<INPUT>){
$count_A=$count_A+($_=~s/A//ig); 
$count_T=$count_T+($_=~s/T//ig); 
$count_G=$count_G+($_=~s/G//ig); 
$count_C=$count_C+($_=~s/C//ig);
}
$total=$count_A+$count_T+$count_G+$count_C;
$GC=($count_G+$count_C)/$total;
print "total count is $total \nGC is $GC";

#total count is 2875363107
#GC is 0.430857696888437
#real   10m23.054s
#user   10m21.412s
#sys    0m1.564s


方法四:
# 将上面最快的方法二写成一行perl代码,方便使用。
[Perl] 纯文本查看 复制代码
perl -ne  '{$count_A=$count_A+($_=~tr/A//);$count_T=$count_T+($_=~tr/T//);$count_G=$count_G+($_=~tr/G//);$count_C=$count_C+($_=~tr/C//)};END{print qq{total count is },$count_A+$count_T+$count_G+$count_C, qq{\nGC%:},($count_G+$count_C)/($count_A+$count_T+$count_G+$count_C),qq{\n} }'  ../4_Removehost/clean.fa

#total count is 2875363107
#GC:0.430857696888437
 #real    1m16.111s
#user    1m14.730s
#sys    0m1.371s


## 以下的代码不会出现将id行的A、T、C、G统计进去,是正确的。
#用于fastq格式文件的碱基数目和GC含量的统计
[Perl] 纯文本查看 复制代码
perl -ne  'if($.%4){chomp;$count_G=$count_G+($_=~tr/G//);$count_C=$count_C+($_=~tr/C//);$cur_length=length($_);$total_length+=$cur_length;}END{print qq{total count is $total_length bp\nGC%:},($count_G+$count_C)/$total_length,qq{\n} }' input.fq


# 用于fasta格式文件的碱基数目和GC含量的统计
[Perl] 纯文本查看 复制代码
perl -ne  'if($.%2){chomp;$count_G=$count_G+($_=~tr/G//);$count_C=$count_C+($_=~tr/C//);$cur_length=length($_);$total_length+=$cur_length;}END{print qq{total count is $total_length bp\nGC%:},($count_G+$count_C)/$total_length,qq{\n} }' input.fa



# 也用于fasta格式文件的碱基数目和GC含量的统计
[Perl] 纯文本查看 复制代码
grep -v '>' input.fa| perl -ne  '{$count_A=$count_A+($_=~tr/A//);$count_T=$count_T+($_=~tr/T//);$count_G=$count_G+($_=~tr/G//);$count_C=$count_C+($_=~tr/C//);$count_N=$count_N+($_=~tr/N//)};END{print qq{total count is },$count_A+$count_T+$count_G+$count_C+$count_N, qq{\nGC%:},($count_G+$count_C)/($count_A+$count_T+$count_G+$count_C+$cont_N),qq{\n} }' 






上一篇:【月亮-Perl练习题1】序列的反向互补
下一篇:小骆驼学习笔记第1章
回复

使用道具 举报

4

主题

29

帖子

119

积分

注册会员

Rank: 2

积分
119
发表于 2016-9-16 22:10:17 | 显示全部楼层
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl

use strict;
use warnings;

my %num;
chomp(my $input = <>);
my @seq = split(//, $input);

for my $word ( @seq ) {
	if($word =~ /\w/) {
	$num{$word}++;
	}
}

for ( sort keys %num ) {
	print "$_ : $num{$_}\n";
}


这样会不会略简单
回复 支持 1 反对 0

使用道具 举报

4

主题

56

帖子

553

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
553
发表于 2016-9-1 11:05:44 | 显示全部楼层
我以为你只有8道题呢,好厉害,你给的题目都是非常适合初学者的!
加油~~~
回复 支持 1 反对 0

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-9-23 20:34:31 | 显示全部楼层
看样子你是学懂了,我以前也玩过用tr来数数,也参照过一些帖子:http://blog.csdn.net/g_r_c/article/details/8020755
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
发表于 2016-10-11 08:12:22 | 显示全部楼层
赞一个,好呀!
苛求远离完美
回复 支持 反对

使用道具 举报

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
发表于 2016-10-11 17:31:10 | 显示全部楼层
楼主,为啥不能用when加continue循环??
苛求远离完美
回复 支持 反对

使用道具 举报

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
发表于 2016-10-11 17:37:10 | 显示全部楼层
我用了when加continue循环,算到和你们公布答案的结果不对!
苛求远离完美
回复 支持 反对

使用道具 举报

58

主题

103

帖子

756

积分

版主

Rank: 7Rank: 7Rank: 7

积分
756
QQ
 楼主| 发表于 2016-10-12 10:59:36 | 显示全部楼层
惠真-市二医院 发表于 2016-10-11 17:37
我用了when加continue循环,算到和你们公布答案的结果不对!

你可以把代码放上来,大家看一下
回复 支持 反对

使用道具 举报

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
发表于 2016-10-12 20:12:06 | 显示全部楼层
Panda姐 发表于 2016-10-12 10:59
你可以把代码放上来,大家看一下
哎呀,我都删了!
类似:
[AppleScript] 纯文本查看 复制代码
for(@array){
when(/A/ig){$count_A++;continue}
when(/T/ig){$count_T++;continue}
when(/C/ig){$count_C++;continue}
when(/G/ig){$count_G++}
}
苛求远离完美
回复 支持 反对

使用道具 举报

1

主题

16

帖子

120

积分

注册会员

Rank: 2

积分
120
发表于 2017-7-19 13:53:16 | 显示全部楼层
if($.%2),当到第二行时,就条件为假了,那fasta不是第二行为碱基序列吗?还有if($.%4),第一行条件也为真,那么索引序列也包括吗?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-17 07:13 , Processed in 0.038182 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.