my code stock.com

Ben Johnson

Ortholog reciprocal blast search script
by Ben Johnson

Brief How To Use…
-Click download snippet and save (change name for quicker use)
-(optional) Open the file in notepad and change the variables $maxEvalue and $minBitScore to suit your search.
-Open cmd, and correct partition (i.e. C:\Users\Ben\Desktop>)
-Type "perl 'example.pl' '(your blast file name.txt/.blast)' 'your blast headers file.txt)' '(name you want to save reciprocal hits.txt)'
-Example run script: perl script.pl Myxo_DK1622.txt Myxo_DK1622_headers.txt Mxanthus.txt

P.S. - Typing perl in the script will depend on operating software and version of perl used.

Snippet options

Download: Download snippet as ortholog-reciprocal-blast-search-script.pl.
Copy snippet: For this you need a free my code stock.com account.
Embed code : You will find the embed code for this snippet at the end of the page, if you want to embed it into a website or a blog!

#!/usr/local/bin/perl
use strict;

print @ARGV;


   my $blastfile = $ARGV[0];
   open(my $in1,  "<",  $blastfile)  or die "Can't open $blastfile: $!";

   my @blast = [[]];
   my $lineNum = 0;
   while ( my $line=<$in1> ) {
      chomp($line);
      my @words = split (/[ \t]+/,$line);

      my $wordNum = 0;
      for my $word (@words) {
         $blast[$lineNum][$wordNum] = $word;
         $wordNum++;
      }
      $lineNum++;
   }   

   my $nlines = $lineNum;
   print "Number of lines read from $blastfile is $nlines";


   my $headersfile = $ARGV[1];
   open(my $in2,  "<",  $headersfile)  or die "Can't open $headersfile: $!";

   my %giToDesc = {};
   $lineNum = 0;
   while ( my $line=<$in2> ) {
      chomp($line);
      my @words = split (/[' ']+/,$line);
      my $giNumber = shift(@words);
      my $description = join(' ',@words);
      #print "GI $giNumber\n";
      #print "DS $description\n";
      $giToDesc{$giNumber} = $description;
      $lineNum++;
   }
   print "Number of lines read from $headersfile is $lineNum";

   #foreach my $giNumber ( keys %giToDesc ) {
   #   print "HASH.$giNumber.$giToDesc{$giNumber}\n";
   #}      
  
  open MYFILE, ">", $ARGV[2];
  
	print MYFILE "Number of lines read from $blastfile is $nlines\n";
	print MYFILE "Number of lines read from $headersfile is $lineNum\n"; 
	chomp;
 
   my $subjNo = 0;
   my $qryNo = 1;
   my $evalNo = 10;
   my $bsNo = 11;

   my $maxEvalue = 0.0000000000000001;
   my $minBitScore = 300;
   
   print MYFILE "Max evalue was $maxEvalue\n";
   print MYFILE "Min bit score was $minBitScore\n";
   chomp;
   
 
   my $numOfHits = 0.0;
   for (my $j=0; $j<$nlines; $j++) {

      #print "$j subj $blast[$j][$subjNo] qry $blast[$j][$qryNo] eval $blast[$j][$evalNo]\n";

      my $subject = $blast[$j][$subjNo];
      my $query = $blast[$j][$qryNo];
      my $evalue =  $blast[$j][$evalNo];
      my $bitscore = $blast[$j][$bsNo];

      my $subDesc = $giToDesc{$subject};
      my $qryDesc = $giToDesc{$query};
      #print "HERE: $subDesc $qryDesc : $subject $query \n";

      if ($subject ne $query) {
         if ($evalue < $maxEvalue )  {
            if ($bitscore > $minBitScore ) {
             
               print MYFILE "---Possible:$j\t$numOfHits\t$subject\t$query\t$evalue\t$bitscore\n";
			   print "---Possible:$j\t$numOfHits\t$subject\t$query\t$evalue\t$bitscore\n";
			   
               for (my $k=0; $k<$nlines; $k++) {
                  if( ($blast[$k][$subjNo] eq $query) && ($blast[$k][$qryNo] eq $subject) ) {
                      my $recipEvalue =  $blast[$k][$evalNo];
                      my $recipBitscore = $blast[$k][$bsNo];
  
                      if ($recipEvalue < $maxEvalue )  {
                         print MYFILE "Reciprocal:$j\t$numOfHits\t$query\t$subject\t$recipEvalue\t$recipBitscore\n";
                         print MYFILE "Reciprocal: Subject: $subDesc\n"; 
                         print MYFILE "Reciprocal: Query: $qryDesc\n";
                      }
                  }
               }

               $numOfHits++;                  
            }
         }
      }   
   }

   print "Max evalue was $maxEvalue\n";
   print "Min bit score was $minBitScore\n";
   print "Number of lines read from $headersfile is $lineNum\n";

Create a free my code stock.com account now.

my code stok.com is a free service, which allows you to save and manage code snippes of any kind and programming language. We provide many advantages for your daily work with code-snippets, also for your teamwork. Give it a try!

Find out more and register now

You can customize the height of iFrame-Codes as needed! You can find more infos in our API Reference for iframe Embeds.