TGICL assembly

• Assembling ~50000 ESTs
  • Took 3.5 hours under 8 core
  • both run and parse scripts are below

#! /usr/bin/perl -w

# *************************************************************
# This Programme does the following in sequence
# -> Makes a single fasta
# -> Runs TGICL Program to assemble Ret files
# -> Contigs and Single
# Author: Rajkumar (r@rajkumar.in)
# Release: SEP 2010

# *************************************************************

# Libraries
	use strict;
	use Bio::Perl;
	use Bio::SeqIO;
	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 $curdir  = "$path/user/rk/castor/unigene";
	my $resdir  = "$curdir/result";
	my $data   = "$curdir/data/castorEST.txt";

	my $tgicl   = "$path/tgicl/tgicl";
	my $cpu     = "8";

	my ($i, $total, $file) = $_;
	my @files       = @_;

#Multi-threading
	my $midstop     = "0"; 		
	my $pm          = new Parallel::ForkManager(1);
	

# Running the TGICL
	if (-e "$resdir") {`rm -fr $resdir`;}
	if (! -e "$resdir") {`mkdir $resdir`;}


# Getiing into the folder


	my $in  = Bio::SeqIO->new(-file => "$data", -format => 'Fasta');
	my ($identifier, $sequence) = $_;
 	while ( my $seq = $in->next_seq() ) {
		$identifier = $seq->id;
  		$sequence = $seq->seq;
		$identifier =~ s/\|/_/g;

		my $fasta = ">".$identifier."\n".$sequence."\n";
			
		open (F, ">>".$resdir."/castorEST.fa");
		print F "$fasta";
		close F;

	}



	# TGILC

	chdir "$resdir";
	`$tgicl $resdir/castorEST.fa -c $cpu -l 40 -p 94 -M`;
			# -l miminum overlap length (default 40)
			# -M ignore lower-case masking in <fasta_db> sequences
			# -p minimum percent identity for overlaps <PID> (default 94)


# 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




#! /usr/bin/perl -w

# *************************************************************
# This Programme does the following in sequence
# -> Open a results folder of TGICL
# -> parse the contigs and singlets
# -> Output as multifasta of single file
# Author: Rajkumar <r@rajkumar.in>
# Release SEP 2010 | For Rajendrakumar
# *************************************************************

# Libraries
	use strict;
	use Bio::Perl;
	use Bio::SeqIO;
	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 $curdir    = "$path/user/rk/castor/unigene";
	my $resdir    = "$curdir/result";
	my $data      = "$resdir/castorEST.fa";
	my $pardir    = "$curdir/parsed";
	my $cdbyank = "$path/tgicl/bin/cdbyank";

	my ($i, $total, $file) = $_;
	my @files       = @_;

#Multi-threading
	my $midstop     = "0"; 		
	my $pm          = new Parallel::ForkManager(4);

# making the parsed folder

	if (-e "$pardir") {`rm -fr $pardir`;}
	if (! -e "$pardir") {`mkdir $pardir`;}


# getting into the results folder

	opendir (DIR, "$resdir") 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 "..");
		next if ($file !~ m/^asm/);		
		#print "$file";
		#next if (! -e "$resdir/$file/contigs");

		$i++;

		# Getting the contigs		
		my $in  = Bio::SeqIO->new(-file => "$resdir/$file/contigs", -format => 'Fasta');
		my ($identifier, $sequence) = $_;
	 	while ( my $seq = $in->next_seq() ) {
			$identifier = $seq->id;
	  		$sequence = $seq->seq;
			
			my $count_seq = length($sequence);
			my $fasta = ">".$identifier."_Size:".$count_seq."bp\n".$sequence."\n";
			
			open (F, ">>".$pardir."/contigs");
			print F "$fasta";
			close F;
		}


	}

	# Fetching singlets
	`$cdbyank $data.cidx <$data.singletons >$pardir/singlets.fa `;


# 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