Co-ordinates blast.pl

From RAJ INFO
Revision as of 17:43, 29 December 2009 by Raj (talk | contribs) (New page: <source lang=perl> #! /usr/bin/perl -w # ************************************************************* # This Programme does the following in sequence # -> Parse the parsed Blast output o...)
(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
# -> Parse the parsed Blast output of scaffold mapping
# -> Get the hit co-ordinates for each scaffold 
# -> The co-ordinate window of EST-scaffolds map that have only 3kb in length
# -> Write the results as scaffold and co-ordinate window
# -> Pick out the selected scaffold sequences only with that co-ordinate window
# -> Some scaffolds have distinct hits at more than 1 places and will written as separate file. They need to cured manyually.
# 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 $datadir   = "$path/euca/blast/result";
	my $coordir   = "$path/euca/blast/coor";
	my $blastdir  = "$path/euca/blast";
	my $scaf_coor = "$blastdir/scaf_coor.stat";
	my $scaf_coor1 = "$blastdir/scaf_big_coor.stat";
	
	my ($i)               = $_;
	my (@files, @files1)  = @_;
	
# Parsing Functions

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

	
	
	opendir (DIR, "$datadir") or die $!;
	@files = readdir(DIR);
	close(DIR);

	foreach my $file(@files) {
		next if ($file eq "."  or $file eq "..");
 
		# Clearing the array
		my (@coor, @coordata)  = ();
		my $scaf_id            = ();
		my $est_id             = ();
		
		foreach my $line (split ("\n", `cat $datadir/$file`)) {
			
			my ($junk1, $junk2, $junk3, $aln_length, $junk5, $junk6, $junk7, $junk8, $start, $end) = split(/\t/, $line);
			push(@coor, "$start", "$end");
			$scaf_id = "$junk2";
			$est_id  = "$junk1";
		}
		
		
		# manipurate the array (arrange)
		@coor = sort({$a <=> $b}@coor);
		my $math = $coor[-1] - $coor[0];
		next if ($math >5000);
		next if ($math <150);
		
		my $tempfile = "$scaf_id\t$coor[0]\t$coor[-1]\n";
		open (FH, ">>".$coordir."/".$scaf_id);
		print FH "$tempfile";
		close FH;
	}
	
# Further filtering the co-ordinates because there are many same scaffolds
	opendir (DIR, "$coordir") or die $!;
	@files1 = readdir(DIR);
	close(DIR);
		
	foreach my $file1(@files1) {
		next if ($file1 eq "."  or $file1 eq "..");
	
		# Clearing the array
		my (@coor1, @coordata1)  = ();
		my $scaf_id1             = ();
		
		
		foreach my $line1 (split ("\n", `cat $coordir/$file1`)) {
			
			my ($junk, $start1, $end1) = split(/\t/, $line1);
			push(@coor1, "$start1", "$end1");
			$scaf_id1 = "$junk";

		}
	
		# manipurate the array (arrange)
		@coor1 = sort({$a <=> $b}@coor1);
		my $math1 = $coor1[-1] - $coor1[0];
		
		if ($math1 <5000) {

			my $tempfile1 = "$scaf_id1\t$coor1[0]\t$coor1[-1]\n";
			open (FH, ">>".$scaf_coor);
			print FH "$tempfile1";
			close FH;
		}

		else {

			my $tempfile2 = "$scaf_id1\t@coor1\n";
			open (FH, ">>".$scaf_coor1);
			print FH "$tempfile2";
			close FH;
		}
	}
	
	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