计算每条染色体上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
网友评论