提取指定基因的fasta序列

简便好用的序列提取的perl脚本

这里,介绍一个非常简便好用的序列提取的perl脚本,用法非常简单。

用法如下:

perl /share/work/huangls/piplines/01.script/get_fa_by_id.pl <id><fa><OUT>

例如:

perl /share/work/huangls/piplines/01.script/get_fa_by_id.pl  id.txt  input.fasta  out.fa

其中 id.txt 为要提取的序列ID,input.fasta 为输入序列文件,out.fa 是输出提取的序列文件。

id.txt 格式如下:

TRINITY_DN116733_c6_g37
TRINITY_DN116733_c6_g70
TRINITY_DN95808_c0_g7
TRINITY_DN104586_c1_g2
TRINITY_DN108413_c2_g23
TRINITY_DN37223_c0_g1
TRINITY_DN107955_c0_g8
TRINITY_DN117047_c0_g2
TRINITY_DN78058_c0_g1

这里是脚本代码:

die "perl $0 <id><fa><OUT>" unless(@ARGV==3);
use Math::BigFloat;
use Bio::SeqIO;
use Bio::Seq;

$in  = Bio::SeqIO->new(-file => "$ARGV[1]" ,
                               -format => 'Fasta');
$out = Bio::SeqIO->new(-file => ">$ARGV[2]" ,
                               -format => 'Fasta');
my%keep;

open IN ,"$ARGV[0]" or die "$!";
while(<IN>){
	chomp;
	next if /^#/;
	#next unless />>/;
	my@tmp=split(/\s+/);
	$keep{$tmp[0]}=1;
}
close(IN);
my$i=0;
while ( my $seq = $in->next_seq() ) {
     my($id,$sequence,$desc)=($seq->id,$seq->seq,$seq->desc);
     if(exists $keep{$id}){
           $out->write_seq($seq);	
     }
}
$in->close();
$out->close();

脚本使用了Bio::SeqIO模块来处理序列文件,简洁而高效,先使用哈希来存储要提取的序列ID,然后利用Bio::SeqIO遍历序列文件,判断每条序列是否是要提取的序列,是的话就输出。

希望对大家有帮助!


  • 发表于 2018-07-27 09:31
  • 阅读 ( 8430 )
  • 分类:perl

0 条评论

请先 登录 后评论
安生水
安生水

335 篇文章

作家榜 »

  1. omicsgene 674 文章
  2. 安生水 335 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. 红橙子 78 文章
  6. CORNERSTONE 72 文章
  7. xun 71 文章
  8. rzx 71 文章