GTH EST-gDNA

#! /usr/bin/perl -w
 
# *************************************************************
# This Programme does the following in sequence
# -> Gene prediction of Tobacco gDNA using all ESTs
# -> makes Multifasta; Runs GTH; converts the XML-to-GFF (gthxmlToGFF.py reuired in the cur.folder)
# -> Multicored
# 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 $datadir = "$path/tobacco/data";
	my $curdir  = "$path/tobacco/annotation/gt";
	my $resdir  = "$curdir/result";
	my $tempdir = "$curdir/temp";
	my $xmldir  = "$resdir/xml";
	my $gffdir  = "$resdir/gff";
	my $dna     = "$datadir/gdna.fa";
	my $est     = "$datadir/unigene.fa";
 
	my $gth     = "$path/gt/bin/gth"; 
	my $bssm    = "$path/gt/bin/bssm/arabidopsis";
 
	my ($i, $i1, $j, $midstop)   = $_;
	my (@files, @files1) = @_;
	my $pm = new Parallel::ForkManager(16);
	
	$midstop = 5000;
 
# cleaning the folders
	if (-e "$tempdir") {`rm -fr $tempdir`;}
	if (! -e "$tempdir") {`mkdir $tempdir`;}
	if (-e "$resdir") {`rm -fr $resdir`;}
	if (! -e "$resdir") {`mkdir $resdir`;}
	if (-e "$xmldir") {`rm -fr $xmldir`;}
	if (! -e "$xmldir") {`mkdir $xmldir`;}
	if (-e "$gffdir") {`rm -fr $gffdir`;}
	if (! -e "$gffdir") {`mkdir $gffdir`;}
 
# making Multifasta
 
	print "\n\tMaking Multi files of gDNA..";
	my $in  = Bio::SeqIO->new(-file => "$dna", -format => 'Fasta');
	my $identifier;
	my $sequence;
	
	$j=0; 
	while ( my $seq = $in->next_seq() ) {
		$identifier = $seq->id;
  		$sequence = $seq->seq;
 
		$j++;  if ($j>$midstop){last;}
		
		my $fasta = $_;
		$fasta = ">$identifier\n$sequence\n";
		open (FASTA, ">".$tempdir."/".$identifier);
		print FASTA $fasta;
		close FASTA;
	}
	print " Done.";
 
# Genome threader RUN
 
	print "\n\tInitiating Genome Threader\n";	
	opendir (DIR, "$tempdir") or die $!;
	@files = readdir(DIR);
	close(DIR);
 
	my $c_files = @files - 2;
 
	$i = 0;
	foreach my $file(@files) {
		next if ($file eq "."  or $file eq "..");
 
		$i++;  if ($i>$midstop){last;}
		my $percent = ($i/$c_files)*100;
		$percent = sprintf("%.2f", $percent);
		
		
		print "\r\t-> $i of $c_files | $percent%| Threading $file";
		$pm->start and next;
 
			`$gth -genomic $tempdir/$file -cdna $est -inverse -force -bssm $bssm -xmlout -o $xmldir/$file.xml`;
 
		$pm->finish;
 
 
	}
	$pm->wait_all_children;
 	print "\n\tDone.";
 
 
 
# Converting XML-to-GFF
 
	print "\n\tConverting XML-to-GFF\n";	
	opendir (DIR, "$xmldir") or die $!;
	@files1 = readdir(DIR);
	close(DIR);
 
 
 	my $c_files1 = @files1 - 2;
 	
	$i1 = 0;
	foreach my $file1(@files1) {
		next if ($file1 eq "."  or $file1 eq "..");
 
 
		my $math = `grep -c "no chain has been computed" $xmldir/$file1`;
		my $stat = `grep -c "general statistics" $xmldir/$file1`;
 
		chomp $math;
		chomp $stat;
		next if ($math ne 0);
		next if ($stat eq 0);
 
 
		my $newfile = $_;
		if ($file1 =~ m/(.+).xml/g){ $newfile = "$1.gff";}
 
 
		$i1++;  if ($i1>$midstop){last;}
		my $percent1 = ($i1/$c_files1)*100;
		$percent1 = sprintf("%.2f", $percent1);
		
		
		chdir $curdir;
		$pm->start and next;
		
		
		print "\r\t-> $i1 of $c_files1 | $percent1%| Converting $file1..";
			`./gthxmlToGFF.py -i $xmldir/$file1 -o $gffdir/$newfile`;

 
		#cleaning the empty files
		my $remove = `grep -c "GenomeThreaderEST" $gffdir/$newfile`;
		chomp $remove;
		if ($remove eq 0) {`rm $gffdir/$newfile`;}
 
		$pm->finish;
 
	
	}	
	$pm->wait_all_children;
 
	print "\n\tDone.";
 
 
# cleaning
 
	print "\n\tTrashing the temporary files and folders..";	
	if (-e "$tempdir") {`rm -fr $tempdir`; }
	if (-e "$xmldir") {`rm -fr $xmldir`; }
 
	print " Done.";
 
	print "\n\tAll Done MAN!\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