#!/usr/bin/env perl

# Autor:	Katrin Grothaus
# Sequenzvergleich zwischen humanen und murinen CGIs mit dem Programm blast2seq

# Eingabe: 	CGI-Datei humane, Nachbarschaftsdatei aus Kapitel 2.2.3
#			CGI-Datei murin, Nachbarschaftsdatei aus Kapitel 2.2.3
#			Sequenz humane CGIs
#			Sequenz murine CGIs

# Aufruf:	perl Sequenzvergleich_human_murin.pl <CGI human> <CGI murin> <Seq human> <Seq murin>

# Hinweis: blast2seq	(http://blast.ncbi.nlm.nih.gov/Blast.cgi) muss lokal auf dem PC installiert sein

use warnings;
use strict;
use Bio::Tools::Run::StandAloneBlast;
use Bio::Seq;
use Bio::SeqIO;
use Bio::Perl;
use Bio::SearchIO;
use List::Util qw[min max];

# Variablen für die Eingabedateien
my $infile_hg19=$ARGV[0];
my $infile_mm10=$ARGV[1];
my $infile_hg19_seq=$ARGV[2];
my $infile_mm10_seq=$ARGV[3];

#Variablen für die Ausgabedateien
my $outfile_same_genes="cgis_same_genes_".$infile_hg19;
my $outfile_same_genes_same_seq="cgis_same_genes_same_seq_".$infile_hg19;

# benötigte Variablen für die unterschiedlichen Berechnungen und Parameter
my @hg19;
my @hg19_seq;
my @hg19_status;
my @mm10;
my @mm10_seq;
my @mm10_status;

my $counter_hg19=0;
my $counter_mm10=0;

my @cgi_hg_parts;
my @cgi_mm_parts;
my @cgi_hg_parts_seq;
my @cgi_mm_parts_seq;


my $cgi_mm_id;
my $cgi_mm_distpre;
my $cgi_mm_genepre;
my $cgi_mm_genepre_up;
my $cgi_mm_genein;
my $cgi_mm_genenext;
my $cgi_mm_genenext_up;
my $cgi_mm_distnext;
my $cgi_mm_org;

my $cgi_mm_nam_seq;
my $cgi_mm_seq_seq;
my $cgi_mm_source_seq;

my $cgi_hg_id;
my $cgi_hg_distpre;
my $cgi_hg_genepre;
my $cgi_hg_genein;
my $cgi_hg_genenext;
my $cgi_hg_distnext;
my $cgi_hg_org;

my $cgi_hg_nam_seq;
my $cgi_hg_seq_seq;
my $cgi_hg_source_seq;

my $cgi_hg_mm_gene=0;

# Eingabedateien jeweils in eine Array einlesen
open (INHG, $infile_hg19)or die "Can't open file: $!\n";
	@hg19= <INHG>;
close INHG;

open (INHGS, $infile_hg19_seq)or die "Can't open file: $!\n";
	@hg19_seq= <INHGS>;
close INHGS;

open (INMMS, $infile_mm10)or die "Can't open file: $!\n";
	@mm10= <INMMS>;
close INMMS;

open (INMM, $infile_mm10_seq)or die "Can't open file: $!\n";
	@mm10_seq= <INMM>;
close INMM;

# Öffnen der Ausgabedateien zum Schreiben
open(OUTS1,">".$outfile_same_genes);
open(OUTS2,">".$outfile_same_genes_same_seq);

