Question

I was stuck when combining the BLAST command into perl script. The problem is that the command line paused when the PART II begin.

PART I is used to crop the fasta sequence. PART II is used to do BLAST with the file generated by PART I. Both the two parts can run well individually, but met the "pause" problem when combining together.

I guess it was because the $ARGV[1] and $ARGV[3] generated by part I cannot be used in part II. I dont know how to fix, though I tried a lot.

Thanks!

#! /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 ");
Was it helpful?

Solution

You changed the way perl reads from 'STDIN' when you set '$/' in this line:

local $/ = ">";

The easiest way to fix this is to add a left bracket right before that line and a right bracket just before the '##### PART II' comment:

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

##### PART II

(I think theoretically, you could put a ">" at the end of the text you input, but that seems strange, so I wouldn't do it)

That will fix your problem. But something that should be addressed is some of the style choices you made. The two big chunks of code in the middle are both identical as far as I can tell and should probably be put into a subroutine and then called twice. This will eliminate duplication and is less error prone.

You should also use the three argument open call to open files.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top