Heterosis blast.pl

#! /usr/bin/perl -w

# *************************************************************
# This Programme does the following in sequence
# -> Performs regular Blast
# -> Multi threads
# -> Parses the out and takes the top hit
# -> Modyfy the infile and outfile and carefully read the variables to change stuff
# Author : Rajkumar (itc@rajkumar.in)
# Release: APR 2010
# *************************************************************

# Libraries
	use strict;
	use Bio::Perl;
	use Bio::SearchIO;
	use Parallel::ForkManager;

# At the top: Time Start
	my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);

# Variables
	my $path    = "/home/raj/bio";
	my $db      = "$path/sorghum/db/sgdna";
	my $curdir  = "$path/user/rk/heterosis";
	my $query   = "$curdir/heterosis_maize.fa";
	my $tophit  = "$curdir/maize-sgdna_tophit.txt";	
	my $tempdir = "$curdir/temp";
	my $resdir  = "$curdir/result";


	my $cpu     = "1";
	my $pm = new Parallel::ForkManager(4); # 4 is the number of cores in your processor

	my (@files, @files_) = @_;
	my ($i)     = $_;
# Functions


	# Having the required dir
	if (-e "$tophit")    {`rm $tophit`;}
	if (-e "$tempdir")   {`rm -fr $tempdir`;}
	if (! -e "$tempdir") {`mkdir $tempdir`;}
	if (-e "$resdir")    {`rm -fr $resdir`;}
	if (! -e "$resdir")  {`mkdir $resdir`;}
	

# Get the MultiFasta

	my $in  = Bio::SeqIO->new(-file => "$query", -format => 'Fasta');
	my $identifier;
	my $sequence;
 
	while ( my $seq = $in->next_seq() ) {
		$identifier = $seq->id;
  		$sequence = $seq->seq;
 
		my $fasta= ">$identifier\n$sequence\n";
	
		open(TEMP, ">>".$tempdir."/".$identifier);
		print TEMP $fasta;
		close TEMP;

	}


# Loop for Blast with multithreading
 
	opendir (DIR, "$tempdir") or die $!;
	@files = readdir(DIR);
	close(DIR);
	
	$i = 0;
	foreach my $file(@files) {
		next if ($file eq "."  or $file eq "..");
		
		$i++;
		$pm->start and next;
		print "\n\tInitiated Blast on $file .. ";
		
		my $outfile = $resdir."/".$file.".blast.out";
              	my $infile  = $tempdir."/".$file;
              	
                `blastall -p blastn -d $db -i $infile -e 10 -m 8 -o $outfile -a $cpu`;
                
                print "  Done.";

		$pm->finish;
	}
	$pm->wait_all_children;
 

			# Blast parameters
			# -p program
			# -d database
			# -i query file
			# -e evalue
			# -m 8 tab delimited out
			# -o outfile
			# -a number of CPUs
			# -q  Penalty for a nucleotide mismatch (blastn only)   default = -3
			# -v  Number of database sequences to show one-line
			# -K  Number of best hits from a region to keep

	# cleaning the temp dir
	if (-e "$tempdir") {`rm -fr $tempdir`;}
		

# Parsing the Blast out to get Top Hit
	
 	opendir (DIR, "$resdir") or die $!;
	@files_ = readdir(DIR);
	close(DIR);

	foreach my $file_(@files_) {
		next if ($file_ eq "."  or $file_ eq "..");
 		$i = 0;
		foreach my $line_ (split ("\n", `cat $resdir/$file_`)) {
 
			# Filters only the top hit EST 
			
			$i++;
			next if ($i>1);

			open (F, ">>".$tophit);
			print F "$line_\n";
			close F;
		}
	}
	

	# Cleaning the Res dir
	if (-e "$resdir") {`rm -fr $resdir`;}
	
	print "\nAll Done!\n\n";

# At the End: Time Calculation
	my ($sec_,$min_,$hour_,$mday_,$mon_,$year_,$wday_,$yday_,$isdst_) = localtime(time);
	my $days = $mday_-$mday;
	my $hours = $hour_-$hour;
	my $mins = $min_-$min;
	my $secs = $sec_-$sec;
	print "\n------------------------\nDuration: ";
	if ($days>0) { print "$days days";}
	if ($hours>0) { print ": $hours hours";}
	if ($mins>0) { print ": $mins mins";}
	if ($secs>0) { print ": $secs secs";}
	print "\n------------------------\n";

# END