#!/usr/local/bin/perl
# ===========================================================================
# CAPTURA D'SNPs humans i calcul de l'ancestralitat
# Clara Balleste
# Núria Conde
# ===========================================================================
use strict;
use warnings;
use LWP::Simple; #nom del modul

####################variables glabals#############################################################
my $utils = "http://www.ncbi.nlm.nih.gov/entrez/eutils";
my $query;#variable on guardarem el rs que ens entrin per pantalla
# = undef; seria interessant asegurar-nos una no definicio en cas que canvièssim el while infinit per el do while.
my $db = "snp";
my $tall; #= undef;variable que recollira el nombre de nucleotids que volem up and downstring.
my $defecte='all'; 
my @entrada; #vector on guardem els rs demanats un a un a cada posicio del vector 
my @claus; #vector on guardem les claus que tenim per el minihash; podriem incicialitzar-lo ara pero hem usat la funcio keys
my@hies; #vector que usa la subrutina minicrida per recorre la informacio rebuda des de l'ncbi
my$r=4;#sabem sempre que sera aquesta poscicio pel format.
my $j;
my %maxihash;
my$n=0;#contador per saber per quin rs dels demanats estem passant
my $seq1 = ""; #variables on hiha la seqüència upstring al snp de la cadena forward
my $seq2 = ""; #variable on hi ha la seqüència downstring al snp de la cadena reverse


#########control de les dades d'entrada el rs query i el nombre de nucleotids#########
while (1) { 
    print "Introdueix l'RS a buscar (màxim 90): ";
    $query = <STDIN>;
    if (($query =~ m/(rs\d+,)+(rs\d+)\n/)||($query =~ m/^rs\d+$/oi)){ #expressions regulars per saber si es un rs o una llista de molts
                                                                       # separats per comes i sense espai al final.
	chomp($query);
	@entrada=split(/,/,$query);
	if (scalar(@entrada)<90){ 
	    last;
	}
    };
    print "### ERROR:Cal un identificador valid (i.e.'rs123')o si es mes d'un separat per comes i sense espais. Mai mes de 90.\n\n";
};


while (1) { 
    print "Nombre de nucleotids al voltant: ";
    $tall  = <STDIN>;
    chomp($tall);
    if(($tall =~ m/^\d+$/oi)||($tall eq $defecte)){ #expressio regualar per assegurarnos que el que entra es un nombre o l'estrig all.
        last;
    };
    print "### ERROR: necessito un nombre sencer\n\n";
};



###########################################programa principal que corre per cada rs##########################################
for ($j=0;$j<(scalar(@entrada));$j=$j+1){
    if ((minicrida($entrada[$j],$tall))==1){ 
	$n=$n+1;
	$maxihash{$entrada[$j]}=crida($entrada[$j],$tall,\@claus);
	pagweb($tall,$seq1,$seq2,$n,\@claus,\%maxihash);
    }
    else{
	print"$entrada[$j] aquest nombre rs no consta a la base de dades\n";
    }
}


#***************************************************************************
#SUBRUTINES*****************************************************************
#*************************************************************************** 

########crida de sols el rs per controlar si crides a qui no esta penjat########

sub minicrida{
#http://www.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=rs455 #direccio que anem a buscar

    my$query=$_[0];
    my$tall=$_[1];
    my $resposta=0;
    my $esearch = "$utils/esearch.fcgi?" .
              "db=$db&term=";
    my $esearch_result = get($esearch . $query);
   # print "$esearch_result"; si la descomentem mostra el que rebem del ncbi.
    @hies=split(/>/,$esearch_result);
    if($hies[$r]eq "1</Count"){ # sabem que r=4 es la que volem pel format del esearch_result (comentari en variables globals)
	$resposta=1;
    }
    return $resposta;
}

