从参考基因组文件中获得指定区间内的序列

实现这个需求的方法很多，最快捷的是使用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";
        }
}