Montagem e anotação de um genoma de referência humano Ashkenazi

Para a criação do primeiro genoma de referência humano a ser montado a partir de um único indivíduo, escolhemos HG002, um indivíduo Ashkenazi que faz parte do Projeto Genoma Pessoal (PGP). O PGP usa o Open Consent Model, a primeira plataforma de acesso verdadeiramente aberto para compartilhar genoma humano individual, fenótipo e dados médicos. O processo de consentimento educa potenciais participantes sobre as implicações e riscos do compartilhamento de dados genômicos, e sobre o que eles podem esperar de sua participação. O consentimento aberto permitiu a criação do primeiro material de referência do genoma humano do mundo (HG002 é NIST Reference Material 8391) do Genome In A Bottle (GIAB), que está sendo usado para calibração, desenvolvimento de métodos de montagem do genoma e medições de desempenho de laboratório. Todos os dados de seqüência bruta para este projeto foram obtidos do GIAB, onde está disponível gratuitamente ao público .

Nós montamos o genoma HG002 a partir de uma combinação de três conjuntos de dados de cobertura profunda: 249-bp Illumina lê, Oxford Nanopore (ONT) lê em média mais de 33 kbp de comprimento, e PacBio “HiFi” de alta qualidade lê em média 9567 bp (Tabela 1).

Criamos inicialmente duas montagens: uma usando as leituras Illumina e ONT, e uma segunda usando todos os três conjuntos de dados, incluindo as leituras do PacBio HiFi. A adição dos dados do PacBio HiFi resultou em uma seqüência ligeiramente mais total no conjunto (2.99 Gb vs. 2.88 Gb) e um contig N50 substancialmente maior (18.2 Mb vs. 4.9 Mb). Este conjunto, designado Ash1 v0.5, foi a base para todos os refinamentos subsequentes.

Maptando o conjunto para cromossomas

Para criar atribuições cromossómicas para o conjunto Ash1 v0.5, usamos alinhamentos com GRCh38 para mapear a maioria dos andaimes em cromossomas. Os passos descritos na seção “Métodos” geraram uma série de conjuntos cromossômicos gradualmente melhorados, resultando em Ash1 v1.7. Ash1 v1.7 tem maior contiguidade e lacunas menores que GRCh38, como mostrado na Tabela 2. Note que no processo de construção desses cromossomos, uma pequena quantidade da seqüência GRCh38 (58,3 Mb, 2% do genoma) foi usada para preencher lacunas em Ash1. Estas regiões incluem algumas regiões difíceis de montar que foram curadas manualmente para o GRCh38. No total, o tamanho estimado de todas as lacunas em Ash1 é 82.9 Mbp, contra 84.7 Mbp em GRCh38.p13.

Como parte do processo de melhoria da montagem, pesquisamos uma das montagens preliminares de Ash1 (v1.1) para as 12.745 variantes estruturais isoladas de alta qualidade (inserções e exclusões ≥ 50 bp) que Zook et al. identificaram ao comparar os dados do trio Ashkenazi com o GRCh37 . Esse estudo utilizou quatro tecnologias diferentes de seqüenciamento e vários chamadores de variantes para identificar variantes e filtrar os falsos positivos. Destes 12.745 SVs, 5807 são homozigotos e 6938 são heterozigotos. Esperávamos que o conjunto Ash1 concordasse com quase todas as variantes homozigotos. Como Ash1 captura apenas um haplótipo, esperávamos que ele concordasse com aproximadamente metade das SVs heterozigotas, assumindo que o algoritmo de montagem escolhesse aleatoriamente entre os haplótipos ao decidir qual variante incluir no consenso final. Das 5807 variantes homozigotas, 5284 (91%) estavam presentes usando nosso critério de combinação (ver seção “Métodos”), e 3922 (56,5%) das 6938 variantes heterozigotas estavam presentes. Todas as variantes foram encontradas no local correto.

Também fizemos pequenas chamadas de variantes (≤ 5 bp) em Ash1 v1.1 e comparamos com as variantes de referência HG002 v4.0 do GIAB, que usamos para corrigir inúmeros erros de substituição e indel (ver a seção “Métodos”), produzindo Ash 1 v1.2. Em seguida, realinhamos o conjunto Ash1 ao GRCh38, recalculamos as variantes e comparamos essas variantes com o novo conjunto de referência GIAB v4.1 desenvolvido. Das variantes dentro das regiões de referência v4.1, as variantes Ash1 corresponderam a 1.256.458 homozigotos e 1.041.476 heterozigotos SNP, e 187.227 homozigotos e 193.524 indels heterozigotos. Após excluir chamadas de variantes dentro de 30 bp de uma variante verdadeira, restaram 79.269 SNPs e 17.439 indels, o que (assumindo que todos são erros em Ash1) corresponde a um valor de qualidade (QV) de aproximadamente Q45 para erros de substituição. A maioria dessas variantes (52.191 SNPs e 4629 indels) cai em duplicações segmentares, possivelmente representando duplicações ausentes em Ash1 ou polimento imperfeito por leituras curtas. Em resumo, a qualidade do conjunto Ash1 é muito alta, com um valor estimado de qualidade de substituição de 62 e uma taxa de erro indel de 2 por milhão de bp após excluir duplicações segmentares conhecidas, repetições tandem e homopolímeros.

