DomStratStats

estadísticas estratificadas para dominios de proteínas

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

Computa valores q y FDRs locales para dominios, y valores q con gradas (combinando informacion de dominios y sequencias de HMMER3)

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 DomStratStats

El propósito de este programa es traer mejores estadísticas al análisis de secuencias de proteina. Las estrategias normales están basadas en valores \(p\) y valores \(E\), pero aquí introducimos los valores \(q\) y el lFDR (siglas en inglés de la Tasa de Falso Descubrimiento local) para los dominios de proteinas. La clave para que los valores \(q\) y el lFDR trabajen mejor es estratificando, en este caso cada familia de dominio por separado, lo cual balancea los errores a través de las familias de dominio.

Prácticamente, la estrategia del "valor \(q\) de dominios" va a ser la más útil. Teóricamente, hemos probado en nuestro trabajo que el lFDR es óptimo; sin embargo, descubrimos en nuestras pruebas que las estimaciones de los lFDRs no funcionan también como los valores \(q\) porque algunos valores \(p\) de HMMER3 son malos (para ciertas familias de dominio repititivas, en particular). Encontramos que los valores \(q\) son mucho más robustos que los lFDRs dados los valores \(p\) que HMMER3 nos da. Nuestro código provée los lFDRs de todos modos para propósitos de investigación, y con la esperanza de que se puedan usar mejor en el futuro, cuando los valores \(p\) de dominio mejoren.

Finalmente, incluímos un modo aún más experimental, los "valores \(q\) con gradas". Este análisis con gradas combina los valores \(p\) de secuencia y de dominio de HMMER en una decisión de significancia combinada basada en el FDR. Por el lado bueno, produce aún más predicciones por FDR combinado que las otras estrategias. Sin embargo, el FDR teórico desvía demasiado del FDR empírico, lo cual significa que usted debería ejercitar precaución escogiendo el contraste de hipótesis.

Código fuente

/dl.png DomStratStats-1.04.tgz (50 KB)

Baje nuestro código fuente Perl, en un archivo tar comprimido con gzip. Tiene 8 paquetes y 5 guiones para un total de 1366 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 DomStratStats-1.03.tgz (47 KB) y manual 1.03.

Código DomStratStats-1.02.tgz (46 KB) y manual 1.02.

Código DomStratStats-1.01.tgz (46 KB) y manual 1.01.

