es-MXen-US

Laboratorio de Aprendizaje e Investigación en Cómputo Biológico

LAICBIO

 

 

Regresa al curso de evolución molecular

 

Esta página contiene ejercicios y tutoriales. Los ejercicios están pensados para mostrar aspectos del funcionamiento de los programas. Por otro lado, los tutoriales pretenden mostrar cómo se usan los programas para hacer un análisis determinado. Ambos, ejericios y tutoriales, forman parte del curso de evolución molecular.

 

Contenido

 

Ejercicios

Ejercicio 1: Búsqueda de secuencias homólogas con BLAST

Ejercicio 2: Búsqueda de homólogos de la proteína S (spike)

Ejercicio 3: Alineación múltiple de secuencias (proteína S)

Ejercicio 4: Infiere una filogenia

Tutoriales

Tutorial 1: Búsqueda de secuencias homólogas con BLAST

Tutorial 2: Alineación múltiple de proteínas

Tutorial 3: Infiere una filogenia utilizando una alineación múltiple de proteínas

Tutorial 4: Infiere una filogenia utilizando una alineación múltiple de genes

(en construcción)

 

 

 

Ejercicio 1: Búsqueda de secuencias homólogas con BLAST

 

Una de las tareas más comunes en bioinformática y evolución molecular es buscar genes o proteínas homólogas. Para ello, los biólogos utilizamos comunmente el software BLAST. En este ejercicio mostraremos cómo podemos buscar secuencias homólogas con BLAST. Además de ello, veremos cómo podemos obtener mayor provecho de BLAST utilizando otras herramientas bioinformáticas complementarias. La metodología que mostramos en este tutorial la ulilizamos para publicar este artículo.

 

Durante los ejercicios, utilizaremos la fuente verdana para indicar que se trata de instrucciones que tienes que escribir en la ventana del sistema en Linux. Para las instrucciones en R para RStudio utilizaremos la fuente verdana. Además, el código de los scripts de Perl que utilizaremos te los mostraremos con la fuente Courier New. Estos scripts los puedes copiar y pegar en un editor de texto simple para usarlos. Sólo asegúrate de que la terminación de los scripts que copies y guardes con un editor de texto sea .pl y no .txt.

 

Paso 1 BLAST

Para hacer el tutorial utilizaremos la siguiente proteína:

 

>XP_018228376

MIKLDHFFFF VRLLCFIVAI IWIFLLPFDY FNRKVYISEN AILPGSVNMY FGGSDSNVLE 

AYRDEIEFIN QQSKNVRFDS LHDIFRKLGL RTAIQDYKIQ FKGKNHSGTN FYAIFDTPRS 

DGTEALVLSA AWKNMNGELN KGGVSLVLAL ARYFKRWSLW SKDIIFLIPE EKEVGVQVWV 

DAYHKIDYEE GVSRLIMKGG EIQAVVDIDF VSEFREFDTI ELLYDGINGQ LSNLDLLNTV 

NHIIQSKSGI KIKMQNVLQG KNSYFRKLLT MINGMISQCF GVLSGTQSCF IPYKIDAITL 

KVHSKKNSEY DDILLGKIIE SLFRSLNNLL EHLHQSFFFY LMLSEKRFVS IATYLPSAIL 

IASSFTIDAF SIWLSTFSNF YKNLSSSKEK LYSSKEKKNK LKYTGMDSEL YFKKQFILSI 

STVIVIHIIA YIFLKILYNS VNYKGLFKDK KAYLMFLLFL QILQHIISIL IQKILKKIIK 

SQNMSFFYNY RLNKVFSQAI LGVELSTLAI VNFSLSLFIG IVTFFLSFVR PAKTKFQQIL 

ISIVLELINP LNWIFLFSYY KNKEVKDLMH LFAFGWKVWY LWTPTVLWLL WWPIWIMTKI 

ILCMPIKDHN E

 

 

Esta proteína pertenece a un hongo ascomiceto llamado Pneumocystis jirovecii RU7. La proteína es parte de un complejo transamidasa que participa en la biosíntesis del GPI (glicosil-fostatidil-inositol). Puedes aprender un poco más de esta proteína aquí y también puedes leer por qué esta proteína es interesante desde el punto de vista de la evolución molecular aquí.

 

Lo primero que tienes que hacer es ingresar a la página de BLAST en el NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi) y seleccionar la opción Protein BLAST (Figura 1.1). Una vez ahí, copia y pega la secuencia de proteína en formato FASTA en el recuadro blanco (Figura 1.2) y ajusta los parámetros de BLAST para que obtengas un máximo de 5000 secuencias homólogas con un límite máximo de E-value de (Expected threshold < 0.01) (Figura 1.3). Observa que el BLAST lo estamos realizando contra la base de datos 'Non-redundant protein sequences (nr)'. Después de ello, da click en el botón de BLAST.

 

 

 

Figura 1.1. Página web de BLAST en NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

 

 

 

Figura 1.2. Copia y pega la secuencia FASTA en el recuardo blanco.

 

 

 

Figura 1.3. Ajusta los parámetros de BLAST para que obtengas un máximo de 5000 secuencias (hits) y con un límite máximo de E-value < 0.01.

 

 

Una vez que tengas el resultado de BLAST, deberás familiarizarte con ellos (Figura 1.4). Observa que la página principal, te proporciona información básica de tu secuencia 'Query' que enviaste con BLAST, tal como su descripción y su longitud. Además, te permite filtrar los resultados por porcentaje de identidad, E-value o cobertura de la secuencia 'Query' por las secuencias que encontró en la base de datos. Además, la página principal tiene cuatro pestañas: 'Descriptions', 'Graphic Summary', 'Alignments' y 'Taxonomy'. Vamos a revisar estas pestañas. 

 

 

 

Figura 1.4. Resultado de BLAST del 'Query' XP_018228376 contra la base de datos nr.

 

En primer lugar, la pestaña de 'Descriptions' contiene una lista con las secuencias que BLAST encontró en la base de datos que cumplieron con un E-value < 0.01. Para cada secuencia se muestra el bit Score (Max y Total Score),  la cobertura de la secuencia 'query' en la alineación pareada (Query Cover), el E-value, el porcentaje de identidad en la región que BLAST logró alinear y el identificador de las proteínas 'subject' que BLAST encontró.

 

En segundo lugar, si miras la pestaña 'Graphic Summary' verás una imagen que representa a las alineaciones pareadas de cada una de las secuencias 'subject' con respecto a la secuencia 'query' (Figura 1.5). El color de cada alineación (línea) revela la significancia estadística en términos del bit score. También puedes observar que esta proteína tiene un dominio llamado: 'Zinc_peptidase_like'. 

 

 

 

Figura 1.5. Resúmen gráfico del resultado de BLAST.

 

En tercer lugar, revisa la pestaña que contiene las alineaciones pareadas. Para cada alineación pareada, se muestra el bit score y el E-value además de otros estadísticos. Además, puedes ver los detalles de la alineación pareada (Figura 1.6).

 

 

 

Figura 1.6. Alineación pareada de BLAST entre la proteína 'query' y una de las proteínas 'subject'.

 

Finalmente revisa la pestaña de 'Taxonomy'. Podrás ver la clasificación taxonómica de las especies de donde se obtuvieron las secuencias.

 

A continuación, regresa a la pestaña 'Descriptions', selecciona todas las secuencias palomeando la casilla 'select all' y en el menú que se despliega con la opción 'Download' selecciona la opción 'Hit Table (CSV)'. Cambia el nombre al archivo que obtienes a: Alignment-HitTable-vs-nr.csv.

 

Ahora, vuelve a hacer el BLAST, pero en lugar de hacerlo contra la base de datos 'nr', realízalo contra 'Ascomycota'. Para hacer nuevamente el BLAST, simplemente selecciona la opción Edit Search y selecciona el grupo biológico (Figura 1.7). Recuerda mantener el resto de los parámetros iguales a los de la búsqueda anterior.

 

 

Figura 1.7. Realiza nuevamente el BLAST pero ahora contra Ascomycota (taxid:4890).

 

Nuevamente, selecciona todas las proteínas y descarga el archivo 'Hit Table (CSV)'. Cambia el nombre al archivo que obtienes a: Alignment-HitTable-vs-Ascomycota.csv.

 

Los archivos 'Hit Table (CSV)' contienen las siguientes columnas:

query

subject

% de identidad

longitud de la alineación

mismatches

apertura de gaps

inicio de la alineación en la secuencia query

fin de la alineación en la secuencia query

inicio de la alineación en la secuencia subject

fin de la alineación en la secuencia subject

E-value

bit score

% de positivos

 

A continuación, vamos a utilizar un script en R para graficar los bit score de ambas búsquedas (Figura 1.8). También graficamos el E-value (Figura 1.9).

 

------------------------------------------------------------------------------------------

tabla_Nr = read.table(file="Alignment-HitTable-vs-nr.csv", header=FALSE, sep=",")

tabla_As = read.table(file="Alignment-HitTable-vs-Ascomycota.csv", header=FALSE, sep=",")

boxplot(tabla_Nr$V12, tabla_As$V12, names = c('nr-database','Ascomycota'), main = 'bit scores')
boxplot(log(tabla_Nr$V11), log(tabla_As$V11), names = c('nr-database','Ascomycota'), main = 'E-value')
------------------------------------------------------------------------------------------

 

 

Figura 1.8. Diferencias en los valores de bit socre entre las búsquedas hechas con la misma secuencia pero bases de datos de distintos tamaños.

 

 

 

Figura 1.9. Diferencias en los valores de log(E-value) entre las búsquedas hechas con la misma secuencia pero a bases de datos de distintos tamaños. 

 

¿A qué se debe la diferencia entre los valores que se muestran en las Figuras 1.8 y 1.9?

 

 

Ejercicio 2: Búsqueda de homólogos de la proteína S (spike)

 

Vas a realizar un BLAST de la proteína S del SARS-CoV-2, del genoma de referencia. Puedes acceder a la versión RefSeq del genoma aquí. La proteína S tiene el identificador YP_009724390.1, pero vamos a hacer la búsqueda con el gen cuyo número de identificador es 43740568. La secuencia del gen la tienes a continuación.

 

>MN908947.3:21563-25384 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome

ATGTTTGTTTTTCTTGTTTTATTGCCACTAGTCTCTAGTCAGTGTGTTAATCTTACAACCAGAACTCAAT

TACCCCCTGCATACACTAATTCTTTCACACGTGGTGTTTATTACCCTGACAAAGTTTTCAGATCCTCAGT

TTTACATTCAACTCAGGACTTGTTCTTACCTTTCTTTTCCAATGTTACTTGGTTCCATGCTATACATGTC

TCTGGGACCAATGGTACTAAGAGGTTTGATAACCCTGTCCTACCATTTAATGATGGTGTTTATTTTGCTT

CCACTGAGAAGTCTAACATAATAAGAGGCTGGATTTTTGGTACTACTTTAGATTCGAAGACCCAGTCCCT

ACTTATTGTTAATAACGCTACTAATGTTGTTATTAAAGTCTGTGAATTTCAATTTTGTAATGATCCATTT

TTGGGTGTTTATTACCACAAAAACAACAAAAGTTGGATGGAAAGTGAGTTCAGAGTTTATTCTAGTGCGA

ATAATTGCACTTTTGAATATGTCTCTCAGCCTTTTCTTATGGACCTTGAAGGAAAACAGGGTAATTTCAA

AAATCTTAGGGAATTTGTGTTTAAGAATATTGATGGTTATTTTAAAATATATTCTAAGCACACGCCTATT

AATTTAGTGCGTGATCTCCCTCAGGGTTTTTCGGCTTTAGAACCATTGGTAGATTTGCCAATAGGTATTA

ACATCACTAGGTTTCAAACTTTACTTGCTTTACATAGAAGTTATTTGACTCCTGGTGATTCTTCTTCAGG

TTGGACAGCTGGTGCTGCAGCTTATTATGTGGGTTATCTTCAACCTAGGACTTTTCTATTAAAATATAAT

GAAAATGGAACCATTACAGATGCTGTAGACTGTGCACTTGACCCTCTCTCAGAAACAAAGTGTACGTTGA

AATCCTTCACTGTAGAAAAAGGAATCTATCAAACTTCTAACTTTAGAGTCCAACCAACAGAATCTATTGT

TAGATTTCCTAATATTACAAACTTGTGCCCTTTTGGTGAAGTTTTTAACGCCACCAGATTTGCATCTGTT

TATGCTTGGAACAGGAAGAGAATCAGCAACTGTGTTGCTGATTATTCTGTCCTATATAATTCCGCATCAT

TTTCCACTTTTAAGTGTTATGGAGTGTCTCCTACTAAATTAAATGATCTCTGCTTTACTAATGTCTATGC

AGATTCATTTGTAATTAGAGGTGATGAAGTCAGACAAATCGCTCCAGGGCAAACTGGAAAGATTGCTGAT

TATAATTATAAATTACCAGATGATTTTACAGGCTGCGTTATAGCTTGGAATTCTAACAATCTTGATTCTA

AGGTTGGTGGTAATTATAATTACCTGTATAGATTGTTTAGGAAGTCTAATCTCAAACCTTTTGAGAGAGA

TATTTCAACTGAAATCTATCAGGCCGGTAGCACACCTTGTAATGGTGTTGAAGGTTTTAATTGTTACTTT

CCTTTACAATCATATGGTTTCCAACCCACTAATGGTGTTGGTTACCAACCATACAGAGTAGTAGTACTTT

CTTTTGAACTTCTACATGCACCAGCAACTGTTTGTGGACCTAAAAAGTCTACTAATTTGGTTAAAAACAA

ATGTGTCAATTTCAACTTCAATGGTTTAACAGGCACAGGTGTTCTTACTGAGTCTAACAAAAAGTTTCTG

CCTTTCCAACAATTTGGCAGAGACATTGCTGACACTACTGATGCTGTCCGTGATCCACAGACACTTGAGA

TTCTTGACATTACACCATGTTCTTTTGGTGGTGTCAGTGTTATAACACCAGGAACAAATACTTCTAACCA

