Um pacote de fácil utilização, autoun para análises de metilação de DNA
Para completar a análise de dados de metilação de DNA de forma mais conveniente, empacotamos todas as funções em um pacote de fácil utilização, autoun para análise de metilação de DNA. A Figura 1 mostra as principais características do BatMeth2: 1) O BatMeth2 tem um desempenho de alinhamento eficiente e preciso. 2) O BatMeth2 pode calcular o nível de metilação de DNA (ML) de locais individuais de citosina ou quaisquer regiões funcionais, tais como cromossomos inteiros, regiões gênicas, elementos transponíveis (TEs), etc. 3) Após a integração de diferentes algoritmos estatísticos, o BatMeth2 pode realizar análise diferencial de metilação de DNA para qualquer região, qualquer número de amostras de entrada e requisitos do usuário. 4) Ao integrar a visualização dos dados do BS-Seq (distribuição da metilação do DNA em cromossomos e genes) e a anotação da metilação diferencial, o BatMeth2 pode visualizar os dados da metilação do DNA com mais clareza. Durante a execução da ferramenta BatMeth2, um relatório html é gerado para as estatísticas da amostra. Detalhes do relatório html da amostra são mostrados em http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.
O fluxo de trabalho do BatMeth2. As duas grandes setas significam arquivos de entrada ou saída
BatMeth2 tem melhor desempenho no mapeamento de dados simulados do BS-Seq
Primeiro avaliamos todos os alinhadores usando conjuntos de dados simulados (sem indels) consistindo em leituras com 75 pares de bases (bp), 100 bp e 150 bp e com diferentes taxas de conversão de bissulfito (variando de 0 a 100% com o passo 10%). Estes conjuntos de dados foram simulados a partir do genoma humano (UCSC hg19) usando as ferramentas FASTX-mutate, wgsim (v0.3.0) e o simulador em SAMtools (v1.1) , o que permite 0,03% de indels, uma taxa de erro base de 1% em todo o genoma e um máximo de dois desajustes por leitura. Mapeamos as leituras simuladas para o genoma de referência, permitindo no máximo dois desajustes. Como as posições originais das leituras simuladas eram conhecidas, pudemos avaliar a precisão de todos os programas, comparando suas saídas de mapeamento com as posições originais.
Para comparar as performances dos diferentes softwares, uma leitura sequencial com indels foi considerada corretamente mapeada se as seguintes condições fossem verdadeiras: 1) a leitura foi unicamente mapeada para a mesma vertente da leitura simulada e a qualidade do mapeamento foi maior que 0; 2) a posição inicial relatada da leitura alinhada estava dentro de dez pares de bases da posição inicial original da leitura simulada; 3) os resultados do mapeamento tinham indels ou desajustes semelhantes à leitura simulada. Se alguma dessas condições fosse violada, a leitura era considerada erroneamente mapeada. Como o BatMeth2 permite uma lacuna na região das sementes, ele pode encontrar locais de sementes incorporando indels com alta precisão e pode evitar locais desalinhados, o que faria com que as leituras incorporando indels fossem desalinhadas. Os resultados na Fig. 2 mostram que o BatMeth2 atingiu o maior número de leituras corretamente alinhadas e o menor número de leituras incorretamente alinhadas em todos os conjuntos de dados de teste com diferentes taxas de conversão de bissulfito.
Avaliação de todos os alinhadores BS-Seq usando conjuntos de dados simulados com diferentes comprimentos de leitura de FASTX e wgsim. Os dados simulados com diferentes taxas de conversão de bissulfito são mostrados em diferentes formas. Resultados de alinhadores diferentes são mostrados com cores diferentes dos símbolos. Os resultados perto do canto superior esquerdo em cada painel mostram que o software conseguiu mais leituras mapeadas corretamente e o menor número de leituras mapeadas incorretamente. Os resultados do nosso alinhador BatMeth2 são os melhores nos diferentes conjuntos de dados simulados de bissulfito
Em resumo, os resultados dos conjuntos de dados indel-aberrantes simulados de wgsim mostram que o BatMeth2 tem melhor desempenho (1~2% melhor do que o segundo alinhador superior) do que os outros métodos ao alinhar leituras BS simuladas gerais contendo uma mistura de desajustes e indels. Podemos ver que com o aumento da taxa de conversão do BS, a precisão do alinhamento de todo o software é reduzida. Nessas diferentes condições, o BatMeth2 tem um melhor desempenho.
BatMeth2 tem melhor desempenho no mapeamento de dados reais do BS-Seq
Para testar o desempenho do BatMeth2 em conjuntos de dados reais do BS-Seq, nós baixamos conjuntos de dados de BS-Seq em ponta de pizza e extraímos aleatoriamente 1 milhão de leituras em ponta de pizza de 2 × 90 bp do SRA SRR847318, 1 milhão 2 × 101 bp paired-end lê de SRA SRR1035722 e 1 milhão 2 × 125 bp paired-end lê de SRA SRR3503136 para fins de avaliação. Como estes conjuntos de dados são de linhas celulares ou tecidos saudáveis, espera-se que contenham um baixo número de variações estruturais. Portanto, alinhamos esses dados reais usando leituras de ponta única a partir dos conjuntos de dados de ponta de tela e avaliamos as taxas de mapeamento concordantes e discordantes dos alinhamentos pareados para estimar as taxas de alinhamento corretas e incorretas. Como o tamanho de inserção das leituras em ponta de par era de aproximadamente 500 bp, um par de leituras de parceiro poderia ser considerado concordante se fosse mapeado dentro de uma distância nominal de 500 bp; caso contrário, um par de leituras de parceiro poderia ser considerado discordante. Similar aos nossos resultados com os dados simulados, BatMeth2 relatou mais concordância e menos alinhamentos discordantes nos conjuntos de dados reais em uma grande variedade de escores de qualidade de mapas, como mostrado na Fig. 3.
Concordância e taxas de discordância de alinhamentos em leituras reais de ponta de pizza a partir de alinhadores diferentes. Contagens cumulativas de alinhamentos concordantes e discordantes de alta a baixa qualidade de mapeamento para leituras reais de seqüenciamento de bissulfito. Há apenas um ponto para o BSmap e os alinhadores baseados em bowtie separadamente, uma vez que esses alinhadores não têm pontuação de qualidade de mapa. Bismark-bowtie2L15 significa alinhamento bowtie2 com comprimento de semente 15
Além disso, a Tabela 1 mostra os tempos relativos de execução dos programas. BatMeth2 com as configurações padrão correu mais rápido do que a maioria dos alinhadores publicados e foi comparável ao BWA-meth e BatMeth. Bismark2 (com Bowtie2 como método de mapeamento fundamental), BS Seeker2 e BSmooth requerem tempos de execução mais longos.
Chamada de metilação de DNA
Para avaliar a precisão da chamada de metilação de DNA entre diferentes softwares, nós baixamos dados do chip de 450 K de esferas da linha de células IMR90 da ENCODE (Enciclopédia de Elementos de DNA). Também baixamos dados de seqüenciamento de bissulfito de genoma inteiro (WGBS-Seq) da linha de células IMR90 da ENCODE (42,6 Gbases). Para cada software, alinhamos as leituras do WGBS-Seq e calculamos o nível de metilação do DNA. Em seguida, comparamos os resultados com os MLs nos mesmos locais nos dados do chip de 450 K Bead Chip. Quando a diferença entre o DNA ML dos dados do WGBS-Seq pelo software e o do chip de esferas de 450 K era inferior a 0,2, o resultado da chamada foi definido como correto; caso contrário, foi considerado incorreto.
Os resultados são mostrados na Tabela 2. A sobreposição entre os resultados corretos de todo o software é mostrada no arquivo Adicional 1: Figura S2. Podemos ver que BatMeth2 e Biscuit têm desempenhos semelhantes, que são melhores do que os do outro software. Em conclusão, BatMeth2 melhora a precisão tanto do alinhamento de leitura de BS quanto da chamada de DNA ML.
BatMeth2 alinha as leituras de BS enquanto permite indels de comprimento variável
Cancer contém uma proporção notavelmente maior de indels do que as células saudáveis. Portanto, para verificar se o BatMeth2 pode alinhar as leituras de BS com indels de comprimentos diferentes, nós baixamos dados WGBS (75 Gbases) e dados de 450 K Bead Chip do HepG2 (carcinoma hepatocelular hepático, uma linha de células cancerígenas) do ENCODE. Verificamos a distribuição do comprimento do indel nas leituras após o alinhamento dos dados do HepG2 WGBS-Seq. Arquivo adicional 1: A figura S3A mostra que os comprimentos dos indels detectados foram distribuídos principalmente na faixa de 1 bp~ 5 bp, e o indel mais longo foi de 40 bp de comprimento. De acordo com nossas estatísticas, 2,3% das leituras de alinhamento continham os indels. A partir destes resultados, sabemos que o BatMeth2 pode alinhar leituras com indels de diferentes comprimentos.
Próximo, testamos o efeito da detecção do indel na chamada de metilação do DNA. Para o BatMeth2, executamos duas opções nos dados do HepG2: com e sem detecção de indel (ou seja, definir o parâmetro -I no BatMeth2). Também executamos Bismark nos dados WGBS-Seq do HepG2 como referência para a chamada de metilação de DNA com detecção de indel, porque Bismark não tem uma função de chamada de indel. Comparamos a chamada de metilação de DNA no BatMeth2 e Bismark com a chamada dos dados do chip de 450 K Bead Chip. Os resultados são mostrados no arquivo adicional 1: Figura S3B, onde “BatMeth2-noIndel” corresponde ao BatMeth2 sem detecção de indel. Podemos ver que, na ausência de detecção de indel, o resultado do BatMeth2 foi apenas ligeiramente melhor do que o do Bismark (com Bowtie1 como método de mapeamento fundamental). O resultado do BatMeth2 com detecção de indel foi significativamente melhor. Além disso, podemos ver que o BatMeth2 pode detectar mais locais de metilação de DNA do que o BatMeth2-noIndel e o Bismark (Bowtie 1). Para entender porque o desempenho do BatMeth2 com detecção de indel é melhor, definimos os locais de metilação chamados pelo BatMeth2 como Resultado A, enquanto os locais de metilação chamados pelo BatMeth2-noIndel e Bismark foram definidos como Resultado B. Então, deixamos o mclA ser os locais de metilação que aparecem no Resultado A mas não no Resultado B. Observamos que o mclA incluía 23.853 locais de metilação de DNA e 15.048 (63%) dos 23.853 locais cobertos pelos alinhamentos de indel chamados pelo BatMeth2 com detecção de indel (ver arquivo adicional 1: Figura S3C). Além disso, descobrimos que as taxas de indel no Resultado A e no Resultado B foram de apenas 5 e 0%, respectivamente. Assim, concluímos que a detecção precisa do indel pode melhorar a chamada de metilação do DNA.
Visualização dos dados de metilação do DNA
BatMeth2 fornece ferramentas para visualizar os dados de metilação. Para ilustrar os recursos de visualização do BatMeth2, nós baixamos (1) 117 Gbases de leituras de ponta única da linha de células humanas H9, (2) 105.2 Gbases de leituras de ponta única da linha de células humanas IMR90 e (3) 12.6 Gbases de leituras de ponta de arroz tipo selvagem. Primeiro, o BatMeth2 pode visualizar a densidade de metilação da citosina a nível cromossômico. Os pontos na Fig. 4a representam uma janela deslizante de 100 kb com um passo de 50 kb. Para permitir a visualização do ML em sites individuais CpG ou não-CpG num navegador do genoma, também fornecemos arquivos nos formatos cama e bigWig (Fig. 4b). Comparando com a densidade dos genes e TE, observamos que o ML estava correlacionado com a densidade do TE e era anticorrelacionado com a densidade do gene (Fig. 4c). Esta tendência foi previamente observada no arroz .
Visualização dos níveis de metilação na escala cromossômica. a A densidade de metil-citosina no cromossomo humano 10. Os pontos representam os níveis de metilação em janelas de correr de 100Kb com um degrau de 50Kb. Os pontos vermelhos referem-se aos níveis de metilação no cordão mais, e os pontos azuis referem-se aos níveis de metilação no cordão menos. b Um exemplo sobre as distribuições dos níveis de metilação do DNA e regiões diferencialmente metiladas (DMRs) entre as linhas de células H9 e IMR90 no cromossoma humano 10. c A densidade dos genes, elementos transpositores (TEs) e o nível de metilação do DNA em todo o genoma do arroz. Painel A são os resultados gerados a partir do Batmeth2. Painel B é o resultado de visualização do navegador UCSC, com os arquivos BED do Batmeth2
Segundo, BatMeth2 pode visualizar os MLs dos genes. Mais precisamente, o BatMeth2 pode visualizar os MLs 2 kb a montante do gene, no local de início da transcrição (TSS), no corpo do gene, no local final da transcrição (TES) e 2 kb a jusante do corpo do gene. Comparando as regiões a montante, corpo e jusante, a Fig. 5a mostra que o DNA ML do corpo do gene é maior do que o da região promotora. Comparando as cinco regiões, existe obviamente um vale na região do SST (Fig. 5b). BatMeth2 também pode calcular os perfis do ML em torno de introns, exons, regiões intergênicas e ETs (arquivo adicional 1: Figura S4). Além disso, o BatMeth2 pode fornecer um mapa de calor de múltiplos genes por região de genes para uma comparação conveniente do gene ML geral de diferentes amostras (Fig. 5c).
Visualização da metilação do DNA sob diferentes contextos. a Os níveis de metilação do DNA em regiões 2Kb a montante dos genes, corpos de genes, 2Kb a jusante dos corpos de genes. b O perfil de agregação dos níveis de metilação do DNA através dos genes. c O mapa de calor de todos os genes em regiões 2Kb a montante dos genes, corpos de genes, 2Kb a jusante dos corpos de genes
Terceiro, BatMeth2 pode visualizar a distribuição da metilação do DNA. Arquivo adicional 1: A figura S5A mostra a distribuição da metilação do DNA nas linhas de células H9 e IMR90. Na figura, o DNA ML é dividido em cinco categorias: metilado (M: > 80%), intermediário entre parcialmente metilado e metilado (Mh: 60-80%), parcialmente metilado (H: 40-60%), intermediário entre não metilado e parcialmente metilado (hU: 20-40%), e não metilado (U: < 20%). Como mostrado no arquivo adicional 1: Figura S5A, o ML foi maior na linha de célula H9 na categoria M do que na linha de célula IMR90, especialmente no contexto CpG. No contexto da seqüência CH, a metilação CpG é a forma predominante, mas uma fração significativa das citosinas metiladas é encontrada nos locais CpA, enquanto o ML é inferior a 40%, particularmente na linha de células H9 (arquivo adicional 1: Figura S5B).
Fourth, BatMeth2 pode analisar a correlação entre o nível de expressão gênica e o DNA promotor do gene ML. Nós ilustramos este recurso usando as linhas de células H9 e IMR90. Os níveis de expressão dos genes em H9 ou IMR90 foram divididos em diferentes categorias. Como mostrado no arquivo adicional 1: Figura S5C, os genes altamente expressos exibiam MLs inferiores nas suas regiões promotoras. Além disso, dividimos os MLs dos promotores dos genes em cinco categorias. O resultado no Arquivo Adicional 1: A Figura S5D mostra que os genes com promotores com valores de ML mais altos exibiam níveis de expressão mais baixos. A correlação negativa entre a expressão dos genes de mamíferos e a metilação do DNA promotor é conhecida . Esta análise indica ainda mais a precisão do BatMeth2.
De encontrar citosinas e regiões diferencialmente metiladas (DMCs/DMRs)
A identificação de citosinas diferencialmente metiladas (DMCs) e regiões diferencialmente metiladas (DMRs) é um dos principais objetivos na análise de dados de metilação. Embora os pesquisadores estejam ocasionalmente interessados em correlacionar locais de citosina simples a um fenótipo , as DMRs são características muito importantes .
Early BS-Seq studies profilled cells sem coleta de réplicas. Para tais conjuntos de dados, usamos o teste exato de Fisher para discernir as citosinas metiladas diferencialmente (DMCs). Para os conjuntos de dados BS-Seq com réplicas, o modelo estatístico mais natural para chamar DMCs é a distribuição beta-binomial . Sabemos que vários programas de software podem realizar análise diferencial de dados de metilação de DNA, como o MetilKit (um programa de análise diferencial que requer réplicas biológicas) e o Methy-Pipe (um programa de análise diferencial sem duplicação biológica). Entretanto, não há um pacote abrangente que inclua tanto o mapeamento quanto a análise diferencial de metilação disponível. Assim, nós desenvolvemos um pacote que integra o mapeamento com a análise diferencial. Para facilitar a identificação de DMRs a partir de dados de bissulfito sem réplicas, integramos o teste exato de Fisher para realizar um teste de hipóteses. Quando uma amostra tem duas ou mais réplicas, utilizamos a distribuição beta-binomial para realizar a análise diferencial de metilação. Nós também fornecemos arquivos de cama ou bigWig para a lista de DMRs. Os DMRs podem ser visualizados em um navegador de genoma (Fig. 4b) com o leito gerado ou arquivos bigWig.
Como ilustração, Fig. 6a mostra os números de DMCs e regiões na linha de células IMR90 e na linha de células H9, como detectado pelo BatMeth2 (valor de p< 0,05, meth.diff > = 0,6). O BatMeth2 pode visualizar se CpGs e DMCs são enriquecidos em algumas regiões, como gene, CDS, intron, intergênico, UTR, TE, LTR, LINE e regiões SINE. A Figura 6b visualiza as proporções dos CpGs em diferentes regiões genômicas. Além das regiões intergênicas, não observamos enriquecimento de DMC em nenhuma região.
Análise da metilação diferencial. a Análise de resultados de regiões diferencialmente metiladas (DMRs), genes diferencialmente metilados (DMGs) e promotores diferencialmente metilados (DMPs) entre linhas de células H9 e IMR90. b Anotação de Citosinas Diferencialmente Metiladas (DMC) contra diferentes propriedades genômicas e elementos repetidos. c Os DMPs contêm H9 ou IMR90-indels específicos (laranja) ocupam uma proporção substancial em todos os DMPs (DNA Methylation differential Promoters)
Uma proporção substancial de promotores diferencialmente metilados (DMPs) contém indels
Sabemos que os indels e a metilação do DNA desempenham um papel importante no desenvolvimento dos tecidos e das doenças . Aqui, nós examinamos a relação entre promotores diferencialmente metilados (DMPs) e indels. Realizamos este estudo usando as leituras do BS-Seq em linhas de células IMR90 e H9. Primeiro alinhamos as leituras do BS-Seq usando BatMeth2; depois, os indels foram chamados usando as ferramentas BisSNP e GATK. Posteriormente, definimos os indels que ocorrem apenas em H9 ou IMR90 como indels específicos de linha de célula.
Então, detectamos 1384 DMPs entre H9 e IMR90 pelo BatMeth2 (valor de p< 0.05, meth.diff > = 0.6). Um total de 236 (17%) entre todos os DMPs acima contém indels, como mostrado na Fig. 6c. Em resumo, uma proporção substancial dos DMPs contém indels. Portanto, o alinhamento preciso da leitura do BS-Seq perto desses indels é muito importante para a pesquisa e exploração da metilação do DNA.