Perl学习18之生信简单运用（二）

"pythonic生物人"的第31篇分享。

---

原创不易，点个“赞“或"在看"鼓励下呗！

摘要

Perl计算sam文件每个染色体匹配reads数，GC碱基数，GC含量。

正文开始啦

---

gccount_reads.pl

#!/usr/bin/perl
use 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

chr7 108663 2212904 40.73%

chr20 44747 990732 44.28%

chr22 24434 586887 48.04%

chr14 63145 1293962 40.98%

chrY 1444 29104 40.31%

chr19 36744 881130 47.96%

chr8 102992 2080246 40.40%

chr1 162445 3403360 41.90%

chr11 94802 1984799 41.87%

。。。。。。。。。。。。。。。

---

同系列文章

---

原创不易，点个"点赞"、"在看"鼓励下呗