Raj: Created page with 'Category: Bioinformatics #! /usr/bin/perl -w # ************************************************************* # This Programme does the following in seque…'

2010-04-15T13:40:00Z

Created page with 'Category: Bioinformatics <source lang=perl> #! /usr/bin/perl -w # ************************************************************* # This Programme does the following in seque…'

New page

[[Category: Bioinformatics]]

<source lang=perl>
#! /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

</source>

Heterosis blast.pl - Revision history

Raj: Created page with 'Category: Bioinformatics #! /usr/bin/perl -w # ************************************************************* # This Programme does the following in seque…'