Un paquete autorun fácil de usar para los análisis de metilación del ADN
Para completar el análisis de datos de metilación del ADN de forma más conveniente, empaquetamos todas las funciones en un paquete autorun fácil de usar para el análisis de metilación del ADN. La figura 1 muestra las principales características de BatMeth2: 1) BatMeth2 tiene un rendimiento de alineación eficiente y preciso. 2) BatMeth2 puede calcular el nivel de metilación del ADN (ML) de sitios individuales de citosina o de cualquier región funcional, como cromosomas enteros, regiones de genes, elementos transponibles (TEs), etc. 3) Tras la integración de diferentes algoritmos estadísticos, BatMeth2 puede realizar análisis diferenciales de metilación del ADN para cualquier región, cualquier número de muestras de entrada y requisitos del usuario. 4) Al integrar la visualización de los datos de BS-Seq (distribución de la metilación del ADN en cromosomas y genes) y la anotación de metilación diferencial, BatMeth2 puede visualizar los datos de metilación del ADN con mayor claridad. Durante la ejecución de la herramienta BatMeth2, se genera un informe html para las estadísticas de la muestra. Los detalles del informe html de la muestra se muestran en http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.
El flujo de trabajo de BatMeth2. Las dos flechas grandes significan archivos de entrada o de salida
BatMeth2 tiene un mejor rendimiento de mapeo en datos simulados de BS-Seq
En primer lugar, evaluamos todos los alineadores utilizando conjuntos de datos simulados (sin indels) que consistían en lecturas con 75 pares de bases (bp), 100 bp y 150 bp y con diferentes tasas de conversión de bisulfitos (que iban de 0 a 100% con un paso del 10%). Estos conjuntos de datos se simularon a partir del genoma humano (UCSC hg19) utilizando FASTX-mutate-tools , wgsim (v0.3.0) y el simulador en SAMtools (v1.1) , que permite un 0,03% de indels, una tasa de error de base del 1% en todo el genoma y un máximo de dos mismatches por lectura. Hemos mapeado las lecturas simuladas con el genoma de referencia, permitiendo como máximo dos desajustes. Como las posiciones originales de las lecturas simuladas eran conocidas, pudimos evaluar la precisión de todos los programas comparando sus resultados de mapeo con las posiciones originales.
Para comparar el rendimiento de los distintos programas, se consideró que una lectura de secuenciación con indels se había mapeado correctamente si se cumplían las siguientes condiciones: 1) la lectura se mapeaba de forma única en la misma hebra de la que se había simulado y la calidad del mapeo era superior a 0; 2) la posición inicial informada de la lectura alineada estaba dentro de los diez pares de bases de la posición inicial original de la lectura simulada; 3) los resultados del mapeo tenían indels o desajustes similares a los de la lectura simulada. Si se incumple alguna de estas condiciones, se considera que la lectura está mal mapeada. Dado que BatMeth2 permite una brecha en la región de siembra, puede encontrar ubicaciones de siembra que incorporan indels con gran precisión y puede evitar ubicaciones desajustadas, que harían que las lecturas que incorporan indels estuvieran mal alineadas. Los resultados de la Fig. 2 muestran que BatMeth2 logró el mayor número de lecturas correctamente alineadas y el menor número de lecturas incorrectamente alineadas en todos los conjuntos de datos de prueba a diferentes tasas de conversión de bisulfitos.
Evaluación de todos los alineadores BS-Seq utilizando conjuntos de datos simulados con diferentes longitudes de lectura de FASTX y wgsim. Los datos simulados con diferentes tasas de conversión de bisulfito se muestran con diferentes formas. Los resultados de los diferentes alineadores se muestran con diferentes colores de los símbolos. Los resultados cerca de la esquina superior izquierda de cada panel muestran que el software consiguió un mayor número de lecturas correctamente mapeadas y un menor número de lecturas incorrectamente mapeadas. Los resultados de nuestro alineador BatMeth2 son los mejores en los diferentes conjuntos de datos de bisulfito simulados
En resumen, los resultados de los conjuntos de datos de indel-aberrante simulados por wgsim muestran que BatMeth2 tiene un mejor rendimiento (1~2% mejor que el segundo alineador superior) que los otros métodos cuando se alinean lecturas generales de BS simuladas que contienen una mezcla de mismatches e indels. Podemos ver que con el aumento de la tasa de conversión de BS, la precisión de la alineación de todos los programas se reduce. En estas condiciones diferentes, BatMeth2 se comporta mejor.
BatMeth2 tiene un mejor rendimiento de mapeo en datos reales de BS-Seq
Para probar el rendimiento de BatMeth2 en conjuntos de datos reales de BS-Seq, descargamos conjuntos de datos de BS-Seq de extremo emparejado y extrajimos aleatoriamente 1 millón de lecturas de extremo emparejado de 2 × 90 pb de SRA SRR847318, 1 millón de lecturas de 2 × 101 pb de SRA SRR1035722 y 1 millón de lecturas de 2 × 125 pb de SRA SRR3503136 para su evaluación. Dado que estos conjuntos de datos proceden de líneas celulares o tejidos sanos, se espera que contengan un bajo número de variaciones estructurales. Por lo tanto, alineamos estos datos reales utilizando lecturas de extremo único de los conjuntos de datos de extremo emparejado y evaluamos las tasas de mapeo concordantes y discordantes de las alineaciones emparejadas para estimar las tasas de alineación correcta e incorrecta. Dado que el tamaño de la inserción de las lecturas de extremo emparejado era de aproximadamente 500 pb, un par de lecturas asociadas podía considerarse concordante si se mapeaba dentro de una distancia nominal de 500 pb; de lo contrario, un par de lecturas asociadas podía considerarse discordante. Al igual que nuestros resultados con los datos simulados, BatMeth2 reportó más alineaciones concordantes y menos discordantes en los conjuntos de datos reales en un amplio rango de puntuaciones de calidad del mapa, como se muestra en la Fig. 3.
Tasas de concordancia y discordancia de los alineamientos en las lecturas reales de extremos pareados de diferentes alineadores. Recuentos acumulativos de alineaciones concordantes y discordantes de alta a baja calidad de mapeo para lecturas reales de secuenciación por bisulfito. Sólo hay un punto para BSmap y los alineadores basados en bowtie por separado, ya que estos alineadores no tienen puntuación de calidad de mapa. Bismark-bowtie2L15 significa alineación bowtie2 con longitud de semilla 15
Además, la Tabla 1 muestra los tiempos de ejecución relativos de los programas. BatMeth2, con la configuración por defecto, funcionó más rápido que la mayoría de los alineadores publicados y fue comparable a BWA-meth y BatMeth. Bismark2 (con Bowtie2 como método fundamental de mapeo), BS Seeker2 y BSmooth requieren tiempos de ejecución más largos.
Llamada a la metilación del ADN
Para evaluar la precisión de la llamada a la metilación del ADN entre los diferentes softwares, descargamos datos de 450 K de chips de cuentas de la línea celular IMR90 de ENCODE (Enciclopedia de Elementos de ADN). También descargamos los datos de la secuenciación de bisulfitos del genoma completo (WGBS-Seq) de la línea celular IMR90 de ENCODE (42,6 Gbases). Para cada software, alineamos las lecturas WGBS-Seq y calculamos el nivel de metilación del ADN. A continuación, comparamos los resultados con los ML en los mismos sitios en los datos de 450 K Bead Chip. Cuando la diferencia entre el ML de ADN de los datos de WGBS-Seq por el software y el de los 450 K Bead Chip era inferior a 0,2, el resultado de la llamada se definía como correcto; en caso contrario, se consideraba incorrecto.
Los resultados se muestran en la Tabla 2. El solapamiento entre los resultados correctos de todos los programas informáticos se muestra en el archivo adicional 1: Figura S2. Podemos ver que BatMeth2 y Biscuit tienen rendimientos similares, que son mejores que los del otro software. En conclusión, BatMeth2 mejora la precisión tanto de la alineación de las lecturas BS como de la llamada de ADN ML.
BatMeth2 alinea las lecturas BS al tiempo que permite indels de longitud variable
El cáncer contiene una proporción notablemente mayor de indels que las células sanas. Por lo tanto, para verificar si BatMeth2 puede alinear lecturas BS con indels de diferentes longitudes, descargamos los datos de WGBS (75 Gbases) y los datos de 450 K Bead Chip de HepG2 (carcinoma hepatocelular de hígado, una línea celular de cáncer) de ENCODE. Comprobamos la distribución de la longitud de los indeles en las lecturas tras la alineación de los datos WGBS-Seq de HepG2. Archivo adicional 1: La Figura S3A muestra que las longitudes de los indels detectados se distribuyeron principalmente en el rango de 1 pb~ 5 pb, y el indel más largo fue de 40 pb de longitud. Según nuestras estadísticas, el 2,3% de las lecturas de alineación contenían indels. A partir de estos resultados, sabemos que BatMeth2 puede alinear lecturas con indels de diferentes longitudes.
A continuación, probamos el efecto de la detección de indels en la llamada de metilación del ADN. Para BatMeth2, ejecutamos dos opciones en los datos de HepG2: con y sin detección de indels (es decir, establecer el parámetro -I en BatMeth2). También ejecutamos Bismark en los datos de WGBS-Seq de HepG2 como referencia para la llamada de metilación del ADN con detección de indel, porque Bismark no tiene una función de llamada de indel. Comparamos la llamada de metilación del ADN en BatMeth2 y Bismark con la llamada de los datos de 450 K Bead Chip. Los resultados se muestran en el archivo adicional 1: Figura S3B, donde «BatMeth2-noIndel» corresponde a BatMeth2 sin detección de indel. Podemos ver que, en ausencia de detección de indels, el resultado de BatMeth2 fue sólo ligeramente mejor que el de Bismark (con Bowtie1 como método de mapeo fundamental). El resultado de BatMeth2 con detección de indels fue significativamente mejor. Además, podemos ver que BatMeth2 puede detectar más sitios de metilación del ADN que BatMeth2-noIndel y Bismark (Bowtie 1). Para entender por qué el rendimiento de BatMeth2 con la detección de indels es mejor, definimos los sitios de metilación llamados por BatMeth2 como Resultado A, mientras que los sitios de metilación llamados por BatMeth2-noIndel y Bismark fueron definidos como Resultado B. Entonces, dejamos que mclA fuera los sitios de metilación que aparecían en el Resultado A pero no en el Resultado B. Observamos que mclA incluía 23.853 sitios de metilación del ADN y 15.048 (63%) de los 23.853 sitios cubiertos por las alineaciones de las lecturas de indel llamadas por BatMeth2 con detección de indel (ver Archivo Adicional 1: Figura S3C). Además, descubrimos que las tasas de indel en el Resultado A y el Resultado B eran sólo del 5 y el 0%, respectivamente. Por lo tanto, concluimos que la detección precisa de indels puede mejorar la llamada de metilación del ADN.
Visualización de los datos de metilación del ADN
BatMeth2 proporciona herramientas para visualizar los datos de metilación. Para ilustrar las características de visualización de BatMeth2, hemos descargado (1) 117 Gbases de lecturas de extremo único de la línea celular humana H9, (2) 105,2 Gbases de lecturas de extremo único de la línea celular humana IMR90 y (3) 12,6 Gbases de lecturas de extremo de par de arroz de tipo salvaje. En primer lugar, BatMeth2 puede visualizar la densidad de metilación de la citosina a nivel cromosómico. Los puntos en la Fig. 4a representan una ventana deslizante de 100 kb con un paso de 50 kb. Para permitir la visualización de la ML en sitios individuales CpG o no CpG en un navegador del genoma, también proporcionamos archivos en formatos bed y bigWig (Fig. 4b). Comparando con la densidad de genes y TEs, observamos que el ML estaba correlacionado con la densidad de TEs y estaba anticorrelacionado con la densidad de genes (Fig. 4c). Esta tendencia se ha observado previamente en el arroz.
Visualización de los niveles de metilación en la escala cromosómica. a La densidad de metil-citosina en el cromosoma 10 humano. Los puntos representan los niveles de metilación en ventanas deslizantes de 100Kb con un paso de 50Kb. Los puntos rojos se refieren a los niveles de metilación en la cadena positiva, y los puntos azules se refieren a los niveles de metilación en la cadena negativa. b Un ejemplo sobre las distribuciones de los niveles de metilación del ADN y las regiones diferencialmente metiladas (DMRs) entre las líneas celulares H9 e IMR90 en el cromosoma 10 humano. c La densidad de genes, elementos transposones (TEs) y el nivel de metilación del ADN en todo el genoma del arroz. El panel A son los resultados generados a partir de Batmeth2. El panel B son los resultados de la visualización del navegador UCSC, con los archivos BED de Batmeth2
En segundo lugar, BatMeth2 puede visualizar los MLs de los genes. Más concretamente, BatMeth2 puede visualizar los MLs 2 kb aguas arriba del gen, en el sitio de inicio de la transcripción (TSS), en el cuerpo del gen, en el sitio final de la transcripción (TES) y 2 kb aguas abajo del cuerpo del gen. Al comparar las regiones aguas arriba, cuerpo y aguas abajo, la Fig. 5a muestra que el ML del ADN del cuerpo del gen es mayor que el de la región promotora. Comparando las cinco regiones, es evidente que hay un valle en la región TSS (Fig. 5b). BatMeth2 también puede calcular los perfiles de ML alrededor de intrones, exones, regiones intergénicas y TEs (Archivo adicional 1: Figura S4). Además, BatMeth2 puede proporcionar un mapa de calor de múltiples genes por región génica para una conveniente comparación de los ML generales de genes de diferentes muestras (Fig. 5c).
Visualización de la metilación del ADN bajo diferentes contextos. a Los niveles de metilación del ADN en regiones de 2Kb aguas arriba de los genes, cuerpos de los genes, 2Kb aguas abajo de los cuerpos de los genes. b El perfil de agregación de los niveles de metilación del ADN a través de los genes. c El mapa de calor de todos los genes en regiones de 2Kb aguas arriba de los genes, cuerpos de los genes, 2Kb aguas abajo de los cuerpos de los genes
En tercer lugar, BatMeth2 puede visualizar la distribución de la metilación del ADN. Archivo adicional 1: Figura S5A muestra las distribuciones de metilación del ADN en las líneas celulares H9 e IMR90. En la figura, el ADN ML se divide en cinco categorías: metilado (M: > 80%), intermedio entre parcialmente metilado y metilado (Mh: 60-80%), parcialmente metilado (H: 40-60%), intermedio entre no metilado y parcialmente metilado (hU: 20-40%), y no metilado (U: < 20%). Como se muestra en el archivo adicional 1: Figura S5A, el ML fue mayor en la línea celular H9 en la categoría M que en la línea celular IMR90, especialmente en el contexto CpG. En el contexto de la secuencia CH, la metilación CpG es la forma predominante, pero una fracción significativa de citosinas metiladas se encuentra en sitios CpA, mientras que el ML es inferior al 40%, especialmente en la línea celular H9 (archivo adicional 1: Figura S5B).
En cuarto lugar, BatMeth2 puede analizar la correlación entre el nivel de expresión del gen y el ML del ADN promotor del gen. Ilustramos esta característica utilizando las líneas celulares H9 e IMR90. Los niveles de expresión de los genes en H9 o IMR90 se dividieron en diferentes categorías. Como se muestra en el archivo adicional 1 Figura S5C, los genes altamente expresados mostraron MLs más bajos en sus regiones promotoras. Además, dividimos los ML de los promotores de los genes en cinco categorías. El resultado en el archivo adicional 1: Figura S5D muestra que los genes con promotores que tienen valores de ML más altos exhibieron niveles de expresión más bajos. La correlación negativa entre la expresión de los genes de mamíferos y la metilación del ADN del promotor es conocida. Este análisis indica además la precisión de BatMeth2.
Encontrar citosinas y regiones metiladas diferencialmente (DMCs/DMRs)
La identificación de citosinas metiladas diferencialmente (DMCs) y regiones metiladas diferencialmente (DMRs) es uno de los principales objetivos en el análisis de datos de metilación. Aunque los investigadores están ocasionalmente interesados en correlacionar sitios de citosina individuales con un fenotipo, las DMRs son características muy importantes.
Los primeros estudios de BS-Seq perfilaban las células sin recoger réplicas. Para tales conjuntos de datos, utilizamos la prueba exacta de Fisher para discernir las citosinas metiladas diferencialmente (DMC). Para los conjuntos de datos de BS-Seq con réplicas, el modelo estadístico más natural para llamar a las DMC es la distribución beta-binomial. Sabemos que varios programas de software pueden realizar análisis de datos de metilación diferencial del ADN, como methylKit (un programa de análisis diferencial que requiere réplicas biológicas) y Methy-Pipe (un programa de análisis diferencial sin duplicación biológica). Sin embargo, no existe ningún paquete completo que incluya tanto el mapeo como el análisis diferencial de metilación. Por ello, desarrollamos un paquete que integra el mapeo con el análisis diferencial. Para facilitar la identificación de DMRs a partir de datos de bisulfito sin réplicas, integramos el test exacto de Fisher para realizar una prueba de hipótesis. Cuando una muestra tiene dos o más réplicas, utilizamos la distribución beta-binomial para realizar el análisis de metilación diferencial. También proporcionamos archivos bed o bigWig para la lista de DMRs. Las DMRs pueden visualizarse en un navegador del genoma (Fig. 4b) con los archivos bed o bigWig generados.
Como ilustración, la Fig. 6a muestra el número de DMCs y regiones en la línea celular IMR90 y en la línea celular H9, tal y como se detectó con BatMeth2 (valor p< 0,05, meth.diff > = 0,6). BatMeth2 puede visualizar si los CpGs y DMCs están enriquecidos en algunas regiones, como las regiones génicas, CDS, intrónicas, intergénicas, UTR, TE, LTR, LINE y SINE. La figura 6b visualiza las proporciones de DMCs en diferentes regiones genómicas. Aparte de las regiones intergénicas, no observamos enriquecimiento de DMC en ninguna región.
Análisis de metilación diferencial. a Resultados del análisis de las regiones metiladas diferencialmente (DMRs), genes metilados diferencialmente (DMGs) y promotores metilados diferencialmente (DMPs) entre las líneas celulares H9 e IMR90. b Anotación de las citosinas diferencialmente metiladas (DMC) frente a diferentes propiedades genómicas y elementos de repetición. c Los DMPs que contienen indels específicos de H9 o IMR90 (naranja) ocupan una proporción sustancial en todos los DMPs (Promotores diferencialmente metilados del ADN)
Una proporción sustancial de promotores diferencialmente metilados (DMPs) contienen indels
Sabemos que los indels y la metilación del ADN juegan un papel importante en el desarrollo de los tejidos y las enfermedades . Aquí, examinamos la relación entre los promotores metilados diferencialmente (DMPs) y los indels. Realizamos este estudio utilizando las lecturas de BS-Seq en las líneas celulares IMR90 y H9. Primero alineamos las lecturas de BS-Seq usando BatMeth2; luego, los indels fueron llamados usando las herramientas BisSNP y GATK. Posteriormente, definimos los indels que ocurren sólo en H9 o IMR90 como indels específicos de la línea celular.
Entonces, detectamos 1384 DMPs entre H9 e IMR90 por BatMeth2 (valor p< 0,05, meth.diff > = 0,6). Un total de 236 (17%) entre todos los DMPs anteriores contienen indels, como se muestra en la Fig. 6c. En resumen, una proporción sustancial de los DMPs contienen indels. Por lo tanto, la alineación precisa de las lecturas de BS-Seq cerca de estos indels es muy importante para la investigación y la exploración de la metilación del ADN.