#+title:      Projet pseudogene
#+date:       [2024-07-16 mar. 10:49]
#+filetags:   :projet:pseudogene:auragen:
#+identifier: 20240716T104901


Ou régions répétées

* Contexte
- on couvre 98% (régions codantes) mais les 2% sont vraiment un "trou"
- sera un des défauts principaux d'auragen après curagen v2
- proposer une annotation dans la limite de ce qu'on peut faire sur ce genre de région
- liste de pseudogene connue +++

* Objectif
:PROPERTIES:
:CUSTOM_ID: h:5e2833d6-77eb-4b94-b82f-0b8f36472af7
:END:
- Porter [[https://github.com/Genome-Bioinformatics-RadboudUMC/Chameleolyser][Chameleolyser]]  [cite:@steyaert2023systematic] en
  - hg19 -> hg38  [[https://github.com/Genome-Bioinformatics-RadboudUMC/ChameleolyserBEDs][fichiers BEDs]]
  - +/- recoder en python
* Notes
[cite:@steyaert2023systematic] définit les régions homologues en fusionnant
- les gènes codant ayant des pseudogènes
- les régions bien couvertes avec une qualité d'alignement = 0 en utilisant 250 exomes

Données 1000 genomes 30x https://www.internationalgenome.org/data-portal/data-collection/30x-grch38 = 3200 échantillons

** Problèmes avec [cite:@steyaert2023systematic] :
- en hg19 et ne donne pas le code pour les générer
- les filtres après l'appel de variant utilisent la cohorte
- les régions homologues sont définies avec la cohorte (régions bien couvertes mais avec une mauvaise qualité d'alignement)
- la détection de délétions ou gène conversion utilise toute la cohorte
- exome seul (en génome on aura beaucoup de régions intergéniques...)

** Essai liftover BED
#+begin_src sh
git clone https://github.com/Genome-Bioinformatics-RadboudUMC/ChameleolyserBEDs ~/code/ChameleolyserBEDs
#+end_src

#+RESULTS:
On utilise l'outil en ligne
https://genome.ucsc.edu/cgi-bin/hgLiftOver
On prend tous les gènes (pas que les OMIM) avec la convention «chr»

#+begin_src sh

wget https://hgdownload.cse.ucsc.edu/goldenpath/hg19/liftOver/hg19ToHg38.over.chain.gz
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/liftOver
gunzip hg19ToHg38.over.chain.gz
chmod +x liftOver
#+end_src

On compare les multiples
#+begin_src julia :dir ~/code/pseudogenes
using CSV, DataFramesMeta, DataFrames

for e in ["forextraction", "formasking", "forvarcall", "homologousexons"]
    f = "All.$e.noalt.chr"
    run(`./liftOver $f.bed hg19ToHg38.over.chain $f.hg38.bed $f.hg38.err`)
    run(`./liftOver -multiple $f.bed hg19ToHg38.over.chain $f.hg38.multiple.bed $f.hg38.multiple.err`)
end
#+end_src

Sans l'option -multiple (en lifte pas les régions qui s'alignent à plusieurs endroits), moins de résultats

Si on essaie de lifter par exemple
chr1	13606679	13609013 [[https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr1%3A13606679%2D13609013&hgsid=2264446990_5w3JSGMGACw9wraSAvbEKw2GTryr][UCSC]]

est coupé en 2 sur ALT (PRAMEF8) car ce gène n'est que dans le ALT
chr1:13281531-13282617 (46.6% of bases, 46.6% of span)
chr1:12920254-12921000 (32.0% of bases, 32.0% of span)
https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr1%3A13281531%2D13282617&hgsid=2328173942_Ut5M2Px0e241mNfbnnOQMKRQ1MSk
https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr1%3A12920254%2D12921000&hgsid=2328173952_fxB39QAXCwFtUBm8czqGaUI8HaIh

Plus de problèmes en T2T
*** Problème 1 : quid des régions qui s'alignent à plusieurs endroits
Pseudogène ?

#+begin_src julia
function ucsc(pos)
    `google-chrome https://genome-euro.ucsc.edu/cgi-bin/hgTracks\?db=hg19\&lastVirtModeType=default\&lastVirtModeExtraState=\&virtModeType=default\&virtMode=0\&nonVirtPosition=\&position=$pos`
end
function readErr(f)
    extr = CSV.read(f, DataFrame; comment="#",header=false)
    @rtransform extr :pos = :Column1 * ":" * string(:Column2) * "-" * string(:Column3)
end
extr = readErr("All.forextraction.noalt.chr.hg38.err")
extrM = readErr("All.forextraction.noalt.chr.hg38.multiple.err")
# extrM = CSV.read("All.forextraction.noalt.chr.hg38.multiple.err", DataFrame; comment="#",header=false)
# antijoin(extr, extrM; on=[:Column1, :Column2, :Column3])
d = antijoin(extr, extrM; on=:pos)
for x in eachrow(d)
    run(ucsc(x.pos))
end
#+end_src
On vérifie dans le fichier de pseudogene : non présent
Faut-il les accepter ??
*** Problème 2 : quid des régions qui ont échoué au liftover

* Pistes
- liftover directement
- se limiter au pseudogènes ENCODE ?
* Données fournies
Supplementary Data 23 : régions où est fait l'appel de variant
Supplementary Data 24 : régions d'où sont extraites les reads
Supplementary Data 25 : régions à masquer dans le génome
Supplementary Data 26 : sous-régions pour délétion et gene conversion
* [[denote:20240716T104934][Bibliographie pseudogène]]
* Annotation pseudogene
** Genecode
[[https://www.gencodegenes.org/human/][Genecode]] propose une annotation par HAVANA des pseudogene dans l'annotation principale (manuel donc idéal)
En bonus, les pseudogene prédits par 2 pipeline mais pas dans havannah
** Pseudogene.org
Contient en plus l'activité
* Validation
** GIAB
[cite:@steyaert2023systematic]

** Patient test Auragen
*** MR-24000544 IKBKG chrX g.154560563_154560564del
*** MR-2202491:  chr1:155235252A>G  sur /GBA/ : probablement recombinaison avec pseudogène /GBAP1.
séquence très proche entre les 2 gènes donc alignement difficile.
Plusieurs reads ont donc une mauvaise qualité -> VAF sous-estimée
*** MR-1900206 ?
variant sur IKBKG cas index (D. Sanlaville) mais on ne peut conclure sur le status en mosaïque : de novo + pseudogène (IKBKGP), garçon klinefelter
***  MR-2300984 ? double délétion STRC
Pseudogène STRCP1
*** MR-2303627 TUBB2B confirmé à Lyon par technique complémentaire
Lyon (Louis Januel)
chr6:g.3225154G>A