Comparação da chamada de variantes usando Ash1 versus GRCh38

Uma das motivações para a criação de novos genomas de referência é que eles fornecem uma melhor estrutura para a análise de dados de seqüências humanas quando se procura por variantes genéticas ligadas a doenças. Por exemplo, um estudo dos judeus Ashkenazi que coletaram dados de caçadeira de genoma inteiro (WGS) deve usar um genoma de referência Ashkenazi em vez de GRCh38. Como a base genética é semelhante, menos variantes devem ser encontradas ao pesquisar contra Ash1.

Para testar essa expectativa, coletamos dados do WGS de um participante masculino no Projeto Genoma Pessoal, PGP17 (hu34D5B9). Este indivíduo é estimado em 66% de Ashkenazi de acordo com o banco de dados PGP, que foi a maior fração estimada disponível de indivíduos já sequenciados do PGP. Baixamos 983.220.918.100-bp leituras (aproximadamente 33x cobertura) e as alinhamos a Ash1 e GRCh38 usando Bowtie2 . Uma fração ligeiramente maior de leituras (3.901.270, 0,5%) alinhadas com Ash1 do que com GRCh38.

Examinamos então todas as variantes de um nucleotídeo (SNVs, veja os “Métodos”) entre o PGP17 e cada um dos dois genomas de referência. Para simplificar a análise, consideramos apenas os locais onde o PGP17 era homozigotos. Em nossas comparações com Ash1, primeiro identificamos todas as SNVs e depois examinamos os dados lidos no Ash1 original para determinar se, para cada uma dessas SNVs, o genoma Ash1 continha um alelo diferente que correspondia ao PGP17.

No total, o número de locais homozigotos no PGP17 que discordaram do Ash1 foi de 1.333.345, contra 1.700.364 quando comparamos locais homozigotos no PGP17 com o GRCh38 (arquivo adicional 1: Tabela S1). Em seguida, analisamos os dados de Ash1 subjacentes para os 1,33 milhões de sites da SNV que inicialmente não correspondiam, e descobrimos que para aproximadamente metade deles, o genoma Ash1 era heterozigoto, com um alelo correspondente ao PGP17. Se restringimos a SNV a locais onde PGP17 e Ash1 são ambos homozigotos (mais um número muito pequeno de locais onde Ash1 contém duas variantes que diferem do PGP17), isso reduziu ainda mais o número total de SNVs, para 707.756. Assim, encontramos pouco menos de 1 milhão de SNVs homozigotos quando usamos Ash1 como referência para o PGP17. Note que, em vez de alinhar com Ash1, pode-se alinhar as leituras ao GRCh38 e depois remover as variantes conhecidas específicas da população. Esse processo em duas etapas, embora mais complexo, pode produzir resultados semelhantes, exceto para regiões de Ash1 que simplesmente faltam no GRCh38.

Comparação contra variantes comuns de Ashkenazi

Para examinar até que ponto Ash1 contém variantes Ashkenazi conhecidas e comuns (em relação ao GRCh38), examinamos SNVs relatadas em alta freqüência em uma população Ashkenazi do Banco de Dados de Agregação de Genoma (gnomAD) . O GnomAD v3.0 contém chamadas SNV de dados do genoma inteiro de 1662 indivíduos Ashkenazi de leitura curta. Como algumas variantes foram chamadas apenas em um subconjunto desses indivíduos, consideramos apenas sites de variantes que foram relatados em um mínimo de 200 pessoas. Em seguida, coletamos SNVs de alelos principais, exigindo que a freqüência dos alelos fosse acima de 0,5 na população amostrada. Restringimos ainda nossa análise a substituições de base única. Isto nos deu 2.008.397 sites gnomAD SNV onde o alelo principal de Ashkenazi diferiu do GRCh38.

