#!/usr/local/bin/perl -w
use strict;
use lib '/software/pathogen/external/apps/usr/bin/perl';
use Bio::SeqIO;
use Getopt::Std;
use File::Basename;
use Data::Dumper;

#Author: np1@sanger.ac.uk

sub usage {

die <<EOF;

This script takes embl files, runs orthomcl and produces new .orthomcl embl files containing 
ortholog and paralog information.  The orthomclxxxxxx_clusters.txt is a tab delimited file
containing the breakdown of clusters along with product information intended for viewing in 
a spreadsheet.

Usage:
$0: [-p <pvalue cutoff>] [-i <percentage ID cutoff>] [-m <mcl inflation>] [-s produce sybil format output]  <embl files>

  Defaults:  -p   1e-5
             -i   off
             -m   1.5
	     -s   off

EOF


}


=pod
NB/
If using this script with a new installation of orthomcl it will require a few modifications.
orthomcl was altered:
*to write results to a directory in the CWD named orthomcl${$}
*to pass the pid back to this script by printing the line "TEST <pid>"
*to not look for data in the sample_data directory
*to use bsub to run the all vs all blast.

=cut

my %options;
getopts( 'si:p:m:o:', \%options );

# -o = a debugging option - provide the mcl file as the arg

my $pvalue = "";
my $id = "";
my $inflation = "";
#store defaults for output
my $p = 0;
my $i = 0;
my $m = 1.5;
my $sybil = 0;
my $mcl;
my $pid;	   
my %count_by_organism;   
if(defined $options{p}){
    $pvalue = "--pv_cutoff=$options{p}";
    $p = $options{p};
}
if(defined $options{i}){
    $id = "--pi_cutoff=$options{i}";
    $i = $options{i};
}
if(defined $options{m}){
    $inflation = "--inflation=$options{m}";
    $m = $options{m};
}

defined $options{o} and $mcl = $options{o};


exists $options{s} and $sybil = 1;


unless (scalar(@ARGV) > 1){
    warn "Need two or more files for comparison";
    usage;
}

foreach (@ARGV){

    unless (-e $_) {
	warn "You must specify the embl file locations: $_ does not exist\n";
	usage;
    }
}


my @tmp_files;
my %products;


print STDERR "Producing protein files...\n";

#write files to protein fasta 
	foreach my $file (@ARGV){
    
   	 open (FASTATMP, ">$file.tmp") or die "Could not write to $file.tmp: $!";
   	 push (@tmp_files, "$file.tmp");
   	 my $stream = Bio::SeqIO -> new ( -file => $file, -format => 'embl');

	 while (my $seq = $stream->next_seq()){
	
		#include psudo genes but add pseudo to name
		foreach my $feature ($seq-> get_SeqFeatures()){
		
			#print STDERR $feature->location->to_FTstring . "\n";
	    		if( $feature -> primary_tag() eq 'CDS'){
			
			#print STDERR $feature->location->to_FTstring . "\n";
			
			
				my $id;
	            	
				if ($feature ->has_tag('systematic_id')){
                
				    ($id) = $feature -> get_tag_values('systematic_id');
				}
				elsif ($feature ->has_tag('temporary_systematic_id')){
                
				    ($id) = $feature -> get_tag_values('temporary_systematic_id');
				}
				elsif ($feature ->has_tag('locus_tag')){
                		
				    ($id) = $feature -> get_tag_values('locus_tag');
				}
				else {
				    warn"Couldn't find an id for the CDS\n";
				    next;
				}
		
				if ($feature ->has_tag('pseudo')){

		
				    $id = 'pseudo'. $id;
				}
		
				
				if (length ($feature-> spliced_seq -> translate -> seq) > 10){
					if ($feature -> has_tag('product')){
					    my @products = $feature->get_tag_values ('product');
					    my $products = join (',',  @products);
					    $products{$id} = $products;
					}
	
					my $seq = $feature -> spliced_seq -> translate -> seq;
			
					$seq =~ s/\*//g;
		
					print FASTATMP ">$id\n$seq\n";
					$count_by_organism{$file}++;
				}
				else {
				
					print STDERR "Excluding $id.... too short\n";
				}
			
    
	    		}

		}	
    	}

    close FASTATMP;

}

unless ($mcl){

print STDERR "Running orthomcl...\n";
#perform orthomcl


	my $files = join (',', @tmp_files);
	print "/software/pathogen/external/applications/OrthoMCL_v1.4/orthomcl.pl --mode=1 --fa_files=$files $pvalue $id $inflation\n";
	open (ORTH, "/software/pathogen/external/applications/OrthoMCL_v1.4/orthomcl.pl --mode=1 --fa_files=$files $pvalue $id $inflation |");
	
	while (<ORTH>){
	    
	if (/TEST (\d+)/){
	
		$pid = $1;
	}	
}

}
else {
	print STDERR "MCL provided...skipping orthorun\n";
	$pid = ${$};
}


$mcl or $mcl = "orthomcl$pid/all_orthomcl.out";

