#!/usr/bin/perl -w

use strict;

#Amb aquest programa contem en quantes especies esta conservat cada nucleotid dels 2000 que hi ha a la seq humana. Llavors ens quedarem amb els intervals de nucleotids que estan conservats en un minim de 4 especies. 

#Declarem un vector per a cada especie, menys la humana, ja que no hem fet aliniaments locals de la seq humana contra ella mateixa.

my $longitud = 2000; #Aquesta variable la utilitzem com a constant per a poder posar 2000 posicions a 0 en cada un dels vectors de les especies. 

my @bos;
my @danio;
my @tetraodon;
my @mus;
my @fugu;
my @canis;
my @rattus;
my @monodelphis;
my @macaca;
my @xenopus;
my @pan;
my @gallus;


my $m = 0; #Aquesta variable ens permet anar creant totes les posicions que ens             interessen en els vectors anteriors. Tenint en compte que la seq                humana te 2000 parells de bases, volem que tots els vectors tinguin             2000 posicions per a poder posar en cada una si en lespecie concreta            el nucleotid de la seq humana esta conservat. Un 1 significa que el             nucleotid esta conservat i un 0 que no ho esta. 

my $inici; #Es la variable on enregistrem la coordenada de homo sapiens on                 comença cada aliniament local.

my $final; #Es la variable on enregistrem la coordenada de homo sapiens on acaba           cada aliniament local.

my $codi; #Es la variable on enregistrem el codi de cada especie que aliniem               amb homo sapiens en cada un dels aliniaments locals.

my $i; #Serveix per recorrer totes les posicions de cada vector entre la               coordenada dinici i la de final per anar entrant els 1 en les posicions         que toquen.

my %vectors = ("ENSBTAG00000000031"=> \@bos,
	       "ENSDARG00000018242"=> \@danio,
	       "GSTENG00029525001"=>\@tetraodon,
	       "ENSMUSG00000014158"=>\@mus,
	       "NEWSINFRUG00000145973"=>\@fugu,
	       "ENSCAFG00000010796"=>\@canis,
	       "ENSRNOG00000001195"=>\@rattus,
	       "ENSMODG00000004471"=>\@monodelphis,
	       "ENSMMUG00000003942"=>\@macaca,
	       "ENSXETG00000018609"=>\@xenopus,
	       "ENSPTRG00000005430"=>\@pan,
	       "ENSGALG00000005152"=>\@gallus);

#Aquest hash ens permet accedir amb cada codi despecie al vector que li correspon.Com que en les linies dentrada hi ha la referencia de cada especie podem utilitzarho per recuperar el vector corresponent. 

my $j = 0; #Aquesta variable ens permet recorrer les 2000 posicions de cada                vector per veure si contenen un 1 o un 0. .

my @totsv = ("ENSBTAG00000000031","ENSDARG00000018242","GSTENG00029525001","ENSMUSG00000014158","NEWSINFRUG00000145973","ENSCAFG00000010796","ENSRNOG00000001195","ENSMODG00000004471","ENSMMUG00000003942","ENSXETG00000018609","ENSPTRG00000005430","ENSGALG00000005152"); 
#Aquest vector conte totes les especies amb que treballem. ens permet per a cada posicio dels vectors que contenen la conservacio de les 2000 posicions sumarhi el valor que correspon a cada especie. Ens cal un vector pq ens interessa poderho fer ordenadament i el hash no te ordre. 

my $st; #A aquesta variable li assignarem el valor de lscalar del vector @totsv.        Aixi podrem recorrer el vector assegurantnos que sumem per a cada               posicio totes les especies. 

my $f = 0; #Aquesta varible ens permet recorrerel vector @totsv. 

my @resultats; #Aquest es el vector que ens permet enregistrar els resultats.                  Posar en cada posicio del vector el nombre despecies que conserva               el nucleotid.

my $n; #Aquesta variable la utilitzem per posar el minim nombre despecies en que       ha destar conservat un nucleotid perque ho considerem significatiu.  

my $flag; #Aquesta variable es una etiqueta que tan sols pendra els vaors 0 i 1.          Tindran el valor 1 quan ens trobem a dins dun interval que esta                 conservat en el nombre despecies que ens ionteressa. Tindra el valor 0          si esta conservat en un nombre despoecies menor del que ens ionteress

my $score; #En aquesta variabel enregidtrem el minim despecies en que es troba             consevada cada regio. 

my @info; #En aquest vector hi posem la informacio de cada interval que ens               interessa:coordenada dinici, de final,i nombre minim despecies en que           esta conservat en format gff, per tal de poder fer posteriorment una            representacio grafica. 

my $h; #Ens serveix per a recorrer @info i posar en cada posicio la informacio         corresponent a una de les regions conservades. 

my $homo = ''; #Es la variavle on enregistrem tota la sequenscia humana per tal d              obtenir el format gff. 
my $homo2 = '';

my $gffFile = "filtrat.gff"; #Aquest es el niom del fitxer on quedara tota la                                informacio que surt del programa. 

my @hsapiens = ();

my @nucleotids = ();

my $c = 0;

open (HOMOS,"<homosapiens.fst");

while (<HOMOS>) {
	chomp;
	if($_=~ /^[actgACTG]+/){
	    $homo = $homo."##$_\n";
	    $homo2 = $homo2."$_";
	} 
   }	

close HOMOS;

