Parse misa.pl: Difference between revisions

From RAJ INFO
Jump to navigation Jump to search
m New page: <source lang=perl> #! /usr/bin/perl -w # ************************************************************* # This Programme does the following in sequence # -> Open the misa results # -> Par...
 
No edit summary
 
(2 intermediate revisions by the same user not shown)
Line 1: Line 1:
<source lang=perl>
<source lang="perl">
 
#! /usr/bin/perl -w
#! /usr/bin/perl -w


Line 83: Line 84:
}
}


print "\nDone All Man\n\n";
print "\n\tDone All Man\nOutfile written into $resfile\n\n";
#Removing Duplicates
 


# END
# END
</source>
</source>

Latest revision as of 11:12, 20 May 2016

#! /usr/bin/perl -w


# *************************************************************
# This Programme does the following in sequence
# -> Open the misa results
# -> Parse the files
# -> compare the original EST file and pick the sequence that have repeats 
# -> Outputs as fasta
# Author: Rajkumar (itc@rajkumar.in)

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

# Libraries
	use strict;
	use Bio::SeqIO;
	use Bio::Perl;


# Variables
	my $path    = "/home/raj/bio";
	my $ssrdir  = "$path/euca/ssr";
	my $datadir = "$ssrdir/data";
	my $resdir  = "$ssrdir/result";
	
	my $est     = "$datadir/globulus_est_trim.fa";
	my $repeat  = "$resdir/globulus_est_trim.fa.misa";
	my $resfile = "$datadir/globulus_est_repeat.fa";
	
	my ($i, $ssr_id, $total, $math ) = $_;
	my (@repeat_id, @identifier) = @_;

	
# Function
	
	# Lets Clean it
	if (-e "$resfile") {`rm $resfile`;}

	# start your hack

	open(CSV, "$repeat") or die "Idiot! Where is that Misa Outfile ? 'CSV': $!\n";
	<CSV>; # skip header line
	while (my $line = <CSV>) {
	
		my ($ssr_id) = split('\t', $line);
		next if ($ssr_id =~ m/^ID/);
		
		push (@repeat_id, "$ssr_id");
	}	

	# Making them Unique
	my %seen;
	my @unique_repeat_id = grep !$seen{$_}++, @repeat_id;
	
	# Getting some statistics 
	my $count_unique_repeat_id = @unique_repeat_id;
	$i = 0;
	
	# Lets Grab the sequence	
	foreach my $repeat_id(@unique_repeat_id) {
	
		$i++;
		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;
			
			next if ($repeat_id ne $identifier);
			
			if ($repeat_id eq $identifier) {
			
				my $fasta  = ">$identifier\n$sequence\n";
				open (FH, ">>".$resfile);
				print FH $fasta ;
				close FH ;
			print "\r\tWritten $identifier. Still $i of $count_unique_repeat_id to be written!";
			}
  		}
	}

	print "\n\tDone All Man\nOutfile written into $resfile\n\n";

# END