0

将 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 ");
4

1 回答 1

2

当您在此行中设置 '$/' 时,您更改了 perl 从 'STDIN' 读取的方式:

local $/ = ">";

解决此问题的最简单方法是在该行之前添加一个左括号,并在 '##### PART II' 注释之前添加一个右括号:

{
    local $/ = ">";
    ...
    ...
}

##### PART II

(我认为理论上,你可以在你输入的文本的末尾加上一个“>”,但这看起来很奇怪,所以我不会这样做)

这将解决您的问题。但是应该解决的问题是您所做的一些样式选择。据我所知,中间的两大块代码都是相同的,可能应该放入一个子例程中,然后调用两次。这将消除重复并且不易出错。

您还应该使用三个参数 open 调用来打开文件。

于 2014-01-05T06:05:03.713 回答