PARiS: Asignación y redistribución probabilística de secuencias de isomiR: un método basado en datos para eliminar el ruido de los datos de recuento de lecturas de isomiR.

Los microARN (miARN) son una subclase de pequeñas moléculas de ARN no codificantes, de cadena sencilla, de aproximadamente 18 a 24 nucleótidos de longitud, con importantes funciones biológicas. Los miARN regulan la expresión génica (^[1]), controlan la diferenciación celular (^[2]) e incluso se utilizan para la comunicación intercelular (^[3]). La expresión de los miARN debe estar estrictamente controlada para el crecimiento y desarrollo saludables de las células (^[1]). Una expresión aberrante de los miARN a menudo da como resultado células enfermas. Por ejemplo, los cambios en la expresión de los miARN se han asociado con muchos tipos diferentes de cáncer (^[6], ^[7]), enfermedades autoinmunes (^[10]), enfermedades cardiovasculares (^[2]) y enfermedades neurológicas, como la enfermedad de Huntington (^[3], ^[13]).

Los perfiles de expresión de los miARN se estudian típicamente utilizando la secuenciación de ARN pequeño (sRNA-seq), que combina tecnologías de secuenciación de alto rendimiento con un paso de selección de tamaño para generar lecturas de secuencias de miARN en diferentes condiciones experimentales. La secuenciación de alto rendimiento ha reemplazado en gran medida la cuantificación basada en microarrays y qPCR debido a su capacidad para identificar nuevos miARN y diferenciar entre secuencias a nivel de un solo nucleótido. El uso de sRNA-seq condujo a la detección de nuevas isoformas de miARN, llamadas isomiRs. Aunque inicialmente se pensó que estos isomiRs eran el resultado de errores técnicos en la secuenciación y/o la alineación, estudios posteriores demostraron que los isomiRs eran demasiado abundantes en las muestras como para ser producidos exclusivamente por errores de secuenciación, lo que sugiere que los isomiRs deben ser productos de vías de biosíntesis en varios organismos (^[4]). Desde este descubrimiento, los investigadores han demostrado que los isomiRs son moléculas biológicamente funcionales (^[5]). Además, los estudios han relacionado la expresión no solo de los miARN, sino de isomiRs específicos con diversas enfermedades, como la obesidad, la diabetes, la enfermedad de Alzheimer y el cáncer (^[6], ^[7]).

Los isomiRs se clasifican en función de cómo difieren de una secuencia de miARN canónica o de referencia informada en una base de datos de miARN, como miRBGeneDB (^[8]) o miRBase (^[9]). Se desarrolló una nomenclatura detallada para describir los isomiRs por el miRNA Transcriptomic Open Project (miRTOP) (^[10]). La adición o eliminación de nucleótidos del extremo 5' o 3' de la secuencia de referencia produce un isomiR 5' o un isomiR 3', respectivamente. En este manuscrito, nos referimos a los isomiRs que resultan de tales cambios como isomiRs de longitud variante. La sustitución de nucleótidos dentro de la secuencia de referencia produce un isomiR polimórfico. Nos referiremos a los isomiRs polimórficos en este manuscrito como isomiRs de secuencia variante. Finalmente, los isomiRs de tipo mixto son el resultado de diferencias tanto de longitud como de secuencia (^[11]).

Aunque los estudios han demostrado que algunas secuencias de isomiRs se producen biológicamente y son funcionales, las secuencias de isomiRs no se producen exclusivamente de forma biológica. Durante el proceso de secuenciación, los nucleótidos pueden agregarse, eliminarse o sustituirse erróneamente de una secuencia verdadera (^[19]). Estos errores aleatorios producen isomiRs técnicos no deseados en los datos de secuenciación. Además, estas variantes técnicas son indistinguibles de los isomiRs producidos biológicamente. No eliminar estos isomiRs técnicos de los datos de secuenciación de miARN antes del análisis puede afectar el cálculo de los niveles de expresión de isomiRs y la estimación de la expresión diferencial (^[12]). Esto sugiere la necesidad de eliminar el ruido de los datos de isomiR antes del análisis.

A pesar de la importancia de los isomiRs, el enfoque más utilizado para analizar los datos de secuenciación de miARN es agregar los datos de recuento a nivel de isomiR al nivel de miARN después de la alineación. La agregación se realiza sumando los recuentos de lecturas para todas las secuencias que se asignan al mismo miARN. Después de la agregación, la expresión diferencial de miARN se estima utilizando métodos desarrollados inicialmente para datos de secuenciación de ARN mensajero (ARNm), como edgeR (^[13]) o DESeq2 (^[14]). Al agregar tanto los isomiRs biológicos como los técnicos, este enfoque reduce sustancialmente el ruido no deseado introducido por los isomiRs técnicos. Sin embargo, la agregación de secuencias de isomiR al nivel de miARN da como resultado una pérdida de información y potencialmente oscurece importantes diferencias funcionales. Además, estudios recientes han demostrado que el uso de datos de expresión a nivel de isomiR para inferir la expresión diferencial a nivel de miARN puede mejorar la inferencia (^[23], ^[15]).

Un enfoque alternativo para manejar la variación técnica no deseada en los datos de secuenciación a nivel de isomiR se implementa en miREC, una herramienta de rectificación de errores de miARN (^[12]). MiREC utiliza un enfoque de red k-mer para identificar y eliminar las secuencias de error de los datos. Aunque el modelo está equipado para manejar tanto los errores de sustitución, lo que da como resultado isomiRs de secuencia variante técnicos, como los errores de adición/eliminación, lo que da como resultado isomiRs de longitud variante técnicos, la evaluación del método se centró en su capacidad para eliminar los errores de sustitución (^[12]). Actualmente, no existen estudios comparativos que demuestren que la inclusión de miREC como primer paso en una canalización analítica para la expresión diferencial de miARN mejora el análisis posterior. Además, miREC requiere la selección de hiperparámetros que pueden alterar sustancialmente el rendimiento.

La necesidad de eliminar el ruido de las lecturas de secuenciación no es exclusiva de los datos de secuenciación de miARN. En la investigación del microbioma, la secuenciación de amplicones se utiliza para cuantificar las variantes de un gen marcador específico, a menudo el gen de ARNr 16S. Las lecturas de la secuencia del gen marcador se agrupan en unidades taxonómicas operativas (OTU) como primer paso en una canalización analítica para los análisis de abundancia diferencial. La agregación de las lecturas de la secuencia del gen marcador al nivel de OTU permite la eliminación de la variación técnica no deseada de los datos, pero impide el análisis de los datos a una resolución más fina que las OTU. Para resolver esto, se desarrollaron varios métodos para eliminar el ruido de los datos del microbioma: Deblur (^[16]), DADA2 (^[17]) y UNOISE2 (^[18]). Todos estos métodos infieren un conjunto de secuencias verdaderas del gen marcador y recuentos de lecturas asociados, aunque la forma en que los métodos determinan la pertenencia al conjunto de secuencias verdaderas del gen marcador difiere. En el contexto de la secuenciación de miARN, un miARN dado sería análogo a un gen marcador, y los isomiRs de ese miARN serían análogos a las OTU. Sin embargo, ninguno de los métodos desarrollados para los datos de secuenciación de amplicones microbianos se puede aplicar directamente a los datos de secuenciación de miARN porque se centran únicamente en las variantes de secuencia. De hecho, Deblur, UNOISE2 y DADA2 implementan un paso de filtrado o recorte para eliminar las variantes de longitud de los datos. Se necesita desarrollar un enfoque novedoso y basado en datos para eliminar el ruido de las secuencias de isomiR que pueda manejar todos los tipos de secuencias de isomiR.

Para mejorar el análisis de los datos de secuenciación de miARN, proponemos PARiS, un algoritmo para la asignación probabilística y el re-particionamiento de secuencias de isomiR. PARiS modela el proceso de generación de lecturas de isomiR de error a partir de secuencias verdaderas de isomiR en función de los datos observados y lo utiliza para definir una partición de secuencias para cada miARN. Para un miARN dado, cada elemento de la partición contiene una única secuencia de isomiR verdadera inferida y el conjunto de todas las secuencias de isomiR de error producidas por la variación técnica alrededor del isomiR verdadero. Los recuentos de lecturas de isomiR eliminados del ruido se producen sumando las lecturas dentro de cada elemento de una partición. Demostramos la eficacia de PARiS utilizando una combinación de simulaciones, un conjunto de datos de referencia sintéticos experimentales y tres líneas de células de adenocarcinoma de colon.

Materiales y métodos

Un experimento de perfilado de miARN consiste típicamente en dos o más condiciones experimentales, con múltiples muestras dentro de cada condición experimental. Sea referirse a las muestras (biológicas) únicas en un experimento de perfilado de miARN. Cada muestra consiste en secuencias de isomiR, , con refiriéndose al recuento de lecturas observado de la secuencia de isomiR en la muestra . Cada secuencia se ha asignado a un miARN , donde cada . Asumimos que hay una o más secuencias de isomiR que se asignan a cada miARN . La asignación de la secuencia de isomiR al miARN se da por , donde la función está determinada por la elección del software de alineación y la base de datos de referencia. Organizamos los recuentos de lecturas de la secuencia de isomiR en una matriz -dimensional, . Aquí tiene la misma definición que arriba y es el número total de secuencias de isomiR únicas en el experimento de perfilado. Organizamos la información adicional de la alineación de las secuencias de isomiR a la base de datos de referencia en otra tabla de datos. La tabla de datos con la información de alineación tiene filas y hasta 3 columnas. La primera columna es para cada secuencia de isomiR única y la segunda columna es para la secuencia de miARN a la que se ha asignado la secuencia por la función . Algunos algoritmos de alineación devuelven el tipo de coincidencia realizada entre una secuencia y la secuencia canónica para el miARN en la base de datos de referencia. Por ejemplo, miRge indica una coincidencia exacta entre la secuencia de isomiR y el miARN canónico, independientemente de la longitud, o un isomiR si existe una variante de secuencia (^[19]). El modelo de error que presentamos utiliza y la información en la tabla de datos de alineación.

Modelo de error.

Una secuencia de isomiR de error , también conocida como secuencia de isomiR de error técnico, no aparece en el sistema biológico real del que se tomó la muestra. Más bien, la variación no deseada en el proceso de secuenciación resultó en errores a lo largo de la lectura de la secuencia de isomiR biológico , produciendo la secuencia . Supongamos que tenemos una secuencia de isomiR de error que se asigna al miARN en la muestra de un experimento de perfilado de miARN con el recuento de lecturas asociado . También se asigna al miARN en la muestra una secuencia de isomiR biológico verdadero, , con el recuento de lecturas asociado . El modelo de error, incluida la estimación de parámetros y la prueba de hipótesis, se basa en una alineación por pares entre el isomiR biológico y la secuencia . Aquí, utilizamos la alineación por pares para significar los resultados del uso del algoritmo de Needleman-Wunsch para obtener una alineación global de dos secuencias de nucleótidos (^[30]).