GGTTGCTGTTCTTTATCAGGATGTTAACTGCACAGAAGTCCCTGTTGCTATTCATGCAGATCAACTTACT

CCTACTTGGCGTGTTTATTCTACAGGTTCTAATGTTTTTCAAACACGTGCAGGCTGTTTAATAGGGGCTG

AACATGTCAACAACTCATATGAGTGTGACATACCCATTGGTGCAGGTATATGCGCTAGTTATCAGACTCA

GACTAATTCTCCTCGGCGGGCACGTAGTGTAGCTAGTCAATCCATCATTGCCTACACTATGTCACTTGGT

GCAGAAAATTCAGTTGCTTACTCTAATAACTCTATTGCCATACCCACAAATTTTACTATTAGTGTTACCA

CAGAAATTCTACCAGTGTCTATGACCAAGACATCAGTAGATTGTACAATGTACATTTGTGGTGATTCAAC

TGAATGCAGCAATCTTTTGTTGCAATATGGCAGTTTTTGTACACAATTAAACCGTGCTTTAACTGGAATA

GCTGTTGAACAAGACAAAAACACCCAAGAAGTTTTTGCACAAGTCAAACAAATTTACAAAACACCACCAA

TTAAAGATTTTGGTGGTTTTAATTTTTCACAAATATTACCAGATCCATCAAAACCAAGCAAGAGGTCATT

TATTGAAGATCTACTTTTCAACAAAGTGACACTTGCAGATGCTGGCTTCATCAAACAATATGGTGATTGC

CTTGGTGATATTGCTGCTAGAGACCTCATTTGTGCACAAAAGTTTAACGGCCTTACTGTTTTGCCACCTT

TGCTCACAGATGAAATGATTGCTCAATACACTTCTGCACTGTTAGCGGGTACAATCACTTCTGGTTGGAC

CTTTGGTGCAGGTGCTGCATTACAAATACCATTTGCTATGCAAATGGCTTATAGGTTTAATGGTATTGGA

GTTACACAGAATGTTCTCTATGAGAACCAAAAATTGATTGCCAACCAATTTAATAGTGCTATTGGCAAAA

TTCAAGACTCACTTTCTTCCACAGCAAGTGCACTTGGAAAACTTCAAGATGTGGTCAACCAAAATGCACA

AGCTTTAAACACGCTTGTTAAACAACTTAGCTCCAATTTTGGTGCAATTTCAAGTGTTTTAAATGATATC

CTTTCACGTCTTGACAAAGTTGAGGCTGAAGTGCAAATTGATAGGTTGATCACAGGCAGACTTCAAAGTT

TGCAGACATATGTGACTCAACAATTAATTAGAGCTGCAGAAATCAGAGCTTCTGCTAATCTTGCTGCTAC

TAAAATGTCAGAGTGTGTACTTGGACAATCAAAAAGAGTTGATTTTTGTGGAAAGGGCTATCATCTTATG

TCCTTCCCTCAGTCAGCACCTCATGGTGTAGTCTTCTTGCATGTGACTTATGTCCCTGCACAAGAAAAGA

ACTTCACAACTGCTCCTGCCATTTGTCATGATGGAAAAGCACACTTTCCTCGTGAAGGTGTCTTTGTTTC

AAATGGCACACACTGGTTTGTAACACAAAGGAATTTTTATGAACCACAAATCATTACTACAGACAACACA

TTTGTGTCTGGTAACTGTGATGTTGTAATAGGAATTGTCAACAACACAGTTTATGATCCTTTGCAACCTG

AATTAGACTCATTCAAGGAGGAGTTAGATAAATATTTTAAGAATCATACATCACCAGATGTTGATTTAGG

TGACATCTCTGGCATTAATGCTTCAGTTGTAAACATTCAAAAAGAAATTGACCGCCTCAATGAGGTTGCC

AAGAATTTAAATGAATCTCTCATCGATCTCCAAGAACTTGGAAAGTATGAGCAGTATATAAAATGGCCAT

GGTACATTTGGCTAGGTTTTATAGCTGGCTTGATTGCCATAGTAATGGTGACAATTATGCTTTGCTGTAT

GACCAGTTGCTGTAGTTGTCTCAAGGGCTGTTGTTCTTGTGGATCCTGCTGCAAATTTGATGAAGACGAC

TCTGAGCCAGTGCTCAAAGGAGTCAAATTACATTACACATAA

 

Para hacer el BLAST, dirígete a la página del programa en el NCBI (Figura 2.1) y selecciona la opción 'Search Betacoronavirus Database'.

 

 

Figura 2.1. Selecciona la opción 'Search Betacoronavirus Database'.

 

Una vez en la página de BLAST, coloca la secuencia del gen de la proteína S en el recuadro blanco y cambia uno de los parámetros para permitir que puedas obtener hasta 5000 hits (Figura 2.2).

 

 

 

Figura 2.2. Cambia el número máximo de hits a 5000.

 

Una vez que obtengas los resultados, selecciona todas las secuencias (revisa que todas ellas tiengan un Evalue significativo) y descarga tanto el archivo 'FASTA aligned sequences' como el 'HIT Table (csv)' (Figura 2.3). Revisa también la pestaña de 'Taxonomy' para que veas los organismos de donde se obtuvieron los distintos coronavirus. Trabajaremos con estas secuencias.

 

 

 

Figura 2.3. Descarga las secuencias alineadas y la tabla csv.

 

 

Ejercicio 3: Alineación múltiple de secuencias (proteína S)

 

Vamos ahora a hacer una alineación múltiple de los genes que codifican para la proteína S y que encontraste con la búsqueda de BLAST. Antes de continuar, es necesario que sepas que nos enfrentamos a un problema que no es trivial. La búsqueda de BLAST arrojó una gran cantidad de secuencias homólogas (en mi caso, fueron más de 5000), muchas de ellas muy similares. Ahora bien, los programas que vamos a utilizar, no estan diseñados para funcionar de forma eficiente con tantas secuencias. Esta dificultad se podría resolver trabajando con la línea de comandos de Linux y programando algunos scrips simples en R o Perl. Sin embargo, dado que este curso se está impartiendo en línea, por simplicidad nos ajustaremos a las herramientas que se ofrecen en Internet. No es la manera más eficiente ni tampoco la que usaría un bioinformático, pero para los fines del curso, funciona.

 

Nuestro objetivo final, es inferir la evolución del gen que codifica para la proteína S del SARS-CoV-2. Para ello, necesitamos elegir las secuencias que vamos a estudiar. Si observas el resultado de BLAST (por ejemplo, en el archivo 'HIT Table (csv)' o si realizas nuevamente el BLAST del ejercicio 2) encontrarás que muchas de las secuencias son prácticamente idénticas a la secuencia query. Estas secuencias idénticas no nos aportan nueva información filogenética. Lo primero que tendrás que hacer es elegir aquellas secuencias de SARS-CoV-2 que sean distintas a tu secuencia query y que no sean de otros coronavirus. Puedes hacerlo manualmente en los archivos que ya descargaste o bien, realizar nuevamente el BLAST y filtrar los resultados. Para filtrar los resultados, puedes escribir en la casilla superior derecha en 'Organism' SARS-CoV2 (taxid:2697049) y en 'Percent identity' '1 to 99' (Figura 3.1).

 

 

 

Figura 3.1. Es posible filtrar los resultados de BLAST para identificar y seleccionar más fácilmente secuencias.

 

Selecciona unas cuantas secuencias (pueden ser unas 10) y descárgalas como 'FASTA (aligned sequences)'. Es buena idea que cuando las selecciones, veas las alineaciones pareadas de las secuencias que vas a descargar en 'Alignments' y también revises la pestaña de 'Graphic Summary'. Ello con el fin de seleccionar secuencias que alinen bien con tu secuencia query y que no vayan a tener muchos nucleótidos indefinidos (es decir 'N'). En mi caso, he seleccionado simplemente las 20 primeras secuencias. El archivo que descargué, lo renombré como SARS-CoV-2.fasta y le agregué la secuencia query.

 

A continuación es necesario descargar secuencias de otros coronavirus distintos a SARS-CoV-2. Para ello, borra los porcentajes de identidad que habías usado nuevamente y simplemente escribe de nuevo en 'Organism' SARS-CoV2 (taxid:2697049) y palomea la opción 'exclude'. De esta forma se desplegarán en la pantalla las secuencias homólogas de otros coronavirus distintos a SARS-CoV-2. Entre estas secuencias están los co-descendientes más cercanos al SARS-CoV-2. Deberás de hacer una selección de secuencias para tratar de encontrar el origen de este virus. Por supuesto, las más similares tienen mayor probabilidad de ser las más cercanas evolutivamente. Haz tu selección de secuencias y descárgalas como 'FASTA (aligned sequences)'. En mi caso, el archivo que descargué lo renombré como SARS-like.fasta. Mi recomendación es que intentes una selección juiciosa de las secuencias.

 

Finalmente, combina los dos archivos (SARS-CoV-2.fasta y SARS-like.fasta) en uno solo (por ejemplo coronavirus-1.fasta). Ahora bien, los encabezados de las secuencias son muy largos y sería recomendable que los abreviaras. De esta forma se podrán apreciar mejor en la filogenia cuando la hagamos. Lo ideal es que conserved el identificador de la secuencia (por ejemplo, el de la secuencia query es MN908947.3). También puedes agregar algún texto que te de más información (por ejemplo SARS-CoV-2 o pangolin). Pero deberás usar solamente números y letras sin acentos. Idealmente los espacios sustitúyelos por guiones bajos '_'. Ello evitará que algunos programas marquen errores. También considera que algunos programas truncan el nombre de los OTUs a 10 caracteres. 

 

Si tienes sistema Linux/MacOS o tienes Windows con ActivePerl y sabes cómo acceder a una terminal, aquí tienes un simple script en Perl que te puede ayudar a simplificar los encabezados de las secuencias.

 

Se usa así:

 

$ perl nombredelscript.pl coronavirus-1.fasta > coronavirus-1e.fasta

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


use strict;


my $file = $ARGV[0];


open (MIA, "$file") or die ("No puedo abrir $file\n");

while (my $linea = <MIA>){

chomp ($linea);

if ($linea =~ />/){

if ($linea =~ />([\d\w\.]+):/){

print (">$1\n");

} else {

die ("el parton fallo: $linea\n");

}

} else {

print ("$linea\n");

}

}

close (MIA);

------------------------------------------------------------------------------------------
 
Además, deberás de tener cuidado en identificar aquellos casos en los que BLAST reporte más de una alineación pareada por secuencia. Dado que BLAST es un algoritmo de alineación local, puede reportar más de una alineación pareada para dos secuencias. En ese caso, tendrás que editar el archivo que descargues de BLAST para que, si se da el caso, fusionar manualmente las distintas regiones de una misma secuencia. La clave para darte cuenta de si tienes esta situación, es que en el archivo de BLAST tendrás distintas secuencias con el mismo identificador.
 

Para hacer la alineación múltiple utilizaremos los servicios que provee el servidor NGPhylogeny.fr. Este servidor ofrece todas las herramientas básicas para inferir una filogenia y nos será muy útil estos días que estamos en cuarentena y sin acceso a computadoras potentes. Puedes leer el artículo que describe los servicios de este servidor aquí (Lemoine et al. 2019).

  

Para comenzar, vamos a realizar una alineación múltiple de las secuencias que encontramos con BLAST. En la página principal de NGPhylogeny.fr selecciona la opción Tools (Figura 3.2).

 

 

Figura 3.2. Selecciona la opción Tools en la página principal de NGPhylogeny.fr.

 

Una vez que entres a Tools, busca el programa MAFFT (también es posible usar otros programas). Ya en la página de MAFFT, encontrarás una opción que te permite subir al servidor el archivo con las secuencias en formato FASTA que encontraste con BLAST (coronavirus-1e.fasta). Rivisa las opciones que te ofrece el programa (no es necesario cambiar ninguna) y da click en 'submit'. 

 

Una vez que termine la alineación múltiple, utiliza el símbolo de (+) asociado al renglón "2. Mafft alignment" para agregar el resultado a tu Workspace (Figura 3.3). Posteriormente, descarga la alineación múltiple con la opción de descarga en el mismo renglón.

 

 

Figura 3.3. Utiliza el símbolo de (+) asociado al renglón "2. Mafft alignment" para agregar a alineación múltiple a tu Workspace.

 

A continuación, en el renglón "2. Mafft alignment" selecciona MSAViewer para ver la alineación múltiple. Es importante que inspecciones la alineación múltiple para que veas que todo está en orden (por ejemplo, que no tienes dos o más veces la misma secuencia). Posteriormente, descarga la alineación con el boton correspondiente (en el boton Export -> Export alignment (FASTA)). Por consistencia, renombra el archivo que descargaste a coronavirus-1e-mafft.fasta.

 

Vamos a visualizar la alineación múltiple utilizando un visualizador más poderoso que se llama JalView. Puedes leer el artículo que lo describe (Waterhouse et al. 2009). Cuando lo hayas abierto, simplemente sigue la ruta 'File -> Input Alignment -> From File' y selecciona el archivo (deberás de tener la opción FASTA seleccionada). Cuando abras el archivo, puedes colorear los nucleótidos con 'Colour -> Nucleotide'. Y bueno ¿en qué fijarnos? sobre todo, queremos que no haya secuencias demasiado pequeñas con respecto a las demas (Figura 3.4). 

 

 

Figura 3.4. Alineación múltiple de genes que codifican para la proteína S. Visualización con JalView.

 

Este ejercicio lo puedes repetir cuantas veces sea necesario hasta que cuentes con una alineación múltiple que te satisfaga. Es decir, que hayas hecho un buen muestreo de la diversidad de secuencias que encontraste con BLAST; que las secuencias tengan relativamente pocas 'N'; que no haya secuencias demasiado pequeñas en relación al resto; y que no haya errores (como por ejemplo, dos secuencias con el mismo nombre). La originalidad de tu trabajo final, va a depender de lo bien que hayas hecho este muestreo.

 

 

Ejercicio 4: Infiere una filogenia

 

