#!/usr/bin/perl -w use strict; ########## ### ### ########## ### ### ### ######### ########## ### ### ########## #### ### ### ########## ### ### ### ### ### ##### ### ### ### ### ### ### ### ### ### ###### ### ### ### ### ###### #### ### ### ### ### ### ### ### ### ###### #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###### ### ### ### ### ### ### ### ### ### ##### ### ### ### ########## ### ### ########## ### #### ### ########## ########## ### ### ########## ### ### ### ######### #Pedro Melgar Rojas #Alex Brenchat Barbera ############ ############ ############ ENTRADA DE DADES ############ ############ ############ #si el nombre d'arguments es menor que 8 el programa ens dirà els arguments que em d'introduir. if (scalar(@ARGV) < 8){ print "\n"; print "escribid: ./exonid.pl fitxero cutoff_acceptors cutoff_donor longitud_minima iniciales internos terminales single-genes\n\n"; print "fitxero: nombre del fichero con la secuencia en formato FASTA\n"; print "cutoff_acceptors: 0-1 (valor del 0 al 1)\n"; print "cutoff_donors: 0-1\n"; print "longitud_minima: tamaño mínimo de los exones predichos\n"; print "iniciales: 0 o 1 (1 para mostra ese tipo de exones y 0 para no mostrarlo)\n"; print "internos: 0 o 1\n"; print "terminaless: 0 o 1\n"; print "single-genes: 0 o 1\n"; exit (1); } #assigno els noms del vector d'arguments (parametres) amb noms significatius a unes variables. my $fitxerseq = $ARGV[0]; my $cutoffacc = $ARGV[1]; my $cutoffdonn = $ARGV[2]; my $long = $ARGV[3]; my $inicials = $ARGV[4]; my $internals = $ARGV[5]; my $terminals = $ARGV[6]; my $single = $ARGV[7]; # llegim les proporcions de codons del fitxer i les enregistrem # al hash %pcodons my %pcodons; open (PROPCODONS,") { # a la variable anonima tenim la linia llegida i l'encaixem amb # la seguent expressio regular m/(\w+) (\d+\.\d+)/; # l'expressio regular anterior encaixa amb dues paraules la primera # de les quals es alfabetica i la segona un numero amb decimals. # les dues paraules estan agrupades amb parentesis lo qual significa # que ens podem referir a la primera com $1 i a la segona com $2 $pcodons{$1}=$2; } close(PROPCODONS); # llegim una sequencia de DNA en format FASTA i la fiquem en un vector open (SEQUENCIA,"<$fitxerseq"); my $iden = ; # llegim primera linia del format FASTA ">XXXX" my @seql = ; # llegim les linies de la sequencia my $seq0 = ""; # declarem una variable on anira a para tota la # sequencia seguida, i la inicialitzem amb la # paraula buida my $i = 0; while ($i < scalar(@seql)) { chomp($seql[$i]); $seq0 = $seq0 . $seql[$i]; $i = $i + 1; } #convertim les matrius de pesos en hashes my %matriu_donors = &hash ("matdon.txt"); my %matriu_acceptors = &hash ("matacc.txt"); $seq0 = uc($seq0); # passem totes les bases a majuscules. Aquesta seqüència correspon a la pauta 0 #treiem 1 o 2 nucleòtids de la seqüència original per obtenir els diferents tipus d'ORFs my $seq1=substr ($seq0,1); #pauta 1 my $seq2=substr ($seq0,2); #pauta 2 # dividim la seqüència en codons i la fiquem en diferents posicions d'un vector my @vseq0 =($seq0 =~ m/.../g); #pauta 0 my @vseq1 =($seq1 =~ m/.../g); #pauta 1 my @vseq2 =($seq2 =~ m/.../g); #pauta 2 #Calculem les diferents subseqüències (ORFS) de cada una de les pautes. #posem cada ORF a una posició d'un vector #@vsubseq: matriu amb dos files: fila 0: vector amb un ORF en cada posicio; fila 1: posicions inicial de cada ORFrespecte la sequencia original. my @vsubseq0 = &subseq (\@vseq0, 0); #pauta 0 my @vsubseq1 = &subseq (\@vseq1, 1); #pauta 1 my @vsubseq2 = &subseq (\@vseq2, 2); #pauta 2 #@matriu amb tots els ORFs de la mateixa pauta, cada fila conté un ORF, cada columna conté els nucleotids de l'ORF my @matriu0 = &matriu (@vsubseq0); my @matriu1 = &matriu (@vsubseq1); my @matriu2 = &matriu (@vsubseq2); ############# ########## ############# EXONS INTERNS ########## ############# ########## #amb la funcio finestra obtenim matrius que contenen: #columna 0: numero del ORF #columna 1: posició del acceptor o del donor #columna 2: puntuació del acceptor o del donor #columna 3: posició inicial del ORF my @acceptors0 = &finestra (\%matriu_acceptors, \@matriu0, \@vsubseq0, 0, $cutoffacc); my @acceptors1 = &finestra (\%matriu_acceptors, \@matriu1, \@vsubseq1, 1, $cutoffacc); my @acceptors2 = &finestra (\%matriu_acceptors, \@matriu2, \@vsubseq2, 2, $cutoffacc); my @donors0 = &finestra (\%matriu_donors, \@matriu0, \@vsubseq0, 0, $cutoffdonn); my @donors1 = &finestra (\%matriu_donors, \@matriu1, \@vsubseq1, 1, $cutoffdonn); my @donors2 = &finestra (\%matriu_donors, \@matriu2, \@vsubseq2, 2, $cutoffdonn); #amb la funcio parelles obtenim matrius que contenen: #columna 0: numero del ORF #columna 1: posició del acceptor #columna 2: puntuació del acceptor #columna 3: posició del donor #columna 4: puntuació del donor #columna 5: seqüència de la posicio d'acceptor a la de donor #columna 6: posicio inicial del ORF #columna 7: pauta del exo my @parelles0 = &parelles (\@acceptors0, \@donors0, $seq0); my @parelles1 = &parelles (\@acceptors1, \@donors1, $seq0); my @parelles2 = &parelles (\@acceptors2, \@donors2, $seq0); #obtenim matrius amb dos files: fila 0: amb els biaixos de cada exo en les diferents columnes; fila 1: amb la pauta de cada exo en les diferents columnes my @loglkh0 = &calcula_loglkh ($seq0, \@parelles0, \%pcodons, 1, \@vsubseq0); my @loglkh1 = &calcula_loglkh ($seq0, \@parelles1, \%pcodons, 1, \@vsubseq1); my @loglkh2 = &calcula_loglkh ($seq0, \@parelles2, \%pcodons, 1, \@vsubseq2); #amb la funcio suma_puntuacions obtenim matrius que contenen: #columna 0: numero del ORF #columna 1: posició del acceptor #columna 2: posició del donor #columna 3: puntuació final (que es la suma de la puntuacio del acceptor+donor+biaix codificant) #columna 4: pauta #columna 5: seqüència de la posicio d'acceptor a la de donor my @puntuacio_internal0= &suma_puntuacio (\@parelles0, \@loglkh0, 1, \@vsubseq0, $seq0, $long); my @puntuacio_internal1= &suma_puntuacio (\@parelles1, \@loglkh1, 1, \@vsubseq1, $seq0, $long); my @puntuacio_internal2= &suma_puntuacio (\@parelles2, \@loglkh2, 1, \@vsubseq2, $seq0, $long); ######### ############# ######### EXONS INICIALS ############# ######### ############# #amb la funcio finestra_atg obtenim matrius que contenen: #columna 0: numero del ORF #columna 1: posició del atg my @atg0 = &finestra_atg (\@vsubseq0, 0); my @atg1 = &finestra_atg (\@vsubseq1, 1); my @atg2 = &finestra_atg (\@vsubseq2, 2); #amb la funcio parelles_atg obtenim matrius que contenen: #columna 0: numero del ORF #columna 1: posició del atg #columna 2: valor de 0 (per poder després reutilitzar un funció posterior #columna 3: posició del donor #columna 4: puntuació del donor #columna 5: seqüència de la posicio d'atg a la de donor my @parelles_atg0 = &parelles_atg(\@donors0, \@atg0, $seq0); my @parelles_atg1 = &parelles_atg(\@donors1, \@atg1, $seq0); my @parelles_atg2 = &parelles_atg(\@donors2, \@atg2, $seq0); #obtenim matrius amb dos files: fila 0: amb els biaixos de cada exo en les diferents columnes; fila 1: amb la pauta de cada exo en les diferents columnes, que en el cas dels exons inicial sempre serà igual a 0 my @loglkh_atg0 = &calcula_loglkh ($seq0, \@parelles_atg0, \%pcodons, 2, \@vsubseq0); my @loglkh_atg1 = &calcula_loglkh ($seq0, \@parelles_atg1, \%pcodons, 2, \@vsubseq1); my @loglkh_atg2 = &calcula_loglkh ($seq0, \@parelles_atg2, \%pcodons, 2, \@vsubseq2); #amb la funcio suma_puntuacions obtenim matrius que contenen: #columna 0: numero del ORF #columna 1: posició del atg #columna 2: posició del donor #columna 3: puntuació final (que es la suma de la puntuacio del donor+biaix codificant) #columna 4: pauta #columna 5: seqüència de la posicio d'atg a la de donor my @puntuacio_atg0= &suma_puntuacio (\@parelles_atg0, \@loglkh_atg0, 1, \@vsubseq0, $seq0, $long); my @puntuacio_atg1= &suma_puntuacio (\@parelles_atg1, \@loglkh_atg1, 1, \@vsubseq0, $seq0, $long); my @puntuacio_atg2= &suma_puntuacio (\@parelles_atg2, \@loglkh_atg2, 1, \@vsubseq0, $seq0, $long); ########### ############## ########### EXONS TERMINALS ############## ########### ############## #obtenim matrius amb dos files: fila 0: amb els biaixos de cada exo en les diferents columnes; fila 1: amb la pauta de cada exo en les diferents columnes my @loglkh_terminal0 = &calcula_loglkh ($seq0, \@acceptors0, \%pcodons, 3, \@vsubseq0); my @loglkh_terminal1 = &calcula_loglkh ($seq0, \@acceptors1, \%pcodons, 3, \@vsubseq1); my @loglkh_terminal2 = &calcula_loglkh ($seq0, \@acceptors2, \%pcodons, 3, \@vsubseq2); #amb la funcio suma_puntuacions obtenim matrius que contenen: #columna 0: numero del ORF #columna 1: posició del acceptor #columna 2: posició final de l'ORF #columna 3: puntuació final (que es la suma de la puntuacio del acceptor+donor+biaix codificant) #columna 4: pauta del exo predit #columna 5: seqüència de la posició d'acceptor fins a l'stop codon. my @puntuacio_terminal0= &suma_puntuacio (\@acceptors0, \@loglkh_terminal0, 2, \@vsubseq0, $seq0, $long); my @puntuacio_terminal1= &suma_puntuacio (\@acceptors1, \@loglkh_terminal1, 2, \@vsubseq1, $seq0, $long); my @puntuacio_terminal2= &suma_puntuacio (\@acceptors2, \@loglkh_terminal2, 2, \@vsubseq2, $seq0, $long); ############# ########### ############# SINGLE-GENES ########### ############# ########### #obtenim matrius amb dos files: fila 0: amb els biaixos de cada exo en les diferents columnes; fila 1: amb la pauta de cada exo en les diferents columnes my @loglkh_single0 = &calcula_loglkh ($seq0, \@atg0, \%pcodons, 4, \@vsubseq0); my @loglkh_single1 = &calcula_loglkh ($seq0, \@atg1, \%pcodons, 4, \@vsubseq1); my @loglkh_single2 = &calcula_loglkh ($seq0, \@atg2, \%pcodons, 4, \@vsubseq2); #amb la funcio suma_puntuacions obtenim matrius que contenen: #columna 0: numero del ORF #columna 1: posició del atg #columna 2: posició final de l'ORF #columna 3: puntuació final (del biaix codificant) #columna 4: pauta del exo predit, que serà sempre 0. #columna 5: seqüència de la posició del atg fins a l'stop codon. my @puntuacio_single0= &suma_puntuacio (\@atg0, \@loglkh_single0, 3, \@vsubseq0, $seq0, $long); my @puntuacio_single1= &suma_puntuacio (\@atg1, \@loglkh_single1, 3, \@vsubseq1, $seq0, $long); my @puntuacio_single2= &suma_puntuacio (\@atg2, \@loglkh_single2, 3, \@vsubseq2, $seq0, $long); ############ ############## ############ SORTIDA DE DADES ############## ############ ############## #amb la funcio unio obtenim una matriu que conté les 3 matrius corresponent als 3 ORFs. #columna 0: posicio inicial del exó #columna 1: posició final del exó #columna 3: puntuació final (del biaix codificant) #columna 4: pauta del exó my @internals = &unio (\@puntuacio_internal0, \@puntuacio_internal1, \@puntuacio_internal2); my @inicials = &unio (\@puntuacio_atg0, \@puntuacio_atg1, \@puntuacio_atg2); my @terminals = &unio (\@puntuacio_terminal0, \@puntuacio_terminal1, \@puntuacio_terminal2); my @single = &unio (\@puntuacio_single0, \@puntuacio_single1, \@puntuacio_single2); #ordena els exons de les matrius anteriors amb les posicions en ordre creixent my @internal_ordenat = &ordre (@internals); my @inicial_ordenat = &ordre (@inicials); my @terminal_ordenat = &ordre (@terminals); my @single_ordenat = &ordre (@single); #imprimim els resultats obtinguts en format gff print "VALORES DE ENTRADA:\n\n"; print "Cutoff donors: $cutoffdonn\n"; print "Cutoff acceptors: $cutoffacc\n"; print "Longitud mínima: $long pb\n\n"; print "SECUENCIA: $iden\n"; &gff(\@inicial_ordenat, "iniciales", $inicials); &gff(\@internal_ordenat, "internos", $internals); &gff(\@terminal_ordenat, "terminales", $terminals); &gff(\@single_ordenat, "single_genes", $single); ############################################################# ############################################################# ##################### FUNCIONS ################## ############################################################# ############################################################# #NOM: hash #PROPOSIT: convertir les matrius de donors i acceptors en hashes #PARAMETRES: primer - matriu de pesos #RETORNA: hash amb el contingut de les matrius de pesos sub hash () { open (MATRIU,"< $_[0]"); my %matpesos; # hash on enregistrarem la matriu de pesos my @nucleotids; # vector on enregistrarem els nucleotids my $llegint_matriu=0; # variable indicadora de quan llegim la matriu my $posicio=0; # posicio de la matriu que s'esta llegint while () { if ($llegint_matriu == 0 && m/\AP0\s+(\w+)\s+(\w+)\s+(\w+)\s+(\w+)/) { $nucleotids[0] = $1; $nucleotids[1] = $2; $nucleotids[2] = $3; $nucleotids[3] = $4; $llegint_matriu = 1; } if ($llegint_matriu == 1 && m/\A\d+\s+([\-\d\.]+)\s+([\-\d\.]+)\s+([\-\d\.]+)\s+([\-\d\.]+)/) { $matpesos{$nucleotids[0]}[$posicio] = $1; $matpesos{$nucleotids[1]}[$posicio] = $2; $matpesos{$nucleotids[2]}[$posicio] = $3; $matpesos{$nucleotids[3]}[$posicio] = $4; $posicio = $posicio + 1; } if ($llegint_matriu == 1 && m/\AXX\s+(\d+)/) { $matpesos{P} = $1; $llegint_matriu = 0; } } return %matpesos; } #NOM: subseq #PROPOSIT: Calcular les diferents subseqüències (ORFS) de cada # una de les pautes. #PARAMETRES: primer - vector amb els codons en les diferents posicions # del vector # segon - valor de la pauta de l'ORF #RETORNA: matriu amb dos files: fila 0: vector amb un ORF en cada # posicio; fila 1: posicions inicial de cada ORFrespecte la # sequencia original. sub subseq () { my @vseq= @{$_[0]}; my $pauta= $_[1]; my @vsubseq; #matriu que ens retorna my $i=0; $vsubseq[0][$i] = ""; $vsubseq[1][0]= $pauta; my $e=0; while ($e < scalar(@vseq)){ $vsubseq[0][$i] = $vsubseq[0][$i] . $vseq[$e]; if ($vseq[$e] eq 'TAA' || $vseq[$e] eq 'TGA' || $vseq[$e] eq 'TAG') { $i = $i + 1; $vsubseq[1][$i] = ($e * 3) + 3 + $pauta; #posicio d'inici d'orf $vsubseq[0][$i]=""; } $e = $e+1; } return @vsubseq; } #NOM: matriu #PROPOSIT: tenir cada nucleotid de cada orf separat en una # posicio diferent per després poder identifica les senyals # de donors i acceptors #PARAMETRES: primer - matriu amb dos files: fila 0: vector amb un ORF # en cada posicio; fila 1: posicions inicial de cada ORF # respecte la sequencia original #RETORNA: matriu amb tots els ORFs de la mateixa pauta, cada fila # conté un ORF, cada columna conté els nucleotids de l'ORF sub matriu (){ my @vsubseq=@_; my @matriu; #matriu que ens retorna my $i=0; while ($i= $posicio) { while ($i < (scalar(@{$matriu[$a]}) - $posicio + 1 )) { #recorrent cada orf my $k = 0; #posicio del nucleotid en l'orf (s'ha d'anar guardant) my $p = 0; #probabilitat de la posicio while ($k < $posicio) { $p = $p + $mpes{$matriu[$a][$i+$k]}[$k]; $k = $k + 1; } if ($p > $cutoff){ $vector[0] = $x; #posicions dels senyals respecte tot l'orf $vector[1] = $p; #probabilitat de la posicio $matrix[$b][0] = $a; #orf $matrix[$b][1] = $vector[0]+ $f; #posicio d'inici o final d'exo $matrix[$b][2] = $vector[1]; #puntuacio del senyal de donor o acceptor $matrix[$b][3] = $vsubseq[1][$a]; #posicio inicial de l'orf $b = $b + 1; } $i = $i + 1; $x = $x + 1; } $x = $x + $posicio-1 ; } else { $x=$x+scalar(@{$matriu[$a]}); } $a = $a + 1; } return @matrix; } #NOM: parelles #PROPOSIT: obtenir una matriu que aparelli les senyals d'acceptors amb les # de donors sempre i quan aquestes senyals pertanyin a un mateix # orf i la posicio de donors sigui posterior a la d'acceptors. #PARAMETRES: primer - matriu amb les senyals d'acceptors # segon - matriu amb les senyals de donors # tercer - seqüència original #RETORNA: matriu que conté: # columna 0: numero del ORF # columna 1: posició del acceptor # columna 2: puntuació del acceptor # columna 3: posició del donor # columna 4: puntuació del donor # columna 5: seqüència de la posicio d'acceptor a la de donor # columna 6: posicio inicial del ORF # columna 7: pauta del exo sub parelles (){ my @macc = @{$_[0]}; #matriu amb posicions i puntuacions d'acceptors my @mdon = @{$_[1]}; #matriu amb posicions i puntuacions de donors my $seq0 =$_[2]; my @seq = split(//, $seq0); my $a = 0; #files de la matriu acceptors my $b = 0; #files de la matriu donors my $i = 0; #files de la matriu @parelles que te 4 columnes my @parelles; #matriu @parelles while ($a < scalar (@macc)) { while ($b < scalar(@mdon)) { if ($mdon[$b][1] > $macc[$a][1]+10 && $mdon[$b][0] == $macc[$a][0]) { my $c = $macc[$a][1]; my $x =""; while ($c <= $mdon[$b][1]) { $x = $x . $seq[$c]; $c = $c + 1; } $parelles[$i][0]= $mdon[$b][0]; #numero de l'orf $parelles[$i][1]= $macc[$a][1]; #posicio d'inici d'exo $parelles[$i][2]= $macc[$a][2]; #puntuacio del senyal acceptor $parelles[$i][3]= $mdon[$b][1]; #posicio de final d'exo $parelles[$i][4]= $mdon[$b][2]; #puntuacio del senyal de donor $parelles[$i][5]= $x; #orf $parelles[$i][6]= $mdon[$b][3]; #posicio inicial absoluta de l'orf $i=$i+1; } $b=$b+1; } $b=0; $a=$a+1; } return @parelles; } #NOM: finestra_atg #PROPOSIT: obtenir un matriu amb les posicions dels codons inicials (atg) #PARAMETRES: primer - matriu amb dos files: fila 1: vector amb un orf # en cada posicio del vector. fila 2: posicions inicial # de l'orf respecte la sequencia original. # segon - valor de la pauta de lectura que servirà per tenir # les posicions respecte la seqüència original #RETORNA: matriu que conté: # columna 0: numero del ORF # columna 1: posició del atg sub finestra_atg (){ my @vsubseq = @{$_[0]}; my $f=$_[1]; my $s=0; my $x=0; #posicio absoluta del atg my @v; my $a=0; my @atg; #matriu que retorna while ($s < scalar(@{$vsubseq[0]})) { #estem a cada orf $i = 0; @v = ($vsubseq[0][$s] =~ m/.../g); while ($i < scalar(@v)){ if ($v[$i] eq "ATG"){ $atg[$a][0]=$s; #numero de l'orf $atg[$a][1]=$x+$f; #posicio del atg respecte la sequencia original $x=$x+3; $a=$a+1; } else { $x=$x+3; } $i=$i+1; } $s=$s+1; } return @atg; #matriu amb dos columnes:numero de subsequencia, posicio de l'atg } #NOM: parelles_atg #PROPOSIT: obtenir una matriu que aparelli les posicions dels codons # inicials amb les posicions dels senyals de donors, sempre i # quan pertanyin al mateix orf i la posició del senyal de donor # sigui posterior a la posició del codó d'inici. #PARAMETRES: primer - matriu amb les senyals de donors # segon - matriu amb les senyals de atgs # tercer - sequencia original #RETORNA: matriu que conté: # columna 0: numero del ORF # columna 1: posició del atg # columna 2: valor de 0 (per poder després reutilitzar un # funció posterior # columna 3: posició del donor # columna 4: puntuació del donor # columna 5: seqüència de la posicio d'atg a la de donor sub parelles_atg (){ my @mdon = @{$_[0]}; #matriu amb posicions i puntuacions de donors my @atg=@{$_[1]}; #matriu amb posicions i puntuacions de atgs my $seq0 =$_[2]; my @seq = split(//, $seq0); my $a = 0; #files de la matriu acceptors my $b = 0; #files de la matriu donors my $i = 0; #files de la matriu @parelles_atg que te 4 columnes my @parelles_atg; #matriu atgs-donors while ($a < scalar (@atg)) { while ($b < scalar(@mdon)) { if ($mdon[$b][1] > $atg[$a][1]+10 && $mdon[$b][0] == $atg[$a][0]) { my $c = $atg[$a][1]; my $x =""; while ($c <= $mdon[$b][1] ) { $x = $x . $seq[$c]; $c = $c + 1; } $parelles_atg[$i][0] = $mdon[$b][0]; #numero de l'orf $parelles_atg[$i][1] = $atg[$a][1]; #posicio atg $parelles_atg[$i][2] = 0; #valor sempre igual a 0 $parelles_atg[$i][3] = $mdon[$b][1]; #posicio del senyal de donor $parelles_atg[$i][4] = $mdon[$b][2]; #puntuacio del senyal de donors $parelles_atg[$i][5]= $x; #orf $i=$i+1; } $b=$b+1; } $b=0; $a=$a+1; } return @parelles_atg; #matriu amb 5 columnes: subsequencia, posicio atg, numero 0, posicio donor, puntuacio. } #NOM: calcula_loglkh #PROPOSIT: calcular el biaix codificant dels exons predits #PARAMETRES: primer - seqüència original # segon - matriu amb els parells de senyals aparellats en # cada tipus d'exó: inicials (atg/donor), interns # (acceptor/donor),terminals (acceptor/stop), # single-genes (atg/stop) # tercer - hash amb els diferents codons i els seus valors de # probabilitat # quart - valor utilitzat per tractar diferents aquesta funció # depenent del tipus d'exó que estiguem analitzant # cinquè - matriu amb dos files: fila 0: vector amb un ORF en # cada posicio; fila 1: posicions inicial de cada ORF # respecte la sequencia original. #RETORNA: matrius amb dos files: fila 0: amb els biaixos de cada exo en # les diferents columnes; fila 1: amb la pauta de cada exo en # les diferents columnes sub calcula_loglkh (){ # ara solament ens queda ja calcular la rao de versemblanc,a en escala # logaritmica. la rao de versemblanc,a es la probabilitat conjunta dels # codons segons la taula donada, respecte (dividit per) la probabilitat # conjunta uniforme dels codons (la probabilitat uniforme d'un unic codo # es 1/64). aixo ho calcularem en escala logaritmica per tal de poder # manipular els nombres molt petits que es poden produir al multiplicar # moltes xifres que tenen valors entre 0 i 1. el logaritme del producte # es la suma de logaritmes i el logaritme de la divisio es la resta de # logaritmes. my @loglkh; #matriu que ens retorna my $seq = $_[0]; my @parelles =@{$_[1]}; my %pcodons= %{$_[2]}; my $n= $_[3]; my @vsubseq= @{$_[4]}; my $x; #variable util per calcular la pauta my $pauta; #valor de la pauta que pot ser 0, 1, 2 my $inici; #posicio a partir de la que calcularem el biaix my $longitud; #longitud de la sequecia en que calcularem el biaix my $a=0; my $q; while ($a< scalar(@parelles)) { if ($n==1) { #pels exons interns $x=$parelles[$a][1]-$parelles[$a][6]; if ($x%3==0){ $pauta=0; } if ($x%3==1){ $pauta=2; } if ($x%3==2){ $pauta=1; } $loglkh[1][$a]=$pauta; $inici= $parelles[$a][1] + $pauta; $longitud=$parelles[$a][3] - $inici +1; } if ($n==2){ #pels exons inicials $inici = $parelles[$a][1]; $loglkh[1][$a]=0; $longitud=$parelles[$a][3] - $inici +1; } if ($n==3){ #pels exons terminals $x=$parelles[$a][1]-$parelles[$a][3]; if ($x%3==0){ $pauta=0; } if ($x%3==1){ $pauta=2; } if ($x%3==2){ $pauta=1; } if ($parelles[$a][0] < scalar (@{$vsubseq[0]})-1) { $inici = $parelles[$a][1]; $loglkh[1][$a]=$pauta; $q = $parelles[$a][0]; $longitud= $vsubseq[1][$q+1] - $inici; #-1+1=0 } } if ($n==4){ #single-genes if ($parelles[$a][0] < scalar (@{$vsubseq[0]})-1) { $q = $parelles[$a][0]; $inici = $parelles[$a][1]; $loglkh[1][$a]=0; $longitud = $vsubseq[1][$q+1] - $inici; } } my $seqfinestra =substr ($seq, $inici, $longitud); my @codons = ($seqfinestra =~m/.../g); my $i=0; my $r=0; while ($i < scalar(@codons)) { $r = $r + log($pcodons{$codons[$i]}) - log(1/64); $i = $i + 1; } $loglkh[0][$a]=$r; $a=$a+1; } return @loglkh; #matriu amb dos files: fila 1 té els logaritmes i la fila dos la pauta de cada logaritme. } #NOM: suma_puntuacio #PROPOSIT: obtenir una matriu que contingui les posicions dels predits # amb les seves puntuacions finals després de sumar el biaix # codificant. #PARAMETRES: primer - matriu que conté les posicions dels exons predits # segon - matriu que conté el biaix de cada exó predit # tercer - valor que permet tracta aquesta funció diferent en # funció del tipus d'exó. # quart - matriu que conté els orf amb les seves posicions inicial # cinquè - seqüència original # sisè - longitud mínima dels exons predits que considerem #RETORNA: matrius que contenen: # columna 0: numero del ORF # columna 1: posició del acceptor # columna 2: posició del donor # columna 3: puntuació final (que es la suma de la puntuacio del # acceptor+donor+biaix codificant) # columna 4: pauta # columna 5: seqüència de la posicio d'acceptor a la de donor sub suma_puntuacio () { my @puntuacio; #matriu que ens retorna my @parelles =@{$_[0]}; my @loglkh =@{$_[1]}; my $r=$_[2]; my @vsubseq=@{$_[3]}; my $seq= $_[4]; my $long=$_[5]; my $i=0; my $b=0; my $q; my $sequencia; while ($i$long){ $puntuacio[$b][0]=$parelles[$i][0]; $puntuacio[$b][1]=$parelles[$i][1]+1; $puntuacio[$b][4]=$loglkh[1][$i]; $puntuacio[$b][2]=$parelles[$i][3]+1; $puntuacio[$b][3]=$parelles[$i][2]+$parelles[$i][4]+$loglkh[0][$i]; $puntuacio[$b][5]=$parelles[$i][5]; $b=$b+1; } } if ($r==2) { #per exons terminals if ($parelles[$i][0] < scalar (@{$vsubseq[0]})-1) { $q = $parelles[$i][0]; if ($vsubseq[1][$q+1]-$parelles[$i][1]>$long){ $puntuacio[$b][0]=$parelles[$i][0]; $puntuacio[$b][1]=$parelles[$i][1]+1; $puntuacio[$b][4]=$loglkh[1][$i]; $puntuacio[$b][2]= $vsubseq[1][$q+1]; $puntuacio[$b][3]=$parelles[$i][2]+$loglkh[0][$i]; $sequencia=substr($seq, $parelles[$i][1], $vsubseq[1][$q+1]-$parelles[$i][1]); $puntuacio[$b][5]= $sequencia; $b=$b+1; } } } if ($r==3){ #single-genes if ($parelles[$i][0] < scalar (@{$vsubseq[0]})-1) { $q = $parelles[$i][0]; if ($vsubseq[1][$q+1]-1-$parelles[$i][1]>$long){ $puntuacio[$b][0]=$parelles[$i][0]; $puntuacio[$b][1]=$parelles[$i][1]+1; $puntuacio[$b][4]=$loglkh[1][$i]; $puntuacio[$b][2]= $vsubseq[1][$q+1]; $puntuacio[$b][3]=$loglkh[0][$i]; $sequencia=substr($seq, $parelles[$i][1], $vsubseq[1][$q+1]-$parelles[$i][1]); $puntuacio[$b][5]= $sequencia; $b=$b+1; } } } $i=$i+1; } return @puntuacio; } #NOM: unio #PROPOSIT: fusionar les 3 matriu que hem obtingut de les 3 pautes per # obtenir una única matriu #PARAMETRES: primer - matriu amb les posicions d'inici i final d'exons i # les seves puntuacions per la pauta 0 # segon - matriu amb les posicions d'inici i final d'exons i # les seves puntuacions per la pauta 1 # segon - matriu amb les posicions d'inici i final d'exons i # les seves puntuacions per la pauta 2 #RETORNA: una matriu que conté les 3 matrius corresponent als 3 ORFs. # columna 0: posicio inicial del exó # columna 1: posició final del exó # columna 3: puntuació final (del biaix codificant) # columna 4: pauta del exó sub unio (){ my @unio; #matriu que ens retorna my @pauta0 =@{$_[0]}; my @pauta1 =@{$_[1]}; my @pauta2 =@{$_[2]}; my $n=0; my $i=0; my $j=0; my $a=0; while ($i ${$b}[0] ? ${$a}[0] <=> ${$b}[0] : ${$a}[1] <=> ${$b}[1]}(@matriu); return @matriu; } #NOM: gff #PROPOSIT: imprimim els resultats obtinguts en format gff #PARAMETRES: primer - matriu amb les posicions i puntuacions dels # diferents exons # segon - nom que utilitzem per visualitzar en pantalla # tercer - variable que pot pendre el valor de 0 o 1 i que # permet visualitzar només aquells exons que ens # interessin #RETORNA: els resultats en format gff sub gff () { my @punt =@{$_[0]}; my $tipus =$_[1]; my $i=0; if ($_[2]==1) { print "Exones $tipus(+) predecidos en la secuencia $iden\n"; while ($i