Foro - Perl en Español

Hola a todos, estoy teniendo un problema con un script para tratar de obtener los promotores de Arabidopsis en formato gff3.

En principio los obtuve simplemente contando -500 nt río arriba del inicio de los genes, pero luego me di cuenta que muchos de ellos se solapan con genes contiguos.

De manera que solo me quiero quedar con los que no solapen con genes.

Para eso escribí un script para que me encuentre coordenadas y si me solapan las exporte a un archivo X y si no a otro Y.

Parto de dos archivos de entrada, el de los genes en formato gff3 y el de los promotores (genes -500 nt) también en gff3.



Sé que mis conocimientos no son muy amplios y que debe haber una mejor manera de hacer esto, pero en general uso esta estructura para varios de mis scripts en donde tengo que hacer coincidir coordenadas.

Mis problemas son:

1. Tuve que poner dos if() porque cuando uso solo uno con todas las condiciones no me reconoce el ($chr1 eq $chr2), no entiendo por qué ya que varias veces usé algo similar y me ha funcionado
2. Los archivos de salida están bien, pero revisándolos encontré algunos promotores en el archivos "sin_solapar" que también aparecen en el "solapados" y de hecho sí se solapan, con lo cual algo no está funcionando. ¿Es el if() o el else{}? En teoría si un ID aparece en un archivo en el otro no debería aparecer, ¿no es así?
3. Otro fallo que tiene es que en el archivo "sin_solapar" me aparecen mas de una vez las lineas, es evidente que tengo algo en mi bucle

Si alguien puede decirme en qué está fallando mi script se los agradecería.

Puedes ahorrar mucho código, en

Y en

Como ves hemos quitado las variables $tair2, $gene2, $column55, $strand2 y $column77 porque no se usan luego. Y que el chomp() lo hacemos antes de la partición, para evitar que el último campo contenga el carácter de nueva línea.

Y en


También, es muy poco eficiente el tener que leer el segundo archivo, continuamente, por cada registro del primer archivo. Lo inteligente es leerlo solo una vez, pasando su información a un array, o mejor aún, a un hash, cuyas claves serán aquél campo que es único por cada registro del segundo archivo (¿cuál es? ¿$chr2?)

Lo del if(), viendo el código, me parece lógico: haces dos operaciones distintas según sean $chr1 igual o distinto a $chr2.

En el tema de la aparición de promotores en el archivo, cuando no deberían, todo depende de la condición puesta en el if(), y recordar que estás revisando todos los registros del segundo archivo con todos los registros del primero, por lo que si una misma condición se repite en todas esas condiciones, se repetirán los print() que sacan la información hacia afuera.

Para resolver esto, a mi me haría falta saber qué campos son únicos en cada archivo (si los hay), o qué hacer en caso de que se detecte una repetición de la condición (¿filtrar o no filtrar?).

Me falta información más detallada de este problema.

explorer, muchas gracias otra vez por tu respuesta.
Voy a empezar a implementar el split() como escribiste para ahorrar código.

¿A qué te refieres con campos únicos?

Los archivos tienen exactamente la misma estructura, con los mismos campos.

La primera condición if ($chr1 eq $chr2) es porque necesito que solape las coordenadas correspondientes al mismo cromosoma.

Los archivos tienen en el primer campo todos los cromosomas de Arabidopsis (Chr1, Chr2, Chr3, Chr4, etc). Para cada $chr igual en ambos archivos, tengo que hacer la segunda parte: (($coord1 > $coord3) and ($coord1 < $coord4)) or (($coord2 > $coord4) and ($coord2 > $coord3)), en donde las $coord1 y $coord2 corresponden al promotor y las $coord3 y $coord4 al gen.

Lo que intenté hacer con estas condiciones es que, si se solapan, las guarde en un archivo, y si no cumple esta condición, que guarde las coordenadas del promotor (primer archivo) en otro archivo distinto (OUTPUTFILE2).

El array lo hago para leer las líneas del segundo archivo una sola vez, guardándolas allí.

Gracias por los aportes.

dzavallo escribiste:¿A qué te refieres con campos únicos?

