根据ID列表提取fasta序列

最新推荐文章于 2024-08-13 04:24:17 发布

原创最新推荐文章于 2024-08-13 04:24:17 发布 · 1.2w 阅读

20 ·

CC 4.0 BY-SA版权

文章标签：

#perl

perl 专栏收录该内容

14 篇文章

订阅专栏

该脚本用于根据指定的基因列表从大型FASTA文件中快速提取相应的DNA或蛋白质序列。使用Perl编写，支持压缩文件输入。

根据基因列表，从总的fasta文件中提取相关的序列，是经常遇到的问题，这个脚本很好地帮助我实现这个动作，很速度。

用法：perl $0 gene.list *fa > out

#! /usr/bin/perl -w
use strict;
die "perl $0 <lst><fa>\n" unless  @ARGV==2;
my ($lst,$fa)=@ARGV;
open IN,$lst||die;
my %ha;
map{chomp;$ha{(split)[0]}=1}<IN>;
close IN;

$fa=~/gz$/?(open IN,"gzip -cd $fa|"||die):(open IN,$fa||die);
$/=">";<IN>;$/="\n";
my %out;
while(<IN>){
	my $info=$1 if(/^(\S+)/);
	$/=">";
	my $seq=<IN>;
	$/="\n";
	$seq=~s/>|\r|\*//g;
	print ">$info\n$seq" if(exists $ha{$info} && ! exists $out{$info});
	$out{$info}=1;
}
close IN;