Código DomStratStats-1.00.tgz (35 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.

También necesitará HMMER3 (sólo los ejecutables hmmpress y hmmscan son necesarios). Nuestro código también asume que gzip está disponible en su sistema.

Corriendo los guiones de DomStratStats

Recomendaciones

Los guiones no tienen contrastes de hipótesis estándar para los valores \(q\) o lFDRs. Valores de contraste differentes le permitirán balances diferentes de positivos verdaderos y falsos. Si quiere predicciones comparables en calidad a las de Pfam, nuestras pruebas sugieren un contraste de valor \(q\) de 4e-4 (más recomendado), o un contraste de lFDR de 2.5e-4 (los lFDRs son menos recomendados en general), o un contraste para valores \(q\) con gradas de 1e-4 (para cada grada, pero el FDR será impreciso en este modo). Pero en algunas maneras Pfam es demasiado conservador, así que si prefiere podrá usar contrastes menos severos. Por ejemplo, \(q \le \mbox{1e-3}\) debería dar un falso positivo por cada 1000 dominios (aunque el FDR verdadero será más alto si sus secuencias contienen hélices superenrolladas, dominios transmembranales, u otras secuencias de baja complejidad, ¡así que tengan cuidado!).

Los FDRs y lFDRs que calculamos se vuelven más precisos cuando las bases de datos son más grandes. En práctica, encontramos que las proteínas de un sólo organismo no son sufucientes para dar resultados óptimos, pero agrupara a varios organismos relacionados funciona mejor (por ejemplo, usualmente trabajo con 15 organismos del filo Apicomplexa). Los ejemplos de abajo rompen esta regla, sólo usan a Plasmodium falciparum para sólo lidiar con archivos pequeños. En general, es mejor empezar con una colección más grande de proteinas que representen la diversidad de los dominios contenidos en el organismo que usted estudie (como los Apicomplexa en el caso de P. falciparum, como lo hice en este otro proyecto). Los lFDRs son mucho más sensibles que los valores \(q\) al tamaño de los datos de entrada.

Bases de datos de HMM

Si escoge a Pfam, que es lo que usamos en los ejemplos, debería bajar Pfam-A.hmm.gz del sitio FTP de Pfam. Será necesario usar hmmpress de HMMER3 en esta base de datos de HMM antes de que hmmscan la pueda usar. Los siguientes comandos logran eso.

# descomprima la base de datos HMM primero 
gunzip Pfam-A.hmm.gz 

# prepare la base de datos HMM para búsquedas con HMMER3 (la cual provée hmmpress) 
# cuatro archivos diferentes se generarán, con nombres Pfam-A.hmm.h3{f,i,m,p} 
hmmpress Pfam-A.hmm 

# (opcional) recomprima el archivo HMM original (HMMER3 no lo usa) 
gzip Pfam-A.hmm

Para Pfam, es opcional pero altamente recomendado bajar Pfam-A.hmm.dat.gz también, ayuda a nuestro código reconocer dominios anidados. Este archivo puede permanecer comprimido, ¡DomStratStats lo leera correctamente!

Otra base de datos popular es Superfamily. Tendrá que registrarse para obtener el archivo hmmlib_1.75.gz. Se usa de la misma manera exactamente que Pfam-A.hmm.gz, pero también tendrá que correr hmmpress con él. No hay un equivalente en Superfamily de Pfam-A.hmm.dat.gz, así que desafortunadamente los dominios anidados no se les tratará correctamente. Hay algunas maneras importantes en que los HMMs de Superfamily difieren de Pfam; sin embargo, nuestro código no trata a Superfamily de manera especial. En todos casos, estratificamos el análisis de valores \(q\) y lFDRs por HMM.

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.3MB (comprimido).

Todos los archivos de salida para estas proteinas están listados abajo usando Pfam 31 y Superfamily 1.75, y la version de 64 bits de HMMER 3.1b2. Ningún archivo de muestra está incluído con nuestro código fuente. Podrán haber pequeñas diferencias en archivos de salida tan sólo usando una arquitectura de 32 bits, para su información. Adicionalmente, HMMER 3.1b2 fechea al archivo de salida, así que cada corrida parecera generar archivos diferentes si se comparan con funciones de resumen criprográfico como SHA o MD5, así que no las use. Es mejor comparar cada linea con una herramienta como zdiff, que es como diff pero para archivos comprimidos.

Descripción Pfam 31 Superfamily 1.75
Salida de HMMER3 sin procesar PLAF7.31.p.txt.gz PLAF7.s.txt.gz
Salida de HMMER3 sin superposiciones PLAF7.31.po.txt.gz PLAF7.so.txt.gz
Salida de domStratStats sin contrastes PLAF7.31.pod.txt.gz PLAF7.sod.txt.gz
Salida de domStratStats, \(q \le \mbox{4e-4}\) PLAF7.31.pod.q4e-4.txt.gz PLAF7.sod.q4e-4.txt.gz
Salida de domStratStats, \(\mbox{lFDR} \le \mbox{2.5e-2}\) PLAF7.31.pod.l2.5e-2.txt.gz PLAF7.sod.l2.5e-2.txt.gz
Salida de domStratStats, \(\mbox{FDR|lFDR} \le \mbox{4e-4}\) PLAF7.31.pod.fl4e-4.txt.gz PLAF7.sod.fl4e-4.txt.gz
Salida de tieredStratQ, \(q_{seq},q_{dom|seq} \le \mbox{1e-4}\) PLAF7.31.pot.1e-4.txt.gz PLAF7.sot.1e-4.txt.gz
Salida de tieredStratQ, \(q_{seq} \le \mbox{1e-4}\) (sin contraste de \(q_{dom|seq}\)) PLAF7.31.pot.1e-4,1.txt.gz PLAF7.sot.1e-4,1.txt.gz

La única excepción son los datos de muestra para el guión 4allManyOrgs.pl, el cual son datos mucho más grandes para otro proyecto, el cual sólo corrí con Pfam 27.

Sinopsis de guiones

Los siguientes comandos pueden correrse en los archivos de ejemplo provistos arriba, mas los dos archivos de Pfam, todos colocados en el mismo directorio que el código y llamados de ese mismo lugar. Todos los archivos de entrada pueden estar comprimidos con gzip y también pueden ser especificados con o sin la extensión GZ. Los archivos de salida serán comprimidos 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.

# produzca predicciones de dominios con HMMER3 (provée hmmscan) con filtros débiles. 
# este es el paso más lento 
perl -w 0runHmmscan.pl hmmscan Pfam-A.hmm PLAF7.fa PLAF7.p.txt 
# (opcional) comprima el archivo de salida, nuestros guiones lo leerán de todos modos 
gzip PLAF7.p.txt 

# remueva las superposiciones, ordenando por valor p, creando un archivo intermedio 
perl -w 1noOvs.pl PLAF7.p.txt PLAF7.po.txt Pfam-A.hmm.dat 

# añada columnas de valores q de dominio, FDR locales, y FDR|lFDR (sin filtrar) 
perl -w 2domStratStats.pl PLAF7.fa PLAF7.po.txt PLAF7.pod.txt 
#... y filtre con un contraste de valor q 
perl -w 2domStratStats.pl PLAF7.fa PLAF7.po.txt PLAF7.pod.q4e-4.txt 4e-4 
#... o filtre con un contraste de FDR local 
perl -w 2domStratStats.pl PLAF7.fa PLAF7.po.txt PLAF7.pod.l2.5e-2.txt 1 2.5e-2 
#... or filtre con un contraste de FDR a través de FDR locales iguales 
perl -w 2domStratStats.pl PLAF7.fa PLAF7.po.txt PLAF7.pod.fl4e-4.txt 1 1 4e-4 

# valores q con gradas (para valores p de secuencia y dominio), con contraste obligatorio 
perl -w 3tieredStratQ.pl PLAF7.fa PLAF7.po.txt PLAF7.pot.1e-4.txt 1e-4 
#... o filtro de contraste para los valores q de secuencia solamente 
perl -w 3tieredStratQ.pl PLAF7.fa PLAF7.po.txt PLAF7.pot.1e-4,1.txt 1e-4 1 

# corre todo el proceso con varios organismos, combinando las estadísticas 
# pase los prefijos ORG en lugar de ORG.fa u ORG.fa.gz 
# las salidas siguen patrones (vease abajo) 
perl -w 4allManyOrgs.pl hmmscan Pfam-A.hmm Pfam-A.hmm.dat \ 
	 Pf Pv Pk Py Pb Pc Bb Ta Tp Tg Nc Et Ch Cp Cm

Todos los ejemplos también fueron corridos con Superfamily 1.75 (archivos salida son PLAF7.s* en lugar de PLAF7.31.p*), y todos los archivos de salida están disponibles arriba. La única excepción son los datos para 4allManyOrgs.pl (vease arriba).

0runHmmscan.pl 1.01 - Adquiera predicciones de dominios para sus secuencias de proteína

perl -w 0runHmmscan.pl <hmmscan> <Pfam-A.hmm> <FASTA input> <output table>
Las entradas requeridas son

Este guión no corre simplemente a hmmscan con parámetros estándares, de hecho cambia varias opciones que son importantes para que DomStratStats funcione. En partcular, se asegura que las salidas reporten valores \(p\) en lugar de valores \(E\), y relaja los filtros heurísticos de valores \(p\) para que pasen más predicciones. Si se desean predicciones con \(q \le \mbox{4e-4}\) en UniRef50, encontramos que debemos conservar valores \(p\) tan grandes como 1e-5 (sin embargo, para la mayoría de las familias el contraste de valor \(p\) efectivo es mucho menor). Este guión fija el contraste de valor \(p\) más estricto a 1e-4 para asegurar que nada falte para el análisis de valores \(q\) estratificados para dominios.

Nótese que los filtros heurísticos de HMMER3 no fueron optimizados en el contexto de valores \(q\) con gradas, es posible que se puedan modificar para ese uso en particular (por ejemplo, no fijar el contraste final de valores \(p\) para dominios en HMMER3, que actualmente es igual al contraste de valores \(p\) de secuencia de 1e-4).

1noOvs.pl 1.02 - Remueve los dominios superimpuestos ordenando por valor \(p\)

perl -w 1noOvs.pl <input table> <output table> [<Pfam-A.hmm.dat>]
Las entradas requeridas son

Esta es una entrada opcional adicional, que es específica a Pfam.

Encontramos que la mejor manera de calcular FDRs y lFDRs correctos es primero removiendo las superposiciones entre dominios, ordenando por valor \(p\) (en la ausencia de una estadística mejor de comienzo, esto funciona muy bien). Asumimos que las superposiciones se removeran últimamente entre las mejores predicciones, como es práctica común entre los dominios del mismo clan de Pfam o superfamilia. Este guión remueve todas las superposiciones entre todas las familias, sin contar con relaciones conocidas entre los HMMs. Se desempeña excelentemente si el contraste final es razonable (por ejemplo, \(q \le \mbox{4e-4}\)), pero este paso llevará a estadísticas demasiado conservadoras (particularmente los lFDRs parecerán mucho peores de lo que son en realidad) para las predicciones más débiles.

Este guión usa una definición permisiva de superposición, en que sólo las superposiciones que excedan 40 aminoácidos o el 50% del largo del dominio más pequeño serán removidas. Estos parámetros están fijos en el código fuente, ¡modifíquelos si sabe lo que hace!

2domStratStats.pl 1.00 - Calcule y añada los valores \(q\), lFDRs y FDR|lFDR para dominios

perl -w 2domStratStats.pl <FASTA input> <input table> <output table> \
  [<q-value> <local FDR> <FDR|lFDR>]
Las entradas requeridas son

Los siguientes filtros de contraste opcionales se pueden fijar. Cada uno o todos se pueden usar simultaneamente, lo cual está permitido aunque sería raro y difícil de interpretar, así que evite usarlo así. Use 1 para cada filtro de contraste que no se quiera usar pero deba listar (por ejemplo, para solo filtrar el contraste de FDR|lFDR, indique 1 para los contrastes de valor \(q\) y lFDR).

Esta es la estrategia más básica para controlar el FDR. Si no está investigando la predicción de dominios, y sólo quiere predicciones con un FDR igual al de Pfam, fije un contraste de valor \(q\) (\(q \le \mbox{4e-4}\) será mejor). Considere a los lFDRs y FDR|lFDR como funcionalidades experimentales que merecen mayor investigación, al menos por el momento.

Nótese que el FDR|lFDR es monotónico con los lFDRs estratificados de dominio, pero no con los valores \(q\) estratificados para dominios. Cuando los lFDRs son precisos (¡es una condición importante!), provéen la manera óptima de controlar el FDR combinado a través del FDR|lFDR.

3tieredStratQ.pl 1.00 - Calcule y añada valores \(q\) para secuencias y dominios

perl -w 3tieredStratQ.pl <FASTA input> <input table> <output table> \
  <seq q-value> [<dom q-value>]
Las entradas requeridas son

Puede opcionalmente especificar el contraste para valores \(q\) de dominio que sea diferente del contraste para valores \(q\) de secuencia.

Esta estrategia "con gradas" consiste en usar un contraste para los valores \(q\) estratificados de secuencias, luego calculando valores \(q\) en el resto de los dominios (estos valores \(q\) estiman el FDR condicionalmente con el primer filtro de contraste de FDR de secuencias) y usando otro contraste. De esta manera, la información de los dominios que se repiten se puede usar, como Pfam lo hace excepto que aquí los contrastes se escojen automáticamente en lugar de ser escogidos por expertos. Sin embargo, como mencioné arriba, esta estrategia no controla el FDR con exactitud, así que los contrastes de valores \(q\) deberán ser escogidos usando evaluaciones de referencia.

4allManyOrgs.pl 1.00 - Obtenga las predicciones finales de dominios a partir de varios archivos de secuencias

perl -w 4allManyOrgs.pl <hmmscan> <Pfam-A.hmm> <Pfam-A.hmm.dat> <orgs>...
Las entradas requeridas son

Este guión complejo corre todo el proceso con varios organismos, produciendo tres tipos de salidas que uno podría querer comparar, específicamente el Pfam Normal (que usa los contrastes GA y otros pasos), las Estadisticas Estratificadas de Dominio (sin contrastes), y los Valores \(q\) Estratificados con Gradas (con un contraste fijo de 1e-4 por grada).

Para cada prefijo ORG indicado, el archivo de entrada se asume ser ORG.fa o ORG.fa.gz (ambos se leeran correctamente), y los archivos de salida serán (todos comprimidos con gzip automáticamente, la extensión GZ omitida abajo)

Nótese que, mientras la mayoría de los pasos se corren independientemente por organismo, todas las estadísticas estratificadas son estimadas combinando la distribución de valores \(p\) a través de los organismos (siempre estratificadas por familia de dominio), las cuales usualmente son mejores que estimaciones en un solo organismo. Además, las salidas permaneceran separadas por organismo, lo cual podría serle más conveniente que simplemente juntar todas las proteinas en un solo archivo de secuencias (lo cual produciría un solo archivo de salida para todos los organismos).

Este es el guión más inflexible de todos estos, pero demuestra el uso del API, el cual está sin documentar por el momento, particularmente el uso de múltiples archivos de entrada. Espero que usted lo modifique para sus necesidades, por ejemplo para cambiar los patrones de las salidas, o si sólo quiere las predicciones de dominios con un contraste de lFDR, o si quiere que los archivos intermedios se borren cuando el proceso esté completo.

Compatibilidad

Este código fue probado con HMMER versiones 3.0, 3.1b1, 3.1b2, contra las bases de datos HMM de Pfam 25, 27-31, y Superfamily 1.75.

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.

Cita

2015-11-17. Alejandro Ochoa, John 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

2017-03-13 - Retested for Pfam 31 (code unchanged)

Rehice los ejemplos incluidos en esta página usando este Pfam más nuevo, confirmando que mi cógido sirve.

2016-08-10 - Retested for Pfam 30 (code unchanged)

Rehice los ejemplos incluidos en esta página usando este Pfam más nuevo, confirmando que mi cógido sirve.

2015-12-23 - Retested for Pfam 29 and HMMER 3.1b2 (code unchanged)

Pfam 29 salió ayer. Rehice los ejemplos incluidos en esta página usando este Pfam más nuevo y también el HMMER más nuevo. Confirmé que mi cógido sirve para ellos.

2015-07-28 - DomStratStats 1.04

El código fue probado con Pfam 28, el cual salió en Mayo 2015. Todos los archivos muestra proveidos ahora son de esta versión. Cuando lo probé, encontré un error en el archivo Pfam-A.hmm.dat, cuyas asociaciones de anidación entre familias de dominio estaban en blanco, y esto causaba que mi lector de Pfam-A.hmm.dat se quejara mucho a STDERR. Ajusté mi código para que esto ya no suceda. Desafortunadamente, para el Pfam-A.hmm.dat actual de Pfam 28, no hay información de anidación disponible para usarse, aunque contacté al equipo de Pfam y puede que el error se arregle en el futuro (y mi código funcionará correctamente cuando eso suceda).

El código del FDR local ahora tiene una correción para un "caso de borde", en que datos sin observaciones ahora obtienen un lFDR de cero en lugar de morir por dividir por cero.

2014-10-24 - DomStratStats 1.03

(Sí, el mismo día que 1.02)

Corregí el código en Hmmer3ScanTab.pm que inserta las columnas de DomStratStats a la salida de HMMER3 para manejar casos de desplazamiento severo cuando los identificadores de proteína son demasiado largos. La salida del código corregido permanecerá desplazada pero las columnas insertadas van a aparecer correctas en cada fila, sin corromper las columnas originales. El código viejo asumía posiciones fijas de columnas excepto por desplazamientos muy pequeños, lo cual podía cometer errores de inserción. En algunos casos, el código viejo abortaba cuando no se identificaba la posición de inserción correcta. Gracias otra vez a Sesh A. Sundararaman por reportar este error.

2014-10-24 - DomStratStats 1.02

Corregí el guión 1noOvs.pl para que cargue Domains.pm en lugar de Hits.pm, el cual era el nombre viejo de ese paquete. Este error había prevenido la compliación del guión, así que no podía correr para nada. Gracias a Sesh A. Sundararaman por reportar este error.

2014-09-18 - DomStratStats 1.01

Añadí el guión 4allManyOrgs.pl

3tieredStratStats.pl se renombró a 3tieredStratQ.pl, ya que es una descripción más precisa (este guión no usa lFDRs, sólo valores \(q\)).

Reorganicé los paquetes internos, incluyendo crear y renombrar paquetes y funciones, para claridad adicional para cualquiera interesado en usar el API que está sin documentar. También añadí funciones no usadas por DomStratStats, para asegurar la compatibilidad de paquetes con proyectos relacionados (particularmente dPUC 2).

Corregí errores en la documentación. También hice más conciso a este manual.

Las salidas no cambian, dado que ningún algoritmo cambió.

2014-06-23 - DomStratStats 1.00

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

VIIIA

Historial