Ciertos aspectos de nuestro modelo de error están inspirados en el modelo DADA2 para eliminar el ruido de las variantes de secuenciación de amplicones (^[17]). Primero, asumimos que los errores de lectura de la secuencia ocurren de forma independiente dentro de una lectura de una secuencia de isomiR y entre lecturas. A partir de estas suposiciones, se deduce que la probabilidad de producir una lectura de a partir de debido a errores técnicos se puede expresar como el producto de las probabilidades de las transiciones de nucleótidos individuales a lo largo de la longitud de la alineación por pares de a . Sea denotar la probabilidad de producir una lectura de la secuencia de error técnico a partir de la secuencia . Sea denotar la probabilidad de que el nucleótido en la posición de de la alineación por pares se lea como el nucleótido en la posición de de la alineación por pares. Finalmente, sea la longitud de la alineación por pares entre y . Entonces, a partir de nuestras suposiciones, tenemos:

La probabilidad de transición de un carácter a otro carácter se define por la identidad de esos caracteres. Nos referimos a estos caracteres individuales de la alineación por pares entre las secuencias y como caracteres y no solo como nucleótidos porque permitimos que estos caracteres sean nucleótidos o un espacio: . Esto nos permite manejar las transiciones de nucleótido a nucleótido y de espacio a nucleótido o de nucleótido a espacio simultáneamente.

Cada uno de estos diferentes tipos de transiciones que permitimos representa un tipo diferente de isomiR. Las transiciones de nucleótido a nucleótido representan errores de sustitución. Un alineamiento por pares entre y que contenga solo transiciones de nucleótido a nucleótido indica que es una variante de secuencia isomiR de . De manera similar, las transiciones de brecha a nucleótido representan errores de adición y las transiciones de nucleótido a brecha representan errores de deleción. Un alineamiento por pares entre y que contenga solo transiciones de brecha a nucleótido o de nucleótido a brecha indica que es una variante de longitud isomiR de . Las variantes de longitud, particularmente en el extremo 3′, son comunes en la secuenciación de miARN (miRNA-seq) debido a la heterogeneidad del corte 5′/3′ y a la adición de una cola 3′ (^[20]); en la secuenciación de amplicones de microbiomas, los cebadores conservados, el recorte de calidad y la fusión de lecturas de doble extremo producen longitudes de amplicones casi uniformes (^[17]). Finalmente, una secuencia con los tres tipos de transiciones es un isomiR de tipo mixto. Organizamos estas probabilidades de transición en una matriz, . Las probabilidades de transición pueden estar preespecificadas o estimadas directamente a partir de los datos. Los detalles sobre la metodología utilizada para la estimación de parámetros se encuentran en los Materiales Suplementarios.

A continuación, denotemos una variable aleatoria que representa el recuento de lecturas de la secuencia de error en la muestra . Sea la abundancia real de la secuencia biológica de isomiR en la muestra . El valor es una variable aleatoria latente. Usamos para denotar su recuento de lecturas observado:

En cada muestra , cada molécula de la secuencia se selecciona para la secuenciación con probabilidad . Por lo tanto, el número de lecturas observadas que se originan en sigue la siguiente distribución:

Condicionado a que se origine en , una lectura se observa como secuencia de error con probabilidad , determinado por el modelo de error Eq. Eq. (1). Por la propiedad de adelgazamiento de Poisson, el recuento de lecturas asignadas a sigue:

o, equivalentemente, condicionado a :

Por lo tanto, consideramos el proceso de secuenciación para como proveniente de dos adelgazamientos de Poisson: 1) selección de lecturas del recuento real basado en la probabilidad de muestreo , y 2) asignación errónea de estas lecturas a la secuencia de error con probabilidad .

Prueba de hipótesis.

Utilizamos el modelo de error descrito anteriormente para probar la hipótesis nula de que la secuencia que se asigna al miARN es una secuencia de isomiR de error de la secuencia biológica de isomiR . Específicamente, comparamos el recuento de lecturas observado de la secuencia en la muestra con lo que esperaríamos según el modelo de error de la Ecuación 5. Si es mucho más abundante de lo que esperaríamos según el modelo de error, rechazamos la hipótesis nula e inferimos que es otra secuencia biológica de isomiR del miARN . Formalizamos la idea de ser "mucho más abundante de lo que esperaríamos" tomando prestado el concepto de un valor p de abundancia del algoritmo DADA2 (^[17]). El valor p de abundancia se define matemáticamente como:

donde representa la función de densidad de Poisson. El valor p de abundancia es la probabilidad de observar un recuento de lecturas para la secuencia en la muestra tan extremo o más extremo que , asumiendo que es una secuencia de isomiR de error de la secuencia biológica de isomiR . El criterio de decisión para rechazar la hipótesis nula compara el valor p de abundancia, , con un umbral definido por el usuario, funciona de manera similar al nivel en un marco de prueba de hipótesis frecuentista tradicional y se establece comúnmente en 0,05.

En la práctica, calculamos un valor p de abundancia para cada secuencia de isomiR tal que y lo comparamos con . Para controlar el número de falsos descubrimientos que cometemos, corregimos los valores p de abundancia para pruebas múltiples utilizando el procedimiento de Benjamini-Hochberg (^[31]). Usaremos para representar el valor p de abundancia ajustado de Benjamini-Hochberg.

PARiS: Asignación y Repartición Probabilística de Secuencias de IsomiR.

El algoritmo PARiS toma, como entrada, una matriz de dimensión de recuentos de lecturas a nivel de isomiR, . El primer paso del algoritmo PARiS es estimar las probabilidades de transición almacenadas en utilizando los métodos de estimación de parámetros descritos en los Materiales Suplementarios. PARiS es un algoritmo iterativo que particiona el conjunto de secuencias que se asignan al miARN en la muestra a través de un proceso de refinamiento secuencial. Eliminamos la notación para y en lo siguiente. Definimos la partición inicial por:

donde para algún miARN .

Usamos para indexar las iteraciones del algoritmo PARiS. En la iteración , sea una partición de en subconjuntos disjuntos (no vacíos), es decir:

Cada bloque tiene una secuencia central designada , que representa un isomiR biológico.

En la iteración , obtenemos una partición más fina dividiendo el bloque . Los pasos son los siguientes: i) aplicar el modelo de error de la Ecuación 5 a cada , tomando para calcular el valor p de abundancia ; ii) después de calcular , corregimos para pruebas múltiples utilizando el procedimiento de Benjamini-Hochberg (^[31]) para generar ; iii) comparamos con , el umbral de significancia definido por el usuario; iv) identificamos una nueva secuencia central, , por:

Después de identificar , el conjunto se da por:

Este proceso iterativo se lleva a cabo de forma similar a una cascada, donde las secuencias "fluyen" de a , y así sucesivamente, hasta que no se puedan formar más subconjuntos o se alcance un número máximo de iteraciones preespecificado.

Después de que se ha asignado cada , hemos definido una partición de secuencias de isomiR que se asignan al miARN en la muestra . Nos referiremos a como la partición a nivel de muestra del miARN en la muestra o , omitiendo la dependencia del miARN cuando quede claro por el contexto.

Formando una partición de consenso y eliminando el ruido de los recuentos.

Después de obtener particiones por muestra para el miARN , las consolidamos en una única partición coherente definiendo primero un conjunto de consenso de secuencias centrales. Sea denotar los conjuntos de secuencias centrales identificadas para el miARN en cada una de las muestras. Definimos el conjunto de consenso de secuencias centrales para el miARN como .

Intuitivamente, contiene las secuencias centrales que aparecen en una proporción suficiente de muestras. Dos extremos ilustran esta idea:

Tomar la unión es permisivo y conlleva el riesgo de falsos positivos (incluyendo isomiRs de error como centros), mientras que tomar la intersección es estricto y conlleva el riesgo de falsos negativos (excluyendo isomiRs biológicos verdaderos).

Para equilibrar estos extremos, introducimos un umbral de prevalencia : una secuencia central candidata se retiene en si aparece en al menos una fracción de las muestras. Establecer produce la unión, mientras que produce la intersección. El conjunto de consenso resultante define los bloques de la partición de consenso, que luego utilizamos para eliminar el ruido de los recuentos en las muestras.

Después de identificar el conjunto central de consenso , reasignamos cada secuencia no central de la siguiente manera. Para cada centro , calcule como en la Ecuación 1, y asigne al bloque cuyo centro obtiene el valor más grande. El número total de subconjuntos para se da por la cardinalidad de .

Después de construir la partición de consenso , eliminamos el ruido de los recuentos agregando lecturas dentro de los bloques. Sea denotar el centro del bloque . Para la muestra , el recuento de ruido eliminado para la secuencia de isomiR biológico putativo es:

La composición inferida para el miARN se da por sus centros de consenso , con abundancias de ruido eliminado en las muestras.

Simulación de datos de la distribución de errores.

Comenzamos evaluando la capacidad de PARiS para identificar y eliminar las secuencias de isomiR de error técnico sin introducir nuevos errores. Para hacer esto, necesitamos un conjunto de datos donde, para cada secuencia de isomiR , sepamos si es un isomiR biológico o un isomiR de error. Generamos nuestros propios datos de secuenciación de miARN con recuentos de lecturas asociados y secuencias de isomiR de error verdaderas a partir de un conjunto de datos de referencia experimental que secuenciaba los miARNs del ratón. El conjunto de datos de referencia experimental que utilizamos consistió en 421.096 secuencias que se asignaban a 759 miARNs distintos después de la alineación a miRGeneDB (^[8]) utilizando sRNAbench (^[21]). Filtramos los miARNs de baja expresión con recuentos medianos por millón menores que 5. Después de filtrar, quedaron 435 miARNs. Para cada miARN, identificamos , la secuencia más abundante que se asigna al miARN en todas las muestras. Luego, para cada miARN, simulamos la creación de adiciones, deleciones o sustituciones en la cadena que representa para producir secuencias de isomiR a partir de la secuencia central. Luego, simulamos recuentos de lecturas a partir del modelo de error dado en la Ecuación 5 para generar recuentos de lecturas de las secuencias de isomiR simuladas. Los detalles de la simulación de datos se pueden encontrar en los Materiales Suplementarios.

Variamos el número de secuencias de isomiR biológico verdaderas que se asignan a cada miARN seleccionando el valor del conjunto {1, 4, 9}. Luego, para un número verdadero dado de secuencias de isomiR biológico, repetimos el proceso de generación de secuencias de isomiR simuladas y simulamos recuentos de lecturas a partir del modelo de error para generar un conjunto de datos simulado 10 veces.

Conjuntos de datos simulados de monocitos.

Para demostrar que PARiS mejora el análisis de expresión diferencial a nivel de miARN, necesitamos un conjunto de datos con varias características clave. Primero, el conjunto de datos debe tener recuentos de lecturas de secuencias de isomiR. En segundo lugar, el conjunto de datos debe tener valores verdaderos para la cantidad de expresión diferencial para cada uno de los miARNs en el conjunto de datos. La expresión diferencial se mide típicamente mediante el cambio logarítmico en la expresión (logFC) entre las condiciones experimentales, por lo que debemos conocer el logFC verdadero para cada miARN o poder calcularlo.