#######################subrutina crida (molt llarga dividida conceptualment en diverses parts)##############
sub crida{ 
##################crida principal a SNPdb del ncbi usant el modul###############################
    my$query=$_[0];
    my$tall=$_[1];
    my $esearch = "$utils/esearch.fcgi?" .
	"db=$db&retmax=1&usehistory=y&term=";
    
    my $esearch_result = get($esearch .$query);
    
    $esearch_result =~
	m|<Count>(\d+)</Count>.*<QueryKey>(\d+)</QueryKey>.*<WebEnv>(\S+)</WebEnv>|s;
    
    #$Count    = $1;
    my$QueryKey = $2;
    my $WebEnv   = $3;
    
    
    my $efetch = "$utils/efetch.fcgi?" .
	"retmode=text&db=$db&query_key=$QueryKey&WebEnv=$WebEnv";

    my$efetch_result = get($efetch);
   #  print $efetch_result; si la descomentes pots veure el que retorna el ncbi al cridar el snp
    
###################fi de crida###################################################
###################variables locals de la funcio crida##########################
    my $i;
    my $concatenat=" ";#variable per guardar les dues cadenes unides
    my $concatenatara; #variable per comparar les dues cadenes unides
    my $seq5=' '; # guarda la seq upstring al snp de la cadena forward cal pk si no del $1 encreuabem les coses
    my $seq3=' '; #guarda la seq downstring al snp de la cadena forward cal pk si no del $1 encreuabem les coses
    my $rev5=' '; #guarda la seq upstring al snp de la cadena reverse cal pk si no del $1 encreuabem les coses
    my $rev3=' ';#guarda la seq downstring al snp de la cadena reverse cal pk si no del $1 encreuabem les coses
    my $ja=0;
    my $avere; 
    #my $seq1 = "";tretes pk interessa k siguin globals per usarles al printar a cada bucle
    #my $seq2 = "";
    my $rev1=""; #guardem la sequencia definitiva (la k es mes gran)
    my $rev2=""; #guardem la sequencia definitiva (la k es mes gran)
    my $rs = ""; #variable que sera el valor de la clau.
    my $esp = "";#variable que sera el valor de la clau.
    my $seq = "";#variable que sera el valor de la clau.
    my $snpf = "";#variable que sera el valor de la clau.
    my $snpr = "";#variable que sera el valor de la clau.
    my @linies; #vector on cada posicio conte la informacio que ens ha retornat el ncbi separat per comes
    my %minihash =("acessionumber" => "$rs", #noms posats estratègicament per a que quedin ordenats amb un sort
		   "especie" => "$esp",
		   "forwarda" => "$seq1",
		   "forwardsnp"=>"$snpf",
                   "forwardz" => "$seq2",
                   "reversea" => "$rev1",
		   "reversesnp"=>"$snpr",
		   "reversez" => "$rev2" );

    @linies = split(/,/,$efetch_result);
    

#####obrim el recorregut del array i omplim el nostre minihash###################
    for ($i=0;$i<(scalar(@linies));$i=$i+1){
	
	if(($linies[$i]=~(m/(rs\d+)/))==1){ #expresio regular per agafar el acession number
	    $rs=$1;
	    $minihash{acessionumber}=$rs; #guardem al minihash
	}
    
  
	if($linies[$i]=~ m/(\[.+\])/){#expresio regular per el agafar l'especie
	    $esp=$1;
	    $minihash{especie}=$esp;#guardem al minihash
	}
	
	if(($ja==0)&&($linies[$i]=~ m/(orient reverse)/)){ #expresio regular per saber si ja hem arribat a l'orient reverse
	    $avere=$1;#no caldria pero recollim l'estring de la posicio del vector que conte orient reverse
	    $ja=1; #si hem passat per orient reverse serà 1 (es un indicador de pas)
	    $concatenat=" "; #per que al comparar quines de les cadenes reverse son mes llargues comenci de zero i no amb el guardat de
                             # les cadenes forward
	}  
  
#per la cadena forward:     
	if(($ja!=1)&&($linies[$i] =~ m/seq5[\s\"]+([acgt\s]+)/i)){#expresio regular per agafar la sequencia upstring
	    $seq5=$1;
	    $seq5=~s/\n//g; #treiem els espais de la sequencia
	    if($linies[$i+2] =~ m/seq3[\s\"]+([acgt\s]+)/i){ #en ppi SEMPRE coincidira que 2 posicions mes aball hi hagi la seq.reverse
		$seq3=$1;
		$seq3=~s/\n//g;
                 $concatenatara = $seq5.$seq3;
		if((max($concatenatara,$concatenat))==1){ #agafem el grup seq 5 i seq3 que doni mes informacio pero sense desaparellar
		    $seq1=$seq5;
		    $seq2=$seq3;
		    $concatenat=$concatenatara;
		    $minihash{forwarda}=$seq1;#guardem al minihash
		    $minihash{forwardz}=$seq2;#guardem al minihash
		}
	    }
	}
    
#podriem haver posat en l'if anterior que ens guardes linies[i+1]que SEMPRE hagues
# quadrat tal com per seq 3; hem decidit però mantenir-ho en un if apart. En canvi per
# la cadena reversa l'hem integrat
 
########if separat per agafar el snp de la cadena forward###################
	if(($ja!=1)&&($linies[$i] =~ m/observed \"(.+)\"/)){ 
	    $snpf=$1;
	    $minihash{forwardsnp}=$snpf;#guardem al minihash
	}


 #####per la cadena reverse:(igual que la forward pero amb el ja==1);
   
   
	if(($ja==1)&& ($linies[$i] =~ m/seq5[\s\"]+([acgt\s]+)/i)){
	    $rev5=$1;
	    $rev5=~s/\n//g;
	    if($linies[$i+1] =~ m/observed \"(.+)\"/){
		$snpr=$1;
		$minihash{reversesnp}=$snpr;
	    }
	    if($linies[$i+2] =~ m/seq3[\s\"]+([acgt\s]+)/i){ #en ppi SEMPRE
		$rev3=$1;
		$rev3=~s/\n//g;
		$concatenatara = $rev5.$rev3;
		if((max($concatenatara,$concatenat))==1){
		    $rev1=$rev5;
		    $rev2=$rev3;
		    $concatenat=$concatenatara;
		    $minihash{reversea}=$rev1;#guardem al minihash
		    $minihash{reversez}=$rev2;#guardem al minihash
		}
	    }
	}
  
}#tanquem el recorregut de l'array linies.

    @claus= sort keys(%minihash);#guardem les claus del minihash ordenades per usar-lo la funcio pagweb
    return \%minihash; #retornem la referencia al hash

} #tanquem la subrutina de la crida

####################subrutina max#################################################
#retorna un 1 si hem de fer un canvi pk el nou es major i un zero quan no cal canvi.
    sub max {
	my $sino=0;
	my $nou;
	my $seq;
	$nou=$_[0];
	$seq=$_[1];
	$seq=~s/\n//g;
	$nou=~s/\n//g;
	if((length $nou) > (length $seq)){
	    $sino=1; 
	}
	return $sino;
    }

################################surutina mostratallat#######################################    
#decidir per decidir; hem decidit de guardar tota la seq mes gran i mostrar sols el k el usuari demani. si el nombre de nucleotids ( que no sobrepassi el max) es major del k dona ncbi llavors dona els k hi ha i treu una nota.

    sub mostratallat{
	my $vull;
	my $seq;
	my $posicio=0;
	my $longitut;
	my @base;
	my $i=(-1);
	my $quina;
	$vull=$_[0];
	$seq=$_[1];
	$quina=$_[2];
	$longitut=(length $seq);#aixi evitem usar 2 cops la funcio length
	    if ($vull >$longitut){
		print"L'annotacio actual del ncbi no disposa de tants nucleotids\n";
		return $seq;
	    }
	else{
	    if($quina==1){ #seq5 s'ha de començar a comptar pel final.
		#@seq=explode $seq;#pag85 no se si es millor opcio...
		#$base=substr($seq,$posicio,1); #pag93 del manual
		$i=$longitut-$vull; #no se per que pero no he pogut integrar la operacio a la condidio del for:$i
		for($posicio=($longitut-1); $posicio>($i-1); $posicio=$posicio-1){
		    $vull=$vull-1;
		    $base[$vull]=substr($seq,$posicio,1);
		}	    
		$seq=join('',@base); #pag93 ho unim tot de nou en un string.
		return $seq;
	    }
	    if($quina==2){ #la seq3 s'ha de contar des del principi
		for( $posicio=0; $posicio<($vull);++$posicio){
		    $base[$posicio]=substr($seq,$posicio,1);
		}
		$seq=join('',@base); #pag93
		return $seq;    
	    }
	}
    }
    

########################################subrutina pagweb per printar####################################
    sub pagweb{
	my $k=0;
	my $p=0;
	my $surtweb1;
	my $surtweb2;
	my $tall=$_[0];
	my $seq1=$_[1];
	my $seq2=$_[2];
	my $n=$_[3];
	my $claus=$_[4];#el vector es passat per referència
	my $maxihash=$_[5];#el hash es passat per referència
	$surtweb1=mostratallat($tall,$seq1,1);
	$surtweb2=mostratallat($tall,$seq2,2);
        ########printem el k anira a primer en la pagina web########
	print"\n$surtweb1  ";
	print"$maxihash{$entrada[$n]}{$claus[3]}";
	print"  $surtweb2\n";
	######printem el minihash k toca per n ###########
	print"\n\n$entrada[$n]\n";
	for ($k=0;$k<scalar(@claus);$k=$k+1){
	    print"$claus[$k]= $maxihash{$entrada[$n]}{$claus[$k]}\n";
	}
    }
   
##############################fi####################################################