A continuación, vamos a inferir una filogenia de las secuencias génicas de la proteína S. Para ello, ingresaremos a la pestaña de "Tools" de NGPhylogeny.fr (Figura 4.1). Ahí encontrarás diversos métodos para realizar análisis filogenéticos. 

 

 

Figura 4.1. En la pestaña "Tools" encontrarás diversos métodos para realizar análisis filogenéticos.

 

En mi caso, utilizaré el programa PhyML+SMS. PhyML (http://www.atgc-montpellier.fr/phyml/) es un programa que infiere filogenias mediante máxima verosimilitud. SMS (smart model selection) es otro programa de la familia de PhyML que identifica el modelo de evolución molecular que mejor se ajusta a la alineación múltiple que usamos. Tú puedes elegir usar otro programa para inferir la filogenia o hacer pruebas con distintos programas. 

 

Una vez en la página de PhyML+SMS, es necesario especificar algunos parámetros del programa. El primero y más obvio es que le tenemos que decir cuál es el archivo que contiene la alineación múltiple. En segundo lugar, tenemos que indicar si tenemos ácidos nucleicos o aminoácidos (también podemos dejar que el programa detecte el tipo de molécula). A continuación, el criterio que SMS va a utilizar para identificar el mejor modelo de evolución molecular, algunos trabajos sugieren que es mejor usar BIC (Luo et al. BMC Evol Biol, 2010;10:242). Después, hay que especificar el algoritmo de búsqeda de la topología más verosimil, NNI realiza una exploración más pequeña que SPR del "bosque de árboles". Para problemas difíciles se recomienda usar SPR, pero como nuestras secuencias no son muy divergentes, está bien NNI. Para el soporte estadístico de las ramas, utilizaremos SH-Like aLRT. Este es un método desarrollado específicamente para PhyML, es mucho más veloz que Bootstrap y está bien evaluado (Anisimova et al. 2011, Sys. Biol. 60(5):685-699). Finalmente podemos comenzar el análisis.

 

Una vez que el programa termina, aparece la pantalla que puedes ver en la Figura 4.2. En ella, encontrarás diversos resultados. Por ejemplo, en la opción "SMS Best Model" encontrarás el mejor modelo de evolución molecular, que en mi caso fue el GTR+G. En la opción "PhyML Statistics" encontrarás las estadísticas que caluculó PhyML. Es interesante verlas, pues encontrarás la verosimilitud de la filogenia, el valor de máxima parsimonia, las frecuencias de las bases y los valores del modelo de evolución molecular, entre otras cosas. En la opción "PhyML newick tree" puedes descargar la filogenia en formato newick o bien visualizarla. Para visualizar la filogenia hay dos opciones, la primera es usar el visor que viene integrado en la página, para ello da click en la opción "Viewer" (Figura 4.3). 

 

 

Figura 4.3. Fiologenia de los genes que codifican para la proteína S, inferida con PhyML+SMS. 

 

Una de las primeras cosas que hay que hacer, es colocar la raíz de la filogenia. En este caso, sabemos que las secuencias de coronavirus que infectan pangolin y de murciélago dieron origen a las de SARS-CoV-2. Entonces, una posibilidad es colocar la raíz en el nodo que separa a estos coronavirus del resto. Para colocar la raíz, seleccionas con el botón derecho del ratón el nodo de interés y en el menú que aparece, seleccionas "Reroot on this node" (Figura 4.4). El resultado lo puedes ver en la Figura 4.5.

 

 

Figura 4.4. Para colocar la raíz en una filogenia hay que seleccionar con el botón derecho del ratón el nodo de interés.

 

Figura 4.5. Resultado de colocar una nueva raíz en la filogenia.

 

Posteriormente, es recomendable mostrar los valores de soporte estadístico de los nodos internos (en este caso SH-Like aLRT) seleccionando la opción "Display support values" y alinear los nombres de los OTU para una mayor claridad seleccionando la opción "align text" (Figura 4.6).

 

 

Figura 4.6 Filogenia mostrando los soportes de las ramas y los OTU alineados.

 

En este caso, puedes ver que el SARS-CoV-2 evolucionó de los coronavirus que infectan murciélagos y pangolines. Esto es cierto al menos para los genes que codifican para la proteína S. Puedes leer más al respecto en este artículo: Li X, et al. Emergence of SARS-CoV-2 through Recombination and Strong Purifying Selection. Science Advances 29 May 2020.

 

 

Tutorial 1: Búsqueda de secuencias homólogas con BLAST

 

Una de las tareas más comunes en bioinformática y evolución molecular es buscar genes o proteínas homólogas. Para ello, los biólogos utilizamos comunmente el software BLAST. En este ejercicio mostraremos cómo podemos buscar secuencias homólogas con BLAST. Además de ello, veremos cómo podemos obtener mayor provecho de BLAST utilizando otras herramientas bioinformáticas complementarias. La metodología que mostramos en este tutorial la ulilizamos para publicar este artículo.

 

Durante los ejercicios, utilizaremos la fuente verdana para indicar que se trata de instrucciones que tienes que escribir en la ventana del sistema en Linux. Para las instrucciones en R para RStudio utilizaremos la fuente verdana. Además, el código de los scripts de Perl que utilizaremos te los mostraremos con la fuente Courier New. Estos scripts los puedes copiar y pegar en un editor de texto simple para usarlos. Sólo asegúrate de que la terminación de los scripts que copies y guardes con un editor de texto sea .pl y no .txt.

 

Paso 1 BLAST

Para hacer el tutorial utilizaremos la siguiente proteína:

 

>XP_018228376

MIKLDHFFFF VRLLCFIVAI IWIFLLPFDY FNRKVYISEN AILPGSVNMY FGGSDSNVLE

AYRDEIEFIN QQSKNVRFDS LHDIFRKLGL RTAIQDYKIQ FKGKNHSGTN FYAIFDTPRS

DGTEALVLSA AWKNMNGELN KGGVSLVLAL ARYFKRWSLW SKDIIFLIPE EKEVGVQVWV

DAYHKIDYEE GVSRLIMKGG EIQAVVDIDF VSEFREFDTI ELLYDGINGQ LSNLDLLNTV

NHIIQSKSGI KIKMQNVLQG KNSYFRKLLT MINGMISQCF GVLSGTQSCF IPYKIDAITL

KVHSKKNSEY DDILLGKIIE SLFRSLNNLL EHLHQSFFFY LMLSEKRFVS IATYLPSAIL

IASSFTIDAF SIWLSTFSNF YKNLSSSKEK LYSSKEKKNK LKYTGMDSEL YFKKQFILSI

STVIVIHIIA YIFLKILYNS VNYKGLFKDK KAYLMFLLFL QILQHIISIL IQKILKKIIK

SQNMSFFYNY RLNKVFSQAI LGVELSTLAI VNFSLSLFIG IVTFFLSFVR PAKTKFQQIL

ISIVLELINP LNWIFLFSYY KNKEVKDLMH LFAFGWKVWY LWTPTVLWLL WWPIWIMTKI

ILCMPIKDHN E

 

Esta proteína pertenece a un hongo ascomiceto llamado Pneumocystis jirovecii RU7. La proteína es parte de un complejo transamidasa que participa en la biosíntesis del GPI (glicosil-fostatidil-inositol). Puedes aprender un poco más de esta proteína aquí y también puedes leer por qué esta proteína es interesante desde el punto de vista de la evolución molecular aquí.

 

Lo primero que tienes que hacer es ingresar a la página de BLAST en el NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi) y seleccionar la opción Protein BLAST (Figura 1). Una vez ahí, copia y pega la secuencia de proteína en formato FASTA en el recuadro blanco (Figura 2) y ajusta los parámetros de BLAST para que obtengas un máximo de 5000 secuencias homólogas con un límite máximo de E-value de (Expected threshold < 0.01) (Figura 3). Después de ello, da click en el botón de BLAST.

 

 

 

Figura 1. Página web de BLAST en NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

 

 

 

Figura 2. Copia y pega la secuencia FASTA en el recuardo blanco.

 

 

 

Figura 3. Ajusta los parámetros de BLAST para que obtengas un máximo de 5000 secuencias (hits) y con un límite máximo de E-value < 0.01.

 

 

Una vez que tengas el resultado de BLAST, lo primero en lo que tienes que fijarte es en si nuestra proteína contiene uno o más dominios. En este caso, nuestra proteína contiene un solo dominio que se llama: 'Zinc_peptidase_like' (Figura 4). También deberás fijarte en el número de aminoácidos de nuestra proteína. En este caso la proteína contiene 611 aminoácidos.

 

 

 

Figura 4. La proteína XP_018228376 contiene el dominio Zinc_peptidase_like y 611 aminoácidos.

 

A continuación, selecciona todas las secuencias (Figura 5). En esta ocasión, BLAST encontró 2080 sequencias. Descarga los archivos FASTA (complete sequence) y Hit Table (text) utilizando el botón de Download (Figura 5).

 

 

 

Figura 5. Resultado de BLAST.

 

 

Con ello, descargarás los archivos seqdump.txt y X76GW1WSP014-Alignment.txt. Nota: en el caso del segundo archivo, los caracteres X76GW1WSP014 serán distintos en el archivo que descargaste. Es importante mencionar que todos los archivos que utilizaremos (incluidos los dos anteriores) se deberán de encontrar en una misma carpeta en tu computadora. Por ejemplo, puedes crear la carpeta TutorialBLAST/ y guardarlos ahí. También es recomendable que cambies los nombres de estos archivos por nombres más simples, por ejemplo (utilizando Linux):

 

$ mv seqdump.txt secuencias.fasta

$ mv X76GW1WSP014-Alignment.txt hit-table.txt 

 

Paso 2 Análisis de resultados

Ahora deberás de identificar las secuencias evolutivamente más cercanas del total de secuencias que encontró BLAST. Para ello, utilizarás el siguiente script de Perl (bitscoreandsize.pl) el cual extrae los 'bit score' del archivo hit-table.txt y también el porcentaje de aminoácidos alineados entre nuestra secuencia ('query') y cada una de las secuencias ('subject') que encontró en la base de datos. El script corre directamente en Linux:

 

$ perl bitscoreandsize.pl hit-table.txt 611 > bitscoreandsize.txt

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: bitscoreandsize.pl

# uso: 

# perl bitscoreandsize.pl hit-table.txt N > bitscoreandsize.txt

# N se refiere al numero de aminoacidos que contiene la secuencia query

 

# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left



use strict;


my $file = $ARGV[0];

my $size = $ARGV[1];

my $l = 0;

my $n = 0;


open (MIA, "$file") or die ("No puedoabrir $file\n");

while (my $linea = <MIA>){

        chomp ($linea);

        $l++;

        if ($linea !~ /^\#/ && $linea =~ /\w/){

                my @a = split (/\t/, $linea);

                $a[11] =~ s/\s//g;

                my $s = 100*(($a[7] - $a[6])/$size);

                print ("$a[6]\t$a[7]\t$a[11]\t$s\n");

                $n++;

        }

}

close (MIA);

------------------------------------------------------------------------------------------

 
El script anterior crea una tabla (bitscoreandsize.txt) con cuatro columnas. Las primeras dos columnas indican en que posición comienza y termina la alineación de nuestra secuencia query con cada una de las secuencias subject que encontró BLAST. La tercera columna es el bitscore de cada una de las alineaciones pareadas y la cuarta columna es el porcentaje de aminoácidos que alinean entre nuestra secuencia query y cada una de las secuencias subject.
 
Para analizar los resultados de BLAST utiizarás el siguiente script de R (en RStudio):
 

setwd('/en/donde/estan/los/archivos')

bitscore <- read.table('bitscoreandsize.txt')

summary(bitscore)

#Figura 6

hist(bitscore$V3, breaks = 100, xlab = 'bitscore', main = 'Histogram of bitscores')

segments(330,0,330,200, col = 'red', lty = 2)

# Figura 7

hist(bitscore$V4, breaks = 100, xlab = '% query coverage', main = 'Histogram of % query coverage')

segments(95,0,95,300, col = 'red', lty = 2)

# Figura 8

plot(bitscore$V4 ~ bitscore$V3, xlab = 'bitscore', ylab = '% query coverage', cex = 0.1, main = 'BLAST results')

segments(330,0,330,110, col = 'red', lty = 2)

segments(0,95,1300,95, col = 'red', lty = 2)

 
Debes de sustituir: /en/donde/estan/los/archivos por la dirección real en donde se encuentran los archivos de BLAST (por ejemplo en la carpeta TutorialBLAST/).

 

Como se puede apreciar en la Figura 6, existe un grupo de secuencias que son más similares a la secuencia query XP_018228376 que el resto. Estas secuencias tienen un 'bit-score' > 330. 

 

  

 

Figura 6. Histograma de los 'bit-score'. Nótese que aproximadamente en el valor 330 del eje de los bitscore la distribución sufre un punto de inflexión.

 

Además, si observas la distribución del porcentaje de aminoácidos alineados de la secuencia query con respecto a las secuencias subject, verás que un número importante de secuencias tienen un porcentaje de cobertura mayor al 95% (Figura 7).

 

 

 

Figura 7. Histograma del porcentaje de aminoácidos alineados de la secuencia query con respecto a las secuencias subject.

 

Finalmente, si graficas el bitscore contra el porcentaje de cobertura de cada alieación pareada, verás que estamos seleccionando a un grupo discreto de homólogos muy cercanos evolutivamente (Figura 8). 

 

 

 

Figura 8. Selección de homólogos cercanos con base en el bitscore y el % de cobertura. 

 

