计算全基因组GC content

2020-07-16  本文已影响0人  余绕

计算每条染色体上5000bp 一个bin的GC百分比。

代码如下:

open FA,"$ARGV[0]";

$/=">";

while(<FA>){
    chomp;
    
    my($ID,$seq)=split /\n/,$_,2;
    $seq=~s/\n//g;
    $seq=uc($seq);
    #print "$ID\n$seq\n";
    $hash{$ID}=$seq;
}

open OU,">$ARGV[1]";
foreach (keys %hash){
    $key=$_;
    $seq1=$hash{$key};
    #print  OU "$key\n";
    $len=length($seq1);
    #print OU "$len\n";
    
    $num=int($len/5000);
    #print OU "--- $num --- \n";
    $num=$num+1 if ($len/5000 > $num);
    #print OU "---- $num ---- \n";
    $start=0;
    $end=0;
    
    for ($i=1;$i<= $num;$i+=1){
        #print "$seq1\n";
    
        $seq2=substr($seq1,$start,5000);
        #print OU "$start\n$seq2\n";
        $start+=5000;
        $G=($seq2=~s/G/G/g);
        $C=($seq2=~s/C/C/g);
        $GC=($G+$C)/5000;
        #$HASH{"$key.$i"}=$GC; #可以创建hash,用于后面输出进行排序!
        print OU "$key\_$i\t$GC\n";
        
    }
    
}

#foreach (sort keys %HASH ){
#   
#   print OU "$_\t$HASH{$_}\n";
#   
#}


运行code:perl GC.pl IRGSP-1.0_genome.fasta OUT.GC
结果展示:


image.png
上一篇 下一篇

猜你喜欢

热点阅读