Parse tobea
#! /usr/bin/perl -w
# *************************************************************
# This Programme does the following in sequence
# -> Parse the Blast output of tobea - est contigs
# -> Write the results into a seperate folder with each EST
# -> Write the statistics files (est_scaf_hit.stat, unique_scaf_hit.stat)
# Author : Rajkumar (itc@rajkumar.in)
# Release: 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 $blastdir = "$path/tobacco/blast/tobea";
my $resdir = "$blastdir/result";
my $statfile = "$blastdir/tobea_est_hit.stat";
my $uniquescaf = "$blastdir/tobea_unique_est_hit.stat";
my $blastoutdir = "$blastdir/blastout";
my ($i, $line, $c_sequence) = $_;
my (@files,@esthits, @uniquescaf, @uniquescaf_ ) = @_;
# Parsing Functions
print "\n\tCreating required folders..\n";
if (-e "$resdir") {`rm -fr $resdir`;}
if (! -e "$resdir") {`mkdir $resdir`;}
# Work with the temp file to pick only the top hits with one scaffold
print "\tOpening the separate Blastout files..\n";
opendir (DIR, "$blastoutdir") or die $!;
@files = readdir(DIR);
close(DIR);
# Writing the headers
open (F, ">>".$resdir."/"."blast.parse.out");
print F "TOBEA_ID\tEST_ID\tTOBEA_SEQ_LENGTH\tALN_LENGTH\tALN_PERCENT\tIDNETITY\tMISMATCH\tGAP\tTOBEA_START\tTOBEA_END\tEST_START\tEST_END\tE-VALUE\n";
close F;
foreach my $file(@files) {
next if ($file eq "." or $file eq "..");
$i = 0;
foreach my $line (split ("\n", `cat $blastoutdir/$file`)) {
next if ($line =~ m/^#/);
my ($tobea_id, $est_id, $identity, $aln_length, $mismatches, $gap, $q_start, $q_end, $s_start, $s_end, $e, $bit_score ) = split(/\t/, $line);
$i++;
next if ($i >1 );
if ($file =~ m/[A-Z]+[0-9]+_([0-9]+)/){$c_sequence = $1;}
my $aln_percent= ($aln_length/$c_sequence)*100;
my $aln_percent = sprintf "%.1f",$aln_percent;
my $csv = "$tobea_id\t$est_id\t$c_sequence\t$aln_length\t$aln_percent\t$identity\t$mismatches\t$gap\t$q_start\t$q_end\t$s_start\t$s_end\t$e";
print "\r\tWriting $tobea_id in blast.parse.out";
open (F, ">>".$resdir."/"."blast.parse.out");
print F "$csv\n";
close F;
}
}
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