En este momento cabe que hagamos una reflexión sobre de la identificación de homólogos que hemos hecho. En este caso particular, hemos seleccionado los homólogos más similares a nuestra secuencia query utilizando como criterios el bitscore y el porcentaje de cobertura. Pero es probable que este no sea siempre el objetivo. Por ejemplo, si quisiéramos identificar todos los homólogos que hay en la base de datos a nuestra secuencia query, entonces es probable que debamos de seleccionar todas aquellas secuencias que tienen un e-value menor a 0.01 ó 0.001. Eso haría que seleccionemos a secuencias que muestran un porcentaje de cobertura (de la secuencia query con respecto a la secuencia subject) de un 30% aproximadamente. Aquí tendríamos que corroborar que dichas secuencias contengan también el dominio Zinc_peptidase_like (y para eso hay otros métodos como HMMER y bases de datos como Pfam). También podría ser que quisieramos identificar solo aquellos homólogos que tuviesen los mismos dominios que nuestra secuencia query. Por ejemplo, si nuestra secuencia query tuviera dos dominios, entonces nos interesaría que todos los homólogos tuvieran ambos dominios. En fin, el criterio que debemos usar (que tan estrictos queremos ser) dependerá de la pregunta que se tenga en cada caso. Eso sí, tenemos que tener en cuenta que si queremos a las secuencias para hacer una filogenia, entonces debemos de garantizar que todas ellas sean homólogas, es decir, que compartan un ancestro común.

 

Paso 3 Obtener las secuencias homólogas

Una vez que decidimos cuáles secuencias de proteínas vamos a utilizar para hacer nuestra filogenia (en este caso, todas aquellas que muestran un bitscore >= 330 y un porcentaje de aminoácidos alineados >= 95%), debemos de descargarlas del NCBI y darles un formato adecuado a los encabezados de las secuencias. Para ello usaremos una serie de scripts de Perl. Primero deberemos obtener las accesiones de las secuencias:

 

$ perl dame-acc-bitandsize.pl hit-table.txt 330 95 611 > accesiones.txt

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl 

 

# Script: dame-acc-bitandsize.pl

# uso: 

# perl dame-acc-bitandsize.pl hit-table.txt N1 N2 N3 > accesiones.txt

# N1 se refiere al bitscore minimo

# N2 se refiere al porcentaje de cobertura minimo

# N3 se refiere al numero de aminoacidos que contiene la secuencia query


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left


use strict;

my $file = $ARGV[0];

my $bits = $ARGV[1];

my $cvge = $ARGV[2];

my $size = $ARGV[3];

my $l = 0;

my $n = 0;

my $h = 0;

open (MIA, "$file") or die ("No puedoabrir $file\n");

while (my $linea = <MIA>){

        chomp ($linea);

        $l++;

        if ($linea !~ /^\#/ && $linea =~ /\w/){

                my @a = split (/\t/, $linea);

                $a[11] =~ s/\s+//g;

                $a[1] =~ s/\.\d+//;

                $h++;

                my $s = 100*(($a[7] - $a[6])/$size);

                if ($s >= $cvge && $a[11] >= $bits){

                        print ("$a[1]\n");

                }

        }

}

close (MIA);

------------------------------------------------------------------------------------------

 
De esta forma hemos seleccionado 412 accesiones:
 
$ wc -l accesiones.txt
$      412 accesiones.txt

 

Ahora vamos a utilizar un programa para obtener la clasificación taxonómica de las especies a las cuales pertenecen las 412 secuencias que seleccionamos. Para ello, requerimos también descargar algunos archivos de la base de datos Taxonomy del NCBI en nuestra computadora.

 

Para descargar los archivos de la base de datos de Taxnonomy necesitamos entrar con un explorador web a https://www.ncbi.nlm.nih.gov/taxonomy y ahí a Taxonomy FTP  (ftp://ftp.ncbi.nih.gov/pub/taxonomy/) y descarga el archivo taxdump.tar.gz. Luego entra a ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/ y ahí descarga el archivo prot.accession2taxid.gz. Si tienes un sistema macOS puedes utilizar la siguiente instrucción:

 

curl 'ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz' > taxdump.tar.gz

curl 'ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz' > prot.accession2taxid.gz

 

Y si tienes un sistema Linux lo puedes hacer así:

 

wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz

wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz

 

Coloca estos archivos en una carpeta distina (tal vez llamada TaxonomyDB/) y descomprímelos:

 

gunzip taxdump.tar.gz

tar -xf taxdump.tar

$ gunzip prot.accession2taxid.gz

 

También requerimos descargar un software para clasificar taxonómicamente las 412 secuencias. Ese programa lo puedes descargar de aquí: https://github.com/dib-lab/2018-ncbi-lineages. Descarga el programa, colócalo en un lugar adecuado en tu computadora (por ejemplo en una carpeta que se llame Software/) y descomprímelo:

 

unzip 2018-ncbi-lineages-master.zip

 

A continuación deberás copiar el archivo accesiones.txt a la carpeta 2018-ncbi-lineages-master/ que acabas de descomprimir. Ahí deberás ejecutar las siguientes instrucciones:

 

$ ./make-acc-taxid-mapping.py accesiones.txt /en/donde/esta/la/base/taxonomy/prot.accession2taxid.gz

 

$ ./make-lineage-csv.py /en/donde/esta/la/base/taxonomy/nodes.dmp /en/donde/esta/la/base/taxonomy/names.dmp accesiones.txt.taxid -o accesiones.lineages.txt

 

Estos programas gene los archivos: accesiones.lineages.txt, accesiones.txt.taxid. Es importante que los archivos que descargaste de Taxonomy NCBI sean tan recientes como la búsqueda de BLAST que realizaste. Ello debido a que los archivos de Taxonomy se actualizan constantemente. Si utilizas archivos de Taxonomy más antiguos, es probable que el software make-lineage-csv.py te marque errores y no encuentre la clasifiación taxonómica de todas las secuencias. 

 

Ahora deberás copiar el archivo: accesiones.lineages.txt a la carpeta en donde estas llevando a cabo este tutorial (por ejemplo TuturoialBLAST/).

 

A continuación extraemos las secuencias de proteínas en formato FASTA y a los encabezados de las secuencias les agregamos su clasificación taxonómica. Para ello vamos a usar un script que se llama dameNsecs.pl. A este programa le tenemos que dar el número mínimo de bitscore que usamos previamente (330), el tamaño de la proteína query (611), el porcentaje mínimo de aminoácidos alineados en la proteína query que usamos (95) y el nivel taxonómico que deseamos (ver más abajo):

 

$ perl dameNsecs.pl hit-table.txt secuencias.fasta accesiones.lineages.txt 330 611 95 4

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: dameNsecs.pl

# uso: 

# perl dameNsecs.pl hit-table.txt secuencias.fasta accesiones.lineages.txt N1 N2 N3 N4

# N1 se refiere al bitscore minimo

# N2 se refiere al numero de aminoacidos de la proteina query

# N3 se refiere al porcentaje minimo de aminoacidos alineados en la proteina query

# N4 se refiere al nivel taxonomico a utilizar para clasificar las secuencias


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left


use strict;


my $file = $ARGV[0];

my $secs = $ARGV[1];

my $taxa = $ARGV[2];

my $bits = $ARGV[3]; # N1

my $size = $ARGV[4]; # N2

my $cvge = $ARGV[5]; # N3

my $tlev = $ARGV[6]; # N4


my $l = 0;

my $n = 0;

my $h = 0;

my $m = 0;

my %secnames;

my %linaje;


open (MIA, "$file") or die ("No puedoabrir $file\n");

while (my $linea = <MIA>){

chomp ($linea);

$l++;

if ($linea !~ /^\#/ && $linea =~ /\w/){

my @a = split (/\t/, $linea);

$a[1] =~ s/\.\d+//;

$h++;

my $s = 100*(($a[7] - $a[6])/$size);

if ($s >= $cvge && $a[11] >= $bits){

#print ("$a[1]\n");

$secnames{$a[1]} = 1;

$n++;

}

}

}

close (MIA);

#print ("N: $n\n");


open (MIA, "$taxa") or die ("No puedoabrir $taxa\n");

while (my $linea = <MIA>){

chomp ($linea);

if ($linea =~ /\w/){

my @a = split (/,/, $linea);

if (exists $secnames{$a[0]}){

print ("$a[0]\t$linea\n");

$linaje{$a[0]} = $linea;

$m++;

}

}

}

close (MIA);

#print ("N: $m\n");


my $r = 0;

my $sec;

my $name;

my $spp;

$m = 0;


open (MIA, "$secs") or die ("No puedoabrir $secs\n");

open (ROB, ">outfile.txt") or die ("No puedoabrir outfile.txt\n");

while (my $linea = <MIA>){

chomp ($linea);

if ($linea =~ />/ && $r == 1){

if (exists $linaje{$name}){

my @a = split (/,/, $linaje{$name});

my $taxa = $a[$tlev];

my $namesec = ();

if ($taxa =~ /\w/){

$namesec = $name.'_'.$taxa.'_'.$spp;

} else {

$namesec = $name.'_Notidentified_'.$spp;

}

print ROB (">$namesec\n$sec\n");

$m++;

} else {

#print ("el accession $name no tiene linaje asignado\n");

}

$r = 0;

$name = $sec = ();

}

if ($linea !~ />/ && $r == 1){

$sec = $sec.$linea;

}

if ($linea =~ />/ && $r == 0){

if ($linea =~ />([\w\d_]+)\.[\d]\s/){

$name = $1;

} elsif ($linea =~ />([\w_]+)\s.+,\s/){

$name = $1;

} elsif ($linea =~ />pir\|(.+)\|\s/){

$name = $1;

} elsif ($linea =~ />prf\|\|(.+)\s/){

$name = $1;

} else {

print ("el patron 1 fallo:\n$linea\n");

$name = 'especie_no_identificada';

}

if ($linea =~ /\[(.+)\]/){

$spp = $1;

$spp =~ s/\s/_/g;

$spp =~ s/\.//g;

$spp =~ s/\-/_/g;

$spp =~ s/\'//g;

} elsif ($linea =~ />.+\[(.+)\]/){

$spp = $1;

} elsif ($linea =~ />.+,\s+(.+)/){

$spp = $1;

} else {

print ("el patron 2 fallo:\n$linea\n");

$spp = 'especie_no_identificada';

}

$r = 1;

}

}

if ($r == 1){

if (exists $linaje{$name}){

my @a = split (/,/, $linaje{$name});

my $namesec = $name.'_'.$a[$tlev].'_'.$spp;

print ROB (">$namesec\n$sec\n");

$m++;

} else {

#print ("el accession $name no tiene linaje asignado\n");

}

$r = 0;

$name = $sec = ();

}

close (ROB);

close (MIA);

print ("N: $m\n");

------------------------------------------------------------------------------------------

 

Este programa genera el archivo: outfile.txt

 

Ahora cambiamos el nombre del archivo de salida:

 

$ mv outfile.txt secuencias-330.fasta

 
El script dameNsecs.pl agrega la categoría taxonómica que solicitamos (en este caso la 4) al encabezado de la secuencia FASTA.
 
Por ejemplo en nuestra secuencia:
 
>XP_018228376_Pneumocystidomycetes_Pneumocystis_jirovecii_RU7
 
el script agregó la clase: Pneumocystidomycetes. De esta forma, el encabezado que tenemos está compuesto en primer lugar por: el identificador de la proteína (XP_018228376), en segundo lugar la clase (Pneumocystidomycetes) y en tercer lugar la cepa (Pneumocystis_jirovecii_RU7).
 

En este caso usamos el nivel taxonómico número 4 que corresponde al nivel de 'clase'. Si quisiéramos utilizar otro nivel taxonómico (por ejemplo, 'familia') hay que ver el contenido del archivo accesiones.lineages.txt e identificar el nivel deseado. Los niveles taxonómicos están separados por comas (,) y se enumeran los campos de izquierda a derecha comenzando con el número 0. Además, es importante mencionari que los niveles taxonómicos superiores están a la izquierda de cada renglón. Por ejemplo, en la Tabla 1 mostramos la clasificación taxonómica de Pneumocystis jirovecii. Ahí podremos ver qué otras categorías taxonómicas están a nuestra disposición.

 

Índice Clasificación taxonómica Nivel taxonómico
0 XP_018228376 Accesión
1 1408657 Taxonomy ID
2 Eukaryota Superreino
3 Ascomycota Phylum
4 Pneumocystidomycetes Clase
5 Pneumocystidales Orden
6 Pneumocystidaceae Familia
7 Pneumocystis Género
8 Pneumocystis jirovecii Especie
9 Pneumocystis jirovecii RU7 Cepa

Tabla 1. Clasificación de Pneumocystis jirovecii de acuerdo a NCBI. 

 

 

Tutorial 2: Alineación múltiple de proteínas

 

Paso 4 Alineación múltiple de proteínas

El objetivo de la alineación múltiple de proteínas es identificar los residuos que son homólogos entre las distintas secuencias. Existen varios programas para realizar alineaciones múltiples. En este tutorial vamos a utilizar mafft. Este programa lo puedes descargar de: https://mafft.cbrc.jp/alignment/software/.
 
Utilizaremos la opción --auto, la cual ajusta automáticamente los parámetros del programa dependiendo del número y tamaño de las secuencias a alinear:
 
$ mafft --auto secuencias-330.fasta > secuencias-330.aa.mafft.fasta
 
Una vez alineadas las secuencias, es importante que revises la alineación múltiple. Esto sirve para ver si no hay alguna secuencia que este mal alineada. Por ejemplo, si tuviéramos secuencias incompletas, estas introducen un gran número de indels (contracción de inserciones deleciones) o gaps
 
