OUTPUT FILE

FarTech has asked for the wisdom of the Perl Monks concerning the following question:

Hi guys.

Thanks for all the people who helped me

I am trying to find the exhaustive homologous match of a given sequence against a given library.

I run the fasta34.exe from my perl script,input the sequence & library file name which is stored in the perl folder and other input values .i got a output in the fasta format, which I parse thru and get the best sequence match above a given threshold value and store the match each at a time in a file. Now use this extracted sequence file to run the fasta again against the same library and inputs and try to call the fasta34 again and repeat the same procedure till u get the exhaustive homologous match. If u find any new sequence comparded to the original one, append to the original fasta file.

Using my code I am able to extract the sequence from the orginal file which fits the criteria above a threshold using the system command in the perl with the command line i can avoid reentring the inputs

The main problem the system command with the command line is not storing in a file.

has anybody experienced in using the command line in fasta as well as comparing the original file with new file generated and appending the new contents to the original file

Following is my code

#!usr/bin/perl -w

print"\n********Running Fasta34**********  \n \n";
print"\n Please enter the file sequence  \n";
my $seqfile = <STDIN>;
print "\n Please Enter the library file name \n";
my $library = <STDIN>;

#my $library = defined $ARGV[1] ? $ARGV[1] : "$default_library"; 
# get library from command line, or use default
my $output = "c:/perl/sam/out.fasta";
open(FASTA,"$output") or die "cant open the output file \n";
$result = system("c:/perl/sam/fasta34.exe -Q $seqfile $library -o $out
+put");
#print"\n\n*****enter the output file name you have given***\n";
#$output = <STDIN>;

my $seqname;
my $iteration;
my $orig = "c:/perl/sam/orginal.fasta";
my $out2in = "c:/perl/sam/output1.aa";
my $compare = "c:/perl/sam/compare.fasta";
print "Enter your cut off percentage";
$cut = <STDIN>;
while (<FASTA>)
{
    #compare the line of matching sequence 
    if(m/>>(.{4,6})(.*)/)
    {
        $seqname = $1;
        $laterhalf = $2;
    }
    #print "SKIP A LINE iF A MATCH \n";
    next if /^ini/;
    if (m/(\d+\.\d+)% identity/) 
    {
        $per = $1;
        #check if match is above the cutoff percentage  
        if ($per > $cut)
        {
            #print "\n\n$seqname$laterhalf \n";\
            #print "identity match $per % \n";
            #store the first line of the input file 
            open(ORG,">>$orig");
            open(OUT2IN,">$out2in");
            print OUT2IN "$seqname$laterhalf \n";
            print ORG "$seqname$laterhalf \n";
            while(<FASTA>)
            {
                if (m/^($seqname)(.*)/)
                {
                    #store the match sequence
                    print OUT2IN "$2 \n";
                    print ORG "$2 \n";
                }
                if(m/>>/)
                {
                    #print "end of given loop";                
                    print "rerunning the fasta program";
                    $iteration = system("c:/perl/sam/fasta34.exe -Q $o
+ut2in $library -O $compare");
                    open(COMPARE,$compare);
                    seek(FASTA,-100,1);
                    last;
                    
                }
            }close(OUT2IN);
        }
    }
}
close (FASTA);
close (COMPARE);
print $result;
#print $iteration;
[download]

Comment on OUTPUT FILE Download Code

Replies are listed 'Best First'.
Re: OUTPUT FILE by graff (Chancellor) on Jan 28, 2005 at 04:41 UTC
Two issues in your code are likely causing trouble (or will cause trouble): In your first system call, "fasta32.exe" is supposed to write its output to "out.fasta", but before running that system call, you open that file name for input with the file handle "FASTA". You should run the system call first -- and check the return value (the exit status from the sub-process) -- then open the file to read it in perl. Since you assign the return value to $return, you would check for bad exit status like this: `if ( $return >> 8 ) { warn "fasta32.exe ended with non-zero exit status\n"; }` [download] The other problem is this embedding of while loops: `while (<FASTA>) { ... if ( /(\d+\.\d+)% identity/ ) { ... if ( $per > $cut ) { ... while (<FASTA>) { ... } ...` [download] The problem there is that once the inner while loop is reached, it will eat up everything in the file that hasn't been read yet. Once the inner while loop finishes, you will have reached the end of the file, and there will be nothing left for the outer loop to read. And then there's some really strange stuff in that inner while loop: you seem to be trying to seek backwards on the FASTA file handle, but how can you be sure that seeking backwards by 100 bytes is the right thing to do? Also, after you run "fasta32.exe" again, sending it's output to a file called "compare.fasta", you open that file, but you never read from it. And then, there may also be some trouble with some of your regexes. Have you tested them against various lines of your data files to make sure they match when you want them to (and do not match when you do not want them to)? For example, `/>>(.{4,6})(.)/` seems a little odd: if only 4, 5 or 6 characters (including spaces and CRLF) follow the double-angle-bracket, this match will succeed, and $2 will be empty. If there are more than 6 characters, everything from the 7th on (including spaces and CRLF) is assigned to $2. Is this your intention? You should start by laying out a description of your task and algorithm in clear, specific terms (e.g. always refer to particular file names, not "a file" or "another file"), and identify a sequence of distinct steps. Try doing the steps manually first, to make sure you have them specified correctly and in the proper order. Create the script in pieces and work with each piece until it does its job the right way for a given step (or some tightly bound subset of steps). If you build the pieces correctly, you'll have more flexibility in assembling them into "higher-order" tools -- not just the one you're working on now. Also, if a particular piece gives you a lot of trouble, it will be easier to post the smaller chunk of code, together with a small sample of the data that is problematic. update:* forgot to mention: if you're going to be reading things like file names and numeric values from STDIN -- and for that matter, if you're going to be manipulating lines of data that you read from a text file -- you need to start using "chomp", to remove the CRLF from the end of each string.	[reply] [d/l] [select]

Replies are listed 'Best First'.

Re: OUTPUT FILE
by graff (Chancellor) on Jan 28, 2005 at 04:41 UTC

In your first system call, "fasta32.exe" is supposed to write its output to "out.fasta", but before running that system call, you open that file name for input with the file handle "FASTA". You should run the system call first -- and check the return value (the exit status from the sub-process) -- then open the file to read it in perl. Since you assign the return value to $return, you would check for bad exit status like this:

    if ( $return >> 8 ) {
       warn "fasta32.exe ended with non-zero exit status\n";
    }
[download]

while (<FASTA>)
{
   ...
   if ( /(\d+\.\d+)% identity/ ) {
      ...
      if ( $per > $cut ) {
         ...
         while (<FASTA>)
         {
            ...
         }
...
[download]

And then there's some really strange stuff in that inner while loop: you seem to be trying to seek backwards on the FASTA file handle, but how can you be sure that seeking backwards by 100 bytes is the right thing to do? Also, after you run "fasta32.exe" again, sending it's output to a file called "compare.fasta", you open that file, but you never read from it.

And then, there may also be some trouble with some of your regexes. Have you tested them against various lines of your data files to make sure they match when you want them to (and do not match when you do not want them to)?

For example, />>(.{4,6})(.*)/ seems a little odd: if only 4, 5 or 6 characters (including spaces and CRLF) follow the double-angle-bracket, this match will succeed, and $2 will be empty. If there are more than 6 characters, everything from the 7th on (including spaces and CRLF) is assigned to $2. Is this your intention?

You should start by laying out a description of your task and algorithm in clear, specific terms (e.g. always refer to particular file names, not "a file" or "another file"), and identify a sequence of distinct steps. Try doing the steps manually first, to make sure you have them specified correctly and in the proper order. Create the script in pieces and work with each piece until it does its job the right way for a given step (or some tightly bound subset of steps).

If you build the pieces correctly, you'll have more flexibility in assembling them into "higher-order" tools -- not just the one you're working on now. Also, if a particular piece gives you a lot of trouble, it will be easier to post the smaller chunk of code, together with a small sample of the data that is problematic.

update: forgot to mention: if you're going to be reading things like file names and numeric values from STDIN -- and for that matter, if you're going to be manipulating lines of data that you read from a text file -- you need to start using "chomp", to remove the CRLF from the end of each string.

[reply]
[d/l]
[select]