while ($m<$longitud){

    $bos[$m] = 0; 
    $danio[$m] = 0;
    $tetraodon[$m] = 0;
    $mus[$m] = 0;
    $fugu[$m] = 0;
    $canis[$m] = 0;
    $rattus[$m] = 0;
    $monodelphis[$m] = 0;
    $macaca[$m] = 0;
    $xenopus[$m] = 0;
    $pan[$m] = 0;
    $gallus[$m] = 0;
    $resultats[$m] = 0;
    
    $m = $m + 1;
}  
#Amb aquest while creem i omplim de 0 les 2000 posicions de cada vector. Llavors posarem 1 a les posicions que estiguin conservades. 

open (TRPV4,"<TRPV4.gff");

while (<TRPV4>){
 
    if($_ =~ /^ENSG00000111199\s+Ens_[a-zA-Z]+\s+Ens_[a-zA-Z]+\s+([0-9]+)\s+([0-9]+)\s+[0-9]+\s+.\s+.\s+([a-zA-Z0-9]+)\s+[0-9]+\s+[0-9]+\s+[0-9]+\s+/){

#Amb aquest if reconeixem del fitxer dentrada totes les linies que segueixen aquest patro, que son les que contenen la informacio que ens interessa de cada aliniament local.El primer parentesi correspon al nucleotid dhomo sapiens en que comenca laliniament local. En el segom hi ha la coordenada en que acaba laliniament local i en el tercer hi ha el nom de lespecie de lalinimanet local. 

	$inici = $1;
	$final = $2;
	$codi = $3; #Assignem a $especie el codi de lespecie. 
	$i = $inici;

#Inicialitzem $i en la posicio del vector que correspon a la coordenada dinici de laliniament local. Posarem 1 fins que arribem a la posicio corresponent al final de laliniament local. 

	while($i<=$final){
	    $vectors{$codi}[$i] = 1;
	    $i = $i + 1;

	}
    }
}

#Ara en cada un dels vectors tenim anotat quins son els nucleotids conservats. Pero ens interessa saber en quantes especies esta conservat cada nucleotid. Per tant, sumem el valor de cada posicio en tots els vectors. 

$j = 0;
$st = scalar (@totsv);
 
while ($j < $longitud){
    $f = 0;
    while ($f < $st){
	$codi = $totsv[$f];
	$resultats[$j] = $resultats[$j] + $vectors{$codi}[$j];

	$f = $f + 1;

    }

    $j = $j + 1;
}
#Ara tenim el @resultats que te 2000 posicions i en cada un el nombre despecies que tenen consevat el nucletotid que ocupa aquella posicio. 

#El seguent pas del programa es filtrar quines son les regions mes conservades. Ens quedarem amb les regions que com a minim estan conservades en quatre de les dotze especies que comparem amb homo sapiens.

$n=4; #Volem que com a minim les regions de nucleotids es trobin conservades en       4 especies.De totes maneres, anirem incrementant $n per tenir les regions       que com a minim estan en 4, les que com a minim estan en 5, les que com a       minim estan en 6...fins a 8. 
$score = 0;
$h = 0;

while ($n<=8){
     
     $flag=0;
     $i = 0;
     while ($i<scalar(@resultats)){

	 if ($resultats[$i]>=$n){
	 
	     if ($flag == 0) {

		 $inici = $i;
		 $flag = 1; 
		 $score=$n;

	     #Si el nucleotid esta conservat en el minim despecies requerit i                $flag es 0, vol dir que es el primer nucleotid de linterval. Per                tant, guardem la seva posicio a $inici i passem $flag a 1 indicant              que ja ens trobem dins de linterval. 
	     }

#Si el nucleotid esta conservat en tantes espec ies com volem pero $flag val1, col dir que ja ens trobavem dins de linterval i, per tant, no ens ionteressa la posicio. 

	 }
	 else {
	     
	     if($flag == 1) {
		 $final = $i - 1;
		 $flag = 0;

		 @hsapiens = split(//,$homo2);
		 $c = $inici;

		 $nucleotids[$h] = "; ";
		 while ($c <= $final) {
		     $nucleotids[$h] = $nucleotids[$h].$hsapiens[$c];
		     $c = $c + 1;
		 }

#Si el nucleotid no esta conservat en el minim despecies i la flag val 1, vol dir que acabem de sortir de linterval. Per tant, guardem loa posicio anterior a la srtida de linterval, ja que ssera la posicio final d'aquest. 

		 $info[$h] ="ENSG00000111199\talig_".$n."_esp\talig_".$n."_esp\t$inici\t$final\t$n\t.\t.\t$nucleotids[$h]";	     
		 $h = $h + 1;
	     }
         
	 
        }

	$i = $i + 1;
     } # end while ($i<scalar(@resultats))
     
     $n=$n+1;
 }

#Ara volem tenir la informacio que hem guardat a @info en el format gff complet. Per tant, cal que entre altres coses imprimim la seq humana.Posem tota la informacio en un fitxer de sortida. 

open(OUT,">$gffFile");


    print OUT "##gff-version 2\n";
    print OUT "##programa2\n";
    print OUT "##date Fri March 10 14:48:20 2006\n";
    print OUT "##DNA ENSG00000111199\n";
    print OUT "$homo\n";

$h = 0;

while($h<scalar@info) {
    
    print OUT "$info[$h]\n"; #no imprimim tot el vector de cop pq ens surtia un                              espai a dalt de tot que no ens interessava.

$h++;

}

close TRPV4;
close OUT;
