RandProt

genere secuencias de proteína aleatorias de Markov de alto orden

por Alejandro Ochoa García, John D Storey, Manuel Llinás, Mona Singh

Implementación rápida que captura sutilezas de secuencias de proteína

Imagen miniatura
VIIIA

es-mx en-us - - Mándame un correo

Aquí puedes bajar nuestro código, y aprender como usarlo. Esta página es el manual del programa.

Acerca de RandProt

El propósito de este programa es generar eficientemente secuencias de proteína aleatorias de un modelo de Markov de alto orden. A grandes rasgos, un modelo de Markov de orden \(k-1\) preserva la distribución de \(k\)-meros de los datos originales, y es así como el modelo está codificado internamente. El código enfatiza la simplicidad y eficiencia computacional. Note, sin embargo, que la mayoría del código es específico a las secuencias de proteína, así que estará difícil reusarlo para generar otras secuencias sin cambios significativos.

La especificación formal de las secuencias aleatorias viene a continuación. La entrada es una colección de secuencias de proteína. La salida tendrá, para cada secuencia de entrada, \(n\) secuencias aleatorias del mismo largo (con la misma identificación más, para \(n > 1\), el sufijo.RAND\(i\) donde \(i\) va de \(0\) a \(n-1\) acolchonado con ceros). La distribución de \(k\)-meros se obtiene de los datos de entrada excluyendo a las metioninas (M) iniciales, si están presentes. Las secuencias aleatorias siempre comienzan con M, seguida por un \(k\)-mero aleatorio sacado de la distribución global, seguido por aminoácidos individuales sacados del modelo Markov (es decir, condicionalmente en los \(k-1\) aminoácidos previos). La secuencia se termina cuando la longitud deseada se haya alcanzado.

Hay más detalles menos importantes. La entrada siempre es "normalizada" para que sólo los 20 aminoácidos canónicos sean considerados. Para lograr esto, las selenocisteínas (U) son reemplazadas por cisteínas (C), y pirrolisina (O) por lisina (K). Adicionalmente, los \(k\)-meros con códigos de aminoácidos ambiguos (B,J,X,Z) siempre son ignorados.

Código fuente

/dl.png RandProt-1.01.tgz (24 KB)

Baje nuestro código fuente Perl, en un archivo tar comprimido con gzip. Es muy pequeño, 4 paquetes y 3 guiones para un total de 468 lineas de código (de acuerdo a cloc).

Todo mi código está publicado bajo la GNU GPLv3 (Licencia Pública General de GNU, versión 3).

Este proyecto ahora está en GitHub, donde usted puede contribuir y mejorar este programa.

Versiones previas

Código RandProt-1.00.tgz (23 KB) y manual 1.00.

Instalación

Nuestro código de perl es autosuficiente, así que sólo Perl 5 es requerido (no usa paquetes externos de perl). Todo lo que necesita hacer es desempacar el código y correr los guiones del directorio que los contiene. En este modo, todos los guiones y paquetes deberían quedarse juntos en el mismo directorio. También puede correr mis guiones desde directorios arbitrarios.

Nuestro código asume que gzip está disponible en su sistema.

Corriendo los guiones de RandProt

Recomendaciones

A parte del archivo de entrada de secuencias, el análisis principal requiere los parámetros \(k\) y \(n\) que son enteros positivos. El valor de \(n\) provee el número de muestras aleatorias por secuencia real. En teoría, mientras más grande mejor, pero la complejidad del análisis posterior impondrá límites superiores para \(n\). Yo he usado personalmente \(n=100\) y \(n=20\) para proteomas pequeños, y \(n=1\) para una base de datos de proteínas mucho más grande.

El valor de \(k\) controla cuanto las secuencias aleatorias se parecerán a las secuencias reales, con valores pequeños de \(k\) dando proteinas aleatorias menos realísticas. Sin embargo, un valor de \(k\) que sea demasiado grande podría ser malo, ya que podría recapitular propiedades de sus proteinas reales que no desee tener en las secuencias aleatorias (por ejemplo, en mi investigación sin publicar encontré que algunos dominios de proteínas se les puede predecir muy bien por la presencia de \(5\)-meros que representan bloques altamente conservados, así que aún \(k=5\) preserva alguna señal biológicamente relevante que tal vez deba ser excluida de sus secuencias aleatorias). En mi trabajo favorecí a \(k=3\).

Archivo de secuencia de muestra

Para tratar estos ejemplos, cualquier archivo de sequencias de proteina formato FASTA se puede usar. Si gusta, puede usar este archivo, PLAF7.fa.gz, que es el proteoma de Plasmodium falciparum de PlasmoDB versión 11.0, con pseudogenes removidos, contiene 5441 proteínas y un tamaño de 2.3 MB (comprimido).

Sinopsis de guiones

Los siguientes comandos pueden correrse en el archivo de ejemplo provisto arriba colocado en el mismo directorio que el código y llamados de ese mismo lugar. El archivo de entrada puede estar comprimido con gzip y también puede ser especificados con o sin la extensión GZ. El archivo de salida será comprimido con gzip, ya sea si la ruta de salida indique la extensión GZ o no. Así que los comandos así como se muestran abajo producirán y funcionarán completamente con archivos comprimidos, sin una sola extension GZ indicada.

# obtenga un estimado rápido de la k máxima a considerar 
perl -w kMax.pl PLAF7.fa 

# trate varios valores de k hasta que una cobertura decente sea lograda 
perl -w kCov.pl PLAF7.fa 5 
perl -w kCov.pl PLAF7.fa 4 
perl -w kCov.pl PLAF7.fa 3 

# ¡genere su proteoma aleatorio deseado! 
perl -w randProt.pl PLAF7.fa PLAF7.k3.n100.fa 3 100

