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
# -> 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/est.fa";
 
	my $gth     = "$path/gt/bin/gth"; 
	my $bssm    = "$path/gt/bin/bssm/arabidopsis";
 
	my ($i, $total)   = $_;
	my (@files, @files1) = @_;
	my $pm = new Parallel::ForkManager(16);

# cleaning the folders
	if (  -e "$tempdir") {`rm $tempdir`;}
	if (! -e "$tempdir") {`mkdir $tempdir`;}
	if (  -e "$resdir") {`rm $resdir`;}
	if (! -e "$resdir") {`mkdir $resdir`;}
	if (  -e "$xmldir") {`rm $xmldir`;}
	if (! -e "$xmldir") {`mkdir $xmldir`;}
	if (  -e "$gffdir") {`rm $gffdir`;}
	if (! -e "$gffdir") {`mkdir $gffdir`;}

# making Multifasta

	print "\n\tMaking Multi files of ESTs..";
	my $in  = Bio::SeqIO->new(-file => "$est", -format => 'Fasta');
	my $identifier;
	my $sequence;

	while ( my $seq = $in->next_seq() ) {
		$identifier = $seq->id;
  		$sequence = $seq->seq;
		
		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..";	
	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++;
	
		print "\n\tThreading $file - ($i of $c_files)";
		$pm->start and next;

			`$gth -genomic $dna -cdna $tempdir/$file -force -bssm $bssm -xmlout -o $xmldir/$file.xml`;
		
		print "\r..OK.";
		$pm->finish;
		

	}
	$pm->wait_all_children;
	
	


# Converting XML-to-GFF
	
	print "\n\tConverting XML-to-GFF..";	
	opendir (DIR, "$xmldir") or die $!;
	@files1 = readdir(DIR);
	close(DIR);

	foreach my $file1(@files1) {
		next if ($file1 eq "."  or $file1 eq "..");
		
		
		my $math = `grep -c "no chain has been computed" $xmldir/$file`;
		my $stat = `grep -c "general statistics" $xmldir/$file`;
		
		chomp $math;
		chomp $stat;
		next if ($math ne 0);
		next if ($stat eq 0);
		
		
		my $newfile = $_;
		if ($file =~ m/(.+).xml/g){ $newfile = "$1.gff";}

		chdir $curdir;
		$pm->start and next;
			`./gthxmlToGFF.py -i $xmldir/$file -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 " Done.";
	

# 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