open (MCL, $mcl) or die "could not open MCL file $mcl:$!";

#read in the orthomcl file


-e "orthomcl$pid" or mkdir("orthomcl$pid") or die "couldn't make dir :$!" ;

open (OUT, ">orthomcl$pid/clusters.txt") or die "Could not open clusters file orthomcl$pid/clusters.txt for write: $!";

if ($sybil){
	open (SYBIL, ">orthomcl$pid/sybil.out") or die "Could not open sybil output file for write: $!";
}

my %results;
my @clusters;
my %sources;

print STDERR "Parsing results file...\n";
while  (my $line = <MCL>){
    
    my %cluster;
    
    my ($cluster, $members ) = split(/:/, $line);

    my ($name) = $cluster =~ /(\S+?)\(/;
    
    $members =~ s/^\s+//;

    my @members = split( /\s+/, $members);

    foreach my $member (@members){
		
	my ($id, $source) = $member =~ /(\S+?)\((\S+?)\)/;

	#$sybil and print SYBIL $source . ":" . $id . "\t";
	
	unless (defined $products{$id}){
		$products{$id} = 'unknown';
	}
	
	#add to cluster hash with product info 
	push (@{$cluster{$source}}, $id);
	
	$sources{$source} = 1;

	foreach my $compare (@members){
	    
	    my ($id2, $source2) = $compare =~ /(\S+?)\((\S+?)\)/;

	    #print "$id $source - $id2 $source2\n";

	    next if $id eq $id2 && $source eq $source2;

	    my $xlog = 'paralog';
	   
	    $xlog = 'ortholog' if $source ne $source2;
	    
	    push (@{$results{$id}}, ([$xlog, $source2, $id2]));
	}
    }
    
    #$sybil and print SYBIL "\n";
    
    #push cluster into clusters array
    push (@clusters, \%cluster);
}

#print STDERR Dumper(\@clusters) . "\n";

foreach my $file (@ARGV){
    
    my $stream = Bio::SeqIO -> new ( -file => $file, -format => 'embl');
    my $base = basename ($file);
    my $out_stream = Bio::SeqIO -> new ( -file => ">$base.orthomcl" , -format => 'embl');
    while (my $seq = $stream->next_seq()){
	
	#include psudo genes but add pseudo to name
	foreach my $feature ($seq-> get_SeqFeatures()){
	    if( $feature -> primary_tag() eq 'CDS'){
		
		my $id;
            
		if ($feature ->has_tag('systematic_id')){
                
			($id) = $feature -> get_tag_values('systematic_id');
		}
		elsif ($feature ->has_tag('temporary_systematic_id')){
               		($id) = $feature -> get_tag_values('temporary_systematic_id');
		}
		elsif ($feature ->has_tag('locus_tag')){
                		
			($id) = $feature -> get_tag_values('locus_tag');
		}
		else {
		    #warn"Couldn't find an id for the CDS\n";
		    next;
		}
		
		if (defined $results{$id}){
		    
		    foreach (@{$results{$id}}){
			
			my ($xlog, $xsource, $xid)  = @{$_};
		       $feature -> add_tag_value( $xlog, "$xsource:$xid;predicted_by_orthomcl" );
		    }
		}
		
		
	    }
	}#features
	$out_stream -> write_seq($seq);
	
    }#seq
}#file

my $count = 0;    





print OUT "Parameters:\npvalue $p\n%id $i\nmcl inflation $m\n\nProcessed:\n";

foreach (sort (keys %count_by_organism)){
	print OUT "$_ : $count_by_organism{$_}\n";
}

print OUT "\n\t";
foreach (sort (keys %sources)){

    print OUT $_ . "\t";
}   

print OUT "\n";
 
foreach my $hash (@clusters){
    
    print OUT "Cluster$count \t";
    
    my %hash = %{$hash};
    
    my $represented_sources = 0;
    my $printed = 0;
    
    
    #print out number in clusters
    foreach (sort (keys %sources)){
	
	if(defined $hash{$_}){
	    print OUT scalar(@{$hash{$_}}) . "\t";
	    $represented_sources++;
	}
	else{
	    print OUT "0\t";
	}
    }
    
    foreach my $source (sort (keys %sources)){
	
	if(defined $hash{$source}){
	
	    if ($sybil && $represented_sources > 1){
	        $printed and print SYBIL "\t";
	    	print SYBIL join ("\t", map ( "$source:$_" , @{$hash{$source}} ) );
		$printed = 1;
	    }
	
	    my $line = join (',',  map ( "$_($products{$_})" ,@{$hash{$source}}));
	    print OUT "$line\t";
	}
	else{
	    print "\t";
	}
	
    }
    $count++;
    print OUT "\n"; 
    $sybil && $printed and print SYBIL "\n";

}

print STDERR "Writing Venn file...\n";
#create file for venn diagram
system ("/nfs/team81/np1/workspace/my_scripts/ortho_to_venn.pl orthomcl$pid/clusters.txt > orthomcl$pid/clusters.venn") and die "Couldn't run Venn script";