Para visualizar la alineación múltiple podemos utilizar un programa que se llama Jalview (http://www.jalview.org ). Una vez que hemos instalado Jalview podemos cargar nuestras secuencias en Jalview siguiendo la siguiente ruta en el menú: File -> Input Alignment -> From File. Luego, en el menú de la alineación múltiple, podrás colorear las columnas de acuerdo a su nivel de conservación con: Colour -> % identity. A continuación, podrás visualizar toda la alineación múltiple con: View -> Overview Window. 
 
 
 
Figura 9. Alineación múltiple de proteínas hecha con mafft. Cada renglón representa a una secuencia primaria de una proteína. Las proteínas están coloreadas en blanco y azul, entanto que en gris se muestran los indels o gaps. El color azul denota las columnas más conservadas entre las distintas proteínas alineadas. Las barras amarillas del inferior de la figura denotan el nivel de conservación de cada una de las columnas.
 
Ahora bien, es recomendable eliminar aquellas secuencias pequeñas que introducen muchos gaps en la alineación múltiple. Esto lo podemos hacer utilizando el programa trimAl (http://trimal.cgenomics.org).
 
A trimAl hay que proporcionarle dos parámetros -resoverlap y -seqoverlap. El primero de ellos, -resoverlap determina la fracción mínima de empalme que debe de presentar una posición determinada (un residuo) con respecto a otras posiciones (otros residuos) en la misma columna para que dicha posición pueda ser considerada una 'buena posición'. El segundo parámetro, -seqoverlap determina el porcentaje mínimo de 'buenas posiciones' que debe de contener una secuencia para que no sea eliminada de la alineación múltiple. De esta forma, si queremos eliminar a las secuencias que no tengan al menos un 90% de sus residuos empalmados en columnas que presenten al menos un 0.7 de empalme: 
 
$ trimalx -in secuencias-330.aa.mafft.fasta -out secuencias-330.aa.mafft.trim.fasta -resoverlap 0.7 -seqoverlap 90
 
Podemos ver que con trimalx eliminamos 90 secuencias:
 
$ grep -c '>' secuencias-340.aa.mafft 
$ 412
$ grep -c '>' secuencias-340.aa.trim.mafft
$ 322
 
Pero ¿cómo sabemos que valores debemos darle a -resoverlap y -seqoverlap? Una posibilidad es hacer un barrido de distintos valores y graficarlo. Para ello, utilizaremos el siguiente script el Perl:
 
$ perl analizaalineamiento.pl secuencias-330.aa.mafft.fasta
 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: analizaalineamiento.pl

# uso: 

# perl analizaalineamiento.pl secuencias-330.aa.mafft.fasta

# out: outfile.txt


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left

 

# Nota: para la version de macOS deberas sustituir trimal por trimalx en la linea 23. La version de linux corre como trimal.


use strict;


my $file = $ARGV[0];


for (my $i = 100; $i > 0; $i = $i-10){

my $outfile = 'trim.'.$i.'.txt';

open (ROB, ">$outfile") or die ("No puedo abrir $outfile\n");

for (my $j = 1; $j >= 0; $j = $j-0.1){

system ("trimalx -in $file -out borrame.fasta -resoverlap $j -seqoverlap $i");

if ($j >= 0.1){

print     ("$i\t$j\t");

print ROB ("$i\t$j\t");

} else {

print     ("$i\t0\t");

print ROB ("$i\t0\t");

}

system ("grep -c '>' borrame.fasta");

system ("grep -c '>' borrame.fasta > borrame");

open (MIA, "borrame") or die ("No puedo abrir borrame\n");

while (my $linea = <MIA>){

print ROB ("$linea");

}

close (MIA);

system ("rm borrame.fasta");

}

close (ROB);

}

 

------------------------------------------------------------------------------------------

 
Y después, usaremos el siguiente script en R:
 
setwd('/en/donde/estan/los/archivos')
trim100 <- read.table('trim.100.txt')
trim90 <- read.table('trim.90.txt')
trim80 <- read.table('trim.80.txt')
plot(trim100$V3 ~ trim100$V2, col = 'darkorange4', type = 'b', lty = 4, lwd = 2, xlab = 'resoverlap', ylab = 'No of sequences', main = 'TrimAl analysis', ylim = c(0,450) )
points(trim90$V3 ~ trim90$V2, col = 'darkorange3', type = 'b', lty = 4, lwd = 2)
points(trim80$V3 ~ trim80$V2, col = 'darkorange', type = 'b', lty = 4, lwd = 2)
legend(0.85, 350, legend = c('100', '90', '80'), col = c("darkorange4", "darkorange3", "darkorange"), lty = 1, cex = 0.8, lwd = 2, title = 'seqoverlap')
 
Debes de sustituir: /en/donde/estan/los/archivos por la dirección real en donde se encuentra la alineación múltiple (por ejemplo en la carpeta TutorialBLAST/). Reconozco que hay mejores formas de programar en R, pero por lo pronto, esta funciona.
 
Podrás ver que para distintos valores de -resoverlap, utilizar -seqoverlap 100 resulta ser demasiado estricto, en tanto que -seqoverlap 80 es muy laxo (Figura 10) . En cambio, la combinación -resoverlap 0.7 -seqoverlap 90 parece razonable. 
 
 
 
Figura 10. Número de secuencias retenidas para distintos valores de -resoverlap y -seqoverlap utilizando el software trimAl.
 
Si revisamos la alineación múltiple después de haber utilizado trimAl (-resoverlap 0.7 -seqoverlap 90 ) con Jalview, encontraremos que el número de columnas coloreadas de azul (color que indica el nivel de conservación) se incrementó (Figura 11).
 
 
Figura 11. El porcentaje de columnas azules se incrementó después de utilizar trimAl con los parámetros -resoverlap 0.7 y -resoverlap 90.
 
De cualquier forma, es buena práctica revisar nuestra alineación múltiple a ojo antes y después de aplicar trimAl. El objetivo es maximizar el número de posiciones sin indels o gaps manteniendo el mayor número de secuencias posibles.
 
En este caso hemos utilizado el software trimAl, pero es válido también eliminar secuencias manualmente. Algunas veces una sola secuencia fragmentaria (por que se secuenció solo un fragmento del gen) introduce una gran cantidad de gaps. Esta secuencia la podemos detectar fácilmente con Jalview y la deberemos de eliminar manualmente. 
 
En esta ocasión hemos usado trimAl para eliminar sequencias que no alinean suficientemente bien con el resto de las secuencias en la alineación múltiple. Pero trimAl también puede eliminar columnas individuales de la alineación múltiple. Sin embargo, a la fecha no es recomendable (en términos generales) eliminar columnas utilizando trimAl u otros programas similares. Para saber más de este tema, te recomendamos leer el artículo publicado por Tan et al (2015)
 
Después de eliminar secuencias con trimAl es bueno revisar que nuestra secuencia query siga en la alineación múltiple:
 
$ grep 'Pneumocystis_jirovecii' secuencias-330.aa.mafft.trim.fasta 
>XP_018228376_Pneumocystidomycetes_Pneumocystis_jirovecii_RU7 1490 bp
 
Finalmente es importante mencionar que al momento de eliminar secuencias no nos podemos olvidar de la pregunta biológica que queremos responder. Por ejemplo, si queremos estudiar la evolución de una secuencia fragmentaria, no la deberemos de eliminar de la alineación múltiple (evidentemente) por muchos gaps que esta introduzca. En fin, lo que queremos decir es que la pregunta biológica deberá de guiar el análisis bioinformático.
 
 

Tutorial 3: Infiere una filogenia utilizando una alineación múltiple de proteínas

 
Existen distintos métodos para inferir filogenias. Los más populares son: máxima verosimilitud, inferencia bayesiana, métodos basados en matrices de distancia y métodos de máxima parsimonia. De cualquier forma, la inferencia filogenética es un tema amplio. Un buen resumen de los métodos actuales lo puedes encontrar aquí: Yang and Rannala, (2012). También te recomendamos leer un buen libro introductorio: Page and Holmes (1991). En este tutorial utilizaremos el método de máxima verosimilitud implementado en el programa PhyML el cual puedes descargar de (http://www.atgc-montpellier.fr/phyml/). 
 
Una vez que contamos con la alineación múltiple de proteínas, el siguiente paso para realizar una filogenia es identificar el modelo de evolución molecular que mejor se ajusta a nuestros datos (en el caso de los análisis que utilizan máxima verosimilitud). Sin embargo, antes de proceder a hacer el análisis filogenético, reduciremos el número de secuencias que tenemos en la alineación múltiple por fines didácticos . Lo hacemos así para que los análisis siguientes no ocupen demasiado tiempo de cómputo.
 
Usaremos un script de Perl que selecciona N=5 secuencias de cada clase taxonómica:
 
$ perl seleccionaconbaseentaxa.pl secuencias-330.aa.mafft.trim.fasta 5
 
------------------------------------------------------------------------------------------
#!/usr/bin/perl

# Script: seleccionaconbaseentaxa.pl
# uso: 
# perl seleccionaconbaseentaxa.pl secuencias-330.aa.mafft.trim.fasta N
# out: outfile.txt
# N representa el numero maximo de secuencias de cada taxa

# Luis Jose Delaye Arredondo
# luis.delaye@cinvestav.mx
# http://www.ira.cinvestav.mx/evolutionary.genomics
# Copy left

use strict;

my $file = $ARGV[0];
my $max  = $ARGV[1];
my $r = 0;
my $sec;
my $name;

my %hash;
my %sele;

open (MIA, "$file") or die ("No puedo abrir $file\n");
while (my $linea = <MIA>){
#chomp ($linea);
if ($linea =~ />/ && $r == 1){
$hash{$name} = $sec;
$r = 0;
$name = $sec = ();
}
if ($linea !~ />/ && $r == 1){
$sec = $sec.$linea;
}
if ($linea =~ />/ && $r == 0){
if ($linea =~ /(>.+)/){
$name = $1;
} else {
die ("el patron 1 fallo: $linea\n");
}
$r = 1;
}
}
if ($r == 1){
$hash{$name} = $sec;
$r = 0;
$name = $sec = ();
}
close (MIA);

my @a = keys (%hash);
open (ROB, ">outfile") or die ("No puedo abrir outfile\n");
for (my $i = 0; $i <= $#a; $i++){
my $taxa = ();
if ($a[$i] =~ />[A-Z]{2}_\d+_([A-Za-z]+)/){
$taxa = $1;
} elsif ($a[$i] =~ />[A-Z\d]+_([A-Za-z]+)/){
$taxa = $1;
} else {
die ("el patron 2 fallo: $a[$i]\n");
}
if (!exists $sele{$taxa}){
$sele{$taxa} += 1;
print ROB ("$a[$i]\n$hash{$a[$i]}\n");
} elsif ($sele{$taxa} < $max){
$sele{$taxa} += 1;
print ROB ("$a[$i]\n$hash{$a[$i]}\n");
}
}
close (ROB);
------------------------------------------------------------------------------------------
 
La salida de este script es un archivo outfile, por lo que hay que cambiarle el nombre:
 
$ mv outfile secuencias-330.aa.mafft.trim.sel.fasta
 
Podemos ver que hemos seleccionado 40 secuencias (máximo 5 de cada clase taxonómica):
 
$ grep -c '>' secuencias-330.aa.mafft.trim.sel.fasta
     40
 
Es importante mencionar que el script seleccionaconbaseentaxa.pl está hecho para seleccionar distintas secuencias cada vez. Así que es casí seguro que las secuencias que tú selecciones sean distintas a las que aquí hemos seleccionado.
 
El programa que vamos a utilizar, PhyML, requiere que las secuencias estén en formato PHYLIP. Una de las particularidades del formato PHYLIP es que no acepta más de 10 caracteres en el encabezado de las secuencias. Por lo que va a ser necesario codificar los nombres de las secuencias para poder trabajar. Vamos a usar el siguiente script para codificar los nombres de las secuencias:
 
$ perl codifica_spp.pl secuencias-330.aa.mafft.trim.sel.fasta
 
------------------------------------------------------------------------------------------
#!/usr/bin/perl

# Script: codifica_spp.pl
# uso: 
# perl codifica_spp.pl secuencias-330.aa.mafft.trim.sel.fasta
# out: outfile.fasta, outfile.txt
# outfile.fasta contiene las secuencias en formato fasta con los encabezados codificados
# outfile.txt contiene la codificacion de los encabezados

# Luis Jose Delaye Arredondo
# luis.delaye@cinvestav.mx
# http://www.ira.cinvestav.mx/evolutionary.genomics
# Copy left

use strict;

my $file = $ARGV[0];

my $c = 0;

open (MIA, "$file") or die ("No puedo abrir $file\n");
open (ROB, ">outfile.fasta") or die ("No puedo abrir outfile.fasta\n");
open (SOL, ">outfile.txt") or die ("No puedo abrir outfile.txt\n");
while (my $linea = <MIA>){
chomp ($linea);
if ($linea =~ />/){
$c++;
my $spp = 'sp_'.$c;
print ROB (">$spp\n");
$linea =~ s/>//;
print SOL ("$spp\t$linea\n");
} else {
print ROB ("$linea\n");
}
}
close (SOL);
close (ROB);
close (MIA);
------------------------------------------------------------------------------------------
 
Este script crea dos archivos: outfile.fasta, outfile.txt. El primero de ellos contiene las secuencias con los encabezados codificados y el segundo de ellos contiene la codificación. Es necesario cambiarles el nombre a estos dos archivos:
 
$ mv outfile.fasta secuencias-330.aa.mafft.trim.sel.cod.fasta
$ mv outfile.txt secuencias-330.aa.mafft.trim.sel.cod.txt
 
Paso 5 Inferencia filogenética
Ahora bien, para identificar el mejor modelo de evolución molecular utilizaremos el software sms (Smart Model Selection) el cual puedes descargar de aquí (http://www.atgc-montpellier.fr/sms/binaries.php). Sin embargo, antes de poder usar sms, debemos de cambiar el formato de las secuencias que tenemos en la alineación múltiple. Nuestras secuencias están actualmente en el formato FASTA, sin embargo sms requiere del formato PHYLIP. El cambio de formato lo podemos hacer con el software readAl, el cual viene junto con trimAl (http://trimal.cgenomics.org).
 
$ readal -in secuencias-330.aa.mafft.trim.sel.cod.fasta -out secuencias-330.aa.mafft.trim.sel.cod.phy -phylip
 

Posteriormente utilizamos sms:

 

$ sms.sh -i secuencias-330.aa.mafft.trim.sel.cod.phy -d aa

 

 El resultado de sms es el siguiente:

 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Selected model : LG +G

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Substitution model : LG

Equilibrium frequencies : Model

Proportion of invariable sites : fixed (0.0)

Number of substitution rate categories : 4

Gamma shape parameter : estimated (0.701)

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 

El programa sms nos dice que el modelo de evolución que mejor se ajusta a nuestra alineación múltiple (usando el criterio de AIC) es el LG+G. Nos dice también que la proporción de sitios invariables es de 0 y que el parámetro alfa de la distribución gamma es de 0.701. También nos dice que las frequencias en equilibrio de los aminoácidos las deberemos de estimar con el modelo de evolución. Con esta información podemos inferir nuestra filogenia con PhyML: 

 

$ PhyML-3.1_macOS-MountainLion -i secuencias-330.aa.mafft.trim.sel.cod.phy -d aa -m LG -v 0 -a 0.701 -c 4 -f m -s BEST -b -4

 

La descripción de los parámetros que utilizamos es la siguiente: para indicar que la alineación múltiple es de aminoácidos usamos -d aa; el modelo lo especificamos con -m LG; para la proporción de sitios invariables igual a 0 usamos -v 0; para especificar el valor alfa de la distribución gamma -a 0.701; el número de categorías para la distribución gamma es -c 4; para estimar las frecuencias de los aminoácidos de acuerdo al modelo usamos -f m; para la búsqueda heurística de la filogenia más verosimil utilizamos -s BEST; y finalmente utilizamos -b -4 para evaluar estadísticamente los nodos internos de la filogenia con una aLRT-SH. Para entender a fondo cómo PhyML va a estimar tu filogenia, deberas leer los artículos y el manual de PhyML. 

 

Una vez que PhyML termina, deberemos de cambiar el nombre a los archivos de salida para indicar que hemos utilizado la opción -b -4 para inferir el soporte de los nodos internos:

 

$ mv secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_stats.txt secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_stats-SH.txt 

 

$ mv secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_tree.txt secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_tree-SH.txt

 

A continuación vamos a descodificar los nombres de las secuencias en la filogenia:

 

$ perl descodifica_spp.pl secuencias-330.aa.mafft.trim.sel.cod.txt secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_tree-SH.txt > secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_tree-SH.dec.tre

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: codifica_spp.pl

# uso: 

# perl descodifica_spp.pl secuencias-330.aa.mafft.trim.sel.cod.txt secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_tree-SH.txt > secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_tree-SH.dec.tre


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left


use strict;


my $file1 = $ARGV[0];

my $file2 = $ARGV[1];

my $file3 = $ARGV[2];


my $name;

my $sp;

my %spp;


open (MIA, "$file1") or die ("No puedo abrir $file1\n");

while (my $linea = <MIA>){

chomp ($linea);

if ($linea =~ /sp_/){

if ($linea =~ /(sp_\d+)\t(.+)/){

#print ("-> $1\n");

$sp = $1;

$name = $2;

$spp{$sp} = $name;

} else {

die ("el patron 1 fallo:\n$linea\n");

}

}

}

close (MIA);


my %acc;


my $tree;

open (MIA, "$file2") or die ("No puedo abrir $file2\n");

while (my $linea = <MIA>){

chomp ($linea);

$tree = $linea;

}


my @keys = sort keys (%spp);

for (my $i = 0; $i <= $#keys; $i++){

$tree =~ s/$keys[$i]:/$spp{$keys[$i]}:/;

}

 

print ("$tree\n");

------------------------------------------------------------------------------------------

 

Luego, para visualizar la filogenia tenemos varias opciones. Una posibilidad es utilizar ETE Toolkit (http://etetoolkit.org/ ). Para ello tendrás que instalar el software siguiendo las instrucciones que vienen en la página web de ETE Toolkit. Una vez que lo tengas listo, puedes visualizar la filogenia utilizando el siguiente script de Python:

 

$ python showtree1.py secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_tree-SH.dec.tre

 

------------------------------------------------------------------------------------------

#!/usr/bin/env python


from ete3 import Tree, NodeStyle, TreeStyle

import sys

import os


# Archivo que contiene la filogenia en formato newick

t1 = Tree(sys.argv[1])

# Colocamos la raiz en el punto medio

R1 = t1.get_midpoint_outgroup()

t1.set_outgroup(R1)

# Muestra la filogenia

print (t1)

 

t1.show()

------------------------------------------------------------------------------------------

 

 El resultado se puede ver en la Figura 12.

 

 

 

Figura 12. Inferencia filogenética por máxima verosimilitud. Los números en los nodos indican el nivel de soporte de acuerdo al estadístico aLRT-SH. La raíz se colocó utilizando el método del punto medio. Nótese que las distintas especies de hongos agrupadas en una misma clase forman grupos monofiléticos.

 

Paso 6 Identificación de la raíz de la filogenia

Recientemente el grupo de una Tal Dagan publicó un método nuevo para identificar la raíz en una filogenia: MAD, Phylogenetic rooting using minimal ancestor deviation (Tria et al. 2017). Como el método se ve prometedor, lo utilizaremos en este tutorial. Puedes descargar el programa aquí (https://www.mikrobio.uni-kiel.de/de/ag-dagan/ressourcen ). Coloca el archivo zip en algún lugar en tu computadora (por ejemplo en Software/) y descomprímelo. Después de ello usa la siguiente instrucción dentro de la carpeta que descomprimiste:

 

$ mad.py secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_tree-SH.dec.tree

 

El archivo de salida tendrá el mismo nombre que el de entrada pero con la extensión 'rooted' añadida. A continuación puedes visualizar nuevamente la filogenia (Figura 13).

 

$ python showtree1.py secuencias-330.aa.mafft.trim.sel.cod.phy_phyml_tree-SH.dec.tre.rooted

 

 

 

Figura 13. Inferencia filogenética por máxima verosimilitud. La raíz se colocó utilizando el método MAD. Nótese que en este caso, la raíz inferida por MAD es la misma que la inferida por el método del punto medio.

 

Aunque en este caso la raíz que infiere MAD es la misma que la que se encontró con el método del punto medio, recomendamos usar MAD siempre que sea posible.

 

 

Tutorial 4: Infiere una filogenia utilizando una alineación múltiple de genes

 

Paso 7 Obtención de genes codificantes

En algunas ocasiones necesitamos tener las secuencias de los genes codificantes ya sea para hacer análisis filogenéticos o para hacer análisis de selección natural, entre otros. Genbank no proporciona fácilmente esta información. Sin embargo, con algunos trucos bioinformáticos es posible descargar de Genbank los genes codificantes a partir de los identificadores de sus proteínas.

 

Lo primero que tenemos que hacer es obtener los identificadores de las proteínas. Esto lo hacemos con el siguiente script:

 

$ perl dameprotacc.pl secuencias-330.aa.mafft.trim.sel.fasta > secuencias-330.aa.mafft.trim.accesiones.txt 

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: dameprotacc.pl

# uso: 

# perl dameprotacc.pl secuencias-330.aa.mafft.trim.sel.fasta > secuencias-330.aa.trim.mafft.accesiones.txt


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left


use strict;


my $file = $ARGV[0];

open (MIA, "$file") or die ("No puedo abrir $file\n");

while (my $linea = <MIA>){

chomp ($linea);

if ($linea =~ />/){

if ($linea =~ />(XP_\d+)_/){

print ("$1\n");

} elsif ($linea =~ />([A-Z0-9]+)_/){

print ("$1\n");

} else {

die ("el patron 1 fallo: $linea\n");

}

}

}

close (MIA);

------------------------------------------------------------------------------------------

 

A continuación obtenemos los archivos Genbank de las proteínas correspondientes (ojo, este script va a generar un archivo Genbank por cada proteína, por lo que tal vez es buena idea direccionar la salida del programa a una carpeta aparte):

 

$ mkdir Genbank_aa

$ perl dameprotgb.pl secuencias-330.aa.mafft.trim.accesiones.txt

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: dameprotgb.pl

# uso: 

# perl dameprotgb.pl secuencias-330.aa.mafft.trim.accesiones.txt


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left


use strict;

my $file = $ARGV[0];

my $n = 0;

open (MIA, "$file") or die ("No puedo abrir $file\n");

while (my $linea = <MIA>){

chomp ($linea);

if ($linea =~ /(\w+)/){

my $acc = $1;

print ("$acc\n");

my $outfile = $acc.'.aa.gbff';

system ("curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=protein&rettype=gb&retmode=text&id=$acc' > Genbank_aa/$outfile");

# system ("wget https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=protein&rettype=gb&retmode=text&id=$acc > Genbank_aa/$outfile");

sleep(20);

$n++;

print ("archivo numero: $n\n");

}

}

close (MIA);

------------------------------------------------------------------------------------------

 

Nota. En este momento el script funciona solo en el sistema operativo macOS. Para usarlo en Linux hay que comentar (con un #) la instrucción 'curl' y quitar el # de la instrucción 'wget'. Ojo, cuando copies y pegues este script en tu computadora, revisa que las instrucciones que comienzan por 'system' se encuentren en una sola línea .

 

Ahora obtenemos los identificadores de los genes correspondientes en los archivos Genbank que acabamos de descargar, también obtenemos los 'locus_tag'. Para ello, deberás de instalar BioPerl (https://bioperl.org/) en tu computadora, pues el siguiente script lo requiere:

 

$ perl damenuctags.pl

$ mv outfile.txt secuencias-330.aa.mafft.trim.accesiones.nuc.txt

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: damenuctags.pl

# uso: 

# perl damenuctags.pl 


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left


use strict;

use Bio::SeqIO;


my @files = glob ("Genbank_aa/*.aa.gbff");


open (ROB, ">outfile.txt") or die ("No puedo abrir outfile.txt\n");

print ROB ("LOCUS_nu\tlocus_tag\tLOCUS_aa\n");

for (my $i = 0; $i <= $#files; $i++){

print ("$files[$i]\n");

my @a = ();

@a = &extraegb($files[$i]);

print ("\n-----\n");

print ("$a[0]\t$a[1]\t$a[2]");

print ("\n-----\n");

if ($a[0] =~/not found/){

print ("No encontro el nombre del archivo Genbank\n");

#my $pausa = <STDIN>;

}

if ($a[1] =~/not found/){

print ("No encontro el nombre del locus_tag\n");

#my $pausa = <STDIN>;

}

if ($a[2] =~/not found/){

print ("No encontro el nombre del protein_id\n");

#my $pausa = <STDIN>;

}

print ROB ("$a[0]\t$a[1]\t$a[2]\n");

}

close (ROB);


sub extraegb {

my $fileGB = $_[0];

my $Gbfile    = 'not found';

my $locus_tag = 'not found';

my $proteinf  = 'not found';

my $n = 0;

# Creo un objeto

my $seqio_object = Bio::SeqIO->new(-file => $fileGB);

#my $pausa = <STDIN>;

# Obtengo cada secuencia dentro del objeto

while (my $seq_object = $seqio_object->next_seq){

my $NC = $seq_object->accession_number;

##print ("\n-----\n");

$proteinf = $NC;

print "$NC\t";

print "length ", $seq_object->length,"\n";

my $anno_collection = $seq_object->annotation;

for my $key ( $anno_collection->get_all_annotation_keys ) {

my @annotations = $anno_collection->get_Annotations($key);

for my $value ( @annotations ) {

print $value->tagname, "\t";

print (" -> ");

print $value->as_text, "\n";

if ($value->as_text =~ /Direct database link to (.+) in database GenBank/){

$Gbfile = $1;

} elsif ($value->as_text =~ /Direct database link to (.+) in database embl/){

$Gbfile = $1;

} elsif ($value->as_text =~ /Direct database link to (.+) in database REFSEQ/){

$Gbfile = $1;

}

}

}

# Obtengo los caracteres de cada secuencia

print ("FEATURES\n");

for my $feat_object ($seq_object->get_SeqFeatures) {  

print "\t",$feat_object->primary_tag, "\n";

my $sequence_string = $feat_object->seq->seq;

##print $sequence_string, "\n";

# Obtengo las notas de cada caracter

for my $tag ($feat_object->get_all_tags) {             

print "\t  tag: ", $tag, "\n";             

for my $value ($feat_object->get_tag_values($tag)) {                

print "\t    value: ", $value, "\n";             

}          

my $GI          = 'not found';

my $PID         = 'not found';

my $gene        = 'not found';

my $product     = 'not found';

my $translation = 'not found';

if ($feat_object->primary_tag eq "CDS"){

$n++;

# Obtengo las notas de cada caracter

for my $tag ($feat_object->get_all_tags) {

if ($tag eq 'db_xref'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /GI:(\d+)/){

$GI = $1;

}

}

} elsif ($tag eq 'protein_id'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$PID = $1;

}

}

} elsif ($tag eq 'gene'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$gene = $1;

}

}

} elsif ($tag eq 'product'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$product = $1;

}

}

} elsif ($tag eq 'translation'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$translation = $1;

}

}

} elsif ($tag eq 'locus_tag'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$locus_tag = $1;

}

}

}

}

#if ($PID ne 'not found' && $GI ne 'not found'){

# print ROB (">GI:$GI $PID\n$sequence_string\n");

# print BEA (">GI:$GI $PID\n$translation\n");

# print MAR (">$PID\t$GI\t$gene\t$product\n");

#}

}

}