kMax.pl 1.00 - Calcule una cota superior débil para \(k\) para análisis de \(k\)-mero

perl -w kMax.pl <input FASTA>
Las entradas requeridas son

El archivo de entrada puede ser comprimido con gzip, y se le puede especificar con o sin la extensión GZ. El análisis es mandado a la salida estándar.

El razonamiento detrás de este guión se describe a continuación. El valor máximo recomendado de \(k\) depende en parte del tamaño de sus secuencias de entrada. Esto es porque conjuntos demasiado pequeños de secuencias puede que no contengan suficientes muestras de todos los \(k\)-meros si \(k\) es demasiado chica. Una secuencia de largo \(l\) tiene a lo mucho \(l-k+1\) de \(k\)-meros únicos, mientras que hay \(20^k\) de \(k\)-meros únicos en total. Por ende, para asegurar que todos los \(k\)-meros puedan aparecer en la entrada, ahora para un proteoma de \(l\) aminoácidos (asumiendo que \(k\) es mucho más pequeña que la mayoría de los largos de proteína) \(k\) deberá satisfacer la siguiente cota superior débil, $$ k \le k_{max} = \left \lfloor \frac{\log(l)}{\log(20)} \right \rfloor, $$ pero usualmente \(k\) debe ser mucho más pequena de lo que esta cota superior sugiere, para asegurar que el espacio esté cubierto con muestras adecuadamente. Para el archivo de entrada de FASTA, este guión manda a la salida los valores de \(l\) y \(k_{max}\), justo como en el ejemplo abajo.

Input: PLAF7.fa.gz
Number of letters: 4139904
kMax: 5
Este es un análisis muy rápido y aproximado a propósito. Para una evaluación más precisa de la covertura de \(k\)-meros para un conjunto de datos en particular, véase el siguiente guión.

kCov.pl 1.00 - Calcule el porcentaje de \(k\)-meros observados en un proteoma

perl -w kCov.pl <input FASTA> <k>
Las entradas requeridas son

El archivo de entrada puede ser comprimido con gzip, y se le puede especificar con o sin la extensión GZ. El análisis es mandado a la salida estándar.

Este guión simplemente obtiene la distribución de \(k\)-meros del archivo de entrada, y compara el número de \(k\)-meros únicos observados con el número de \(k\)-meros únicos posibles en teoría. Para el archivo de muestra, las salidas para \(k=5,4,3\) (corridas por separado) se muestran abajo, lo cual demuestra que \(k=3\) debería ser el valor máximo de \(k\) por considerar para estos datos.

k: 5
Number of theoretical amino acid k-mers: 3200000
Scanning proteome...
Number of observed amino acid k-mers: 1226974
Proportion of k-mers observed: 38.34 %

k: 4
Number of theoretical amino acid k-mers: 160000
Scanning proteome...
Number of observed amino acid k-mers: 148284
Proportion of k-mers observed: 92.68 %

k: 3
Number of theoretical amino acid k-mers: 8000
Scanning proteome...
Number of observed amino acid k-mers: 7999
Proportion of k-mers observed: 99.99 %
Recomiendo áltamente escoger una \(k\) para la cual la mayoría de los \(k\)-meros sean observados, como \(k=3\) en el caso de arriba, pero esta es sólamente una cota superior, y puede que usted desee escover una \(k\) más pequeña usando otro criterio.

Esta es la manera más refinada de escover una cota superior para \(k\), es más precisa que la salida de kMax.pl. Sin embargo, kMax.pl será mucho más rápido con entradas grandes, así que es recomendado correrlo primero para tener una idea aproximada de la gama de \(k\) por probar.

randProt.pl 1.00 - Cree secuencias de proteinas aleatorias con un modelo de Markov de alto orden

perl -w randProt.pl <input FASTA> <output FASTA> <k> <n>
Las entradas requeridas son

El archivo de entrada puede ser comprimido con gzip, y se le puede especificar con o sin la extensión GZ. El archivo de salida será comprimido automáticamente con gzip.

Para el archivo de entrada de muestra, proveo aquí dos archivo de salida aleatorio. El primero usó \(k=3\) y \(n=1\), PLAF7.k3.n1.fa.gz (2.3 MB). El segundo usó \(k=3\) y \(n=100\), PLAF7.k3.n100.fa.gz (227 MB). De este modo pueden ver el formato, pero no va a ser identico a tus archivos de salida porque siempre son diferentes. En mi máquina, me tomó 20 segundos, y 24 minutos, respectivamente, para general estos proteomas aleatorios.

Notas compartidas con otros proyectos

Por favor lea notas adicionales siguiendo el vínculo a continuación.

Mi código público
Mi código público, descripción de mis paquetes originales de Perl
La unión de código a través de las distribuciones de DomStratStats, dPUC, y RandProt más recientes.

Secuencias aleatorias basadas en UniRef50

En nuestro artículo citado abajo, usamos el siguiente archivo de secuencias de proteína aleatorias generado por nuestro programa, usando \(n=1\). ¡Es un archivo grande!

/dl.png UNIREF50.k3.fa.gz (764 MB)

Cita

2015-11-17. Alejandro Ochoa, John D Storey, Manuel Llinás, and Mona Singh. Beyond the E-value: stratified statistics for protein domain prediction. PLoS Comput Biol. 11 e1004509. Pubmed Artículo arXiv 2014-09-23.

Notas de la versión

2015-06-04 - RandProt 1.01

Ahora el caso \(n=1\) genera archivos con identificaciones originales y sin sufijos.

Aumenté README.md, concordante con lo que tengo en GitHub.

FileGz.pm tuvo actualizaciones sin relación a este proyecto, el cual no fue afectado.

2014-09-12 - RandProt 1.00

Esta fue la primera versión pública del código.

VIIIA

Historial