Para generar un conjunto de datos que cumpla con estos requisitos, utilizamos el mismo proceso utilizado por Baran et al. para inyectar artificialmente una señal de expresión diferencial en un conjunto de datos real (^[15]). El conjunto de datos inicial consistió en 39 muestras de monocitos. Después de filtrar los miARNs de baja expresión y eliminar las secuencias que no eran una coincidencia exacta con la base de datos de referencia, el conjunto de datos inicial consistió en 122 miARNs que se asignaban a 3538 secuencias únicas. A partir de este conjunto de datos inicial, se crearon 50 conjuntos de datos simulados utilizando el siguiente proceso: las 39 muestras se dividieron en dos grupos (A y B) al azar para imitar un diseño experimental simple (^[15]). De los 122 miARNs, se seleccionaron 20 para que se sobreexpresaran en el grupo A al azar y se seleccionaron 20 para que se sobreexpresaran en el grupo B al azar. Los recuentos de lecturas a nivel de isomiR de los miARNs seleccionados para que se sobreexpresaran en cada grupo se multiplicaron por valores muestreados de una distribución normal truncada con una media de 2 y una desviación estándar de 1. Más detalles sobre este proceso se pueden encontrar en Baran et al. (^[15]). El proceso de simulación se repitió para generar 50 conjuntos de datos sintéticos, con valores de logFC conocidos para cada miARN y cada secuencia de isomiR en el conjunto de datos. No observamos si cada secuencia de isomiR es una secuencia de isomiR biológico o una secuencia de isomiR de error.

Aplicamos diferentes flujos de trabajo analíticos, que consisten en un método de eliminación de ruido seleccionado del conjunto de {PARiS, miREC, Agregación o Ninguno} y un método apropiado para estimar la expresión diferencial a nivel de miARN dado la resolución de los datos de ruido eliminado. Para el flujo de trabajo analítico que utiliza PARiS, utilizamos como el umbral para identificar valores p de abundancia significativos. Variamos el valor de para crear conjuntos de consenso de secuencias centrales para cada miARN seleccionando del conjunto {0,025, 0,05, 0,10, 0,20, 0,40, 0,60, 0,80, 1,00}. PARiS y miREC producen datos de ruido eliminado a nivel de isomiR, y los datos brutos (la opción de eliminación de ruido Ninguno) también están a nivel de isomiR. Estimamos la expresión diferencial a nivel de miARN utilizando los datos de ruido eliminado de cada uno de estos métodos utilizando miRglmm, un marco de modelado lineal mixto generalizado (^[15]). El método de agregación produce datos de ruido eliminado a nivel de miARN. Estimamos la expresión diferencial a nivel de miARN a partir de los datos agregados utilizando DESeq2 (^[14]) y edgeR (^[13]). En total, aplicamos 11 flujos de trabajo analíticos diferentes a cada uno de los 50 conjuntos de datos simulados.

Conjunto de datos de referencia ERCC.

Evaluamos el rendimiento de PARiS en una secuencia de procesamiento analítico para el análisis de expresión diferencial de miRNA en un conjunto de datos de referencia sintético y experimental. Debido a la naturaleza sintética del conjunto de datos, no hay variación biológica presente. Utilizamos el conjunto de datos de referencia del Consorcio de Comunicación de ARN Extracelular (ERCC) (^[22]). El conjunto de datos ERCC es una secuenciación de ARN pequeño (sRNA-seq) de grupos racionales (A y B) de ARN pequeños sintetizados en proporciones de 10:1 a 1:10. El resultado es un conjunto de datos de expresión de secuencias de isomiR con 15 niveles de expresión diferencial real, que van desde log(0.1) = −2.3 hasta log(^[10]) = 2.3. El conjunto de datos contiene 286 miRNA humanos que se mapean a 8001 secuencias después de eliminar las secuencias que son ≥ 40 nucleótidos o < 16 nucleótidos.

Realizamos una limpieza del conjunto de datos ERCC utilizando un conjunto similar de métodos que se aplicaron a los conjuntos de datos de monocitos simulados. Los métodos de limpieza que aplicamos al conjunto de datos ERCC son los siguientes: miREC, agregación y . Para cada aplicación de PARiS, establecemos . En conjunto, esto da un total de 22 secuencias de procesamiento analítico diferentes. Después de la limpieza en cada secuencia de procesamiento, estimamos la expresión diferencial a nivel de miRNA utilizando el método apropiado dado la resolución de los datos limpios. Para los datos corregidos por miREC y limpios por PARiS, estimamos la expresión diferencial a nivel de miRNA utilizando miRglmm (^[15]). Para los datos agregados, estimamos la expresión diferencial a nivel de miRNA utilizando DESeq2 (^[14]) y edgeR (^[13]).### Comparación de líneas celulares de adenocarcinoma colorrectal.

Aplicamos PARiS como el paso de limpieza en una secuencia de procesamiento para el análisis de expresión diferencial de miRNA aplicada a un conjunto de datos experimental real. Los datos de sRNA-seq (N=9) provienen de 3 líneas celulares de adenocarcinoma colorrectal (DLD-1, DKO-1 y DKS-8), que varían según el estado de KRAS (^[23]). Los datos de secuenciación brutos consisten en 5223 secuencias únicas que se mapean a 143 miRNA después de filtrar las miRNA con baja expresión con lecturas por millón < 5 y mantener solo las secuencias que se alinean exactamente con la base de datos de referencia. Aplicamos PARiS con y y luego estimamos la expresión diferencial a nivel de miRNA utilizando miRglmm (^[15]). También agregamos los recuentos a nivel de isomiR al nivel de miRNA y estimamos la expresión diferencial utilizando DESeq2 con fines de comparación (^[14]). Identificamos las miRNA como diferencialmente expresadas utilizando un nivel de significancia de . Para las miRNA identificadas como diferencialmente expresadas por DESeq2 y no por miRglmm, realizamos una prueba para el uso diferencial de isomiR utilizando una prueba de razón de verosimilitud.### Métricas de evaluación del rendimiento.

Datos nulos simulados.

Para cada miRNA en los conjuntos de datos simulados del modelo de error, sabemos si cada secuencia es una secuencia de isomiR biológica o una secuencia de isomiR de error. Para evaluar el rendimiento del modelo de error en los datos simulados, informamos la tasa de falsos positivos (FPR) y la tasa de verdaderos positivos (TPR):
Tasa de falsos positivos: , muestra la proporción de secuencias de isomiR de error que se identifican erróneamente como secuencias de isomiR biológicas por el algoritmo. Tasa de verdaderos positivos: , muestra la proporción de secuencias de isomiR biológicas que se identifican correctamente como secuencias de isomiR biológicas por el algoritmo.
donde los verdaderos positivos son las secuencias de isomiR biológicas que se identifican correctamente, son las secuencias de isomiR de error que han sido etiquetadas como secuencias de isomiR biológicas por el algoritmo, son las secuencias de isomiR de error que han sido etiquetadas como secuencias de isomiR de error por el algoritmo, y son las secuencias de isomiR biológicas que han sido etiquetadas como secuencias de isomiR de error por el algoritmo. Dado que aplicamos PARiS a 10 conjuntos de datos simulados para cada una de las configuraciones de simulación, informamos la distribución de las tasas de FP y TP en los conjuntos de datos simulados.### Conjuntos de datos simulados y conjunto de datos de referencia experimental con cambio de pliegue logarítmico real.

Tanto para los 50 conjuntos de datos de monocitos simulados como para el conjunto de datos ERCC, evaluamos el rendimiento de un método en un conjunto de datos dado con el error cuadrático medio (MSE) de las estimaciones del cambio de pliegue logarítmico y con la proporción de cobertura de los intervalos de confianza del 95% estimados. Utilizamos las definiciones de ambos términos que se dan comúnmente en la literatura.### Software, estructuras de datos, resultados y reproducibilidad.

El algoritmo PARiS está escrito en el lenguaje de programación R. Utilizamos el paquete ‘SummarizedExperiment’ en el lenguaje de programación R para organizar los datos de recuento, , la tabla de datos con la información del alineamiento definida por la función , y la información adicional a nivel de muestra en un solo objeto. En el lenguaje del paquete SummarizedExperiment, los datos de recuento son un ensayo. La tabla de datos con los datos adicionales producidos por el alineamiento es la ‘rowData’ y la tabla de datos con la información de la covariable a nivel de muestra es la ‘colData’. Muchos de los paquetes R creados para realizar análisis de expresión diferencial, ya sea de datos de miRNA o de datos de RNA-seq a granel, esperan un objeto SummarizedExperiment con estos elementos como entrada. Por lo tanto, PARiS toma un objeto SummarizedExperiment como entrada y devuelve un objeto SummarizedExperiment limpio como salida.

Modelo de error.

Una secuencia de isomiR de error , también conocida como secuencia de isomiR de error técnico, no aparece en el sistema biológico real del que se tomó la muestra. Más bien, la variación no deseada en el proceso de secuenciación resultó en errores a lo largo de la lectura de la secuencia de isomiR biológica , produciendo la secuencia . Supongamos que tenemos una secuencia de isomiR de error que se mapea a miRNA en la muestra de un experimento de perfilado de miRNA con recuento de lecturas asociado . También se mapea a miRNA en la muestra una secuencia de isomiR biológica verdadera, , con recuento de lecturas asociado . El modelo de error, que incluye la estimación de parámetros y las pruebas de hipótesis, se basa en un alineamiento por pares entre el isomiR biológico y la secuencia . Aquí, utilizamos el alineamiento por pares para referirnos a los resultados del uso del algoritmo de Needleman-Wunsch para obtener un alineamiento global de dos secuencias de nucleótidos (^[30]).

Ciertos aspectos de nuestro modelo de error están inspirados en el modelo DADA2 para la limpieza de variantes de secuenciación de amplicones (^[17]). Primero, asumimos que los errores de lectura de la secuencia ocurren de forma independiente dentro de una lectura de una secuencia de isomiR y entre lecturas. A partir de estas suposiciones, se deduce que la probabilidad de producir una lectura de a partir de debido a errores técnicos se puede expresar como el producto de las probabilidades de las transiciones de nucleótidos individuales a lo largo de la longitud del alineamiento por pares de a . Sea denota la probabilidad de producir una lectura de la secuencia de error técnico a partir de la secuencia . Sea denota la probabilidad del nucleótido en la posición de del alineamiento por pares que se lee como el nucleótido en la posición de del alineamiento por pares. Finalmente, sea la longitud del alineamiento por pares entre y . Entonces, a partir de nuestras suposiciones, tenemos:

La probabilidad de transición de un carácter a otro carácter se define por la identidad de esos caracteres. Nos referimos a estos caracteres individuales del alineamiento por pares entre las secuencias y como caracteres y no solo nucleótidos porque permitimos que estos caracteres sean nucleótidos o un espacio: . Esto nos permite manejar las transiciones de nucleótido a nucleótido y de espacio a nucleótido o de nucleótido a espacio simultáneamente.