#my $pausa = <STDIN>;

}

if ($n > 1){

print ("hay mas de 1 CDS\n");

my $pausa = <STDIN>;

}

return ($Gbfile, $locus_tag, $proteinf);

}

------------------------------------------------------------------------------------------

 

 Es importante que revises si el script encontró los identificadores de todas las secuencias:

 

$ grep 'not found' secuencias-330.aa.mafft.trim.accesiones.nuc.txt

   XM_024046157.1    not found     XP_023901925 

   XM_024071215.1    not found     XP_023926983  

 

En este caso no encontró los 'locus_tag' de dos genes. Esta situación no nos afecta pues el siguiente script utiliza la primer columna para encontrar los archivos Genbank de los genes.

 

A continuación obtenemos los archivos Genbank de los genes:

 

$ mkdir Genbank_nu

$ perl damenucgb.pl secuencias-330.aa.mafft.trim.accesiones.nuc.txt

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: damenucgb.pl

# uso: 

# perl damenucgb.pl secuencias-330.aa.mafft.trim.accesiones.nuc.txt


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left


use strict;

my $file = $ARGV[0];

my $n = 0;

my $l = 0;

open (MIA, "$file") or die ("No puedo abrir $file\n");

while (my $linea = <MIA>){

chomp ($linea);

$l++;

if ($linea =~ /\w+/ && $l > 1){

my @a = split (/\t/, $linea);

print ("$a[0]\n");

my $outfile = $a[0].'.nu.gbff';

system ("curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=gb&retmode=text&id=$a[0]' > Genbank_nu/$outfile");

# system ("wget https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=gb&retmode=text&id=$a[0]' > Genbank_nu/$outfile");

sleep(20);

$n++;

print ("listo: $n\n");

#my $pausa = <STDIN>;

}

}