Fomos capazes de mapear precisamente 1.790.688 dos 2.008.397 sites gnomAD do GRCh38 para Ash1 (veja a seção “Métodos”). Comparamos então a base do GRCh38 com a base dos alelos principais de Ashkenazi em cada posição, e também examinamos os alelos alternativos em Ash1 nos locais onde é heterozigoto. Para locais onde os dados lidos mostraram que HG002 era heterozigoto e tinha tanto o alelo principal de Ashkenazi quanto o alelo GRCh38, nós substituímos a base de Ash1, se necessário, para garantir que ela correspondesse ao alelo principal. Após essas substituições, Ash1 continha o alelo principal de Ashkenazi em 88% (1.580.866) dos 1,79 milhões de locais. Nos demais locais, Ash1 ou correspondeu ao alelo GRCh38 porque HG002 é homozigoto para o alelo de referência (204.729 locais), ou continha um terceiro alelo que não corresponde nem ao GRCh38 nem ao alelo gnomAD (5093 locais). Essas modificações devem reduzir ainda mais o número de SNVs reportados ao alinhar um indivíduo Ashkenazi a Ash1.

Olhando que, como a freqüência do alelo principal na população de Ashkenazi gnomAD aumenta, a proporção de locais onde Ash1 coincidiu com o alelo principal também aumenta. Por exemplo, dos SNVs que têm uma frequência de alelos > 0,9 na população gnomAD Ashkenazi, mais de 98% estão representados em Ash1 (Tabela 3).

Annotação

Uma parte vital de qualquer genoma de referência é a anotação: a coleção de todos os genes e outras características encontradas no genoma. Para permitir que Ash1 funcione como um verdadeiro genoma de referência, mapeamos todos os genes conhecidos usados pela comunidade científica em seu sistema de coordenadas, usando os mesmos nomes de genes e identificadores. Vários bancos de dados de anotações diferentes foram criados para o GRCh38, e para Ash1, optamos por usar o banco de dados de genes humanos CHESS porque ele é abrangente, incluindo todos os genes codificadores de proteínas tanto no GENCODE quanto no RefSeq, os dois outros bancos de dados de genes amplamente utilizados, e porque ele retém os identificadores usados nesses catálogos. Os genes não-codificadores diferem entre as três bases de dados, mas o CHESS tem o maior número de loci e isoformas genéticas. Usamos CHESS versão 2.2, que tem 42.167 genes nos cromossomos primários (excluindo os andaimes alternativos GRCh38), dos quais 20.197 são codificadores de proteínas.

Mapping genes de um conjunto para outro é uma tarefa complexa, particularmente para genes que ocorrem em famílias de genes altamente similares, multi-cópia. A tarefa é ainda mais complexa quando os dois conjuntos representam indivíduos diferentes (ao invés de simplesmente conjuntos diferentes do mesmo indivíduo), devido à presença de diferenças de nucleotídeos únicos, inserções, deleções, rearranjos e variações genuínas no número de cópias entre os indivíduos. Precisávamos de um método que fosse robusto diante de todas essas diferenças potenciais.

Para resolver esse problema, usamos a ferramenta de mapeamento Liftoff recentemente desenvolvida, que em nossos experimentos foi a única ferramenta que poderia mapear quase todos os genes humanos de um indivíduo para outro. Liftoff leva todos os genes e transcrições de um genoma e os mapeia, cromossomo a cromossomo, para um genoma diferente. Para todos os genes que não conseguem mapear para o mesmo cromossoma, Liftoff tenta mapeá-los através dos cromossomas. Ao contrário de outras ferramentas, ele não depende de um mapa detalhado construído a partir de um alinhamento de genoma inteiro, mas, em vez disso, mapeia cada gene individualmente. Os genes são alinhados em nível de transcrição, incluindo introns, para que pseudogenes processados não sejam erroneamente identificados como genes.

Tentamos mapear todas as 310.901 transcrições de 42.167 loci de genes nos cromossomos primários do GRCh38 até Ash1. No total, mapeamos com sucesso 309.900 (99,7%) transcrições de 42.070 loci gênicos nos cromossomos principais (arquivo adicional 1: Tabela S2). Consideramos uma transcrição a ser mapeada com sucesso se a sequência de mRNA em Ash1 for pelo menos 50% maior do que a sequência de mRNA no GRCh38. Entretanto, a grande maioria das transcrições excede em muito este limite, com 99% das transcrições mapeadas em uma cobertura maior ou igual a 95% (Arquivo Adicional 2: Figura S2). A identidade de seqüência das transcrições mapeadas é similarmente alta, com 99% das transcrições mapeadas com uma identidade de seqüência maior ou igual a 94% (arquivo adicional 2: Figura S3).

Gens translocados

Desses genes com pelo menos uma isoforma mapeada com sucesso, 42.059 (99,7%) mapeados para os locais correspondentes no mesmo cromossomo em Ash1. Dos 108 genes que inicialmente falharam em mapear, 11 genes mapeados para um cromossomo diferente em 7 blocos distintos (mostrado na Tabela 4), sugerindo uma translocação entre os dois genomas. Curiosamente, 16 dos 22 locais envolvidos nas translocações estavam em regiões subteloméricas, o que ocorreu em 8 pares, onde ambos os locais estavam próximos a telômeros. Isto é consistente com estudos anteriores relatando que rearranjos envolvendo telômeros e subtelômeros podem ser uma forma comum de translocação em humanos .

