用法:
perl getseqbyid.pl <id.list> <database.fasta> <output.fasta>
备注:这是perl脚本,<id.list>是 目标基因名称列表的txt文件,一行一个基因名称。 <database.fasta> 是你的整个fasta文件。<output.fasta> 是输出文件的名称。
代码如下:
[Perl] 纯文本查看 复制代码 #!/usr/bin/perl -w
use strict;
if(@ARGV==0){print "perl $0 <id.list> <database.fasta> <output.fasta>\n"; exit 1}
my $list = $ARGV[0];
my $data = $ARGV[1];
my $outs = $ARGV[2];
my @list = ();
open (LIST,$list) or die "Cannot open file $list: $!\n";
while (my $id = <LIST>) {
chomp $id;
$id =~ s/\s//g;
push @list, $id
}
close LIST;
my %seq = ();
my $sid = ();
open (IN, $data) or die "Cannot open file $data: $!\n";;
while (<IN>) {
if (/^\>(\S+)/) {
$sid = $1;
my @w = split /\|/, $sid;
if (@w > 2) {
$sid = $w[2];
} else {
$sid = $w[0];
}
} else {
$seq{$sid} .= $_;
}
}
close IN;
open (OUT, ">$outs") or die "Cannot create file $outs: $!\n";
foreach my $id (@list) {
print OUT ">$id\n";
print OUT $seq{"$id"};
}
close OUT;
脚本如下附件,使用的时候请把txt拓展名去掉。
|