Parse tobea

From RAJ INFO
Jump to navigation Jump to search


#! /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