Cada uno de estos diferentes tipos de transiciones que permitimos representa un tipo diferente de isomiR. Las transiciones de nucleótido a nucleótido representan errores de sustitución. Un alineamiento por pares entre y que contiene solo transiciones de nucleótido a nucleótido indica que es una secuencia variante de isomiR de . De manera similar, las transiciones de espacio a nucleótido representan errores de adición y las transiciones de nucleótido a espacio representan errores de deleción. Un alineamiento por pares entre y que contiene solo transiciones de espacio a nucleótido o de nucleótido a espacio indica que es una variante de longitud de isomiR de . Las variantes de longitud, particularmente en el extremo 3′, son comunes en la secuenciación de miRNA porque la heterogeneidad de la escisión 5′/3′ y el acoplamiento del extremo 3′ (^[20]); en la secuenciación de amplicones de microbiomas, los cebadores conservados, el recorte de calidad y la fusión de lectura en pares producen longitudes de amplicones casi uniformes (^[17]). Finalmente, una secuencia con los tres tipos de transiciones es una secuencia de isomiR de tipo mixto. Organizamos estas probabilidades de transición en una matriz, . Las probabilidades de transición pueden especificarse previamente o estimarse directamente a partir de los datos. Los detalles sobre la metodología utilizada para la estimación de parámetros se encuentran en los Materiales Suplementarios.

A continuación, denote una variable aleatoria que representa el recuento de lecturas de la secuencia de error en la muestra . Sea representa la abundancia verdadera de la secuencia de isomiR biológica en la muestra . El valor es una variable aleatoria latente. Utilizamos para denotar su recuento de lecturas observado:

En cada muestra , cada molécula de secuencia se selecciona para la secuenciación con probabilidad . Por lo tanto, el número de lecturas observadas que se originan en sigue

Condicionalmente a que se origina en , una lectura se observa como secuencia de error con probabilidad , determinado por el modelo de error Ecuación Ecuación (1). Por la propiedad de adelgazamiento de Poisson, el recuento de lecturas asignadas a sigue:

o, equivalentemente, condicionado a :

Por lo tanto, vemos el proceso de secuenciación para como que surge de dos adelgazamientos de Poisson: 1) selección de lecturas del recuento verdadero basado en la probabilidad de muestreo , y 2) asignación errónea de estas lecturas a la secuencia de error con probabilidad .

Prueba de hipótesis.

Utilizamos el modelo de error descrito anteriormente para probar la hipótesis nula que la secuencia que se mapea a miRNA es una secuencia de isomiR de error de la secuencia de isomiR biológica . Específicamente, comparamos el recuento de lecturas observado de la secuencia en la muestra con lo que esperaríamos bajo el modelo de error de la Ecuación 5. Si es mucho más abundante de lo que esperaríamos bajo el modelo de error, rechazamos la hipótesis nula e inferimos que es otra secuencia de isomiR biológica de miRNA . Formalizamos la idea de ser “mucho más abundante de lo que esperaríamos” tomando prestado el concepto de valor p de abundancia del algoritmo DADA2 (^[17]). El valor p de abundancia se define matemáticamente como:

donde representa la función de densidad de Poisson. El valor p de abundancia es la probabilidad de observar un recuento de lecturas para la secuencia en la muestra tan extremo o más extremo que , asumiendo que es una secuencia de isomiR de error de la secuencia de isomiR biológica . El criterio de decisión para rechazar la hipótesis nula compara el valor p de abundancia, , con un umbral definido por el usuario, funciona de manera similar al nivel en una prueba de hipótesis frecuentista tradicional y se establece comúnmente en 0.05.

En la práctica, calculamos un valor p de abundancia para cada secuencia de isomiR tal que y lo comparamos con . Para controlar el número de falsos descubrimientos que cometemos, corregimos los valores p de abundancia para pruebas múltiples utilizando el procedimiento de Benjamini-Hochberg (^[31]). Utilizaremos para representar el valor p de abundancia ajustado de Benjamini-Hochberg.

PARiS: Asignación y Repartición Probabilística de Secuencias de IsomiR.

El algoritmo PARiS toma, como entrada, una matriz -dimensional de recuentos de isomiR, . El primer paso del algoritmo PARiS es estimar las probabilidades de transición almacenadas en utilizando los métodos de estimación de parámetros descritos en los Materiales Suplementarios. PARiS es un algoritmo iterativo que particiona el conjunto de secuencias que se mapean a miRNA en la muestra a través de un proceso de refinamiento secuencial. Eliminamos la notación para y en lo siguiente. Definimos la partición inicial por:

donde para algún miRNA .

Utilizamos para indexar las iteraciones del algoritmo PARiS. En la iteración , sea una partición de en subconjuntos disjuntos (no vacíos), es decir,

Cada bloque tiene una secuencia central designada , que representa un isomiR biológico.

En la iteración , obtenemos una partición más fina dividiendo el bloque . Los pasos son los siguientes: i) aplicar el modelo de error de la Ecuación 5 a cada , tomando para calcular el valor p de abundancia; ii) después de calcular , ajustamos para la corrección por múltiples pruebas utilizando el procedimiento de Benjamini-Hochberg (^[31]) para generar ; iii) comparamos con , el umbral de significancia definido por el usuario; iv) identificamos una nueva secuencia central, , mediante

Después de identificar , el conjunto se define como

Este proceso iterativo se lleva a cabo de forma secuencial, donde las secuencias “fluyen” de a , y así sucesivamente, hasta que no se puedan formar más subconjuntos o se alcance un número máximo preespecificado de iteraciones.

Después de que se ha asignado cada , hemos definido una partición de las secuencias isomiR que se mapean al miARN en la muestra . Nos referiremos a como la partición a nivel de muestra del miARN en la muestra o , omitiendo la dependencia del miARN cuando quede claro por el contexto.

Formando una partición de consenso y reduciendo el ruido en los recuentos.

Después de obtener las particiones por muestra para el miARN , las consolidamos en una única partición coherente definiendo primero un conjunto de consenso de secuencias centrales. Sea el conjunto de secuencias centrales identificadas para el miARN en cada una de las muestras. Definimos el conjunto de consenso de secuencias centrales para el miARN como .

Intuitivamente, contiene las secuencias centrales que aparecen en una proporción suficiente de muestras. Dos extremos ilustran esta idea:

Tomar la unión es permisivo y conlleva el riesgo de falsos positivos (incluyendo isomiRs de error como centros), mientras que tomar la intersección es estricto y conlleva el riesgo de falsos negativos (excluyendo isomiRs biológicos verdaderos).

Para equilibrar estos extremos, introducimos un umbral de prevalencia : una secuencia central candidata se retiene en si aparece en al menos una fracción de las muestras. Establecer produce la unión, mientras que produce la intersección. El conjunto de consenso resultante define los bloques de la partición de consenso, que luego utilizamos para reducir el ruido en los recuentos en todas las muestras.

Tras la identificación del conjunto de consenso de secuencias centrales , reasignamos cada secuencia no central de la siguiente manera. Para cada secuencia central , calculamos como en la Ecuación 1, y asignamos a el bloque cuya secuencia central alcanza el valor más alto de . El número total de subconjuntos para es dado por la cardinalidad de .

Después de construir la partición de consenso , reducimos el ruido en los recuentos agregando las lecturas dentro de los bloques. Sea el centro del bloque . Para la muestra , el recuento reducido en el ruido para la secuencia isomiR biológica putativa es

La composición inferida para el miARN se da por sus centros de consenso , con abundancias reducidas en el ruido en todas las muestras.

Simulación de datos a partir de la distribución de errores.

Comenzamos evaluando la capacidad de PARiS para identificar y eliminar las secuencias isomiR de error técnico sin introducir nuevos errores. Para ello, necesitamos un conjunto de datos en el que, para cada secuencia isomiR , sepamos si es un isomiR biológico o un isomiR de error. Generamos nuestros propios datos de secuenciación de miARN con recuentos de lecturas asociados y secuencias isomiR de error verdaderas a partir de un conjunto de datos de referencia experimental que secuenciaba los miARN del ratón. El conjunto de datos de referencia experimental que utilizamos consistió en 421.096 secuencias que se mapearon a 759 miARN distintos después del alineamiento a miRGeneDB (^[8]) utilizando sRNAbench (^[21]). Filtramos los miARN de baja expresión con recuentos medianos por millón menores que 5. Después del filtrado, quedaron 435 miARN. Para cada miARN, identificamos , la secuencia más abundante que se mapea al miARN en todas las muestras. Luego, para cada miARN, simulamos la adición, eliminación o sustitución de la cadena que representa para producir secuencias isomiR a partir de la secuencia central. Luego, simulamos los recuentos de lecturas a partir del modelo de error dado en la Ecuación 5 para generar los recuentos de lecturas de las secuencias isomiR simuladas. Los detalles de la simulación de datos se pueden encontrar en los materiales suplementarios.

Variamos el número de secuencias isomiR biológicas verdaderas que se mapean a cada miARN seleccionando el valor del conjunto {1, 4, 9}. Luego, para un número dado de secuencias isomiR biológicas verdaderas, repetimos el proceso de generación de secuencias isomiR simuladas y la simulación de recuentos de lecturas a partir del modelo de error para generar un conjunto de datos simulado 10 veces.

Conjuntos de datos simulados de monocitos.

Para demostrar que PARiS mejora el análisis de expresión diferencial a nivel de miARN, necesitamos un conjunto de datos con varias características clave. En primer lugar, el conjunto de datos debe tener recuentos de lecturas de secuencias isomiR. En segundo lugar, el conjunto de datos debe tener valores verdaderos para la cantidad de expresión diferencial para cada uno de los miARN en el conjunto de datos. La expresión diferencial se mide típicamente mediante el cambio logarítmico (logFC) entre las condiciones experimentales, por lo que debemos conocer el logFC verdadero para cada miARN o ser capaces de calcularlo.

Para generar un conjunto de datos que cumpla con estos requisitos, utilizamos el mismo proceso utilizado por Baran et al. para inyectar artificialmente una señal de expresión diferencial en un conjunto de datos real (^[15]). El conjunto de datos inicial consistió en 39 muestras de monocitos. Después de filtrar los miARN de baja expresión y eliminar las secuencias que no coincidían exactamente con la base de datos de referencia, el conjunto de datos inicial consistió en 122 miARN que se mapearon a 3538 secuencias únicas. A partir de este conjunto de datos inicial, se crearon 50 conjuntos de datos simulados utilizando el siguiente proceso: las 39 muestras se dividieron en dos grupos (A y B) de forma aleatoria para imitar un diseño experimental simple (^[15]). De los 122 miARN, se seleccionaron 20 para que se sobreexpresaran en el grupo A de forma aleatoria y se seleccionaron 20 para que se sobreexpresaran en el grupo B de forma aleatoria. Los recuentos de lecturas a nivel de isomiR de los miARN seleccionados para que se sobreexpresaran en cada grupo se multiplicaron por valores muestreados a partir de una distribución normal truncada con una media de 2 y una desviación estándar de 1. Se pueden encontrar más detalles sobre este proceso en Baran et al. (^[15]). El proceso de simulación se repitió para generar 50 conjuntos de datos sintéticos, con valores de logFC conocidos para cada miARN y cada secuencia isomiR en el conjunto de datos. No observamos si cada secuencia isomiR es una secuencia isomiR biológica o una secuencia isomiR de error.

