实现这个需求的方法很多,最快捷的是使用bedtools。
还有就是写个脚本。
下面是perl简单实现一下。暂时没考虑正负链反向互补序列问题。过段时间再填坑。
还有这位博主使用的python实现的 https://www.jianshu.com/p/4c19389b2f43
,思路是一样的,只是语言换了,功能比较完善。
#!/usr/bin/perl -w
my %hash;
#先打开bed文件。
#chr start end
#bed文件存储为chr作为键,start和end作为值的hash
open FA,"<$ARGV[0]";
while (<FA>){
chomp;
my ($id,$arry) = split(/\t/,$_,2);
#print "$arry\n";
$hash{$id} = $arry;
}
#打开参考基因组文件,
#以>为分隔符。把fasta序列分为,ID和seq两个部分。
#同时把序列seq中所有的空格或者换行替换为空,也就是去掉所有空格和换行符
#接下来,判断fasta文件的染色体是否在上个步骤的hash中
#如果在。那么就按区间分割序列
#分割序列用到了substr这个函数。
#最后就是按自己习惯打印输出序列
open IN,"<$ARGV[1]";
$/ = '>';
<IN>;
while (<IN>){
chomp;
my ($chr,$seq) = split(/\n/,$_,2);
$seq =~ s/\s+//g;
#print "$chr\n$seq\n";
if (exists($hash{$chr})){
my @arry2 = split(/\s+/,$hash{$chr},2);
my $len = $arry2[1]-$arry2[0];
my $res_seq = substr($seq,$arry2[0],$len);
print "$chr:$arry2[0]-$arry2[1]\n$res_seq\n";
}
}