Examinamos a translocação entre os cromossomas 15 e 20, que contém três dos genes da Tabela 4, observando mais de perto o alinhamento entre GRCh38 e Ash1. A translocação está no telômero de ambos os cromossomos, da posição 65.079.275 a 65.109.824 (30.549 bp) de Ash1 chr20 e 101.950.338 a 101.980.928 (30.590 bp) de GRCh39 chr15. Para confirmar a translocação, nós alinhamos um conjunto independente de PacBio muito longo lendo, todos do HG002, para o conjunto Ash1 v1.7 (veja a seção “Métodos”) e avaliamos a região ao redor do ponto de parada em chr20. Estes alinhamentos mostram uma cobertura profunda e consistente que se estende por muitos kilobases em ambos os lados do ponto de quebra, suportando a correcção do conjunto Ash1 (Fig. 1).

Snapshot mostrando alinhamentos de PacBio longo lido para o genoma Ash1, centrado na extremidade esquerda do local no cromossomo 20 (posição 65.079.275) onde ocorreu uma translocação entre os cromossomos 15 (GRCh38) e 20 (Ash1). A parte superior da figura mostra as coordenadas em chr20. Abaixo disso há um histograma de cobertura lida, e a leitura individual preenche a parte inferior da figura. Os indels nas leituras, mostrados como barras coloridas em cada leitura, são devidos à taxa de erro relativamente alta das leituras longas

Genes ausentes

Sixty-two genes falharam completamente em mapear de GRCh38 para Ash1, e outros 32 genes mapeados apenas parcialmente (abaixo do limiar de cobertura de 50%), como mostrado na Tabela 5. Todos os genes que falharam em mapear ou que mapearam parcialmente eram membros de famílias multi-gene, e em cada caso, havia pelo menos uma outra cópia do gene ausente presente em Ash1, a uma identidade média de 98,5%. Assim, não há nenhum caso de um gene que está presente no GRCh38 e que está completamente ausente em Ash1; os genes mostrados na Tabela 5 representam casos onde Ash1 tem menos membros de uma família multi-gene. Três genes adicionais (2 de codificação protéica, 1 lncRNA) mapeados para duas contigs não colocadas, que fornecerão um guia para colocar essas contigs em futuras liberações do conjunto Ash1.

Após o mapeamento dos genes em Ash1, extraímos as sequências codificadoras das transcrições que foram mapeadas completamente (cobertura igual a 100%), as alinhamos com as sequências codificadoras do GRCh38, e chamadas de variantes relativas ao GRCh38 (veja a seção “Métodos”). Dentro dos 35.513.365 bp dessas transcrições codificadoras de proteínas, encontramos 20.864 variantes e indels de nucleotídeos simples. Catorze mil e quatrocentos e noventa e nove dessas variantes caíram dentro das regiões “chamáveis” do GIAB para variantes de alta confiança, embora 3963 delas estivessem em regiões “difíceis” repetitivas do GIAB, para as quais os alinhamentos são frequentemente ambíguos. Das 10.536 variantes que não estavam nessas regiões difíceis, 10.456 (99,2%) estavam de acordo com o conjunto de variantes de alta confiança do GIAB. Nas regiões difíceis, 3804/3963 (96,0%) concordaram com o conjunto GIAB.

Em seguida, anotamos as alterações em aminoácidos causadas por variantes e mapeamento incompleto de todas as sequências de codificação de proteínas. Das 124.238 transcrições de codificação protéica de 20.197 genes, 92.600 (74,5%) têm seqüências proteicas 100% idênticas. Outros 26.566 (21,4%) têm pelo menos uma mudança de aminoácidos mas produzem proteínas com o mesmo comprimento, e 1561 (1,3%) têm mutações frame-preserving que inserem ou apagam um ou mais aminoácidos, deixando o resto da proteína inalterado. A Tabela 6 mostra estatísticas sobre todas as mudanças nas sequências de proteínas. Se uma proteína tivesse mais de 1 variante, contámo-la sob a variante mais conseqüente, ou seja, se uma proteína tivesse uma variante missense e um códon de parada prematura, incluímo-la no grupo “stop gained”.

De particular interesse são aquelas transcrições com variantes que perturbam significativamente a seqüência proteica e podem resultar em perda de função. Estas incluem transcrições afetadas por um framehift (2158), stop loss (58), stop gain (416), start loss (58), ou truncagem devido a um mapeamento incompleto (564). Essas isoformas perturbadas representam 885 loci gênicos; entretanto, 505 desses genes têm pelo menos 1 outra isoforma que não é afetada por uma variante perturbadora. Isso deixa 380 genes nos quais todas as isoformas têm pelo menos uma ruptura; a lista completa é fornecida no arquivo adicional 1: Tabela S1.