Systematic benchmarking of small variant calling pipelines for long-read RNA sequencing data

Resumo Técnico: Avaliação Sistemática de Pipelines de Chamada de Pequenas Variantes para Dados de Sequenciamento de RNA de Leituras Longas

Definição do Problema
O sequenciamento de RNA de leituras longas (lrRNA-seq) permite a detecção de variantes resolvidas em nível de transcrito, vinculando mudanças genéticas a isoformas específicas e facilitando a análise de eventos alílicos e de splicing. No entanto, avaliações sistemáticas e neutras de pipelines de chamada de pequenas variantes (SNVs e indels) para lrRNA-seq permanecem limitadas. O desempenho das ferramentas existentes através de diferentes tecnologias de sequenciamento (Oxford Nanopore Technologies [ONT] e Pacific Biosciences [PacBio]), estratégias de alinhamento, escolhas de chamadores de variantes, contextos genômicos e faseamento de haplótipos a jusante não é totalmente compreendido. Além disso, não está claro se a transformação do arquivo de alinhamento (converter alinhamentos de RNA spliceados em representações contíguas semelhantes ao DNA) continua sendo necessária para chamadores de variantes especificamente otimizados para lrRNA-seq.

Metodologia
Os autores conduziram um benchmark sistemático usando conjuntos de dados Genome in a Bottle (GIB) compreendendo três linhagens celulares (HG002, HG004, HG005) sequenciadas com quatro protocolos de preparação de biblioteca ONT (cDNA R9.4.1, cDNA R10.4.1, dRNA002, dRNA004) e dois protocolos PacBio (Iso-Seq e Mas-Seq/Kinness). O estudo também estendeu a avaliação para conjuntos de dados recentes do LongBench (linhagens celulares de câncer de pulmão H211 e H526).

O fluxo de trabalho de benchmarking envolveu:

1. Alinhamento: As leituras foram alinhadas usando minimap2 e pbmm2 (um wrapper otimizado para PacBio). O estudo avaliou o impacto da transformação de alinhamento (usando SplitNCigarReads para mimetizar alinhamentos de DNA contíguos) versus o uso de alinhamentos de RNA spliceados brutos.
2. Chamada de Variantes: Cinco chamadores específicos para lrRNA-seq foram avaliados: Clair3-RNA, DeepVariant (usando o modelo MASSEQ para PacBio e um modelo ONT de DNA para ONT), isoLASER, longcallR e longcallR-nn. O GATK HaplotypeCaller foi incluído como um baseline de leitura curta.
3. Faseamento de Haplótipos: O desempenho do faseamento foi avaliado usando WhatsHap, LongPhase, HapCUT2, HiPhase e as saídas nativas de faseamento do isoLASER e longcallR.
4. Métricas de Avaliação: A precisão foi quantificada usando scores de precisão, recall e F1 contra chamadas de verdade fundamental (ground-truth) dentro de regiões de alta confiança, estratificadas por profundidade de sequenciamento, tipo de variante (SNV/indel) e contexto genômico (homopolímeros, repetições em tandem, conteúdo GC). O faseamento foi avaliado usando taxa de erro de troca (switch error rate), cobertura de faseamento e N50 de bloco de fase.

Principais Contribuições e Resultados

• Características de Sequenciamento como Determinante Primário: O estudo identifica as características de sequenciamento (química, taxa de erro, throughput e mapeabilidade) como o determinante mais forte do desempenho da chamada de variantes, explicando mais variância nos scores F1 e contagens de verdadeiros positivos do que a escolha do chamador de variantes. Os conjuntos de dados PacBio (particularmente Mas-Seq) superaram consistentemente os conjuntos de dados ONT devido à maior precisão de sequenciamento e a um maior número de bases chamáveis. Entre as químicas ONT, dRNA004 e cDNA R10 tiveram o melhor desempenho, enquanto o dRNA002 apresentou o pior desempenho.
• Desempenho do Chamador de Variantes:
  • Clair3-RNA demonstrou a aplicabilidade mais ampla e robusta para ambas as químicas ONT e PacBio, tanto para SNVs quanto para indels.
  • DeepVariant teve um desempenho forte em dados PacBio, particularmente para chamada de indels, mas teve um desempenho ruim em dados ONT devido à falta de um modelo dedicado de RNA-seq para ONT.
  • longcallR foi competitivo para chamada de SNV em dados ONT sem depender de modelos pré-treinados e forneceu saída de faseamento direta.
  • isoLASER mostrou alta precisão, mas recall mais baixo, particularmente em dados ONT, e foi conservador na chamada de variantes.
  • longcallR-nn mostrou desempenho competitivo em conjuntos de dados específicos, mas teve dificuldades com o dRNA002.