close (MIA);

------------------------------------------------------------------------------------------

 

Es buena idea revisar si el script encontró y descargó todos los archivos Genbank. Para ello, puedes usar la siguiente instrucción:

 

$ ls -l Genbank_nu/

 

Si encuentras que algunos de los archivos están vacíos, puedes volver a correr el script anterior. Para hacerlo más eficiente, haz una copia del archivo secuencias-330.aa.mafft.trim.accesiones.nuc.txt y edíta la copia para dejar solo los renglones que corresponden con los archivos que hicieron falta. Por supuesto, ejecuta el script con la copia que editaste.

 

A continuación obtenemos la secuencia FASTA de los genes. También obtendremos un archivo extra con los nombres (locus_tag) de los genes para los cuales no pudimos obtener sus secuencias FASTA. Este script también requere BioPerl.

 

$ perl damegenes-CDS.pl secuencias-330.aa.mafft.trim.accesiones.nuc.txt

$ mv outfile.fasta secuencias-330.aa.mafft.trim.CDS.fasta

$ mv outfile.txt secuencias-330.aa.mafft.trim.faltan.txt

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: damegenes-CDS.pl

# uso: 

# perl damegenes-CDS.pl secuencias-330.aa.mafft.trim.accesiones.nuc.txt


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left


use strict;


use Bio::SeqIO;


my $file = $ARGV[0];

my $l = 0;


open (MIA, "$file") or die ("No puedo abrir $file\n");

open (ROB, ">outfile.fasta") or die ("No puedo abrir outfile.fasta\n");

open (SOL, ">outfile.txt") or die ("No puedo abrir outfile.txt\n");

while (my $linea = <MIA>){

chomp ($linea);

my @a = split (/\t/, $linea);

$l++;

if ($linea =~ /\w/ && $l > 1){

my $gbfile = $a[0].'.nu.gbff';

my $locus_tag_E = $a[1];

my $pID = $a[2];

my @b = &extraegb($gbfile, $locus_tag_E, $pID);

if ($b[0] =~ /\w/ && $b[1] =~ /\w/ && $b[2] =~ /\w/){

print ROB (">$b[0] $b[1]\n$b[2]\n");

} else {

print SOL ("$gbfile, $locus_tag_E, $pID\n");

}

}

}

close (SOL);

close (ROB);

close (MIA);


sub extraegb {

my $gbfile = $_[0];

my $locus_tag_E = $_[1];

my $pID = $_[2];

my $n = 0;

# Creo un objeto

my $seqio_object = Bio::SeqIO->new(-file => 'Genbank_nu/'.$gbfile);

#my $pausa = <STDIN>;

# Obtengo cada secuencia dentro del objeto

while (my $seq_object = $seqio_object->next_seq){

my $NC = $seq_object->accession_number;

##print ("\n-----\n");

print "$NC\t";

print "length ", $seq_object->length,"\n";

my $anno_collection = $seq_object->annotation;

for my $key ( $anno_collection->get_all_annotation_keys ) {

my @annotations = $anno_collection->get_Annotations($key);

for my $value ( @annotations ) {

print $value->tagname, "\t";

print $value->as_text, "\n";

}

}

# Obtengo los caracteres de cada secuencia

print ("FEATURES\n");

for my $feat_object ($seq_object->get_SeqFeatures) {  

##print "\t",$feat_object->primary_tag, "\n";

my $sequence_string = $feat_object->spliced_seq->seq;

##my $sequence_string = $feat_object->seq->seq;

##print $sequence_string, "\n";

# Obtengo las notas de cada caracter

for my $tag ($feat_object->get_all_tags) {             

##print "\t  tag: ", $tag, "\n";             

for my $value ($feat_object->get_tag_values($tag)) {                

##print "\t    value: ", $value, "\n";             

}          

my $GI          = 'not found';

my $PID         = 'not found';

my $gene        = 'not found';

my $product     = 'not found';

my $translation = 'not found';

my $locus_tag   = 'not found';

if ($feat_object->primary_tag eq "CDS"){

$n++;

# Obtengo las notas de cada caracter

for my $tag ($feat_object->get_all_tags) {

if ($tag eq 'db_xref'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /GI:(\d+)/){

$GI = $1;

}

}

} elsif ($tag eq 'protein_id'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$PID = $1;

}

}

} elsif ($tag eq 'gene'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$gene = $1;

}

}

} elsif ($tag eq 'product'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$product = $1;

}

}

} elsif ($tag eq 'translation'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$translation = $1;

}

}

} elsif ($tag eq 'locus_tag'){

for my $value ($feat_object->get_tag_values($tag)) {

if ($value =~ /(.+)/){

$locus_tag = $1;

}

}

}

}

if ($locus_tag_E eq $locus_tag){

print (">$pID $locus_tag\n$sequence_string\n");

#my $pausa = <STDIN>;

return ($pID, $locus_tag_E, $sequence_string);

}

}

}

#my $pausa = <STDIN>;

}

}

------------------------------------------------------------------------------------------

 
 Como se puede apreciar, el script no encontró todos los genes (CDS) de todas las proteínas:
 
$ grep -c '>' secuencias-330.aa.mafft.trim.CDS.fasta
   30
$ wc -l secuencias-330.aa.mafft.trim.faltan.txt
   10 secuencias-330.aa.mafft.trim.faltan.txt

 

Eso quiere decir que no encontró 10 de las 40 proteínas originales. Si vemos el contenido del archivo secuencias-330.aa.mafft.trim.faltan.txt veremos que uno de los genes que no encontró es el que está en el archivo: Genbank_nu/HG792021.1.nu.gbff. Si buscamos en NCBI el contenido de HG792021.1 veremos que si bien si están las proteínas, el DNA está ausente. Lo que encontraremos en su lugar son los nombres de otros archivos en donde se encuentra el DNA. Si se quiere contar con los genes de cualquier forma, entonces hay que buscar esos archivos en donde se encuentra el DNA y extraer el gen de interés. En este tutorial no buscaremos dichos genes faltantes.

 

Paso 8 Preparar el formato de las secuencias de genes para los análisis filogenéticos

En esta ocasión, usaremos los nombres de las proteínas codificadas por los genes. Usaremos estos nombres porque ya tienen la clasificación taxonómica. Con el siguiente script hacemos la sustitución de nombres:

 

$ perl cambianombre.pl secuencias-330.fasta secuencias-330.aa.mafft.trim.CDS.fasta

$ mv outfile.nu.fasta  secuencias-330.aa.mafft.trim.CDS.nu.fasta

$ mv outfile.aa.fasta  secuencias-330.aa.mafft.trim.CDS.aa.fasta

 

------------------------------------------------------------------------------------------

#!/usr/bin/perl


# Script: cambianombre.pl

# uso: 

# perl cambianombre.pl secuencias-330.fasta secuencias-330.aa.mafft.trim.CDS.fasta


# Luis Jose Delaye Arredondo

# luis.delaye@cinvestav.mx

# http://www.ira.cinvestav.mx/evolutionary.genomics

# Copy left


use strict;


my $file1 = $ARGV[0];

my $file2 = $ARGV[1];


my %hashN;

my %hashS;


my $r = 0;

my $sec;

my $name;


open (MIA, "$file1") or die ("No puedo abrir $file1\n");

while (my $linea = <MIA>){

chomp ($linea);

if ($linea =~ />/ && $r == 1){

$hashS{$name} = $sec;

$r = 0;

$name = $sec = ();

}

if ($linea !~ />/ && $r == 1){

$sec = $sec.$linea;

}

if ($linea =~ />/ && $r == 0){

#$name = $linea;

if ($linea =~ />(XP_\d+)_(.+)/){

#print ("-> $1\n");

$name = $1;

$hashN{$name} = $linea;

} elsif ($linea =~ />([A-Z0-9]+)_(.+)/){

#print ("-> $1\n");

$name = $1;

$hashN{$name} = $linea;

} else {

die ("el patron 1 fallo:\n$linea\n");

}

$r = 1;

}

}

if ($r == 1){

$hashS{$name} = $sec;

$r = 0;

$name = $sec = ();

}

close (MIA);



my $n = 0;

$r = 0;

$name = $sec = ();

open (MIA, "$file2") or die ("No puedo abrir $file2\n");

open (ROB, ">outfile.nu.fasta") or die ("No puedo abrir outfile.nu.fasta\n");

open (SOL, ">outfile.aa.fasta") or die ("No puedo abrir outfile.aa.fasta\n");  

while (my $linea = <MIA>){

chomp ($linea);

if ($linea =~ />/ && $r == 1){

if (exists $hashN{$name}){

print ROB ("$hashN{$name}\n$sec\n");

print SOL ("$hashN{$name}\n$hashS{$name}\n");

$n++;

}

$r = 0;

$name = $sec = ();

}

if ($linea !~ />/ && $r == 1){

$sec = $sec.$linea;

}

if ($linea =~ />/ && $r == 0){

if ($linea =~ />(XP_\d+)\s/){

$name = $1;

} elsif ($linea =~ />([A-Z0-9]+)\s/){

$name = $1;

} else {

die ("el patron 2 fallo:\n$linea\n");

}

$r = 1;

}

}

if ($r == 1){

if (exists $hashN{$name}){

print ROB ("$hashN{$name}\n$sec\n");

print SOL ("$hashN{$name}\n$hashS{$name}\n");

$n++;

}

$r = 0;

$name = $sec = ();

}

close (SOL);

close (ROB);

close (MIA);


 

print ("cambio el nombre a $n genes\n");

------------------------------------------------------------------------------------------
 
A continuación vamos a revisar si las secuencias tienen intrones y si es así, eliminarlos. Para ello, va a ser necesario que instales un programa que se llama exonerate que lo puedes descargar de aquí. Este script utiliza exonerate para alinear las proteínas a sus genes correspondientes.
 
$ perl alineaproteinaDNA.pl secuencias-330.aa.mafft.trim.CDS.aa.fasta secuencias-330.aa.mafft.trim.CDS.nu.fasta
 
------------------------------------------------------------------------------------------
#!/usr/bin/perl

# Script: alineaproteinaDNA.pl
# uso: 
# perl alineaproteinaDNA.pl secuencias-330.aa.mafft.trim.CDS.aa.fasta secuencias-330.aa.mafft.trim.CDS.nu.fasta

# Luis Jose Delaye Arredondo
# luis.delaye@cinvestav.mx
# http://www.ira.cinvestav.mx/evolutionary.genomics
# Copy left

use strict;

my $fileAA = $ARGV[0];
my $fileAN = $ARGV[1];

my %keysMs; 
my %hashAA;
my %hashAN;

my $r = 0;
my $name;
my $sec;
my $naa = 0;
open (MIA, "$fileAA") or die ("No puedo abrir $fileAA\n");
while (my $linea = <MIA>){
chomp ($linea);
if ($linea =~ />/ && $r == 1){
my $gi = ();
if ($name =~ />(.+)/){
$gi = $1;
} else {
die ("el patron I fallo: $name\n");
}
#$name =~ s/\s//g;
$sec =~ s/-//g;
$hashAA{$gi} = $sec;
$keysMs{$gi} += 1;
$name = $sec = ();
$r = 0;
$naa++;
}
if ($linea !~ />/ && $r == 1){
$sec = $sec.$linea;
}
if ($linea =~ />/ && $r == 0){
$name = $linea;
$r = 1;
}
}
if ($r == 1){
my $gi = ();
if ($name =~ />(.+)/){
$gi = $1;
} else {
die ("el patron I fallo: $name\n");
}
#$name =~ s/\s//g;
$sec =~ s/-//g;
$hashAA{$gi} = $sec;
$keysMs{$gi} += 1;
$name = $sec = ();
$r = 0;
$naa++;
}
close (MIA);

$r = 0;
$name = ();
$sec = ();
my $nan = 0;
open (MIA, "$fileAN") or die ("No puedo abrir $fileAN\n");
while (my $linea = <MIA>){
chomp ($linea);
if ($linea =~ />/ && $r == 1){
my $gi = ();
if ($name =~ />(.+)/){
$gi = $1;
} else {
die ("el patron I fallo: $name\n");
}
#$name =~ s/\s//g;
$sec =~ s/-//g;
$hashAN{$gi} = $sec;
$keysMs{$gi} += 1;
$name = $sec = ();
$r = 0;
$nan++;
}
if ($linea !~ />/ && $r == 1){
$sec = $sec.$linea;
}
if ($linea =~ />/ && $r == 0){
$name = $linea;
$r = 1;
}
}
if ($r == 1){
my $gi = ();
if ($name =~ />(.+)/){
$gi = $1;
} else {
die ("el patron I fallo: $name\n");
}
#$name =~ s/\s//g;
$sec =~ s/-//g;
$hashAN{$gi} = $sec;
$keysMs{$gi} += 1;
$name = $sec = ();
$r = 0;
$nan++;
}
close (MIA);

print ("$naa $nan\n");

if ($naa != $nan){
die ("el numero de proteinas no es igual al numero de genes\n");
}

system ("mkdir exodir");

my @keys = sort keys (%keysMs);

for (my $i = 0; $i <= $#keys; $i++){
if ($keysMs{$keys[$i]} == 2){
my $outfile = $keys[$i];
$outfile =~ s/>//;
$outfile = $outfile.'.exo.txt';
open (ROB, ">sequence.faa") or die ("No puedo abrir secuence.faa\n");
print ROB (">GI:$keys[$i]\n$hashAA{$keys[$i]}\n");
close (ROB);
open (ROB, ">sequence.fan") or die ("No puedo abrir secuence.fan\n");
print ROB (">GI:$keys[$i]\n$hashAN{$keys[$i]}\n");
close (ROB);
system ("exonerate --model protein2genome sequence.faa sequence.fan > exodir/$outfile");
system ("rm sequence.faa sequence.fan");
} else {
die ("error: $keys[$i] $keysMs{$keys[$i]}\n");
}
}
------------------------------------------------------------------------------------------
 