Aplicamos diferentes flujos de trabajo analíticos, que consisten en un método de reducción de ruido seleccionado del conjunto {PARiS, miREC, Agregación o Ninguno} y un método apropiado para estimar la expresión diferencial a nivel de miARN dada la resolución de los datos reducidos en el ruido. Para el flujo de trabajo analítico que utiliza PARiS, utilizamos como umbral para identificar los valores p de abundancia significativos. Variamos el valor para crear conjuntos de consenso de secuencias centrales para cada miARN seleccionando de {0,025, 0,05, 0,10, 0,20, 0,40, 0,60, 0,80, 1,00}. PARiS y miREC producen datos reducidos en el ruido a nivel de isomiR, y los datos brutos (la opción de reducción de ruido Ninguno) también están a nivel de isomiR. Estimamos la expresión diferencial a nivel de miARN utilizando los datos reducidos en el ruido de cada uno de estos métodos utilizando miRglmm, un marco de modelado lineal mixto generalizado (^[15]). El método de agregación produce datos reducidos en el ruido a nivel de miARN. Estimamos la expresión diferencial a nivel de miARN a partir de los datos agregados utilizando DESeq2 (^[14]) y edgeR (^[13]). En total, aplicamos 11 flujos de trabajo analíticos diferentes a cada uno de los 50 conjuntos de datos simulados.

Conjunto de datos de referencia ERCC.

Evaluamos el rendimiento de PARiS en un flujo de trabajo analítico para el análisis de expresión diferencial de miARN en un conjunto de datos de referencia sintético y experimental. Debido a la naturaleza sintética del conjunto de datos, no hay variación biológica presente. Utilizamos el conjunto de datos de referencia del Consorcio de Comunicación de ARN Extracelular (ERCC) (^[22]). El conjunto de datos ERCC es sRNA-seq en grupos ratiométricos (A y B) de ARN pequeños sintetizados en proporciones de 10:1 a 1:10. El resultado es un conjunto de datos de expresión de secuencias isomiR con 15 niveles de expresión diferencial verdadera, que van desde log(0,1) = -2,3 hasta log(10) = 2,3. El conjunto de datos contiene 286 miARN humanos que se mapean a 8001 secuencias después de eliminar las secuencias que tienen ≥ 40 nucleótidos o < 16 nucleótidos.

Reducimos el ruido en el conjunto de datos ERCC utilizando un conjunto similar de métodos que se aplicaron a los conjuntos de datos simulados de monocitos. Los métodos de reducción de ruido que aplicamos al conjunto de datos ERCC son los siguientes: miREC, agregación y . Para cada aplicación de PARiS, establecimos . En total, esto da un total de 22 flujos de trabajo analíticos diferentes. Después de la reducción de ruido en cada flujo de trabajo, estimamos la expresión diferencial a nivel de miARN utilizando el método apropiado dado la resolución de los datos reducidos en el ruido. Para los datos corregidos por miREC y los datos reducidos en el ruido por PARiS, estimamos la expresión diferencial a nivel de miARN utilizando miRglmm (^[15]). Para los datos agregados, estimamos la expresión diferencial a nivel de miARN utilizando DESeq2 (^[14]) y edgeR (^[13]).

Comparación de líneas celulares de adenocarcinoma de colon.

Aplicamos PARiS como el paso de reducción de ruido en un flujo de trabajo de análisis de expresión diferencial de miARN aplicado a un conjunto de datos experimental verdadero. Los datos de sRNA-seq (N=9) provienen de 3 líneas celulares de adenocarcinoma de colon (DLD-1, DKO-1 y DKS-8), que varían según el estado de KRAS (^[23]). Los datos de secuenciación brutos consisten en 5223 secuencias únicas que se mapean a 143 miARN después de filtrar los miARN de baja expresión con lecturas por millón < 5 y mantener solo las secuencias que se alinean exactamente con la base de datos de referencia. Aplicamos PARiS con y y luego estimamos la expresión diferencial a nivel de miARN utilizando miRglmm (^[15]). También agregamos los recuentos a nivel de isomiR al nivel de miARN y estimamos la expresión diferencial utilizando DESeq2 para fines de comparación (^[14]). Identificamos los miARN como diferencialmente expresados utilizando un nivel de significancia de . Para los miARN identificados como diferencialmente expresados por DESeq2 y no por miRglmm, probamos la utilización diferencial de isomiR utilizando una prueba de razón de verosimilitud.

Métricas de evaluación del rendimiento.

Datos nulos simulados.

Para cada miARN en los conjuntos de datos simulados a partir del modelo de error, sabemos si cada secuencia es una secuencia isomiR biológica o una secuencia isomiR de error. Para evaluar el rendimiento del modelo de error en los datos simulados, informamos la tasa de falsos positivos (FPR) y la tasa de verdaderos positivos (TPR):
Tasa de falsos positivos: , muestra la proporción de secuencias isomiR de error que se identifican erróneamente como secuencias isomiR biológicas por el algoritmo. Tasa de verdaderos positivos: , muestra la proporción de secuencias isomiR biológicas que se identifican correctamente como secuencias isomiR biológicas por el algoritmo.
donde los verdaderos positivos son las secuencias isomiR biológicas correctamente identificadas, son las secuencias isomiR de error que han sido etiquetadas como secuencias isomiR biológicas por el algoritmo, son las secuencias isomiR de error que han sido etiquetadas como secuencias isomiR de error por el algoritmo, y son las secuencias isomiR biológicas que han sido etiquetadas como secuencias isomiR de error por el algoritmo. Dado que aplicamos PARiS a 10 conjuntos de datos simulados para cada uno de los ajustes de simulación, informamos la distribución de las tasas de FP y TP en los conjuntos de datos simulados.

Conjuntos de datos simulados y conjunto de datos de referencia experimental con cambio logarítmico verdadero.

Para ambos, los 50 conjuntos de datos simulados de monocitos y el conjunto de datos ERCC, evaluamos el rendimiento de un método en un conjunto de datos dado con el error cuadrático medio (MSE) de las estimaciones del cambio logarítmico y con la proporción de cobertura de los intervalos de confianza del 95% estimados. Utilizamos las definiciones de ambos términos que se dan comúnmente en la literatura.

Datos nulos simulados.

Para cada miRNA en los conjuntos de datos simulados a partir del modelo de errores, sabemos si cada secuencia es una secuencia biológica de isomiR o una secuencia de isomiR de error. Para evaluar el rendimiento del modelo de errores en los datos simulados, informamos sobre la tasa de falsos positivos (FPR) y la tasa de verdaderos positivos (TPR):
Tasa de falsos positivos: , muestra la proporción de secuencias de isomiR de error que se identifican erróneamente como secuencias biológicas de isomiR por el algoritmo. Tasa de verdaderos positivos: , muestra la proporción de secuencias biológicas de isomiR que se identifican correctamente como secuencias biológicas de isomiR por el algoritmo.
donde los verdaderos positivos son las secuencias biológicas de isomiR que se identifican correctamente, son las secuencias de isomiR de error que han sido etiquetadas como secuencias biológicas de isomiR por el algoritmo, son las secuencias de isomiR de error que han sido etiquetadas como secuencias de isomiR de error por el algoritmo, y son las secuencias biológicas de isomiR que han sido etiquetadas como secuencias de isomiR de error por el algoritmo. Dado que aplicamos PARiS a 10 conjuntos de datos simulados para cada una de las configuraciones de simulación, informamos sobre la distribución de las tasas de FP y TP en los conjuntos de datos simulados.

Conjuntos de datos simulados y conjunto de datos de referencia experimental con cambio logarítmico real.

Tanto para los 50 conjuntos de datos simulados de monocitos como para el conjunto de datos ERCC, evaluamos el rendimiento de un método en un conjunto de datos dado utilizando el error cuadrático medio (MSE) de las estimaciones del cambio logarítmico y con la proporción de cobertura de los intervalos de confianza del 95% estimados. Utilizamos las definiciones de ambos términos que se dan comúnmente en la literatura.

Software, estructuras de datos, resultados y reproducibilidad.

El algoritmo PARiS está escrito en el lenguaje de programación R. Utilizamos el paquete ‘SummarizedExperiment’ en el lenguaje de programación R para organizar los datos de recuento, , la tabla de datos con la información del alineamiento definida por la función , y la información adicional a nivel de muestra en un único objeto. En el lenguaje del paquete SummarizedExperiment, los datos de recuento son un ensayo. La tabla de datos con los datos adicionales producidos por el alineamiento es la ‘rowData’ y la tabla de datos con la información de la covariable a nivel de muestra es la ‘colData’. Muchos de los paquetes de R creados para realizar análisis de expresión diferencial, ya sea de datos de miRNA o de datos de RNA-seq a granel, esperan un objeto SummarizedExperiment con estos elementos como entrada. Por lo tanto, PARiS toma un objeto SummarizedExperiment como entrada y devuelve un objeto SummarizedExperiment "limpio" como salida.

Resultados

El modelo de errores en PARiS identifica correctamente las secuencias de isomiR de error de variante de longitud técnica y de secuencia en los datos simulados.

Comenzamos con los resultados del uso del modelo de errores propuesto en el que se basa PARiS para identificar las secuencias de isomiR de error de variante de longitud técnica. Permitimos que el número real de secuencias biológicas de isomiR por miRNA tenga un valor diferente del conjunto {1, 4, 9}. Para cada miRNA en el conjunto de datos, calculamos la tasa de falsos positivos según la definición proporcionada en la sección anterior. Cuando hay 1 verdadera secuencia biológica de isomiR por miRNA, la tasa de falsos positivos promedio es de aproximadamente 0,015 (Figura 1, A, superior). Cuando hay 4 verdaderas secuencias biológicas de isomiR por miRNA, la tasa de falsos positivos promedio fue de 0,09 (Figura 1, B, centro). Finalmente, para el escenario en el que hay 9 verdaderas secuencias biológicas de isomiR por miRNA, la tasa de falsos positivos promedio fue de 0,10 (Figura 1, A, inferior).

