将 BLAST 命令组合到 perl 脚本时,我被卡住了。问题是第二部分开始时命令行暂停了。
PART I 用于裁剪 fasta 序列。PART II 用于对 PART I 生成的文件进行 BLAST。两个部分单独运行都可以,但组合在一起时遇到了“暂停”问题。
我猜是因为第一部分生成的 $ARGV[1] 和 $ARGV[3] 不能在第二部分中使用。我不知道如何解决,虽然我尝试了很多。
谢谢!
#! /usr/bin/perl -w
use strict;
#### PART I
die "usage:4files fasta1 out1 fasta2 out2\n" unless @ARGV==4;
open (S, "$ARGV[0]") || die "cannot open FASTA file to read: $!";
open OUT,">$ARGV[1]" || die "no out\n";
open (S2, "$ARGV[2]") || die "cannot open FASTA file to read: $!";
open OUT2,">$ARGV[3]" || die "no out2\n";
my %s;# a hash of arrays, to hold each line of sequence
my %seq; #a hash to hold the AA sequences.
my $key;
print "how long is the N-terminal(give number,e.g. 30. whole length input \"0\") \n";
chomp(my $nl=<STDIN>);
##delete "\n" for seq.
local $/ = ">";
<S>;
while (<S>){ #Read the FASTA file.
chomp;
my @line=split/\n/;
print OUT ">",$line[0],"\n";
splice @line,0,1;
#print OUT join ("",@line),"\n";
#@line = join("",@line);
#print @line,"\n";
if ($nl == 0){ #whole length
my $seq=join("",@line);
my @amac = split(//,$seq);
splice @amac,0,1; # delete first "MM"
#push @{$s{$key}},@amac;
print OUT @amac,"\n";
}
else { # extract inital aa by number ##Guanhua
my $seq=join("",@line);
#print $seq,"\n";
my @amac = split(//,$seq);
splice @amac,0,1; # delete first "MM"
splice @amac,$nl; ##delete from the N to end
#print @amac,"\n";
#push (@{$s{$key}}, @amac);
print OUT @amac,"\n";
}
}
<S2>;
while (<S2>){ #Read the FASTA file.
chomp;
my @line=split/\n/;
print OUT2 ">",$line[0],"\n";
splice @line,0,1;
#print OUT join ("",@line),"\n";
#@line = join("",@line);
#print @line,"\n";
if ($nl == 0){ #whole length
my $seq=join("",@line);
my @amac = split(//,$seq);
splice @amac,0,1; # delete first "MM"
#push @{$s{$key}},@amac;
print OUT2 @amac,"\n";
}
else { # extract inital aa by number ##Guanhua
my $seq=join("",@line);
#print $seq,"\n";
my @amac = split(//,$seq);
splice @amac,0,1; # delete first "MM"
splice @amac,$nl; ##delete from the N to end
#print @amac,"\n";
#push (@{$s{$key}}, @amac);
print OUT2 @amac,"\n";
}
}
##### PART II
print "nucl or prot?\n";
chomp(my $tp = <STDIN>);
system ("makeblastdb -in $ARGV[1] -dbtype prot");
system ("makeblastdb -in $ARGV[3] -dbtype $tp");
print "blast type? (blastp,blastn)\n";
chomp(my $cmd = <STDIN>);
system ("blastp -query $ARGV[1] -db $ARGV[3] -outfmt 6 -evalue 1e-3 -out 12.out ");
system ("$cmd -db $ARGV[1] -query $ARGV[3] -outfmt 6 -evalue 1e-3 -out 21.out ");