其他
Perl学习18之生信简单运用（二）
"pythonic生物人"的第31篇分享
原创不易，点个“赞“或"在看"鼓励下呗
摘要
Perl计算sam文件每个染色体匹配reads数，GC碱基数，GC含量。正文开始啦
gccount_reads.pl
#!/usr/bin/perluse strict;use warnings;
my (%hash,%GC,%ATGC);open IN,"<",$ARGV[0];while(<IN>){ chomp; next if /^@/;#跳过@开头的行 my @line=split; next if $line[2] eq "*";#跳过未匹配上的序列 next if $line[4] < 25;#排除比对质量小于25的序列
#计算ATGC的碱基个数 my $G = ($line[9] =~ s/G/G/g); my $C = ($line[9] =~ s/C/C/g); my $A = ($line[9] =~ s/A/A/g); my $T = ($line[9] =~ s/T/T/g); my $GC = $G + $C; #计算每个染色体匹配到的read数，GC碱基个数，总碱基个数 $hash{$line[2]} ++; $GC{$line[2]} += $GC; $ATGC{$line[2]} += ($GC + $A + $T);}close IN;
foreach (keys %hash){ printf "$_\t$hash{$_}\t$GC{$_}\t%.2f%%\n",($GC{$_}/$ATGC{$_})*100;
}perl gccount_reads.pl *sam