Las secuencias de isomiR de variante de longitud no son el único tipo de secuencia de isomiR que puede estar presente en los datos. Las secuencias de isomiR de variante de secuencia también pueden estar presentes, y es importante que el modelo de errores en PARiS pueda manejarlas también. Simulamos secuencias de isomiR de variante de secuencia a partir de los datos, permitiendo nuevamente que el número real de secuencias biológicas de isomiR tenga diferentes valores del conjunto {1, 4, 9}. Para cada miRNA, informamos sobre la tasa de falsos positivos y la tasa de verdaderos positivos de la identificación de secuencias de isomiR de error de variante de secuencia técnica a partir de los datos. Cuando hay solo 1 verdadera secuencia biológica de isomiR por miRNA, la tasa de falsos positivos promedio fue de aproximadamente 0,04 (Figura 1, B, superior). Cuando hay 4 verdaderas secuencias biológicas de isomiR por miRNA, la tasa de falsos positivos promedio fue de aproximadamente 0,04 (Figura 1, B, centro). Finalmente, cuando hay 9 verdaderas secuencias biológicas de isomiR por miRNA, la tasa de falsos positivos promedio fue de 0,03 (Figura 1, B, inferior).

PARiS disminuye el MSE y aumenta la proporción de cobertura de las estimaciones del cambio logarítmico a nivel de miRNA utilizando datos de expresión a nivel de isomiR en datos simulados.

Comenzamos examinando el MSE de los diferentes flujos de trabajo analíticos utilizados para "limpiar" los datos y estimar la expresión diferencial a nivel de miRNA en el conjunto de 50 conjuntos de datos simulados de monocitos. Un MSE más bajo indica que el flujo de trabajo analítico genera estimaciones de logFC que son, en promedio, más cercanas a los valores de logFC reales. En los 50 conjuntos de datos simulados, los modelos miRglmm ajustados a los datos de recuento a nivel de isomiR "limpios" con PARiS, con y , tuvieron el MSE medio más bajo (MSE = 0,0162) de todos los flujos de trabajo analíticos evaluados (Tabla 1). Sorprendentemente, el MSE promedio en los 50 conjuntos de datos simulados es en realidad mayor para los modelos miRglmm ajustados a los datos corregidos con miREC (media = 0,034) que para los modelos miRglmm ajustados a los datos brutos (media = 0,033).

El flujo de trabajo analítico que "limpió" los datos con PARiS, con y , y estimó la expresión diferencial a nivel de miRNA utilizando miRglmm generó estimaciones de logFC con el MSE promedio más alto. Los flujos de trabajo analíticos que utilizaron PARiS para "limpiar" los datos con y cualquier otro valor de , y luego estimaron la expresión diferencial a nivel de miRNA utilizando miRglmm, generaron estimaciones de logFC con un MSE promedio menor que los datos "limpios" con cualquier otro método.

Además del MSE, también evaluamos los flujos de trabajo analíticos en términos de la proporción de cobertura de los intervalos de confianza del 95% estimados. La estimación de los intervalos de confianza requiere estimaciones puntuales del parámetro de interés y una estimación del error estándar. Debido a que edgeR no proporciona al usuario estimaciones del error estándar, solo informamos sobre el rendimiento en términos de la proporción de cobertura para los flujos de trabajo analíticos que estimaron la expresión diferencial utilizando miRglmm o DESeq2.

De todos los flujos de trabajo analíticos utilizados para estimar la expresión diferencial, los flujos de trabajo analíticos que "limpiaron" los datos con PARiS, con y , y luego estimaron la expresión diferencial utilizando miRglmm, estimaron intervalos de confianza del 95% con la mayor proporción de cobertura promedio (proporción de cobertura promedio = 0,90) (Tabla 2). La siguiente proporción de cobertura más alta provino del flujo de trabajo analítico que "limpió" los datos utilizando PARiS, con y , y luego estimó la expresión diferencial con miRglmm. Independientemente del valor de utilizado, cualquier flujo de trabajo analítico que "limpie" los datos con PARiS y estime la expresión diferencial con miRglmm, estimó intervalos de confianza del 95% con una mayor proporción de cobertura promedio que los métodos de comparación. Cualquier flujo de trabajo analítico que utilice PARiS, independientemente del valor de utilizado, también estimó intervalos de confianza del 95% con una mayor proporción de cobertura mínima y una mayor proporción de cobertura máxima en las 50 simulaciones.

Además de evaluar el rendimiento de los diferentes flujos de trabajo analíticos en términos de MSE y proporción de cobertura en todos los miRNAs, también evaluamos el rendimiento de cada método dentro de cada uno de los grupos de verdad en las simulaciones de monocitos. Cada conjunto de datos en el conjunto de 50 conjuntos de datos simulados de monocitos contiene 3 grupos de verdad. El primer grupo es el grupo de miRNAs con un cambio positivo inducido del grupo al grupo (logFC igual a log(2)). El segundo grupo es el grupo de miRNAs con un cambio negativo inducido del grupo al grupo (logFC igual a log(0,5)). El tercer y último grupo es el grupo de miRNAs sin cambio inducido del grupo al grupo (logFC igual a log(1)). La distribución de los MSE en las 50 simulaciones para cada flujo de trabajo analítico, dentro del grupo de verdad, se muestra en la Figura 2 (A). La distribución de la proporción de cobertura en las 50 simulaciones para cada flujo de trabajo analítico, dentro del grupo de verdad, se muestra en la Figura 2 (B). Nuevamente, no se informan los resultados de la proporción de cobertura para el flujo de trabajo analítico que agrega los recuentos de secuencia a nivel de miRNA y utiliza edgeR para estimar la expresión diferencial. Esto se debe a que edgeR no informa los errores estándar de las estimaciones de logFC.

En los grupos con un efecto inducido, el uso de un flujo de trabajo analítico que utiliza miRglmm para estimar la expresión diferencial después de "limpiar" con PARiS, con , independientemente del valor de seleccionado, produce estimaciones de logFC con un MSE menor o igual al MSE de las estimaciones de logFC producidas por "limpiar" con métodos de comparación (Figura 2, A). En el grupo sin logFC, los flujos de trabajo analíticos que utilizan técnicas de agregación y estiman la expresión diferencial utilizando DESeq2 o edgeR minimizan el MSE (Figura 2, A). Sin embargo, todos los flujos de trabajo analíticos logran un MSE comparablemente bajo en el grupo de verdad sin efecto inducido, excepto el flujo de trabajo que utiliza PARiS con .

De manera similar, en los grupos con un efecto inducido, el uso de un flujo de trabajo analítico que utiliza miRglmm para estimar la expresión diferencial después de "limpiar" con PARiS, con , independientemente del valor de seleccionado, estima intervalos de confianza del 95% con proporciones de cobertura mayores o iguales a las proporciones de cobertura de los intervalos de confianza del 95% estimados de los métodos de comparación (Figura 2, (B)). En el grupo de verdad sin efecto inducido, el flujo de trabajo analítico que agrega los recuentos de secuencia a nivel de miRNA y estima la expresión diferencial utilizando DESeq2 estima intervalos de confianza del 95% con la mayor proporción de cobertura. Si bien el flujo de trabajo analítico de DESeq2 logra la mayor proporción de cobertura, todos los métodos de comparación logran una proporción de cobertura comparablemente alta de los intervalos de confianza del 95% estimados en el grupo de verdad sin cambio inducido.

PARiS identifica y elimina las secuencias de isomiR que miREC no detecta, mejorando así las estimaciones del cambio logarítmico a nivel de miRNA de los datos de expresión a nivel de isomiR cuando no hay variabilidad biológica de la secuencia.

Primero, observamos cómo cambia la distribución del número de isomiRs por miRNA en los datos "limpios" con PARiS frente a los datos corregidos con miREC a medida que cambiamos el valor de utilizado para generar conjuntos de consenso de secuencias centrales. En los datos tratados con miREC, el número promedio de isomiRs por miRNA es 20. Dada la naturaleza sintética de los datos ERCC, el número real de isomiRs por miRNA es 1. "Limpiar" los datos con y cualquier valor de mueve la distribución de isomiRs por miRNA lejos de la distribución de los datos brutos y hacia el valor real (Figura 3, panel A).

A continuación, evaluamos los diferentes flujos de trabajo analíticos en términos del MSE de las estimaciones de logFC que produjo cada flujo de trabajo. Recuerde que cuanto menor sea el MSE, mejor será el rendimiento del método, ya que las estimaciones se acercarán más al valor real. El uso de un flujo de trabajo analítico que utiliza PARiS para "limpiar" los datos, con , independientemente del valor de , produce estimaciones de logFC a nivel de miRNA con un MSE menor que las producidas por los otros flujos de trabajo analíticos (Figura 3, panel B). Además, a medida que aumenta de 0,025 a 1, el MSE de las estimaciones de logFC disminuye. El flujo de trabajo analítico que "limpió" los datos con miREC y estimó la expresión diferencial con miRglmm produjo estimaciones de logFC a nivel de miRNA con el MSE más grande.

Finalmente, utilizamos las proporciones de cobertura de los intervalos de confianza del 95% estimados a partir de los diferentes métodos de análisis para evaluarlos. La definición de la proporción de cobertura es la misma que la utilizada para evaluar los diferentes métodos de análisis en el conjunto de datos simulados de monocitos. Es decir, la proporción de cobertura es la proporción de miRNAs cuyo valor real de logFC se encuentra dentro del intervalo de confianza del 95% estimado. Cuanto mayor sea la proporción de cobertura, mejor será el rendimiento del método. Todos los métodos de análisis estimaron intervalos de confianza del 95% con proporciones de cobertura iguales o superiores al nivel nominal de 0,95 (Figura 3, panel c). Aunque existe variabilidad en la proporción de cobertura de los intervalos de confianza del 95% con respecto al valor utilizado en PARiS para eliminar el ruido de los datos, esta es mínima. El uso de un método de análisis que elimina el ruido de los datos con PARiS, con cualquier valor, estima intervalos de confianza del 95% con proporciones de cobertura superiores al nivel nominal de 0,95. Esto indica que la eliminación de ruido con PARiS mejora las estimaciones puntuales sin sacrificar las estimaciones del error estándar de las estimaciones puntuales.### El uso de PARiS como paso en un método de análisis aplicado a un conjunto de datos experimental real puede evitar que los investigadores confundan el uso diferencial de isomiRs con la expresión diferencial.

Realizamos múltiples comparaciones por pares de tres líneas celulares de adenocarcinoma colorrectal. La primera comparación que realizamos comparó los perfiles de expresión de isomiRs de las líneas celulares DKO-1 y DKS-8, y la segunda comparación que realizamos comparó los perfiles de expresión de isomiRs de las líneas celulares DKO-1 y DLD-1. Utilizamos dos métodos de análisis diferentes para el análisis de expresión diferencial, uno que eliminó el ruido de los datos con PARiS, utilizando y , y estimó la expresión diferencial utilizando miRglmm. El otro método de análisis agrupó los datos a nivel de isomiR al nivel de miRNA y estimó la expresión diferencial utilizando DESeq2. El método de análisis que utiliza PARiS y miRglmm encontró 53 miRNAs que se expresan de forma diferencial entre las líneas celulares DKS-8 y DKO-1 (Figura 4, panel A). El método de análisis que utiliza la agregación y DESeq2 encontró 34 miRNAs que se expresan de forma diferencial entre las líneas celulares DKS-8 y DKO-1. Los dos métodos compartieron 19 miRNAs en común.

