滑窗口统计基因组GC含量的分布

最新推荐文章于 2025-04-30 15:50:23 发布

xflee0608

最新推荐文章于 2025-04-30 15:50:23 发布

阅读量8.2k

点赞数

CC 4.0 BY-SA版权

分类专栏： perl 文章标签： perl

本文链接：https://blog.youkuaiyun.com/hugolee123/article/details/38441927

perl 专栏收录该内容

14 篇文章

订阅专栏

本文提供了一个使用Perl脚本在50k窗口大小下统计香蕉基因组GC含量的方法。该脚本适用于常见的基因组序列文件，通过逐段读取序列并计算每段内的GC碱基数量来得出GC含量比例。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

在基因组学分析中，滑窗口统计GC含量是一个常规动作，笔者在这里提供一个50k窗口统计香蕉基因组GC含量的脚本。输入文件是常见的基因组序列。

#! /usr/bin/perl -w
use strict;
die "#usage:perl $0 <BananaB.chr_V2.1.final.fa>\n" unless @ARGV==1;
my $fa=shift;
my $bin=50000; ##50k
open IN,$fa||die;
$/=">";<IN>;$/="\n";
print "#Chr\tStart\tEnd\tGC_num\tRatio\n";
while(<IN>){
	my $chr=$1 if /^(\S+)/;
	$/=">";
	chomp(my $seq=<IN>);
	$/="\n";
	$seq=~s/\n+//g;
	my $len=length$seq;
	for (my $i=0;$i<$len/$bin;$i++){
		my $loc=$i*$bin;
		my $sub_fa=uc(substr($seq,$loc,$bin));
		my $GC=$sub_fa=~tr/GC//;
		my $ratio=sprintf "%.4f",$GC/$bin;
		my $start=$i*$bin+1;
		my $end=($i+1)*$bin;
		my $out=join "\t",$chr,$start,$end,$GC,$ratio;
		print $out,"\n";
		#print "$out\t$sub_fa\n" unless $GC;
	}	
}
close IN;