Detección de patógenos en datos de RNA-seq con Pathonoia | BMC Bioinformática

En primer lugar, describimos Pathonoia, que tiene como objetivo descubrir secuencias de ARN orgánico a partir de lecturas de secuenciación metagenómica de origen desconocido. Estas lecturas pueden provenir de la parte no alineada de los datos de secuenciación de ARN de, por ejemplo, muestras transcriptómicas humanas. Además, describimos nuestros métodos de evaluación comparativa.

En segundo lugar, describimos una tubería de análisis basada en la salida de Pathonoia. Aquí, el objetivo es distinguir las señales biológicas de la contaminación experimental. Proporcionamos dos estudios ejemplares. Esta canalización está disponible como plantilla de análisis guiado en un Jupyter Notebook en GitHub.

patonoia

Las lecturas que no son de mapeo de host se asumen como entrada para Pathonoia para detectar señales biológicas en datos de secuenciación de ARN ruidosos. SAMherramientas [35] se puede utilizar para montar el correspondiente rápidoQ archivos después de usar cualquier alineador (ver Archivo adicional 1: Secc. S1 para más detalles). Esta entrada se procesa en dos pasos principales: alineación metagenómica con Kraken 2 [11] y la agregación de k-mers de toda la muestra de Pathonoia.

Kraken identifica el ancestro común más bajo (LCA) para cada k-mer en una muestra. En Kraken 2, los desarrolladores optimizaron el tamaño del índice y la velocidad de clasificación mediante un esquema de minimización. Un minimizador es una subsecuencia de un k-mer con longitud yo. Se puede utilizar como representante del k-mer a costa de la precisión de la clasificación. Usamos longitudes minimizadoras y k-mer (l=k=31) para la configuración de mayor precisión de Kraken 2 y el índice Kraken 2 para todos los genomas virales y bacterianos en la base de datos NCBI nt. La precisión es menor para los minimizadores. yo con (l ya que podría haber dos k-mers de diferentes especies con el mismo minimizador. Pathonoia no puede funcionar correctamente con los índices de Kraken que utilizan (l. El LCA para un k-mer es la descripción taxonómica más específica de un conjunto de organismos que comparten esta secuencia. Utilizamos el identificador taxonómico (taxID) para la implementación del algoritmo. El resultado del paso de alineación es el kraken-alinear archivo, que contiene los k-mers clasificados para cada lectura de una muestra. Un formato de ejemplo de una lectura clasificada es:

C K025:418/1 Shamonda o.virus (taxid 15915) 100 0:20 15915:7 0:15 15915:6 1:5 0:2 230658:1 0:10

Los cinco campos son:

1. 1.

   C o tu si la lectura fue clasificada o no.

2. 2.

   El identificador de lectura como se indica en el rápidoQ archivo.

3. 3.

   El nombre taxonómico y taxID con el que se clasificó según el algoritmo de Kraken.

4. 4.

   La longitud de lectura original.

5. 5.

   Secuencia de pares en formato taxID:Y, donde Y es el número de k-mers consecutivos que se clasifican con el mismo identificador taxonómico.

En este ejemplo, no se pudieron identificar 20 k-mers (taxID = 0), siete k-mers pertenecen al Virus Shamonda orthobunya genoma, seguido de 15 k-mers no identificados y así sucesivamente. Al sumar el número de k-mers en una lectura, en este ejemplo (20+7+15+6+5+2+1+10 = 66)el resultado es siempre la longitud de lectura restada por (k+1)aquí (100-(31+1) = 66).

El segundo paso clave de Pathonoia es la interpretación de la salida de alineación de Kraken y separarla de la noción de lecturas. También se describe en el archivo adicional 1: Fig. SF1. Usando un hashmap, cada secuencia identificada de una muestra que es única, se almacena como clave y su ID de impuesto asignado como valor. Llamamos a estas secuencias (teclas) como x-mersya que tienen longitud (x = k + Y) (k: k-mer longitud, Y: número de k-mers consecutivos identificados con el mismo taxID).

A continuación, todas las longitudes X de secuencias de un mismo organismo (taxID) se resumen: (A_O = sum _^ x_i)con (x en X)el conjunto de norte distintos x-mers del organismo O. Solo organismos que superan un umbral de (nuestro valor predeterminado) (A_ = 100) los nucleótidos se consideran para el siguiente paso y el resultado final. Este umbral puede ser ajustado por el usuario. Cuanto más alto es el umbral, menos organismos se detectan. Archivo adicional 1: la Fig. SF2 muestra cómo la puntuación F1 puede aumentar con una disminución de las especies detectadas en un conjunto de datos simulado. Archivo adicional 1: la Fig. SF3 muestra cuánto puede disminuir el número de especies detectadas en una muestra in vitro con un umbral creciente. Para comparaciones justas con los otros algoritmos, seleccionamos un umbral relativamente bajo de (A_ = 100) racionalmente, correspondiente a al menos una lectura de mapeo completo o tres k-mers independientes.

Para aumentar la certeza sobre un organismo específico, la medida de abundancia (A_O) de cada ID de impuesto a nivel de especie se suma con la abundancia de niveles de género y familia, además: (A_O = A_S+A_G+A_F). La intuición detrás de esto es que si un organismo es de hecho parte de la muestra, se muestrean y procesan al azar varias áreas diferentes del mismo. Esto puede incluir áreas que no son específicas del genoma del organismo. Sin embargo, si una especie se puede detectar a nivel de especie, la evidencia se puede aumentar agregando recuentos de nivel más alto que deben provenir de una especie específica en cualquier caso. Finalmente, los organismos pueden clasificarse por su abundancia. (A_O). Sin embargo, todo el potencial se despliega al comparar grupos de muestra entre sí.

Evaluación comparativa de Pathonoia

El objetivo de nuestro algoritmo es la reducción de organismos detectados erróneamente en una muestra metagenómica y, por lo tanto, lograr la mayor precisión posible. La precisión se puede medir evaluando un conjunto de datos simulados, donde se conoce la presencia (y ausencia) de cada organismo. Además, para probar Pathonoia cualitativamente, evaluamos una muestra biológica (muestra GEO GSM1444167fibroblastos infectados por HVV [17]), donde la verdad fundamental de los falsos positivos difícilmente puede conocerse. Comparamos Pathonoia con las técnicas de medición de abundancia basadas en recuento de lecturas y Kraken 2. Queremos enfatizar en comparar la técnica de medir la abundancia y no herramientas específicas. Comparamos las adaptaciones de Kraken 2 para el objetivo de la indicación específica de organismos poco abundantes en una muestra metagenómica ruidosa.

Kraken 2 [11]Centrífuga [12]y helecho [23] representan las medidas de abundancia basadas en el recuento de lecturas. Kraken 2 y un índice de genomas bacterianos y virales con tamaño minimizador y longitud k-mer (l=k=31) son la línea de base para todos los algoritmos en este punto de referencia. Ejecutamos todos los demás algoritmos de medición de abundancia (excepto Centrifuge) sobre la salida de Kraken 2, lo que hace que los resultados de precisión sean comparables. No se puede detectar ningún organismo si Kraken 2 no lo detectó con al menos un k-mer (es decir, hay un valor máximo para TP y FP). Sin embargo, la forma de contar y evaluar las especies difiere en los distintos algoritmos. Nuestro punto de referencia incluye “Kraken 2 con cut-off (ge 5) lee”, lo que significa que un organismo cuenta como detectado cuando Kraken 2 identifica al menos cinco lecturas que se originan en ese mismo organismo. La métrica de abundancia pura de Kraken 2 cuenta un organismo como detectado, si al menos una lectura se clasifica con ella. Bracken es otra herramienta basada en Kraken que corrige las medidas de abundancia, incluida la distribución estadística de los genomas disponibles en la base de datos subyacente. Pathonoia puede detectar organismos con los que no se clasificó ninguna lectura, ya que evalúa en el nivel k-mer.

Para evaluar el rendimiento, utilizamos siete muestras simuladas, que fueron construidas por Ye et al. [16] para la evaluación comparativa de clasificadores de taxonomía en metagenómica. Son muestras que contienen secuencias de ADN de doce a 50 organismos que se encuentran en microbiomas humanos comunes o metagenomas ambientales, por ejemplo, el intestino humano o el hogar.

Nos oponemos a las técnicas comunes de medición de precisión, que indican si las lecturas se identificaron correctamente o no. En cambio, nos enfocamos en la presencia (in)correcta de especies (Fig. 2C) y definimos falsos positivos (FP), verdaderos positivos (TP) y falsos negativos (FN) en consecuencia. Llamando la atención sobre esta clase, medimos la precisión (()), recordar (()) y su media armónica, la puntuación F1 (()). La figura 2D muestra la precisión promedio, la recuperación y la puntuación F1 en las siete muestras (valores detallados en el archivo adicional 1: tablas ST1 y ST2).

Análisis aguas abajo

La aplicación de Pathonoia en varias muestras en un conjunto de datos da como resultado una matriz de datos que contiene para cada muestra las medidas de abundancia (A_O) de todos los organismos encontrados en al menos una de las muestras. Estos datos de abundancia, así como algunos metadatos sobre las muestras, son la entrada para el análisis posterior. El objetivo es identificar las diferencias de abundancia patógena entre grupos de muestras, por ejemplo, entre muestras enfermas y de control. Proporcionamos una secuencia de comandos de análisis de plantilla en línea a la que nos referimos como cuaderno Jupyter de “análisis guiado”. A continuación, demostramos el flujo de trabajo de este análisis ejemplar en dos conjuntos de datos (compárese con la Fig. 3).

Los conjuntos de datos utilizados para el estudio de casos comprenden 48 y 63 muestras. El primero es de un estudio interno de Demencia Frontotemporal (FTD) [25], observando muestras de tejido cerebral divididas en enfermedad y control de acuerdo con la Fig. 3A. El segundo conjunto de datos (Fig. 3E) contiene muestras de hígado de pacientes en diferentes etapas fibróticas del hígado debido a una de cuatro enfermedades: hepatitis autoinmune (AIH), enfermedad del hígado graso no alcohólico (NAFLD), colangitis biliar primaria (PBC) y esclerosante primario Colangitis (PSC).

Se ejecuta un análisis de componentes principales (PCA) como primer paso de la tubería de análisis. Sirve como una “verificación de valores atípicos”, donde se pueden identificar muestras que pueden estar demasiado contaminadas. Además, el gráfico de PCA está coloreado por varios metadatos disponibles (archivo adicional 1: Fig. SF6) para identificar si un determinado entorno experimental puede provocar una contaminación importante. En este paso, las muestras atípicas deben eliminarse del análisis (manualmente). En el estudio FTD extrajimos 3 muestras y en el estudio de fibrosis se excluyeron 12 muestras. (Archivo adicional 1: Fig. SF5) Después de volver a ejecutar el PCA para FTD y colorear según la edad, el sexo, la celda de flujo y la puntuación RIN, no se pudo identificar ningún sesgo en función de estos metadatos.

El paso de comparación por grupos incluye el cálculo de la abundancia media por organismo para todas las muestras de un grupo, por ejemplo, todas las muestras de control. DESeq2 [36] se utiliza para este análisis de abundancia diferencial. Originalmente, esta herramienta fue desarrollada para el análisis diferencial de expresión génica. Sin embargo, puede usarse para nuestros propósitos ya que el mismo modelo estadístico, una distribución binomial negativa, puede asumirse para nuestros datos. Como argumentaron previamente Simon et al. [8] y Rahman et al. [7] la fuente de datos subyacente, los datos de RNA-seq, es la misma en nuestro caso y el gen cuenta. Características metagenómicas, como nuestra (A_O), son como conteos de genes comúnmente caracterizados por una gran cantidad de ceros y pocos aciertos de conteo bajo. Archivo adicional 1: la Fig. SF4 muestra esta distribución en uno de nuestros conjuntos de datos de ejemplo.

Para aumentar la estabilidad de DESeq2, los organismos que tienen abundancia cero en la mayoría de las muestras se excluyen del análisis. Esto generalmente se hace también para el análisis del transcriptoma. [37]. El resultado de este paso es una lista de organismos junto con su valor de cambio de veces log2 entre los grupos de muestra y un valor p (prueba de Wald) ajustado para pruebas múltiples (con Benjamini-Hochberg). Un gráfico de volcanes coloreado por el número de muestras en las que está presente un organismo ofrece una descripción general de estos resultados. (Resultados para el conjunto de datos de FTD y fibrosis en la Fig. 3 B, F y Archivo adicional 1: Tablas ST3 y ST7)

La visualización de la abundancia es el siguiente paso. La abundancia se traza por grupo de muestra para los organismos más significativos. Aquí, se puede identificar si un organismo juega un papel en todo el grupo de muestra, o si la media solo se elevó debido a un caso extremo. El último caso debería ocurrir con menos frecuencia si se eliminaran más valores atípicos en el primer paso. En ambos estudios, observamos que la mayoría de los organismos muestran un aumento constante en las muestras enfermas pero no en los controles (Archivo adicional 1: Figs. SF8 y SF10). El resultado de este paso es la selección de un organismo de interés. En el estudio FTD, seleccionamos Burkholderia stabilis, por ser el organismo más prevalente de los diferencialmente abundantes. Su aparición en las muestras se visualiza en el gráfico PCA en el archivo adicional 1: Fig. SF7.

Un análisis de expresión génica diferencial puede ayudar a comprender si el organismo de interés tiene un origen biológico. Aquí, los datos del transcriptoma primario se utilizan para comparar muestras que contienen el organismo y muestras que no contienen el organismo. Estos grupos de muestra pueden ser subgrupos de otro grupo de muestra. Por ejemplo, en el conjunto de datos de FTD, solo seleccionamos muestras de pacientes para esta comparación, para comprender qué diferencia B. estabilizador podría hacer en el caso enfermo (y dado que no está presente en ningún control, también). Usando DESeq2 nuevamente, en los recuentos de lectura originales, recuperamos un conjunto de genes regulados hacia arriba y hacia abajo. Idealmente, este conjunto concluye como un efecto de la presencia del organismo de interés. En el estudio FTD, este paso resultó en un conjunto de 143 significativamente (valor p-adj. (< 0.05)) genes expresados ​​diferencialmente. (Más detalles en Archivo Adicional 1: Secc. S3)

Realizamos un análisis de ontología de genes para comprender el efecto del organismo de interés, utilizando un análisis de enriquecimiento de sobrerrepresentación (ORA) de los conjuntos de genes con WebGestaltR [24]. (Para conocer la configuración exacta, consulte Archivo adicional 1: Secc. S2.1). Para el estudio FTD, los 34 genes regulados positivamente se compararon con los conjuntos de genes en la base de datos de procesos biológicos (archivo adicional 1: tabla ST6).