#+title: Systematic analysis of paralogous regions in 41,755 exomes uncovers clinically relevant variation
#+date: [2024-07-26 Fri 09:52]
#+filetags: :bib:chameleolyser:pseudogène:
#+identifier: 20240726T095249
#+reference: steyaert2023systematic
Identification en exome sur région homologues (3.5%)
1. SNV,indel
2. CNV et "ectopic gene conversion"
Principe : on prend les reads correspondant à des séquences paralogues sont alignés sur un génome de référence où sont masqués *tous les paralogues sauf un*
* Contexte
:PROPERTIES:
:CUSTOM_ID: h:0637a6fa-b137-434c-94d7-b08c4caba476
:END:
- région paralogues
- 17000 gènes codant ont une séquence fortement identiques à d’autres régions du génome (=paralogue)
- origin : duplication/rétro-transcription dans l’évolution
- divergence avec accumulation de mutation -> parfois perte de fonction -> paralogue non codant (pseudogene) ou codant mais différente fonction
- ectopic gene conversion
- copie d’une séquence homologue d’ADN (donneur) sur une autre (accepteur). Mécanisme = erreur pendant recombinaison homologue durant la méoise
- 1% des gènes impliqués maladies héritées
- problème :
1. gene conversion non indentifié en exome/genome
- les reads sur la région "acceptrice" sont aligné sur la région "donneuse" -> peu de reads sur la région acceptrice -> pas de SNVs appelé
- identifié comme délétion lors de l’appel de CNV
2. les régions identiques à 100% auront une qualité d’alignement de 0 (car multiple endroits possibles)
- Méthodes existantes
- CNV:
- utilisation de profondeur sur des nucléotides uniques : inconvénient : génome (pas un problème pour nous...) et ne fonctionne pas pour les gene conversion
- *rien pour les SNVs/indel sur paralogue identique*
- [cite:@ebbert2019systematic] : caractérisation mais pas de solution concrète d’après les auteurs.
* Algorithme
Sur les régions homologues (3.5%), il y a donc l'appel de SNVs (et indel) d'une part et les CNV/gene conversion
1. extraction des reads (3.5% exomes) affecté par l’homologie de séqueunce
2. ré-aligné sur génome de référence en masquant les régions homologues
3. appels de variant
Pour différencier délétion homozygote et ectopic gene conversion
1. analyse de la couverture dans l’alignement initial :
- les reads du site accepteur seront aligné sur le site donner -> pas de reads à cet pas de reads sur le site et 2x sur le donneur
Attention, des SNV/indel hétérozygotes n'ont pas de position bien définie (car paralogues) -> "variants with ambiguous positions"
* Validation
- Long-read Pacbio sur 20 échantillons
- sur 769 SNV/indel (non du à une gene conversion) -> 678 confirmé
- 8 gene conversion homozygotes -> toutes retrouvées
- 15 délétion homozygotes -> 13/15 confirmée
- 5 patients GIAB : NA12878 HG002-5
- attention: il faut le meme kit, donc ils n’ont regardé que les SNV/indel ne résultant pas d’un gene
- 118SNV/indel non résultat d’un gene conversion -> 98 concordent
** Données disponible
- pacbio maison: accès restreint, à demander
- https://ega-archive.org/datasets/EGAD00001009109 6 trio + 1 proband en HiF 30X
- https://ega-archive.org/datasets/EGAD00001011305 STRPC1 sur cohorte
- giab [[https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-023-42531-9/MediaObjects/41467_2023_42531_MOESM11_ESM.xlsx][Liste des variants]] disponible mais à revoir
* Résultats autres
** Comparaison avec autres appels de variants
- SNV et indel non dépendant gene conversion) vs GTK et Deepvariant
- régions avec une qualité d'alignement à 0 (logique)
- sensibilité de 0 (logique) vs 43% pour chameleolyzer
- région avec aligement unique : sensbilité 86.3% vs 88% chameleo
- délétion homozygotes,: pas assez d'échantillons en long-read pour calculer la spécificité
- comparaison à 2 autres outils d'appel de CNV (ExomeDepth et Confier)
- 3 deletions et 1 gene conversions uniques dont 1 deletion et la gene conversion retrouvé en long-read
- 12 del et 7 conversion retrouvé par exomedepth (marqué comme deletion)
- 201 deletion homozygotes uniquement retrouvé par exome-depth
- précision estimé à 32.5%
** 30-40% variants "ambigus" (= SNV/indel htz) dans région codante (estimation)
- 10% sont soit codant soit avec 2 alternatives dont 1 seule et codante
- 2 approches pour estimer la proportion de VAP codants
1. comparaison avec données long-read : 38% des 65 variants (les 10% ci-dessus) sont codants
2. utilisation du taux de synonymes, faux-sens et perte de fonction
- pipeline d'exome vs chameleolyse r: plus de faux-sens et synonyme dans les régions homologue
- hypthèse 1: pas de pression de sélection de variants synonyme -> la moitié sont dans des régions codantes
- hypèthes 2 : faux/sense/synome et synonyme/LOF semblable entre exome et homologue
** SNV/indel: 14 nouveaux diagnostics
~17k patients sans diagnostic
- filtre
-faux-sens et perte de fonction < 0.5% dans la cohorte
- sur gènes d'intéret (selon demande clinicien)
- + 1 variant synonyme de SMN1 (connu pour conduire à une protéine tronquée)
- 1 071 htz
- 7 dans /STRC/ avec délétion multipe-exon -> MPLA et PCR long-range
- tous les délétions sont dans /STRC/ (et pas dans son pseudogène)
- 4/7 SNV/indel également -> 4 diago
- 57 hmz (position non ambigue)
- 21 gene conversion :
- 10 diag, autres estimé non causal.
- seulement dans 3 gènes : /STRC/, /OTOPA/ /SMN1/
** délétion homozygotes: 11 diag
- Filtre : <= 0.5% dans la cohorte
- 2 sur OTOA, 9 sur SMN1 (toutes confirmées en MLPA)
** Différencier délétions homozygotes de gene conversion
- 58 délétion homozygotes (pipeline maison) confirmées en MLPA
- Chameleolyser : perte homozygote alèlle /STRC/ mais
- 37 = délétion homozygote
- 22 = gene conversion homozygote STRCP1 -> STRC
Pathogénicité identique donc ne sert à rien sur cet exemple
Les auteurs argument que les gene conversion sont bénignes donc important de faire la différence
* Méthodes
** Identification des régions paralogues
Fusion de 2 régions définissant des régions paralogues
I. gènes codant ayant des pseudo-gènes
1. pour chaque pseudo-gènes annoté par Encode, sélection du gène codant correspondant (en utilisant l'identifiant HGNC /NB: les pseudogene se finissent par P ou P$N avec n entier/)
2. alignement entre le gène codant et ses pseuodgene (MAFFT 7.047)
3. filtre : seulement une homologie >= 90% avec le pseudogène et avec exactement 1 pseudegen paralogue
II. régions avec une faible qualité d'alignement (MAPQ < 10)
1. on enlève les régions mal couvertes (profondeur < 10)
2. fusion de ces régions qui sont fragmentées
1. 3 versions (ou listes) : distance < 250bp, 500bp, 5% de la région (/slopping distance/))
- pour 5% de la liste, on supprime les régions < 50bp
2. pour chaque liste, définition d'un ensemble de région paralogues tels que
- toutes les régions d'un ensemble ont un score d'alignement >= 0.9
- aucune ne s'aligne avec un score >= 0.9 en dehors de l'ensemble
Pour chaque liste, les ensembles sans recouvrement exonique et avec > 5 éléments
3. la liste à 500 est fusionnée avec les 2 autres s'il n'y a pas de recouvrement avec les régions déjà fusionnées
(on fait d'abord 250bp puis 5%)
** Génération du génome masqué
Pour un ensemble de régions paralogues, on veut toutes les masquer, sauf 1 -> on ne retient que celles qui ont le moins de recouvrement avec une séquence codante
** Alignment et appel de variants
1. extraction des reads contenant les régions définies ci-dessus
2. conversion en FASTQ
3. suppression des doublons (picard)
4. appel de variant avec lofreq (car la VAF ne sera pas à 50%). Paramètres : no-default-filter, use-orphan, no-baq, no-mq, sig = 1
** Identification des délétions et des conversions
1. Détections des sous-régions
- pour I., intersection avec les exons codant pour les protéines (Gencode) +200bp -> donc exons seuls a priori ?
- pour II., inutile car ces régions sont petites mais on ne garde que les ensembles qui n'ont que 2 paralogues
2. Nombre de reads par sous-régions
3. kernel density estimation (KDE): on utilise la cohorte pour «lisser» les données avec un kernel (pour chaque sous-région probablement). On regarde la distribution sur tous les échantillons. Si un échantillon mal couvert est dans un groupe d'échantillons (pic), ce n'est probablement pas une délétion. Mais s'il est isolé, oui
- NB : je pense qu'on doit travailler avec des données de la forme y=nombre d'échantillons et x=couverture
4. Délétion complète des paralogues dans un ensemble de paralogue : on utilise une KDE mais en supprimant les pics > 10 readt ou > 10% des échantillons (pourquoi)
- NB: x = nb de reads sur cet intervalle donc a priori y = compte
5. Gene conversion ou délétion d'une région sur les 2
Également KDE mais avec un ratio
R = (nb reads s'alignant de manière unique sur X)/(ceux s'alignant sur X ou Y)
Filtre : on enlève les échantillons avec < 30 reads, région où nb moyen de reads sur X ou Y < 60, et pics > 10% des échantillons et /altérations ne recouvrant pas un gène codant pour une protéine/
Donc
- délétion homozygote de X si <= 10 reads sur X et R <= 0.05
- gene conversion si percentile maison <= seuil et si nombre de reads s'alignant sur Y > 3/4 (normalisé par le nombre de reads sur tous les paralogues)
NB : l'idée est de forcer un seuil extrême sur le site donneur s'il y a une couverture très importante sur le site accepteur
- sinon délétion de X
Enfin les délétions et conversion sont fusionnées : par gène puis pour les gènes proches
** Traitement des /shorts variants/
1. Appel de variant sur l'alignement initial et celui masqué
2. on ne garde que ceux
- profondeur >= 60 après masquage
- VAF >= 0.15
- variant absent du VCF initial
3. on ne sait pas sur quelle région sont les variants (sauf homozygotes) -> on calcule toutes les possibilités (x4 mais ce n'est pas intuitif)
4. Autres filtres : sont exclus
- variants > 10% des échantillons
- bonne qualité sur 1 échantillon mais pas sur les autres
- homopolymers
- sous régions avec 5 nucléotides uniques sur <= 10bp (diminue les faux positifs)