Y a continuación ejecutamos el siguiente script. Este script elimina los intrones (en caso de haberlos) de las secuencias de DNA.
 
$ perl eliminaintrones.pl
 
------------------------------------------------------------------------------------------
#!/usr/bin/perl

# Script: eliminaintrones.pl
# uso: 
# perl eliminaintrones.pl secuencias-330.aa.mafft.trim.CDS.aa.fasta secuencias-330.aa.mafft.trim.CDS.nu.fasta

# Luis Jose Delaye Arredondo
# luis.delaye@cinvestav.mx
# http://www.ira.cinvestav.mx/evolutionary.genomics
# Copy left

use strict;

my %genetic_code = (
'TCA' => 'Ser', # Serine
'TCC' => 'Ser', # Serine
'TCG' => 'Ser', # Serine
'TCT' => 'Ser', # Serine
'TTC' => 'Phe', # Phenylalanine
'TTT' => 'Phe', # Phenylalanine
'TTA' => 'Leu', # Leucine
'TTG' => 'Leu', # Leucine
'TAC' => 'Tyr', # Tyrosine
'TAT' => 'Tyr', # Tyrosine
'TAA' => '_', # Stop
'TAG' => '_', # Stop
'TGC' => 'Cys', # Cysteine
'TGT' => 'Cys', # Cysteine
'TGA' => '_', # Stop
'TGG' => 'Trp', # Tryptophan
'CTA' => 'Leu', # Leucine
'CTC' => 'Leu', # Leucine
'CTG' => 'Leu', # Leucine
'CTT' => 'Leu', # Leucine
'CCA' => 'Pro', # Proline
'CAT' => 'His', # Histidine
'CAA' => 'Gln', # Glutamine
'CAG' => 'Gln', # Glutamine
'CGA' => 'Arg', # Arginine
'CGC' => 'Arg', # Arginine
'CGG' => 'Arg', # Arginine
'CGT' => 'Arg', # Arginine
'ATA' => 'Ile', # Isoleucine
'ATC' => 'Ile', # Isoleucine
'ATT' => 'Ile', # Isoleucine
'ATG' => 'Met', # Methionine
'ACA' => 'Thr', # Threonine
'ACC' => 'Thr', # Threonine
'ACG' => 'Thr', # Threonine
'ACT' => 'Thr', # Threonine
'AAC' => 'Asn', # Asparagine
'AAT' => 'Asn', # Asparagine
'AAA' => 'Lys', # Lysine
'AAG' => 'Lys', # Lysine
'AGC' => 'Ser', # Serine
'AGT' => 'Ser', # Serine
'AGA' => 'Arg', # Arginine
'AGG' => 'Arg', # Arginine
'CCC' => 'Pro', # Proline
'CCG' => 'Pro', # Proline
'CCT' => 'Pro', # Proline
'CAC' => 'His', # Histidine
'GTA' => 'Val', # Valine
'GTC' => 'Val', # Valine
'GTG' => 'Val', # Valine
'GTT' => 'Val', # Valine
'GCA' => 'Ala', # Alanine
'GCC' => 'Ala', # Alanine
'GCG' => 'Ala', # Alanine
'GCT' => 'Ala', # Alanine
'GAC' => 'Asp', # Aspartic Acid
'GAT' => 'Asp', # Aspartic Acid
'GAA' => 'Glu', # Glutamic Acid
'GAG' => 'Glu', # Glutamic Acid
'GGA' => 'Gly', # Glycine
'GGC' => 'Gly', # Glycine
'GGG' => 'Gly', # Glycine
'GGT' => 'Gly'  # Glycine
);

my %aminoacid_code = (
'Ser' => 'S', # Serine
'Phe' => 'F', # Phenylalanine
'Leu' => 'L', # Leucine
'Tyr' => 'Y', # Tyrosine
'Cys' => 'C', # Cysteine
'Trp' => 'W', # Tryptophan
'Pro' => 'P', # Proline
'His' => 'H', # Histidine
'Gln' => 'Q', # Glutamine
'Arg' => 'R', # Arginine
'Ile' => 'I', # Isoleucine
'Met' => 'M', # Methionine
'Thr' => 'T', # Threonine
'Asn' => 'N', # Asparagine
'Lys' => 'K', # Lysine
'Val' => 'V', # Valine
'Ala' => 'A', # Alanine
'Asp' => 'D', # Aspartic Acid
'Glu' => 'E', # Glutamic Acid
'Gly' => 'G', # Glycine
);

my $nonpar = 1;
my $protein;
my $nucleic;
my $v = 0;

my @files = glob ("exodir/*.txt");
my @file;

open (ROB, ">outfile.nu.e.fasta") or die ("No puedo abrir outfile.nu.e.fasta\n");
open (SOL, ">outfile.aa.e.fasta") or die ("No puedo abrir outfile.aa.e.fasta\n");
for (my $i = 0; $i <= $#files; $i++){
print ("$files[$i]\n");
open (MIA, "$files[$i]") or die ("No puedo abrir $files[$i]\n");
CICLOA:
while (my $linea = <MIA>){
chomp ($linea);
if ($linea =~ /C4 Alignment:/ && $v == 1){
last (CICLOA);
}
if ($linea =~ /C4 Alignment:/ && $v == 0){
$v = 1;
}
if ($linea =~ /\s*\d+\s+:\s{1}(.+)\s{1}:\s+\d+/ && $v == 1){
my $str = $1;
#print ("$str\n");
if ($nonpar == 1){
$protein = $protein.$str;
$nonpar = 2;
} else {
$nucleic = $nucleic.$str;
$nonpar = 1;
}
}
}
close (MIA);
$protein =~ s/\s/./g;
#print ("\n");
#print ("$protein\n");
#print ("$nucleic\n");
my $lengthaa = length($protein);
my $lengthan = length($nucleic);
die ("las secuencias no son iguales: $lengthaa =! $lengthan") if ($lengthaa != $lengthan);
my @DNA = split (//, $nucleic);
my @Pro = split (//, $protein);
my $gen;
my $pep;
my $icodon;
my $iresid;
my $phase = 0;
for (my $i = 0; $i <= $#DNA; $i = $i +1){
if ($Pro[$i] =~ /[A-Z]/ && $Pro[$i+1] =~ /[a-z]/ && $Pro[$i+2] =~ /[a-z]/){
if ($DNA[$i] =~ /[ATCG]/ && $DNA[$i+1] =~ /[ATCG]/ && $DNA[$i+2] =~ /[ATCG]/){
$phase = 1;
}
}
if ($phase == 1){
my $codon = $DNA[$i].$DNA[$i+1].$DNA[$i+2];
my $resid = $Pro[$i].$Pro[$i+1].$Pro[$i+2];
#print ("$codon\t$genetic_code{$codon}\t$resid\n");
if (exists $genetic_code{$codon}){
if ($genetic_code{$codon} eq $resid){
$gen = $gen.$codon;
if (exists $aminoacid_code{$resid}){
$pep = $pep.$aminoacid_code{$resid};
} else {
die ("el aminoacido no existe: ($resid)\n");
}
} else {
$gen = $gen.'NNN';
$pep = $pep.'X';
#die ("error en el codigo: $genetic_code{$codon} ne ($resid) (1)\n");
print ("error en el codigo: $genetic_code{$codon} ne ($resid) (1)\n");
}
}
$phase = 0;
} else {
if ($DNA[$i] =~ /\{/ && length($icodon) == 0){
my $tres = $DNA[$i].$DNA[$i+1].$DNA[$i+2];
if ($tres =~ /\{(\w)\}/){
$icodon = $1;
} elsif ($tres =~ /\{(\w\w)/){
$icodon = $1;
} else {
die ("error 1 $DNA[$i] $DNA[$i+1] $DNA[$i+2]\n");
}
$tres = $Pro[$i].$Pro[$i+1].$Pro[$i+2];
if ($tres =~ /\{(\w)\}/){
$iresid = $1;
} elsif ($tres =~ /\{(\w\w)/){
$iresid = $1;
} else {
die ("error 2 $Pro[$i] $Pro[$i+1] $Pro[$i+2]\n");
}
} elsif ($DNA[$i] =~ /\{/ && length($icodon) > 0){
my $tres = $DNA[$i].$DNA[$i+1].$DNA[$i+2];
if ($tres =~ /\{(\w)\}/){
$icodon = $icodon.$1;
} elsif ($tres =~ /\{(\w\w)/){
$icodon = $icodon.$1;
} else {
die ("error 3 $DNA[$i] $DNA[$i+1] $DNA[$i+2]\n");
}
$tres = $Pro[$i].$Pro[$i+1].$Pro[$i+2];
if ($tres =~ /\{(\w)\}/){
$iresid = $iresid.$1;
} elsif ($tres =~ /\{(\w\w)/){
$iresid = $iresid.$1;
} else {
die ("error 4 $Pro[$i] $Pro[$i+1] $Pro[$i+2]\n");
}
}
if (length($icodon) == 3){
if (exists $genetic_code{$icodon}){
if ($genetic_code{$icodon} eq $iresid){
#print ("$icodon\t$genetic_code{$icodon}\t$iresid\n");
$gen = $gen.$icodon;
if (exists $aminoacid_code{$iresid}){
$pep = $pep.$aminoacid_code{$iresid};
} else {
die ("el aminoacido no existe: ($iresid)\n");
}
$icodon = ();
} else {
$gen = $gen.'NNN';
$pep = $pep.'X';
#die ("error en el codigo: $genetic_code{$icodon} ne ($iresid) (2)\n");
print ("error en el codigo: $genetic_code{$icodon} ne ($iresid) (2)\n");
$icodon = ();
}
}
} elsif (length($icodon) > 3){
die ("length $icodon > 3\n");
}
}
}
my $id = $files[$i];
$id =~ s/\.exo\.txt//;
$id =~ s/exodir\///;
print ROB (">$id\n$gen\n");
print SOL (">$id\n$pep\n");
#my $pausa = <STDIN>;
}
close (SOL);
close (ROB);

------------------------------------------------------------------------------------------

 
El script eliminaintrones.pl tiene como salida dos archivos: uno con los genes sin intrones (outfile.nu.e1.fasta) y el otro con las proteínas (outfile.aa.e1.fasta)  correspondientes. Hay que cambiarles los nombres
 
$ mv outfile.nu.e.fasta secuencias-330.aa.mafft.trim.CDS.nu.e.fasta
$ mv outfile.aa.e.fasta secuencias-330.aa.mafft.trim.CDS.aa.e.fasta
 
Ahora eliminamos el codon de STOP de cada una de los genes. Para ello usamos el siguiente script:
 
$ perl eliminaSTOP.pl secuencias-330.aa.mafft.trim.CDS.nu.e.fasta
 
------------------------------------------------------------------------------------------
#!/usr/bin/perl

# Este programa elimina los codones de paro al final de las secuencias
# de un archivo en formato fasta

# Script: eliminaSTOP.pl
# uso: 
# perl eliminaSTOP.pl secuencias-330.aa.mafft.trim.CDS.nu.e.fasta
# out: secuencias-330.aa.mafft.trim.CDS.nu.e-e1.fasta

# Luis Jose Delaye Arredondo
# luis.delaye@cinvestav.mx
# http://www.ira.cinvestav.mx/evolutionary.genomics
# Copy left

use strict;

my $file = $ARGV[0];

my $out = $file;
if ($out =~ s/.fasta/-e1.fasta/){

} else {
die ("error, el archivo $file no termina con extension fasta\n");
}

open(MIA,"$file") or die ("No puedo abrir $file\n");
open(ROB,">$out") or die ("No puedo abrir $out\n");
my $nombre;
my $secuencia;
my $r = 0;
while (my $linea = <MIA>) { 
        chomp ($linea);

        if ($linea =~ />/ && $r == 1) {
                # Aqui manipulo las secuencias
                my @a = split (//, $secuencia);
                my $edit = ();
                my $codon = $a[-3].$a[-2].$a[-1];
                if ($codon =~ /TAA/i || $codon =~ /TAG/i || $codon =~ /TGA/i){
               for (my $i = 0; $i <= ($#a -3); $i++){
                $edit = $edit.$a[$i];
}
} else {
for (my $i = 0; $i <= $#a; $i++){
                $edit = $edit.$a[$i];
}
}
my $res = ($#a +1) % 3;
if ($res == 0){
print ROB ("$nombre\n$edit\n");
} else {
print ("secuencia eliminada: $nombre\n");
}
                $nombre = ();
                $secuencia = ();
                $r = 0;
        }
        if ($r == 1) {
                $secuencia = $secuencia.$linea;
        }
        if ($linea =~ />/ && $r == 0){
                $nombre = $linea;
                $r = 1;
        }
}
if ($r == 1) {
        # Aqui manipulo las secuencias
        my @a = split (//, $secuencia);
        my $edit = ();
        my $codon = $a[-3].$a[-2].$a[-1];
        if ($codon =~ /TAA/i || $codon =~ /TAG/i || $codon =~ /TGA/i){
    for (my $i = 0; $i <= ($#a -3); $i++){
        $edit = $edit.$a[$i];
}
} else {
for (my $i = 0; $i <= $#a; $i++){
        $edit = $edit.$a[$i];
}
}
my $res = ($#a +1) % 3;
if ($res == 0){
print ROB ("$nombre\n$edit\n");
} else {
print ("secuencia eliminada: $nombre\n");
}
$nombre = ();
        $secuencia = ();
        $r = 0;
}
close (ROB);
close (MIA);
------------------------------------------------------------------------------------------
 
Este script tiene como salida el archivo: secuencias-330.aa.mafft.trim.CDS.nu.e-e1.fasta.
 
  
¡Estamos trabajando!
 
:-)