#Start der Vergleiche, zeilenweises Durchgehen der <CGI human>
for(my $i=0;$i<@hg19;$i++){
	#Zeile wird aufgesplittet in die einzelnen Parameter
	my $cgi_hg =$hg19[$i];
	chomp($cgi_hg);
	@cgi_hg_parts=split("\t", $cgi_hg);
	$cgi_hg_id=$cgi_hg_parts[0];
	$cgi_hg_distpre=$cgi_hg_parts[1];
	$cgi_hg_genepre=$cgi_hg_parts[2];
	$cgi_hg_genein=$cgi_hg_parts[3];
	$cgi_hg_genenext=$cgi_hg_parts[4];
	$cgi_hg_distnext=$cgi_hg_parts[5];
	$cgi_hg_org=$cgi_hg_parts[6];
	# die dazugehörige Zeile der Sequenzdatei wird ebenfalls aufgesplittet in die einzelnen Parameter
	my $cgi_hg_seq=$hg19_seq[$i];
	chomp($cgi_hg_seq);
	@cgi_hg_parts_seq=split("\t", $cgi_hg_seq);
	$cgi_hg_nam_seq=$cgi_hg_parts_seq[0];
	$cgi_hg_seq_seq=$cgi_hg_parts_seq[1];
	$cgi_hg_source_seq=$cgi_hg_parts_seq[2];
	my $cgi_hg_seq_length=length($cgi_hg_seq_seq);
		
		# Murine Daten werden für den Vergleich auch aufgesplittet 
		for(my $j=0;$j<@mm10;$j++){
			my $cgi_mm =$mm10[$j];
			chomp($cgi_mm);
			@cgi_mm_parts=split("\t", $cgi_mm);
			$cgi_mm_id=$cgi_mm_parts[0];
			$cgi_mm_distpre=$cgi_mm_parts[1];
			$cgi_mm_genepre=$cgi_mm_parts[2];
			$cgi_mm_genein=$cgi_mm_parts[3];
			$cgi_mm_genenext=$cgi_mm_parts[4];
			$cgi_mm_distnext=$cgi_mm_parts[5];
			$cgi_mm_org=$cgi_mm_parts[6];
			if($cgi_mm_genepre ne "none"){
				$cgi_mm_genepre_up=uc($cgi_mm_genepre);
			}
			if($cgi_mm_genenext ne "none"){
				$cgi_mm_genenext_up=uc($cgi_mm_genenext);
			}
			
			my $cgi_mm_seq=$mm10_seq[$j];
			chomp($cgi_mm_seq);
			@cgi_mm_parts_seq=split("\t", $cgi_mm_seq);
			$cgi_mm_nam_seq=$cgi_mm_parts_seq[0];
			$cgi_mm_seq_seq=$cgi_mm_parts_seq[1];
			$cgi_mm_source_seq=$cgi_mm_parts_seq[2];
			my $cgi_mm_seq_length=length($cgi_mm_seq_seq);
				@cgi_mm_parts_seq=split("\t", $cgi_mm_seq);
				$cgi_mm_nam_seq=$cgi_mm_parts_seq[0];
				$cgi_mm_seq_seq=$cgi_mm_parts_seq[1];
				$cgi_mm_source_seq=$cgi_mm_parts_seq[2];
				
				# Überprüfen, ob das CGI die selben Nachbargene hat
				if (($cgi_mm_genepre_up eq $cgi_hg_genepre)&&($cgi_mm_genenext_up eq $cgi_hg_genenext)){
					# Ausgabe
					print OUTS1 $cgi_hg."\t".$cgi_mm."\n";

					# Vergleich der Sequenzen mit Hilfe von blast2seq
						my $input1 = Bio::Seq->new( -seq => $cgi_hg_seq_seq, -id  => $cgi_hg_nam_seq );
						my $input2 = Bio::Seq->new( -seq => $cgi_mm_seq_seq, -id  => $cgi_mm_nam_seq );
						
						my $seq1_length=length($cgi_hg_seq_seq);
						my $seq2_length=length($cgi_mm_seq_seq);
						
						my $min_seq= min($seq1_length,$seq2_length);
						
						# Objekt für blast-Aufruf generieren
						my $factory = Bio::Tools::Run::StandAloneBlast->new('-p' =>'blastn', '-G' =>'-5', '-E' =>'-2', '-W' => '11', '-r' => '2' );

						# Aufruf von blast
						my $blast_report=$factory->bl2seq($input1, $input2);
						
						my $seq1=$input1->seq();
						my $seq2=$input2->seq();
						
						# Ergebnisse von blast durchgehen, e value < 0.05 
						while (my $result = $blast_report->next_result) {
							if ($result->num_hits >0){
								while (my $hit = $result->next_hit) {
									my $hsp = $hit->next_hsp; 
										if ($hsp->evalue < 0.05){
											my $hsp_length=$hsp->length('total');
											my $identity=($hsp_length/$min_seq)*100;
											my $identity_rounded=sprintf("%.2f", $identity);
											if ($identity_rounded > 50){
												# Ausgabe
												print OUTS2 $cgi_hg."\t".$cgi_mm."\t SEQ \n";
												$cgi_hg_mm_gene=1;
												$mm10_status[$j]=1;
												$hg19_status[$i]=1;
											}
										}
								}
							}
						}
				}
		}
}
# Augabedateien schließen
close(OUTS1);
close(OUTS2);