A continuación, comparamos los patrones de expresión de miRNAs entre las líneas celulares DLD-1 y DKO-1. PARiS / miRglmm identificó 45 miRNAs que se expresan de forma diferencial entre las líneas celulares DLD-1 y DKO-1 (Figura 4, panel B). El método de agregación / DESeq2 identificó 3 miRNAs que se expresan de forma diferencial. Dos miRNAs fueron compartidos por ambos enfoques.

Además de estimar las estimaciones de logFC a nivel de miRNA, también podemos utilizar el marco de modelado miRglmm para realizar una prueba de razón de verosimilitud para el uso diferencial de isomiRs (^[15]). Se hipotetiza que, en el uso diferencial de isomiRs, miRglmm devolverá estimaciones de logFC que estén más cerca de 0 que los métodos de agregación. En este escenario, miRglmm a menudo solo encuentra un uso diferencial de isomiRs y no identifica la expresión diferencial a nivel de miRNA. En nuestra comparación de las líneas celulares DKO-1 y DKS-8, examinamos esta hipótesis. Hay 14 miRNAs que se expresan de forma diferencial entre las líneas celulares DKO-1 y DKS-8 por el método de análisis que utiliza DESeq2 pero no PARiS / miRglmm (Figura 4, panel A). De esos 14 miRNAs, la prueba de razón de verosimilitud identificó 6 de ellos como que tienen un uso diferencial de isomiRs (Tabla 3). Finalmente, de los 6 miRNAs que se infiere que tienen un uso diferencial de isomiRs, 5 tenían estimaciones de logFC más pequeñas del marco de modelado miRglmm que las estimaciones de logFC de DESeq2.

El modelo de error en PARiS identifica y elimina con éxito las secuencias de isomiRs con variantes de longitud y secuencia que miREC no detecta en los datos simulados.

Comenzamos con los resultados del uso del modelo de error propuesto en el que se basa PARiS para identificar las secuencias de isomiRs con variantes de longitud y secuencia que representan errores técnicos. Permitimos que el número real de secuencias de isomiRs biológicos por miRNA tenga un valor diferente del conjunto {1, 4, 9}. Para cada miRNA en el conjunto de datos, calculamos la tasa de falsos positivos según la definición proporcionada en la sección anterior. Cuando había 1 isomiR biológico real por miRNA, la tasa promedio de falsos positivos es de aproximadamente 0,015 (Figura 1, A, superior). Cuando hay 4 isomiRs biológicos reales por miRNA, la tasa promedio de falsos positivos fue de 0,09 (Figura 1, B, centro). Finalmente, para el escenario en el que hay 9 secuencias de isomiRs biológicos reales por miRNA, la tasa promedio de falsos positivos fue de 0,10 (Figura 1, A, inferior).

Los isomiRs con variantes de longitud no son el único tipo de secuencia de isomiR que puede estar presente en los datos. Los isomiRs con variantes de secuencia también pueden estar presentes, y es importante que el modelo de error en PARiS también pueda manejarlos. Simulamos isomiRs con variantes de secuencia a partir de los datos, permitiendo nuevamente que el número real de secuencias de isomiRs biológicos tenga diferentes valores del conjunto {1, 4, 9}. Para cada miRNA, informamos la tasa de falsos positivos y la tasa de verdaderos positivos de la identificación de isomiRs con variantes de secuencia que representan errores técnicos a partir de los datos. Cuando solo hay 1 isomiR biológico real por miRNA, la tasa promedio de falsos positivos fue de aproximadamente 0,04 (Figura 1, B, superior). Cuando hay 4 isomiRs biológicos reales por miRNA, la tasa promedio de falsos positivos fue de aproximadamente 0,04 (Figura 1, B, centro). Finalmente, cuando hay 9 isomiRs biológicos reales por miRNA, la tasa promedio de falsos positivos fue de 0,03 (Figura 1, B, inferior).

PARiS disminuye el MSE y aumenta la proporción de cobertura de las estimaciones de log fold change a nivel de miRNA utilizando datos de expresión a nivel de isomiR en datos simulados.

Comenzamos examinando el MSE de los diferentes métodos de análisis utilizados para eliminar el ruido de los datos y estimar la expresión diferencial a nivel de miRNA en el conjunto de 50 conjuntos de datos simulados de monocitos. Un MSE más bajo indica que el método de análisis genera estimaciones de logFC que, en promedio, están más cerca de los valores de logFC reales. En los 50 conjuntos de datos simulados, los modelos miRglmm ajustados a los datos de recuento a nivel de isomiR eliminados con PARiS, con y , tuvieron el MSE promedio más bajo (MSE = 0,0162) de todos los métodos de análisis evaluados (Tabla 1). Sorprendentemente, el MSE promedio en los 50 conjuntos de datos simulados es en realidad mayor para los modelos miRglmm ajustados a los datos corregidos con miREC (media = 0,034) que para los modelos miRglmm ajustados a los datos sin procesar (media = 0,033).

El método de análisis que eliminó el ruido de los datos con PARiS, con y , y estimó la expresión diferencial a nivel de miRNA utilizando miRglmm generó estimaciones de logFC con el MSE promedio más alto. Los métodos de análisis que utilizaron PARiS para eliminar el ruido de los datos con y cualquier otro valor, y luego estimaron la expresión diferencial a nivel de miRNA utilizando miRglmm, generaron estimaciones de logFC con un MSE promedio menor que los datos eliminados con cualquier otro método.

Además del MSE, también evaluamos los métodos de análisis en términos de la proporción de cobertura de los intervalos de confianza del 95% estimados. La estimación de los intervalos de confianza requiere estimaciones puntuales del parámetro de interés y una estimación del error estándar. Debido a que edgeR no proporciona al usuario estimaciones del error estándar, solo informamos el rendimiento en términos de la proporción de cobertura para los métodos de análisis que estimaron la expresión diferencial utilizando miRglmm o DESeq2.

De todos los métodos de análisis utilizados para estimar la expresión diferencial, los métodos de análisis que eliminaron el ruido de los datos con PARiS, con y , y luego estimaron la expresión diferencial utilizando miRglmm, estimaron intervalos de confianza del 95% con la mayor proporción promedio de cobertura (proporción promedio de cobertura = 0,90) (Tabla 2). La siguiente proporción de cobertura más alta provino del método de análisis que eliminó el ruido de los datos utilizando PARiS, con y , y luego estimó la expresión diferencial con miRglmm. Independientemente del valor utilizado, cualquier método de análisis que eliminara el ruido de los datos con PARiS y estimara la expresión diferencial con miRglmm estimó intervalos de confianza del 95% con una mayor proporción promedio de cobertura que los métodos de comparación. Cualquier método de análisis que utilice PARiS, independientemente del valor utilizado, también estimó intervalos de confianza del 95% con una mayor proporción mínima de cobertura y una mayor proporción máxima de cobertura en las 50 simulaciones.

Además de evaluar el rendimiento de los diferentes métodos de análisis en términos de MSE y proporción de cobertura en todos los miRNAs, también evaluamos el rendimiento de cada método dentro de cada uno de los grupos de verdad en las simulaciones de monocitos. Cada conjunto de datos en el conjunto de 50 conjuntos de datos simulados de monocitos contiene 3 grupos de verdad. El primer grupo es el grupo de miRNAs con un cambio positivo inducido del grupo al grupo (logFC igual a log(2)). El segundo grupo es el grupo de miRNAs con un cambio negativo inducido del grupo al grupo (logFC igual a log(0,5)). El tercer y último grupo es el grupo de miRNAs sin un cambio inducido del grupo al grupo (logFC igual a log(1)). La distribución de los MSE en las 50 simulaciones para cada método de análisis, dentro del grupo de verdad, se muestra en la Figura 2 (A). La distribución de la proporción de cobertura en las 50 simulaciones para cada método de análisis, dentro del grupo de verdad, se muestra en la Figura 2 (B). Nuevamente, no se informan los resultados de la proporción de cobertura para el método de análisis que agrega los recuentos de lectura a nivel de secuencia al nivel de miRNA y utiliza edgeR para estimar la expresión diferencial. Esto se debe a que edgeR no informa los errores estándar de las estimaciones de logFC.

En los grupos con un efecto inducido, el uso de un método de análisis que utiliza miRglmm para estimar la expresión diferencial después de eliminar el ruido con PARiS, con , independientemente del valor seleccionado, produce estimaciones de logFC con un MSE menor o igual al MSE de las estimaciones de logFC producidas por la eliminación de ruido con otros métodos (Figura 2, A). En el grupo sin logFC, los métodos de análisis que utilizan técnicas de agregación y estiman la expresión diferencial utilizando DESeq2 o edgeR minimizan el MSE (Figura 2, A). Sin embargo, todos los métodos de análisis logran un MSE comparablemente bajo en el grupo de verdad sin un efecto inducido, excepto el método que utiliza PARiS con y .

De manera similar, en los grupos con un efecto inducido, el uso de un método de análisis que utiliza miRglmm para estimar la expresión diferencial después de eliminar el ruido con PARiS, con , independientemente del valor seleccionado, estima intervalos de confianza del 95% con proporciones de cobertura mayores o iguales a las proporciones de cobertura de los intervalos de confianza del 95% estimados de los métodos de comparación (Figura 2, (B)). En el grupo de verdad sin un efecto inducido, el método de análisis que agrega los recuentos de lectura a nivel de secuencia al nivel de miRNA y estima la expresión diferencial utilizando DESeq2 estima intervalos de confianza del 95% con la mayor proporción de cobertura. Si bien el método de análisis DESeq2 logra la mayor proporción de cobertura, todos los métodos de comparación logran una proporción de cobertura comparablemente alta de los intervalos de confianza del 95% estimados en el grupo de verdad sin un cambio inducido.

PARiS identifica y elimina las secuencias de isomiRs que miREC no detecta, mejorando así las estimaciones de log fold change a nivel de miRNA de los datos de expresión a nivel de isomiR cuando no hay variabilidad biológica de la secuencia.

Primero, observamos cómo cambia la distribución del número de isomiRs por miRNA en los datos eliminados con PARiS en comparación con los datos corregidos con miREC a medida que cambiamos el valor utilizado para generar conjuntos de consenso de secuencias centrales. En los datos tratados con miREC, el número promedio de isomiRs por miRNA es 20. Dada la naturaleza sintética de los datos ERCC, el número real de isomiRs por miRNA es 1. La eliminación del ruido de los datos con y cualquier valor mueve la distribución de isomiRs por miRNA lejos de la distribución de los datos sin procesar y hacia el valor real (Figura 3, panel A).

A continuación, evaluamos los diferentes flujos de trabajo analíticos en términos del error cuadrático medio (MSE) de las estimaciones de logFC que produce cada flujo de trabajo. Recordemos que cuanto menor sea el MSE, mejor será el rendimiento del método, ya que las estimaciones se acercarán más a los valores verdaderos. El uso de un flujo de trabajo analítico que utiliza PARiS para eliminar el ruido de los datos, con , independientemente del valor de , produce estimaciones de logFC a nivel de miRNA con un MSE menor que las producidas por los demás flujos de trabajo analíticos (Figura 3, panel B). Además, a medida que aumenta de 0,025 a 1, el MSE de las estimaciones de logFC disminuye. El flujo de trabajo analítico que eliminó el ruido de los datos con miREC y estimó la expresión diferencial con miRglmm produjo estimaciones de logFC a nivel de miRNA con el MSE más alto.

