Parse blast.pl

#! /usr/bin/perl -w

# *************************************************************
# This Programme does the following in sequence
# -> Parse the Blast output of scaffold mapping
# -> 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: DEC 2009
# *************************************************************

# 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 $infile    = "$path/euca/blast/est_scaffold.txt";
	my $tempdir   = "$path/euca/blast/temp";
	my $resdir    = "$path/euca/blast/result";
	my $statfile  = "$path/euca/blast/est_scaf_hit.stat";
	my $uniquescaf= "$path/euca/blast/unique_scaf_hit.stat";
	
	my ($i)               = $_;
	my (@files,@esthits, @uniquescaf, @uniquescaf_ )  = @_;
	
# Parsing Functions

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

	
	# Let's Parse the file to have files written for indiavidual ESTs
	
	foreach my $line (split ("\n", `cat $infile`)) {
		
		next if ($line !~ m/^CL|gi/);
		next if ($line =~ m/^#/);
		
		
		foreach my $tab (split ("\t", $line)) {
			next if ($tab !~ m/^CL|gi|sca/);
			$tab =~ s/\|/_/g;
			
			if ($tab =~ m/^CL|gi/) {
				
				open (FH, ">>".$tempdir."/".$tab.".gff");
				print FH "$line\n";
				close FH;
			}
		}
        }


	
# Work woth the temp file to pick only the top hits with one scaffold
	opendir (DIR, "$tempdir") or die $!;
	@files = readdir(DIR);
	close(DIR);

	foreach my $file(@files) {
		next if ($file eq "."  or $file eq "..");
 
		# Clearing the array
		my @scaf_ids = ();
		
		foreach my $line_ (split ("\n", `cat $tempdir/$file`)) {
			
			my ($est_id, $scaf_id) = split(/\t/, $line_);
			push(@scaf_ids, "$scaf_id");

			# Filters only the top hit scaffold 
			next if ($scaf_ids[0] ne "$scaf_id");
			#cutoff if the alignment length is less than 70
			next if ($aln_length < 70);
			
			open (F, ">>".$resdir."/".$file);
			print F "$line_\n";
			close F;
		}
	}
	
	# cleaning the temp dir
	if (-e "$tempdir") {`rm -fr $tempdir`;}
	
	
	
# Generating statistics
	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`)) {
			
			$i++;
			my ($est_id, $scaf_id) = split(/\t/, $line__);
			next if ( $i>1 );
			
			push(@esthits, "$est_id\t$scaf_id\n");
			push(@uniquescaf, "$scaf_id\n");
			
		}
	}	
	
	
	# Writing the EST vs scaffolds
	@esthits = sort(@esthits);
	open (H, ">".$statfile);
	print H "@esthits";
	close H;
	
	#Writing the unique scaffold names
	my %seen;
	@uniquescaf_ = grep !$seen{$_}++, @uniquescaf;
	@uniquescaf_ = sort(@uniquescaf_);
	
	open (H1, ">".$uniquescaf);
	print H1 "@uniquescaf_";
	close H1;


	print "\nDone!\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