Blast tobea

From RAJ INFO
Revision as of 04:52, 2 March 2010 by Raj (talk | contribs) (Created page with 'Category: Bioinformatics <source lang=perl> #! /usr/bin/perl -w # ************************************************************* # This Programme does the following in seque…')
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search


#! /usr/bin/perl -w

# *************************************************************
# This Programme does the following in sequence
# -> Performs regular Blast
# -> Modyfy the infile and outfile
# Author : Rajkumar (itc@rajkumar.in)
# lease: MAR 2010
# *************************************************************

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


# 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/tobacco/data/blastdb/est_contig";
	my $query   = "$path/tobacco/data/tobea_unigene.fa";
	my $outfile = "$path/tobacco/blast/tobea/tobea_est.blast.out";
	
	my $blastdir      = "$path/tobacco/blast/tobea";
	my $multifastadir = "$blastdir/multifasta";
	my $blastoutdir   = "$blastdir/blastout";	
	
	my $cpu     = "4";
	my $i= $_;	

# Functions

	print "\n\tCreating required folders..\n";
	# Having the required dir
	if (-e "$multifastadir") {`rm -fr $multifastadir`;}
	if (! -e "$multifastadir") {`mkdir $multifastadir`;}
	if (-e "$blastoutdir") {`rm -fr $blastoutdir`;}
	if (! -e "$blastoutdir") {`mkdir $blastoutdir`;}
	
	
	my $total = `grep -c ">" $query`;
	chomp $total;

	my $in  = Bio::SeqIO->new(-file => "$query", -format => 'Fasta');
	my $identifier;
	my $sequence;
			
			
	$i = 0;		
	while ( my $seq = $in->next_seq() ) {
		$identifier = $seq->id;
  		$sequence   = $seq->seq;

		my $c_sequence = length($sequence);
		$identifier = $identifier."_".$c_sequence;
		
		my $fasta   = ">$identifier\n$sequence\n";
		
		
		$i++;
		my $math = $total - $i;
		

		open (FASTA, ">".$multifastadir."/".$identifier);
		print FASTA $fasta;
		close FASTA;
		
		print "\r\tInitiated Blast on $identifier  [$math to go..]";
		`blastall -p blastn -d $db -i $multifastadir/$identifier -e 10 -m 9 -q -3 -v 3 -o $blastoutdir/$identifier -a $cpu`;
		print "  done.";

			
			# -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 "$multifastadir") {`rm -fr $multifastadir`;}
	
	print "All 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