• Estratégias de Alinhamento:
  • Escolha do Alinhador: O minimap2 geralmente gerou scores F1 mais altos do que o pbmm2, pois o pbmm2 tendeu a introduzir mais falsos positivos, particularmente em regiões de splice, provavelmente devido a diferenças nos presets de parâmetros e à falta de entrada de anotação de junção de splice.
  • Transformação de BAM: Para chamadores específicos de lrRNA-seq, a transformação de alinhamento (divisão de leituras em junções éxon-íntron) não melhorou o desempenho e, em alguns casos, reduziu a precisão ao interromper características de alinhamento específicas de RNA (ex: padrões de CIGAR sensíveis a splice) e informações de haplótipo. A transformação permaneceu benéfica apenas para o DeepVariant em conjuntos de dados PacBio para detecção de indels.
• Contexto Genômico e Edição de RNA:
  • O desempenho degradou em contextos genômicos desafiadores, particularmente em homopolímeros (especialmente A/T) e regiões de GC extremo. A chamada de indels foi mais sensível ao comprimento de homopolímeros do que a chamada de SNV.
  • A edição de RNA (A-para-I) contribuiu para falsos positivos. Chamadores que integraram bancos de dados de edição de RNA (ex: marcação de sítios via REDIportal no Clair3-RNA) ou usaram faseamento iterativo (longcallR) foram melhores em distinguir variantes reais de eventos de edição.
• Faseamento de Haplótipos: Nenhum método de faseamento único foi ideal para todas as métricas. O WhatsHap favoreceu a cobertura de faseamento e contiguidade (blocos de fase mais longos) ao custo de taxas de erro de troca mais altas, enquanto o HapCUT2 priorizou a precisão com uma cobertura ligeiramente menor. O longcallR forneceu uma alternativa forte para chamada de variantes e faseamento conjuntos.
• Generalização: As tendências de desempenho observadas nos conjuntos de dados GIAB foram consistentes com os resultados dos conjuntos de dados de linhagem celular de câncer do LongBench, reforçando a dominância das características de sequenciamento e a hierarquia relativa de desempenho dos chamadores.

Significância
O artigo estabelece um benchmark fundamental para a chamada de variantes de lrRNA-seq, destacando que, embora a seleção da ferramenta importe, os fatores subjacentes de tecnologia de sequenciamento e preparação de biblioteca são os mais críticos para o sucesso. Ele fornece orientações específicas para a construção de pipelines:

1. Clair3-RNA é recomendado como a ferramenta mais versátil para chamada geral de variantes de lrRNA-seq.
2. DeepVariant é uma escolha forte para dados PacBio, especialmente para indels.
3. Transformação de alinhamento deve geralmente ser omitida para chamadores específicos de lrRNA-seq para preservar informações de longo alcance, exceto ao usar DeepVariant em dados PacBio.
4. WhatsHap e HapCUT2 oferecem trade-offs complementares para o faseamento a jusante, enquanto o longcallR é adequado para chamada e faseamento integrados.

O estudo conclui que melhorias futuras na sensibilidade podem exigir estratégias de enriquecimento direcionadas para superar limitações dependentes de expressão, e que mais trabalho é necessário para avaliar a atribuição de variantes em nível de transcrito e aplicações de célula única.