Finalmente, utilizamos las proporciones de cobertura de los intervalos de confianza del 95% estimados de los diferentes flujos de trabajo analíticos para evaluarlos. La definición de la proporción de cobertura es la misma que la definición utilizada para evaluar los diferentes flujos de trabajo analíticos en el conjunto de datos simulados de monocitos. Es decir, la proporción de cobertura es la proporción de miRNAs cuyo valor de logFC verdadero se encuentra dentro del intervalo de confianza del 95% estimado. Cuanto mayor sea la proporción de cobertura, mejor será el rendimiento del método. Todos los flujos de trabajo analíticos estimaron intervalos de confianza del 95% con proporciones de cobertura iguales o superiores al nivel nominal de 0,95 (Figura 3, panel c). Aunque existe variabilidad en la proporción de cobertura de los intervalos de confianza del 95% con respecto al valor de utilizado en PARiS para eliminar el ruido de los datos, esta es mínima. El uso de un flujo de trabajo analítico que elimina el ruido de los datos con PARiS, con y cualquier valor de , estima intervalos de confianza del 95% con proporciones de cobertura superiores al nivel nominal de 0,95. Esto indica que la eliminación de ruido con PARiS mejora las estimaciones puntuales sin sacrificar las estimaciones del error estándar de las estimaciones puntuales.

El uso de PARiS como paso en un flujo de trabajo analítico aplicado a un conjunto de datos experimental real puede evitar que los investigadores confundan el uso diferencial de isomiRs con la expresión diferencial.

Realizamos múltiples comparaciones por pares de tres líneas celulares de adenocarcinoma colorrectal. La primera comparación que hicimos comparó los perfiles de expresión de isomiRs de las líneas celulares DKO-1 y DKS-8, y la segunda comparación que hicimos comparó los perfiles de expresión de isomiRs de las líneas celulares DKO-1 y DLD-1. Utilizamos dos flujos de trabajo analíticos diferentes para el análisis de expresión diferencial, uno que eliminó el ruido de los datos con PARiS, utilizando y , y estimó la expresión diferencial utilizando miRglmm. El otro flujo de trabajo analítico agrupó los datos a nivel de isomiR al nivel de miRNA y estimó la expresión diferencial utilizando DESeq2. El flujo de trabajo analítico que utiliza PARiS y miRglmm encontró 53 miRNAs que se expresan de forma diferencial entre las líneas celulares DKS-8 y DKO-1 (Figura 4, panel A). El flujo de trabajo analítico que utiliza la agregación y DESeq2 encontró 34 miRNAs que se expresan de forma diferencial entre las líneas celulares DKS-8 y DKO-1. Los dos flujos de trabajo compartieron 19 miRNAs en común.

A continuación, comparamos los patrones de expresión de miRNAs entre las líneas celulares DLD-1 y DKO-1. PARiS / miRglmm identificó 45 miRNAs que se expresan de forma diferencial entre las líneas celulares DLD-1 y DKO-1 (Figura 4, panel B). El método de agregación / DESeq2 identificó 3 miRNAs que se expresan de forma diferencial. Dos miRNAs fueron compartidos por ambos enfoques.

Además de estimar las estimaciones de logFC a nivel de miRNA, también podemos utilizar el marco de modelado miRglmm para realizar una prueba de razón de verosimilitud para el uso diferencial de isomiRs (^[15]). Se hipotetiza que, en el uso diferencial de isomiRs, miRglmm devolverá estimaciones de logFC que estén más cerca de 0 que los métodos de agregación. En este escenario, miRglmm a menudo encuentra solo un uso diferencial de isomiRs y no identifica la expresión diferencial a nivel de miRNA. En nuestra comparación de las líneas celulares DKO-1 y DKS-8, examinamos esta hipótesis. Hay 14 miRNAs que se expresan de forma diferencial entre las líneas celulares DKO-1 y DKS-8 por el flujo de trabajo analítico que utiliza DESeq2, pero no PARiS / miRglmm (Figura 4, panel A). De esos 14 miRNAs, la prueba de razón de verosimilitud identificó 6 de ellos como que tienen un uso diferencial de isomiRs (Tabla 3). Finalmente, de los 6 miRNAs que se infiere que tienen un uso diferencial de isomiRs, 5 tenían estimaciones de logFC más pequeñas del marco de modelado miRglmm que las estimaciones de logFC de DESeq2.

Conclusiones

En los experimentos de perfilamiento de miRNAs, hay múltiples secuencias de isomiRs que se asignan a un miRNA, y los estudios han demostrado que estos isomiRs se producen biológicamente y son funcionales. Sin embargo, los isomiRs no se producen exclusivamente a través de vías biosintéticas. La variación técnica resultante de los errores en el proceso de secuenciación de los isomiRs biológicos verdaderos también puede producir secuencias de isomiRs. Para reducir el impacto de estos errores en el análisis de isomiRs, hemos desarrollado PARiS, un método para reasignar las secuencias erróneas a su posible origen.

Demostramos aquí, tanto en un entorno de simulación como en un conjunto de datos de referencia experimental sintético, que la eliminación de ruido con PARiS y la estimación de la expresión diferencial a nivel de miRNA utilizando miRglmm mejoran la inferencia. Las mejoras en la inferencia se reflejan tanto en la disminución del error cuadrático medio (MSE) de las estimaciones de logFC a nivel de miRNA como en el aumento de la proporción de cobertura de los intervalos de confianza del 95% estimados. Comparamos nuestro flujo de trabajo analítico propuesto con dos flujos de trabajo comúnmente aplicados. El primero aplica miREC a los datos para identificar y corregir las lecturas de error, y luego estima la expresión diferencial a nivel de miRNA utilizando miRglmm. El segundo agrega los recuentos a nivel de secuencia al nivel de miRNA y luego estima la expresión diferencial a nivel de miRNA utilizando DESeq2. Idealmente, después de utilizar un método de eliminación de ruido (PARiS, miREC o agregación), podríamos estimar la expresión diferencial a nivel de miRNA utilizando el mismo método. Sin embargo, actualmente no existe tal método. Se implementa un marco de modelado de efectos mixtos en miRglmm y requiere al menos 2 secuencias que se asignen a un miRNA dado para estimar los efectos aleatorios. Por otro lado, DESeq2 requiere que las lecturas se colapsen al nivel de miRNA, lo que equivale a asumir que hay 1 secuencia de isomiR biológico verdadero por miRNA. Un marco de modelado ideal de expresión diferencial a nivel de miRNA sería capaz de manejar ambos escenarios simultáneamente. No obstante, en nuestro análisis de los 50 conjuntos de datos de monocitos simulados, podemos comparar el rendimiento de miRglmm en los datos tratados con miREC y PARiS con el rendimiento de miRglmm en los datos sin procesar para tener una idea de cómo la eliminación de ruido cambia el rendimiento de miRglmm. Suponiendo que estamos utilizando miRglmm para estimar la expresión diferencial a nivel de miRNA, vemos que la eliminación de ruido con PARiS, para y un amplio rango de valores de , mejora la inferencia en comparación con los datos sin procesar.

Finalmente, aplicamos PARiS con y a un análisis real de 3 líneas celulares de adenocarcinoma colorrectal: DLD-1, DKS-8 y DKO-1. Las tres líneas celulares difieren en el estado de desactivación de sus alelos KRAS, donde DKS-8 es de tipo salvaje, DLD-1 es un mutante KRAS heterocigoto y DKO-1 es un mutante KRAS homocigoto. Debido a que la línea celular DLD-1 es heterocigota, su nivel de expresión de KRAS se situaría entre los niveles de expresión de las líneas celulares DKO-1 y DKS-8. Aunque los logFC más pequeños son más difíciles de detectar (como se espera entre DLD-1 y las otras líneas celulares), el enfoque de eliminar el ruido con PARiS y estimar la expresión diferencial con miRglmm puede identificar estos cambios. Por ejemplo, miRglmm identifica que hsa-miR-200c-3p se expresa de forma diferencial entre las líneas celulares DLD-1 y DKO-1; DESeq2 no hace este descubrimiento. Estudios anteriores que examinan la relación entre la expresión de miRNAs y los cánceres colorrectales han identificado una regulación al alza significativa de miR-200c-3p en las células de cáncer colorrectal con una mutación del gen KRAS (^[35]). Por lo tanto, es probable que el enfoque basado en isomiRs PARiS/miRglmm identifique resultados más relevantes desde el punto de vista biológico que el enfoque de agregación/DESeq2.

Un tratamiento a nivel de secuencia de los datos, que se basa en tener una versión de alta calidad y sin ruido de los datos, también evita que los investigadores confundan el uso diferencial de isomiRs con la expresión diferencial a nivel de miRNA. Por ejemplo, en una comparación de los perfiles de expresión de miRNAs entre las líneas celulares DKS-8 y DKO-1, hay 14 miRNAs que se identifican como que se expresan de forma diferencial por DESeq2 y no por miRglmm (Figura 4, 2). Baran et al. señaló que es posible que cuando haya un uso diferencial de isomiRs, miRglmm devolverá estimaciones de logFC que estén más cerca de 0 que los métodos de agregación, no identificará esos miRNAs como que se expresan de forma diferencial, pero identificará esos miRNAs como que tienen un uso diferencial de isomiRs (^[15]). De los 14 miRNAs identificados como que se expresan de forma diferencial, 6 tienen un uso diferencial de isomiRs (Tabla 3). De esos 6, miRglmm devolvió estimaciones de logFC más cercanas a 0 para 5 de ellos que DESeq2. Por lo tanto, analizar los datos a nivel de secuencia permite a los investigadores mejorar la inferencia al identificar potencialmente miRNAs que se expresan de forma diferencial con logFC de magnitudes más pequeñas y no confundir el uso diferencial de isomiRs con la expresión diferencial de miRNAs.

En conclusión, hemos desarrollado PARiS, un algoritmo que elimina el ruido de los datos de secuenciación de isomiRs a través de la asignación y el reordenamiento probabilísticos de las secuencias de isomiRs. En los datos simulados a partir del modelo de error propuesto utilizado por PARiS, el modelo funciona como se espera. En entornos en los que se conoce la expresión diferencial real medida por el logFC, la eliminación de ruido con PARiS y la estimación de la expresión diferencial utilizando miRglmm mejoran la inferencia. La mejora en la inferencia se muestra tanto en la disminución del MSE como en el aumento de la proporción de cobertura de los intervalos de confianza del 95%. Además, un análisis de datos reales que utiliza PARiS mejoró la inferencia al identificar correctamente los miRNAs con logFC de menor magnitud como que se expresan de forma diferencial y no confundir los miRNAs con un uso diferencial de isomiRs como que se expresan de forma diferencial.