A que no sabemos si estos archivos contienen alguno de los campos, en las distintas columnas, que podamos considerar únicos (no se repiten en el resto de las líneas.

dzavallo escribiste:La primera condición if ($chr1 eq $chr2) es porque necesito que solape las coordenadas correspondientes al mismo cromosoma.

Bien, esto ya es una pista. La pregunta anterior entonces se puede replantear de la siguiente manera: ¿puede aparecer, en el archivo primero o en el segundo, más de una línea, con el mismo cromosoma? O... ¿cada cromosoma solo aparece una vez en todo el archivo primero y una sola vez en todo el archivo segundo?

dzavallo escribiste:Los archivos tienen en el primer campo todos los cromosomas de Arabidopsis (Chr1, Chr2, Chr3, Chr4, etc). Para cada $chr igual en ambos archivos, tengo que hacer la segunda parte: (($coord1 > $coord3) and ($coord1 < $coord4)) or (($coord2 > $coord4) and ($coord2 > $coord3)), en donde las $coord1 y $coord2 corresponden al promotor y las $coord3 y $coord4 al gen.

Entonces, si se repiten las condiciones de la condición if(), obtendrás varios resultados, desde luego.

dzavallo escribiste:El array lo hago para leer las líneas del segundo archivo una sola vez, guardándolas allí.

Me temo que no... estás leyendo el segundo archivo cada vez por cada línea del primero...

¿Puedes publicar dos trozos pequeños de ejemplo, de cada archivo de entrada, con los correspondientes archivos de salida, en que aparezcan todas las combinaciones posibles de las condiciones que quieres procesar? Viéndolo, será mucho más fácil buscar una solución.

Bien, esto ya es una pista. La pregunta anterior entonces se puede replantear de la siguiente manera: ¿puede aparecer, en el archivo primero o en el segundo, más de una línea, con el mismo cromosoma? O... ¿cada cromosoma solo aparece una vez en todo el archivo primero y una sola vez en todo el archivo segundo?

Claro, aparecen muchas veces, miles.

Aquí pego 10 líneas del archivo "solo_prom" correspondiente a los promotores:

Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 27999 28499 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 37872 38372 . - . ID=AT1G01060;Note=protein_coding_gene;Name=AT1G01060
Chr1 TAIR10 gene 40945 41445 . - . ID=AT1G01070;Note=protein_coding_gene;Name=AT1G01070
Chr1 TAIR10 gene 44176 44676 . + . ID=AT1G01073;Note=protein_coding_gene;Name=AT1G01073
Chr1 TAIR10 gene 47020 47520 . - . ID=AT1G01080;Note=protein_coding_gene;Name=AT1G01080

Aquí pego 10 líneas del archivo " solo_gene" correspondiente a los genes:

Chr1 TAIR10 gene 3631 5899 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 5928 8737 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 11649 13714 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 28500 28706 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 31170 33153 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 33379 37871 . - . ID=AT1G01060;Note=protein_coding_gene;Name=AT1G01060
Chr1 TAIR10 gene 38752 40944 . - . ID=AT1G01070;Note=protein_coding_gene;Name=AT1G01070
Chr1 TAIR10 gene 44677 44787 . + . ID=AT1G01073;Note=protein_coding_gene;Name=AT1G01073
Chr1 TAIR10 gene 45296 47019 . - . ID=AT1G01080;Note=protein_coding_gene;Name=AT1G01080
Chr1 TAIR10 gene 47485 49286 . - . ID=AT1G01090;Note=protein_coding_gene;Name=AT1G01090

y los archivos de salida son "genes_solapados":

Chr1 33154 33654 Chr1 33379 37871 ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
ID=AT1G01060;Note=protein_coding_gene;Name=AT1G01060
Chr1 47020 47520 Chr1 47485 49286 ID=AT1G01080;Note=protein_coding_gene;Name=AT1G01080
ID=AT1G01090;Note=protein_coding_gene;Name=AT1G01090
Chr1 63812 64312 Chr1 64166 67625 ID=AT1G01130;Note=protein_coding_gene;Name=AT1G01130
ID=AT1G01140;Note=protein_coding_gene;Name=AT1G01140
Chr1 72139 72639 Chr1 72339 74096 ID=AT1G01150;Note=protein_coding_gene;Name=AT1G01150
ID=AT1G01160;Note=protein_coding_gene;Name=AT1G01160
Chr1 71838 72338 Chr1 70115 72138 ID=AT1G01160;Note=protein_coding_gene;Name=AT1G01160
ID=AT1G01150;Note=protein_coding_gene;Name=AT1G01150

y "genes_sin_solapar":

Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 27999 28499 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 27999 28499 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 27999 28499 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050

Remarqué en color uno de los ID como ejemplo. Si observas las coordenadas del promotor AT1G01050 del archivo "solo_prom", se solapa con el gen AT1G01060 del archivo "solo_gene", porque lo que según las condiciones que quiero ajustar, debería archivarse solo en "genes_solapados", como efectivamente aparece, e imprimir los campos que quiero.

Pero no debería aparecer en "genes_sin_solapar" por la condición else{}. Sin embargo también me lo imprime allí (última línea).

Además, como se puede observar, me repite las líneas varias veces, no sé por qué.

Subo un gráfico en donde creo se ve más claramente las combinaciones que quiero hacer.


Aquí se ven las dos posibles maneras en que se solapen el promotor y un gen (si está en la hebra positiva o en la negativa). Es por ello que necesito todas esas condiciones juntas y que son excluyentes entre sí. No me debería aparecer en los dos archivos.



Ahora que tal vez me explique mejor, y con las correcciones que me pasaste, vuelvo a mostrar mi script, salvo que en esta oportunidad, el archivo "sin_solapar" repite la primera línea indefinidamente.



¿Qué estoy haciendo mal con mis condiciones?

Espero que haya sido más claro que mi anterior explicación, y desde ya agradezco mucho la ayuda.

Las condiciones están bien, salvo que has puesto dos veces $coord4 en la segunda parte. Te falta un $coord3.

Pero el principal problema no son las condiciones (la lógica de la condición es correcta). El problema está en que recorres todas las líneas de genes por todas las líneas de los promotores, imprimiendo el resultado por cada línea del archivo de genes, cuando en realidad debería imprimirse al final del bucle de los genes, sabiendo entonces si el promotor está solapado o no con alguno de los genes.

Creo que he encontrado una posible solución. Dados los dos primeros archivos de ejemplo, con este programa:
Primero leemos el archivo de los genes a memoria, pasándolo a un array de arrays (estructura bidimensional).

Luego, hacemos el doble bucle clásico, con las siguientes características:

• línea 38: si no coinciden los nombres de los cromosomas, saltamos a comparar al siguiente (next) gen
• líneas 40 y 41. Estas líneas contienen la misma lógica que tu condición de solapamiento, pero de diferente manera. En lugar de preguntar si el inicio y final del promotor está dentro de las coordenadas de inicio y final del gen, pregunto por la situación contraria: ¿están las coordenadas del promotor fuera del rango de las del gen? En la línea 40 pregunto si el inicio del promotor está más allá del final del gen. Si es así, no está solapado, y paso a analizar el siguiente gen. En la línea 41 pregunto si el final del promotor está antes del inicio del gen. Si es así, tampoco está solapado, y termino el bucle de genes, porque, si supongo que los archivos están ordenados según la posición, los siguientes genes también estarán "por detrás" del final del promotor (y no es posible más solapamientos)
• línea 43: aquí solo llegamos en caso de un solapamiento. Guardamos los datos que nos interesan en la variable (que nos servirá además, de bandera, fuero del bucle) e imprimimos el solapamiento
• línea 52: una vez terminado de buscar por todos los genes, preguntamos si hemos encontrado un solapamiento o no. Si no hemos encontrado ninguno, imprimiremos el promotor, como gen no solapado.

Solo lo he podido probar con los archivos de ejemplo que nos has pasado, y sale esto:

Como hemos dicho antes, estamos haciendo varias suposiciones, como que los archivos están perfectamente ordenados, según por cromosoma y según por la posición de inicio/final de cada gen.

Date cuenta de una cosa: este programa es prácticamente igual al tuyo, salvo por un detalle básico: solo imprimimos el gen promotor, si no está solapado, fuera del bucle de comparación con el resto de los genes, y en caso de que no haya habido ningún solapamiento con ellos.

Una vez más, ¡muchas gracias! Voy a probarlo ya que seguramente corra mucho más rápido que el mio.

Finalmente pude resolverlo agregando una variable (solapado) y la inicio en falso cada vez que lee el solo_prom, y cuando encuentra un solapado imprime en el archivo de salida y le pone verdadero. Al final le puse que si la variable (solapado) es igual a falso, lo imprima en no solapados. Pero me tardó mucho la corrida.


Gracias por la ayuda.

Dinos a qué velocidad va. Y si quieres, se puede acelerar aún más.