AWWCOU4XZXRGWAREDCX5KTC2OLMND45IZNMRYLHLVZCYZS4DUFEAC
OPTRA6QLBQIMDFPEJGFDVSMFRC4CLE6JIMSCA6GKKFIA33W5JCJQC
KNU6ADO757Q42ZLJUEVR6Z4FFSAWZKH6HZNSDNFMVONLEYBQIRIAC
4UOH23PPEC4TH5GALRKZOOGOSLV4FJB6R4GA2N3TCEDCOIZTMJSAC
WL5R6AMW2SA27H5YOFYKRSTM344ZTIDVGBA55TS7QOU5LAKLV7VQC
RJLKEE5RK3UVUB3NVU6BPFMWVUWUPQTOVTHHFG5KL53VHXHOTI5QC
RHWQQAAHNHFO3FLCGVB3SIDKNOUFJGZTDNN57IQVBMXXCWX74MKAC
RR76MPYLRKMP6MMEAVZJNWEACRDPIFP6H7AKV3MROPLWRJIZ3WRQC
3LLSOLDOJ5OSKQN2DKYHYGWTBUJAC5KPRX2SAVNARRUYCPRM44RQC
RUR2HQCDDUOIP7GD2GZV2UPCUF5QS6COVIVDKVIGZF22TVRFBKPAC
FXA3ZBV64FML7W47IPHTAJFJHN3J3XHVHFVNYED47XFSBIGMBKRQC
JH5INW2G6JKUQ3QD2VMSAFEIGUUBCBD5RJGSQOVEVF3I7BKL7SKQC
JJTKCDYDZIY34JHW7BKMYXGNZBFWETAUGXAI3KH6MSTAN7KQ2MWQC
7JUOBG7KBW66RI6SUDOPRLXB7PPDBY7HS4LJDA6ON4Z7ZQ7DBAFQC
KNZKRW5WZVVQL73LUG3GGFR5QHJJWI4LMOLHJUXUU4PHNBCBLJLAC
4FZ6627CHEHJTPLHE7MF6ZVFJKSDUAKOBGMJZ6QBR2HARWWCXOFQC
l#downloadAlt
- alt = séquences alternatives (utilisables)
- fix = patch (correction ou amélioration)
- random = séquence connue sur un chromosome mais non encore utilisée
** Pipelines prêt-à-l’emploi nextflow
Problème : nécessite singularity ou docker (ou conda)
Potentiellement utilisable avec nix...
** Validation : Quelles données de référence ?
Discussion avec Alexis
- Platinum genomes = génome seul
*** [[https://github.com/genome-in-a-bottle/giab_data_indexes][Genome in a bottle]]
- NA12878 :
- Illumina HiSeq Exome : fastq + capture en hg37
- Illumina TruSeq Exome : bam, pas de capture
- Exomes en hg37 https://zenodo.org/record/3597727 avec capture
- HiSeq2000
- NextSeq 500
- HiSeq 2500
- HG002,3,4
- Illumina Whole Exome : bam. le kit de capture est "Agilent SureSelect Human All Exon V5 kit" selon [[https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/README.txt][README]]. On il faut les régions [[https://kb.10xgenomics.com/hc/en-us/articles/115004150923-Where-can-I-find-the-Agilent-Target-BED-files-][selon ce site]]
Un autre fichier est disponible (capture ???)
https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/wex_Agilent_SureSelect_v05_b37.baits.slop50.merged.list
"target region" +/- 50bp
testé sur chr311780-312086 : ok
Autres technologies non adaptées au pipeline (vu avec Alexis)
*** [[https://www.illumina.com/platinumgenomes.html][Platinum genome
]] Que du génome « sequenced to 50x depth on a HiSeq 2000 system”
Genome possible
*** 1000 genomes
- intersection des capture + CCDS [[id:b77e64fa-06a8-4ffa-8b5b-ab3fda684b61][Données brutes exome 1000 Genomes (fastq + capture)]]
- Broad instute : SureSelect human all exon v2 target capture kit : non disponible sur le site d'agilent (V6 ou plus)
*** Zone de capture
GIAB fourni le .bed pour l'exome . INfo : https://support.illumina.com/sequencing/sequencing_kits/nextera-rapid-capture-exome-kit/downloads.html
*** Valider la méthode
- 1000 genomes + SureSelect human all exon v2 target capture kit : non disponible sur le site d'agilent (V6 ou plus)
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-019-2928-9
- GIAB + liftover du fichire de capture en hg38
Ce qui est aussi fait par
https://bcbio-nextgen.readthedocs.io/en/stable/contents/germline_variants.html
Mais avec UCSC liftover
** Centogène
https://www.twistbioscience.com/node/23906
Bed non fourni pour exactement cette capture
On prend https://www.twistbioscience.com/resources/data-files/twist-alliance-vcgs-exome-401mb-bed-files
qui content la majeure partie
* Réunion
** <2023-08-10 Thu> Alexis
Ok pour bloquer le développment d'ici mardi prochain
Dév:
- pipeline jusque VEP en T2T + GRCh38
- ok pour valider spip T2T sur quelques variant => à intégrer au pipeline
- annotation :
- ok pour mobidetails hg38
- +OMIM T2T+ non
- +franklin hg38+ non pour le moment
- métriques (fastq a minima) + rapport multiqc
- optionnel
- reformater la sortie
- on abandonne
- XAMScissors ave indel
- parallélisation haplotype caller
- spliceai à la vollée
- pangolin
Test
- GIAB:
- hg38: ok pour refaire les tests NA12878 avec données cento, sinon ok pour "c'est difficile" sur les 3 fichiers de capture
- T2T: ok pour faire des tests rapides mais probablement pas assez de temps !
- patient de synthèse : variant cento confirém par sanger seuls
Résultats
- ok pour scale up bwa mem et haplotyecaller
Manuscrit
- validation de méthode : laisser tomber la version actuelle et faire comme strasbourg (cf ngs diag) dans la présentatino
- a envoyé le powerponit avec les références des différsences articles
- ok pour robo4 si résultat
- architecture cible = VM : 78 coeurs 54Go RAUM et 1To espace disque
Passage en production : ok pour présentation rapide du code
* Nixpkgs :nix:
** DONE GATK
CLOSED: [2023-05-06 Sat 08:51]
*** DONE [[https://github.com/NixOS/nixpkgs/pull/185819][Binaire]]
CLOSED: [2022-09-10 Sat 23:53] SCHEDULED: <2022-08-10 Wed>
/Entered on/ [2022-08-09 Tue 10:57]
PR submitted
*** KILL Corriger code pour utiliser source
CLOSED: [2022-09-11 Sun 22:05]
*** DONE Corriger PATH pour include java et python
CLOSED: [2022-10-11 Tue 11:46]
https://github.com/NixOS/nixpkgs/pull/191548
Review <2022-10-10 Mon> , corrigé dans la journée
*** DONE Update 4.3.0.0
CLOSED: [2023-04-13 Thu 09:01]
** HOLD Nextflow
*** KILL version script seule
CLOSED: [2023-04-01 Sat 18:29]
Fix pour SGE et nextflow
https://github.com/NixOS/nixpkgs/issues/192396
*** KILL Version avec gradle
CLOSED: [2022-10-09 Sun 22:51]
*** HOLD [[https://github.com/NixOS/nixpkgs/issues/192396][Bug report Version 22.10.6]]
**** Notes
Erreur :
ERROR: Cannot download nextflow required file -- make sure you can connect to the internet
Alternatively you can try to download this file:
https://www.nextflow.io/releases/v22.10.6/nextflow-22.10.6-all.jar
and save it as:
.//nix/store/md2b1ah4d7ivj82k8xxap30dmdci00pa-nextflow-22.10.6/bin/.nextflow-wrapped
Dans la mise à jour, il y a la création d'un environnement virtuel qui casse l'exécution de nextflow (besoin de télécharger)
Fix = désactiver
**** KILL Patch NXF_OFFLINE=true
CLOSED: [2023-07-02 Sun 11:02] SCHEDULED: <2023-06-11 Sun>
** TODO Multiqc
SCHEDULED: <2023-08-13 Sun>
** TODO Mutalyzer
SCHEDULED: <2023-08-13 Sun>
** TODO SPIP T2T
*** DONE PR upstream
CLOSED: [2023-08-12 Sat 18:23] SCHEDULED: <2023-08-12 Sat 18:00>
*** DONE Mail R. Lemann
CLOSED: [2023-08-12 Sat 18:23] SCHEDULED: <2023-08-12 Sat 18:00>
*** TODO Mise à jour packages nix
** TODO VEP :vep:
*** DONE [[https://github.com/NixOS/nixpkgs/pull/185691][BioPerl]]
SCHEDULED: <2022-08-10 Wed>
/Entered on/ [2022-08-09 Tue 10:57]
PR submitted
*** TODO BioDBBBigFile
:PROPERTIES:
:ORDERED: t
:END:
/Entered on/ [2022-08-10 Wed 14:28]
On utilise la dernière version de kent, donc plus de problème.
PRête à être mergé. Rebase faite<2023-07-02 Sun>
**** DONE Version de kent déjà packagée : forcer version 335
CLOSED: [2023-07-02 Sun 11:20]
***** KILL [[https://github.com/NixOS/nixpkgs/pull/206991][Restore building kent 404]]
CLOSED: [2023-05-06 Sat 17:40]
Review faite <2023-03-26 Sun> , atteinte merge]
Relancé <2023-05-06 Sat>
Kent 446 n'a pas ce problème donc PR inutile
***** DONE [[https://github.com/NixOS/nixpkgs/pull/223411][Ajouter les header to package]] (inc folder)
CLOSED: [2023-05-08 Mon 10:18] SCHEDULED: <2023-05-07 Sun>
Review à faire
https://github.com/NixOS/nixpkgs/pull/223411
Corrigé et plus besoin de la PR précédente
***** KILL [[https://github.com/NixOS/nixpkgs/pull/186462][BioDBBBigFile]] avec ces 2 changements
CLOSED: [2023-07-02 Sun 11:20]
**** KILL Version de kent déjà packagée : 404
CLOSED: [2023-03-27 Mon 16:43]
Compile mais les tests de passent pas
**** DONE Modifier selon PR https://github.com/NixOS/nixpkgs/pull/186462
CLOSED: [2023-07-30 Sun 22:01] SCHEDULED: <2023-07-30 Sun 20:00>
:LOGBOOK:
CLOCK: [2023-07-30 Sun 19:13]--[2023-07-30 Sun 20:50] => 1:37
:END:
Modification nécessaire pour kent :
- plus de patch
- suppression d'une boucle dans postPatch
On supprime aussi NIX_BUILD_TOP
*** DONE [[https://github.com/NixOS/nixpkgs/pull/186459][BioDBHTS]]
CLOSED: [2023-05-06 Sat 08:49] SCHEDULED: <2023-04-15 Sat>
/Entered on/ [2022-08-10 Wed 14:28]
Correction pour review faites <2022-10-10 Mon>
*** DONE [[https://github.com/NixOS/nixpkgs/pull/186464][BioExtAlign]]
CLOSED: [2022-10-22 Sat 12:43] SCHEDULED: <2022-08-10 Wed>
/Entered on/ [2022-08-10 Wed 14:28]
Review <2022-10-10 Mon>, correction dans la journée.
Correction 2e passe, attente
Impossible de faire marcher les tests Car il ne trouve pas le module Bio::Tools::Align, qui est dans un dossier ailleurs dans le dépôt. Même en compilant tout le dépôt, cela ne fonctionne pas... On skip les tests.
*** TODO VEP
** WAIT [[https://github.com/N
l#downloadAlt
- alt = séquences alternatives (utilisables)
- fix = patch (correction ou amélioration)
- random = séquence connue sur un chromosome mais non encore utilisée
** Pipelines prêt-à-l’emploi nextflow
Problème : nécessite singularity ou docker (ou conda)
Potentiellement utilisable avec nix...
** Validation : Quelles données de référence ?
Discussion avec Alexis
- Platinum genomes = génome seul
*** [[https://github.com/genome-in-a-bottle/giab_data_indexes][Genome in a bottle]]
- NA12878 :
- Illumina HiSeq Exome : fastq + capture en hg37
- Illumina TruSeq Exome : bam, pas de capture
- Exomes en hg37 https://zenodo.org/record/3597727 avec capture
- HiSeq2000
- NextSeq 500
- HiSeq 2500
- HG002,3,4
- Illumina Whole Exome : bam. le kit de capture est "Agilent SureSelect Human All Exon V5 kit" selon [[https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/README.txt][README]]. On il faut les régions [[https://kb.10xgenomics.com/hc/en-us/articles/115004150923-Where-can-I-find-the-Agilent-Target-BED-files-][selon ce site]]
Un autre fichier est disponible (capture ???)
https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/wex_Agilent_SureSelect_v05_b37.baits.slop50.merged.list
"target region" +/- 50bp
testé sur chr311780-312086 : ok
Autres technologies non adaptées au pipeline (vu avec Alexis)
*** [[https://www.illumina.com/platinumgenomes.html][Platinum genome
]] Que du génome « sequenced to 50x depth on a HiSeq 2000 system”
Genome possible
*** 1000 genomes
- intersection des capture + CCDS [[id:b77e64fa-06a8-4ffa-8b5b-ab3fda684b61][Données brutes exome 1000 Genomes (fastq + capture)]]
- Broad instute : SureSelect human all exon v2 target capture kit : non disponible sur le site d'agilent (V6 ou plus)
*** Zone de capture
GIAB fourni le .bed pour l'exome . INfo : https://support.illumina.com/sequencing/sequencing_kits/nextera-rapid-capture-exome-kit/downloads.html
*** Valider la méthode
- 1000 genomes + SureSelect human all exon v2 target capture kit : non disponible sur le site d'agilent (V6 ou plus)
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-019-2928-9
- GIAB + liftover du fichire de capture en hg38
Ce qui est aussi fait par
https://bcbio-nextgen.readthedocs.io/en/stable/contents/germline_variants.html
Mais avec UCSC liftover
** Centogène
https://www.twistbioscience.com/node/23906
Bed non fourni pour exactement cette capture
On prend https://www.twistbioscience.com/resources/data-files/twist-alliance-vcgs-exome-401mb-bed-files
qui content la majeure partie
* Réunion
** <2023-08-10 Thu> Alexis
Ok pour bloquer le développment d'ici mardi prochain
Dév:
- pipeline jusque VEP en T2T + GRCh38
- ok pour valider spip T2T sur quelques variant => à intégrer au pipeline
- annotation :
- ok pour mobidetails hg38
- +OMIM T2T+ non
- +franklin hg38+ non pour le moment
- métriques (fastq a minima) + rapport multiqc
- optionnel
- reformater la sortie
- on abandonne
- XAMScissors ave indel
- parallélisation haplotype caller
- spliceai à la vollée
- pangolin
Test
- GIAB:
- hg38: ok pour refaire les tests NA12878 avec données cento, sinon ok pour "c'est difficile" sur les 3 fichiers de capture
- T2T: ok pour faire des tests rapides mais probablement pas assez de temps !
- patient de synthèse : variant cento confirém par sanger seuls
Résultats
- ok pour scale up bwa mem et haplotyecaller
Manuscrit
- validation de méthode : laisser tomber la version actuelle et faire comme strasbourg (cf ngs diag) dans la présentatino
- a envoyé le powerponit avec les références des différsences articles
- ok pour robo4 si résultat
- architecture cible = VM : 78 coeurs 54Go RAUM et 1To espace disque
Passage en production : ok pour présentation rapide du code
* Nixpkgs :nix:
** DONE GATK
CLOSED: [2023-05-06 Sat 08:51]
*** DONE [[https://github.com/NixOS/nixpkgs/pull/185819][Binaire]]
CLOSED: [2022-09-10 Sat 23:53] SCHEDULED: <2022-08-10 Wed>
/Entered on/ [2022-08-09 Tue 10:57]
PR submitted
*** KILL Corriger code pour utiliser source
CLOSED: [2022-09-11 Sun 22:05]
*** DONE Corriger PATH pour include java et python
CLOSED: [2022-10-11 Tue 11:46]
https://github.com/NixOS/nixpkgs/pull/191548
Review <2022-10-10 Mon> , corrigé dans la journée
*** DONE Update 4.3.0.0
CLOSED: [2023-04-13 Thu 09:01]
** HOLD Nextflow
*** KILL version script seule
CLOSED: [2023-04-01 Sat 18:29]
Fix pour SGE et nextflow
https://github.com/NixOS/nixpkgs/issues/192396
*** KILL Version avec gradle
CLOSED: [2022-10-09 Sun 22:51]
*** HOLD [[https://github.com/NixOS/nixpkgs/issues/192396][Bug report Version 22.10.6]]
**** Notes
Erreur :
ERROR: Cannot download nextflow required file -- make sure you can connect to the internet
Alternatively you can try to download this file:
https://www.nextflow.io/releases/v22.10.6/nextflow-22.10.6-all.jar
and save it as:
.//nix/store/md2b1ah4d7ivj82k8xxap30dmdci00pa-nextflow-22.10.6/bin/.nextflow-wrapped
Dans la mise à jour, il y a la création d'un environnement virtuel qui casse l'exécution de nextflow (besoin de télécharger)
Fix = désactiver
**** KILL Patch NXF_OFFLINE=true
CLOSED: [2023-07-02 Sun 11:02] SCHEDULED: <2023-06-11 Sun>
** WAIT [[https://github.com/NixOS/nixpkgs/pull/249329][Multiqc]]
HG002,sanger-chr20,data/HG002-sanger-inserted-chr20_1.fq.gz,data/HG002-sanger-inserted-chr20_2.fq.gz
** TODO Mutalyzer
SCHEDULED: <2023-08-13 Sun>
** TODO SPIP T2T
*** DONE PR upstream
CLOSED: [2023-08-12 Sat 18:23] SCHEDULED: <2023-08-12 Sat 18:00>
*** DONE Mail R. Lemann
CLOSED: [2023-08-12 Sat 18:23] SCHEDULED: <2023-08-12 Sat 18:00>
*** TODO Mise à jour packages nix
** TODO VEP :vep:
*** DONE [[https://github.com/NixOS/nixpkgs/pull/185691][BioPerl]]
SCHEDULED: <2022-08-10 Wed>
/Entered on/ [2022-08-09 Tue 10:57]
PR submitted
*** TODO BioDBBBigFile
:PROPERTIES:
:ORDERED: t
:END:
/Entered on/ [2022-08-10 Wed 14:28]
On utilise la dernière version de kent, donc plus de problème.
PRête à être mergé. Rebase faite<2023-07-02 Sun>
**** DONE Version de kent déjà packagée : forcer version 335
CLOSED: [2023-07-02 Sun 11:20]
***** KILL [[https://github.com/NixOS/nixpkgs/pull/206991][Restore building kent 404]]
CLOSED: [2023-05-06 Sat 17:40]
Review faite <2023-03-26 Sun> , atteinte merge]
Relancé <2023-05-06 Sat>
Kent 446 n'a pas ce problème donc PR inutile
***** DONE [[https://github.com/NixOS/nixpkgs/pull/223411][Ajouter les header to package]] (inc folder)
CLOSED: [2023-05-08 Mon 10:18] SCHEDULED: <2023-05-07 Sun>
Review à faire
https://github.com/NixOS/nixpkgs/pull/223411
Corrigé et plus besoin de la PR précédente
***** KILL [[https://github.com/NixOS/nixpkgs/pull/186462][BioDBBBigFile]] avec ces 2 changements
CLOSED: [2023-07-02 Sun 11:20]
**** KILL Version de kent déjà packagée : 404
CLOSED: [2023-03-27 Mon 16:43]
Compile mais les tests de passent pas
**** DONE Modifier selon PR https://github.com/NixOS/nixpkgs/pull/186462
CLOSED: [2023-07-30 Sun 22:01] SCHEDULED: <2023-07-30 Sun 20:00>
:LOGBOOK:
CLOCK: [2023-07-30 Sun 19:13]--[2023-07-30 Sun 20:50] => 1:37
:END:
Modification nécessaire pour kent :
- plus de patch
- suppression d'une boucle dans postPatch
On supprime aussi NIX_BUILD_TOP
*** DONE [[https://github.com/NixOS/nixpkgs/pull/186459][BioDBHTS]]
CLOSED: [2023-05-06 Sat 08:49] SCHEDULED: <2023-04-15 Sat>
/Entered on/ [2022-08-10 Wed 14:28]
Correction pour review faites <2022-10-10 Mon>
*** DONE [[https://github.com/NixOS/nixpkgs/pull/186464][BioExtAlign]]
CLOSED: [2022-10-22 Sat 12:43] SCHEDULED: <2022-08-10 Wed>
/Entered on/ [2022-08-10 Wed 14:28]
Review <2022-10-10 Mon>, correction dans la journée.
Correction 2e passe, attente
Impossible de faire marcher les tests Car il ne trouve pas le module Bio::Tools::Align, qui est dans un dossier ailleurs dans le dépôt. Même en compilant tout le dépôt, cela ne fonctionne pas... On skip les tests.
*** TODO VEP
** WAIT [[https://github.com/N
1.860963 | 2.247273 |
| INDEL | PASS | 549 | 489 | 60 | 899 | 64 | 340 | 8 | 17 | 0.890710 | 0.885510 | 0.378198 | 0.888102 | NaN | NaN | 1.860963 | 2.247273 |
| SNP | ALL | 21973 | 21462 | 511 | 26285 | 563 | 4263 | 68 | 16 | 0.976744 | 0.974435 | 0.162184 | 0.975588 | 3.00711 | 2.784686 | 1.591810 | 1.816145 |
| SNP | PASS | 21973 | 21462 | 511 | 26285 | 563 | 4263 | 68 | 16 | 0.976744 | 0.974435 | 0.162184 | 0.975588 | 3.00711 | 2.784686 | 1.591810 | 1.816145 |
******** Résumé
T2T
| Type | TRUTH.TOTAL | TRUTH.TP | TRUTH.FN | QUERY.TOTAL | QUERY.FP | QUERY.UNK | FP.gt | FP.al | METRIC.Recall | METRIC.Precision |
| INDEL | 413 | 246 | 167 | 751 | 289 | 215 | 2 | 93 | 0.595642 | 0.460821 |
| SNP | 11236 | 10985 | 251 | 23597 | 9771 | 2841 | 26 | 58 | 0.977661 | 0.529245 |
Hg38
| Type | TRUTH.TOTAL | TRUTH.TP | TRUTH.FN | QUERY.TOTAL | QUERY.FP | QUERY.UNK | FP.gt | FP.al | METRIC.Recall | METRIC.Precision |
| INDEL | 549 | 489 | 60 | 899 | 64 | 340 | 8 | 17 | 0.890710 | 0.885510 |
| SNP | 21973 | 21462 | 511 | 26285 | 563 | 4263 | 68 | 16 | 0.976744 | 0.974435 |
****** DONE Interesection des bed: similaire
CLOSED: [2023-07-04 Tue 23:11]
HG38
#+begin_src sh
bedtools intersect -a capture/Agilent_SureSelect_All_Exons_v7_hg38_Regions.bed -b /Work/Groups/bisonex/data/giab/GRCh38/HG001_GRCh38_1_22_v4.2.1_benchmark.bed | wc -l
#+end_src
204280
T2T
#+begin_src sh
bedtools intersect -a /Work/Groups/bisonex/data/giab/T2T/Agilent_SureSelect_All_Exons_v7_hg38_Regions_hg38_T2T.bed -b /Work/Groups/bisonex/data/giab/T2T/HG001_GRCh38_1_22_v4.2.1_benchmark_hg38_T2T.bed | wc -l
#+end_src
204021
****** DONE Vérifier la ligne de commande
CLOSED: [2023-07-04 Tue 23:38]
#+begin_src sh
hap.py \
HG001_GRCh38_
1_22_v4_lifted_merged.vcf.gz \
HG001-SRX11061486_SRR14724513-T2T.vcf.gz \
\
--reference chm13v2.0.fa \
--threads 6 \
\
-T Agilent_SureSelect_All_Exons_v7_hg38_Regions_hg38_T2T.bed \
--false-positives HG001_GRCh38_1_22_v4.2.1_benchmark_hg38_T2T.bed \
\
-o HG001
#+end_src
****** DONE Corriger FILTER : mieux mais toujours trop de négatifs. 3/4 SNP retrouvés
CLOSED: [2023-07-08 Sat 15:19] SCHEDULED: <2023-07-08 Sat>
Type Filter TRUTH.TOTAL TRUTH.TP TRUTH.FN QUERY.TOTAL QUERY.FP QUERY.UNK FP.gt FP.al METRIC.Recall METRIC.Precision METRIC.Frac_NA METRIC.F1_Score TRUTH.TOTAL.TiTv_ratio QUERY.TOTAL.TiTv_ratio TRUTH.TOTAL.het_hom_ratio QUERY.TOTAL.het_hom_ratio
INDEL ALL 413 246 167 751 289 215 2 98 0.595642 0.460821 0.286285 0.519629 NaN NaN 2.428571 2.465116
INDEL PASS 413 246 167 751 289 215 2 98 0.595642 0.460821 0.286285 0.519629 NaN NaN 2.428571 2.465116
SNP ALL 15883 15479 404 23597 5277 2841 46 44 0.974564 0.745760 0.120397 0.844947 3.017198 2.85705 5.560099 2.114633
SNP PASS 15883 15479 404 23597 5277 2841 46 44 0.974564 0.745760 0.120397 0.844947 3.017198 2.85705 5.560099 2.114633
******* DONE Vérifier qu'il ne reste plus de filtre autre que PASS
CLOSED: [2023-07-08 Sat 15:19]
#+begin_src
$ zgrep -c 'PASS' HG001_GRCh38_1_22_v4_lifted_merged.vcf.gz
3730505
$ zgrep -c '^chr' HG001_GRCh38_1_22_v4_lifted_merged.vcf.gz
3730506
#+end_src
****** TODO 1/4 SNP manquant ?
******* DONE Regarder avec Julia si ce sont vraiment des FP: 61/5277 qui ne le sont pas
CLOSED: [2023-07-09 Sun 12:09]
******* DONE Examiner les FP
CLOSED: [2023-07-30 Sun 22:05]
******* DONE Tester un FP
CLOSED: [2023-07-30 Sun 22:05]
2 │ chr1 608765 A G ./.:.:.:.:NOCALL:nocall:. 1/1:FP:.:ti:SNP:homalt:188
liftDown UCSC: rien en GIAB : vrai FP
3 │ chr1 762943 A G ./.:.:.:.:NOCALL:nocall:. 1/1:FP:.:ti:SNP:homalt:287
4 │ chr1 762945 A T ./.:.:.:.:NOCALL:nocall:. 1/1:FP:.:tv:SNP:homalt:287
Remaniements complexes ? Pas dans le gène en HG38
******* DONE La plupart des FP (4705/5566) sont homozygotes: erreur de référence ?
CLOSED: [2023-07-12 Wed 21:10] SCHEDULED: <2023-07-09 Sun>
Sur les 2 premiers variants, ils montrent en fait la différence entre T2T et GRCh38
Erreur à l'alignement ?
******** KILL relancer l'alignement
CLOSED: [2023-07-09 Sun 17:36]
******** DONE vérifier reads identiques hg38 et T2T: oui
CLOSED: [2023-07-09 Sun 16:36]
T2T CHR1608765
38 chr1:1180168-1180168 (
SRR14724513.24448214
SRR14724513.24448214
******* DONE Vérifier quelques variants sur IGV
CLOSED: [2023-07-09 Sun 17:36]
******* KILL Répartition des FP : cluster ?
CLOSED: [2023-07-09 Sun 17:36]
****** DONE Examiner les FP restant après correction selon séquence de référence
CLOSED: [2023-08-12 Sat 15:57]
****** HOLD Examiner les variants supprimé
****** TODO Enlever les FP qui correspondent à un changement dans le génome
SCHEDULED: <2023-08-14 Mon>
******* Condition:
- pas de variation à la position en GRCh38
- variantion homozygote
- la varation en T2T correspond au changement de pair de base GRC38 -> T2T
pour les SNP:
alt_T2T[i] = DNA_GRC38[j]
avec i la position en T2T et j la position en GRCh38
Note: définir un ID n'est pas correct car les variants peuvent être modifié par happy !
******* Idée
- Pour chaque FP, c'est un "faux" FP si
- REF en hg38 == ALT en T2T
- et REF en hg38 != REF en T2T
- et variant homozygote
Comment obtenir les séquences de réferences ?
1. liftover
2. blat sur la séquence autour du variant
3. identifier quelques reads contenant le variant et regarder leur aligneement en hg38
Après discussion avec Alexis: solution 3
******* Algorithme
1. Extraire les coordonnées en T2T des faux positifs *homozygote*
2. Pour chaque faux positif
1. lister 10 reads contenant le variant
2. pour chacun de ces reads, récupérer la séquence en T2T et GRCh38 via le nom du read dans le bam
3. si la séquence en T2T modifiée par le variant est "identique" à celle en GRCh38, alors on ignore ce faux positif
Note: on ignore les reads qui ont changé de chromosome entre les version
******* DONE Résultat préliminaire
CLOSED: [2023-07-23 Sun 14:30]
cf [[file:~/roam/research/bisonex/code/giab/giab-corrected.csv][script julia]]
3498 faux positifs en moins, soit 0.89 sensibilité
julia> tp=15479
julia> fp=5277
julia> tp/(tp+fp)
0.7457602620928888
julia> tp/(tp+(fp-3498))
0.8969173716537258
On est toujours en dessous des 97%
******* HOLD Corriger proprement VCF ou résultats Happy
******* TODO Adapter pour gérer plusieurs variants par read
SCHEDULED: <2023-08-14 Mon>
****** KILL Méthodologie du pangenome
CLOSED: [2023-07-31 Mon 22:29] SCHEDULED: <2023-07-30 Sun>
***** KILL Mail Yannis
CLOSED: [2023-07-08 Sat 10:44]
***** DONE Mail GIAB pour version T2T
CLOSED: [2023-07-07 Fri 18:37]
**** TODO HG002 :hg002:T2T:
**** TODO HG003 :hg003:T2T:
**** TODO HG004 :hg004:T2T:
**** DONE Plot : ashkenazim trio :
1.860963 | 2.247273 |
| INDEL | PASS | 549 | 489 | 60 | 899 | 64 | 340 | 8 | 17 | 0.890710 | 0.885510 | 0.378198 | 0.888102 | NaN | NaN | 1.860963 | 2.247273 |
| SNP | ALL | 21973 | 21462 | 511 | 26285 | 563 | 4263 | 68 | 16 | 0.976744 | 0.974435 | 0.162184 | 0.975588 | 3.00711 | 2.784686 | 1.591810 | 1.816145 |
| SNP | PASS | 21973 | 21462 | 511 | 26285 | 563 | 4263 | 68 | 16 | 0.976744 | 0.974435 | 0.162184 | 0.975588 | 3.00711 | 2.784686 | 1.591810 | 1.816145 |
******** Résumé
T2T
| Type | TRUTH.TOTAL | TRUTH.TP | TRUTH.FN | QUERY.TOTAL | QUERY.FP | QUERY.UNK | FP.gt | FP.al | METRIC.Recall | METRIC.Precision |
| INDEL | 413 | 246 | 167 | 751 | 289 | 215 | 2 | 93 | 0.595642 | 0.460821 |
| SNP | 11236 | 10985 | 251 | 23597 | 9771 | 2841 | 26 | 58 | 0.977661 | 0.529245 |
Hg38
| Type | TRUTH.TOTAL | TRUTH.TP | TRUTH.FN | QUERY.TOTAL | QUERY.FP | QUERY.UNK | FP.gt | FP.al | METRIC.Recall | METRIC.Precision |
| INDEL | 549 | 489 | 60 | 899 | 64 | 340 | 8 | 17 | 0.890710 | 0.885510 |
| SNP | 21973 | 21462 | 511 | 26285 | 563 | 4263 | 68 | 16 | 0.976744 | 0.974435 |
****** DONE Interesection des bed: similaire
CLOSED: [2023-07-04 Tue 23:11]
HG38
#+begin_src sh
bedtools intersect -a capture/Agilent_SureSelect_All_Exons_v7_hg38_Regions.bed -b /Work/Groups/bisonex/data/giab/GRCh38/HG001_GRCh38_1_22_v4.2.1_benchmark.bed | wc -l
#+end_src
204280
T2T
#+begin_src sh
bedtools intersect -a /Work/Groups/bisonex/data/giab/T2T/Agilent_SureSelect_All_Exons_v7_hg38_Regions_hg38_T2T.bed -b /Work/Groups/bisonex/data/giab/T2T/HG001_GRCh38_1_22_v4.2.1_benchmark_hg38_T2T.bed | wc -l
#+end_src
204021
****** DONE Vérifier la ligne de commande
CLOSED: [2023-07-04 Tue 23:38]
#+begin_src sh
hap.py \
HG001_GRCh38_1_22_v4_lifted_merged.vcf.gz \
HG001-SRX11061486_SRR14724513-T2T.vcf.gz \
\
--reference chm13v2.0.fa \
--threads 6 \
\
-T Agilent_SureSelect_All_Exons_v7_hg38_Regions_hg38_T2T.bed \
--false-positives HG001_GRCh38_1_22_v4.2.1_benchmark_hg38_T2T.bed \
\
-o HG001
#+end_src
****** DONE Corriger FILTER : mieux mais toujours trop de négatifs. 3/4 SNP retrouvés
CLOSED: [2023-07-08 Sat 15:19] SCHEDULED: <2023-07-08 Sat>
Type Filter TRUTH.TOTAL TRUTH.TP TRUTH.FN QUERY.TOTAL QUERY.FP QUERY.UNK FP.gt FP.al METRIC.Recall METRIC.Precision METRIC.Frac_NA METRIC.F1_Score TRUTH.TOTAL.TiTv_ratio QUERY.TOTAL.TiTv_ratio TRUTH.TOTAL.het_hom_ratio QUERY.TOTAL.het_hom_ratio
INDEL ALL 413 246 167 751 289 215 2 98 0.595642 0.460821 0.286285 0.519629 NaN NaN 2.428571 2.465116
INDEL PASS 413 246 167 751 289 215 2 98 0.595642 0.460821 0.286285 0.519629 NaN NaN 2.428571 2.465116
SNP ALL 15883 15479 404 23597 5277 2841 46 44 0.974564 0.745760 0.120397 0.844947 3.017198 2.85705 5.560099 2.114633
SNP PASS 15883 15479 404 23597 5277 2841 46 44 0.974564 0.745760 0.120397 0.844947 3.017198 2.85705 5.560099 2.114633
******* DONE Vérifier qu'il ne reste plus de filtre autre que PASS
CLOSED: [2023-07-08 Sat 15:19]
#+begin_src
$ zgrep -c 'PASS' HG001_GRCh38_1_22_v4_lifted_merged.vcf.gz
3730505
$ zgrep -c '^chr' HG001_GRCh38_1_22_v4_lifted_merged.vcf.gz
3730506
#+end_src
****** TODO 1/4 SNP manquant ?
******* DONE Regarder avec Julia si ce sont vraiment des FP: 61/5277 qui ne le sont pas
CLOSED: [2023-07-09 Sun 12:09]
******* DONE Examiner les FP
CLOSED: [2023-07-30 Sun 22:05]
******* DONE Tester un FP
CLOSED: [2023-07-30 Sun 22:05]
2 │ chr1 608765 A G ./.:.:.:.:NOCALL:nocall:. 1/1:FP:.:ti:SNP:homalt:188
liftDown UCSC: rien en GIAB : vrai FP
3 │ chr1 762943 A G ./.:.:.:.:NOCALL:nocall:. 1/1:FP:.:ti:SNP:homalt:287
4 │ chr1 762945 A T ./.:.:.:.:NOCALL:nocall:. 1/1:FP:.:tv:SNP:homalt:287
Remaniements complexes ? Pas dans le gène en HG38
******* DONE La plupart des FP (4705/5566) sont homozygotes: erreur de référence ?
CLOSED: [2023-07-12 Wed 21:10] SCHEDULED: <2023-07-09 Sun>
Sur les 2 premiers variants, ils montrent en fait la différence entre T2T et GRCh38
Erreur à l'alignement ?
******** KILL relancer l'alignement
CLOSED: [2023-07-09 Sun 17:36]
******** DONE vérifier reads identiques hg38 et T2T: oui
CLOSED: [2023-07-09 Sun 16:36]
T2T CHR1608765
38 chr1:1180168-1180168 (
SRR14724513.24448214
SRR14724513.24448214
******* DONE Vérifier quelques variants sur IGV
CLOSED: [2023-07-09 Sun 17:36]
******* KILL Répartition des FP : cluster ?
CLOSED: [2023-07-09 Sun 17:36]
****** DONE Examiner les FP restant après correction selon séquence de référence
CLOSED: [2023-08-12 Sat 15:57]
****** HOLD Examiner les variants supprimé
****** TODO Enlever les FP qui correspondent à un changement dans le génome
SCHEDULED: <2023-08-17 Thu>
******* Condition:
- pas de variation à la position en GRCh38
- variantion homozygote
- la varation en T2T correspond au changement de pair de base GRC38 -> T2T
pour les SNP:
alt_T2T[i] = DNA_GRC38[j]
avec i la position en T2T et j la position en GRCh38
Note: définir un ID n'est pas correct car les variants peuvent être modifié par happy !
******* Idée
- Pour chaque FP, c'est un "faux" FP si
- REF en hg38 == ALT en T2T
- et REF en hg38 != REF en T2T
- et variant homozygote
Comment obtenir les séquences de réferences ?
1. liftover
2. blat sur la séquence autour du variant
3. identifier quelques reads contenant le variant et regarder leur aligneement en hg38
Après discussion avec Alexis: solution 3
******* Algorithme
1. Extraire les coordonnées en T2T des faux positifs *homozygote*
2. Pour chaque faux positif
1. lister 10 reads contenant le variant
2. pour chacun de ces reads, récupérer la séquence en T2T et GRCh38 via le nom du read dans le bam
3. si la séquence en T2T modifiée par le variant est "identique" à celle en GRCh38, alors on ignore ce faux positif
Note: on ignore les reads qui ont changé de chromosome entre les version
******* DONE Résultat préliminaire
CLOSED: [2023-07-23 Sun 14:30]
cf [[file:~/roam/research/bisonex/code/giab/giab-corrected.csv][script julia]]
3498 faux positifs en moins, soit 0.89 sensibilité
julia> tp=15479
julia> fp=5277
julia> tp/(tp+fp)
0.7457602620928888
julia> tp/(tp+(fp-3498))
0.8969173716537258
On est toujours en dessous des 97%
******* HOLD Corriger proprement VCF ou résultats Happy
******* TODO Adapter pour gérer plusieurs variants par read
****** PROJ Méthodologie du pangenome
***** KILL Mail Yannis
CLOSED: [2023-07-08 Sat 10:44]
***** DONE Mail GIAB pour version T2T
CLOSED: [2023-07-07 Fri 18:37]
**** TODO HG002 :hg002:T2T:
**** TODO HG003 :hg003:T2T:
**** TODO HG004 :hg004:T2T:
**** DONE Plot : ashkenazim trio :
anger :xamscissors:
**** DONE Générer la liste
CLOSED: [2023-08-13 Sun 22:26] SCHEDULED: <2023-08-13 Sun 12:00>
**** DONE 2 variants chr20 dans NA12878 avec XAMscissors
CLOSED: [2023-08-15 Tue 00:13] SCHEDULED: <2023-08-13 Sun 13:00>
***** DONE Insertion
CLOSED: [2023-08-13 Sun 23:31] SCHEDULED: <2023-08-13 Sun>
On récupère le BAM
#+begin_src sh :dir /home/alex/code/bisonex/out
scp /Work/Users/apraga/bisonex/out/HG002-151002_7001448_0359_AC7F6GANXX-GRCh38/preprocessing/mapped/HG002-151002_7001448_0359_AC7F6GANXX-GRCh38.bam
#+end_src
On formate les données
#+begin_src sh :dir /home/alex/roam/research/bisonex/code/sanger
julia --project=.. formatXamscissors.jl
#+end_src
On génère les fastq
#+begin_src sh :dir /home/alex/roam/research/bisonex/code/sanger
julia --project=.. xamscissors.jl
#+end_src
On upload
#+begin_src sh :dir /home/alex/roam/research/bisonex/code/sanger
cd out
mv inserted_1.fq.gz HG002-sanger-inserted-chr20_1.fq.gz
mv inserted_2.fq.gz HG002-sanger-inserted-chr20_2.fq.gz
#+end_src
On upload
#+begin_src sh :dir /home/alex/roam/research/bisonex/code/sanger
scp HG002-sanger-inserted-chr20_*.fq.gz meso:/Work/Users/apraga/bisonex/data/
#+end_src
Fichier de configuration
patient,sample,fastq1,fastq2
HG002,sanger-chr20,data/HG002-sanger-inserted-chr20_1.fq.gz,data/HG002-sanger-inserted-chr20_2.fq.gz
On lance la simulation
#+begin_src
nextflow run main.nf -profile standard,helios --input=samples-synthetic.csv --genome=GRCh38 -bg
#+end_src
***** DONE Vérifier post alignement
CLOSED: [2023-08-13 Sun 23:44] SCHEDULED: <2023-08-13 Sun>
***** DONE Vérifier haplotypecaller
CLOSED: [2023-08-15 Tue 00:13] SCHEDULED: <2023-08-13 Sun>
***** DONE Vérifier post-filtre
CLOSED: [2023-08-15 Tue 00:13]
**** TODO Tout insérer dans NA12878 avec XAMscissors
SCHEDULED: <2023-08-14 Mon>
**** TODO Données simuscop 200x
SCHEDULED: <2023-08-22 Tue>
* Résultats
** TODO Speed-up BWA-mem
SCHEDULED: <2023-08-26 Sat>
** TODO Speed-up Hapotypecaller
SCHEDULED: <2023-08-26 Sat>
anger :xamscissors:
**** DONE Générer la liste
CLOSED: [2023-08-13 Sun 22:26] SCHEDULED: <2023-08-13 Sun 12:00>
**** DONE 2 variants chr20 dans HG002 avec XAMscissors: homozygote
CLOSED: [2023-08-15 Tue 00:13] SCHEDULED: <2023-08-13 Sun 13:00>
***** DONE Insertion
CLOSED: [2023-08-13 Sun 23:31] SCHEDULED: <2023-08-13 Sun>
On récupère le BAM
#+begin_src sh :dir /home/alex/code/bisonex/out
scp /Work/Users/apraga/bisonex/out/HG002-151002_7001448_0359_AC7F6GANXX-GRCh38/preprocessing/mapped/HG002-151002_7001448_0359_AC7F6GANXX-GRCh38.bam
#+end_src
On formate les données et on génère les fastq
#+begin_src sh :dir /home/alex/roam/research/bisonex/code/sanger
julia --project=.. xamscissors.jl
#+end_src
On upload
#+begin_src sh :dir /home/alex/roam/research/bisonex/code/sanger
cd out
mv inserted_1.fq.gz HG002-sanger-inserted-chr20_1.fq.gz
mv inserted_2.fq.gz HG002-sanger-inserted-chr20_2.fq.gz
#+end_src
On upload
#+begin_src sh :dir /home/alex/roam/research/bisonex/code/sanger
scp HG002-sanger-inserted-chr20_*.fq.gz meso:/Work/Users/apraga/bisonex/data/
#+end_src
Fichier de configuration
patient,sample,fastq1,fastq2
HG002,sanger-chr20,data/HG002-sanger-inserted-chr20_1.fq.gz,data/HG002-sanger-inserted-chr20_2.fq.gz
On lance la simulation
#+begin_src
nextflow run main.nf -profile standard,helios --input=samples-synthetic.csv --genome=GRCh38 -bg
#+end_src
***** DONE Vérifier post alignement
CLOSED: [2023-08-13 Sun 23:44] SCHEDULED: <2023-08-13 Sun>
***** DONE Vérifier haplotypecaller
CLOSED: [2023-08-15 Tue 00:13] SCHEDULED: <2023-08-13 Sun>
***** DONE Vérifier post-filtre
CLOSED: [2023-08-15 Tue 00:13]
**** DONE 2 variant chr20 dans HG002 heterozygote
CLOSED: [2023-08-15 Tue 13:45] SCHEDULED: <2023-08-15 Tue>
La zygosité n'était pas correcte...
**** TODO [#B] Tout insérer dans NA12878 avec XAMscissors
SCHEDULED: <2023-08-14 Mon>
***** Insertion
On récupère le BAM
#+begin_src sh :dir /home/alex/code/bisonex/out
scp /Work/Users/apraga/bisonex/out/HG002-151002_7001448_0359_AC7F6GANXX-GRCh38/preprocessing/mapped/HG002-151002_7001448_0359_AC7F6GANXX-GRCh38.bam
#+end_src
#+begin_src sh :dir /home/alex/roam/research/bisonex/code/sanger
julia --project=.. xamscissors.jl
cd out-all
mv inserted_1.fq.gz HG002-sanger-inserted-_1.fq.gz
mv inserted_2.fq.gz HG002-sanger-inserted-_2.fq.gz
scp HG002-sanger-inserted-all_*.fq.gz meso:/Work/Users/apraga/bisonex/data/
#+end_src
Fichier de configuration
patient,sample,fastq1,fastq2
HG002,sanger-all,data/HG002-sanger-inserted-all_1.fq.gz,data/HG002-sanger-inserted-all_2.fq.gz
On lance la simulation
#+begin_src
nextflow run main.nf -profile standard,helios --input=samples-synthetic.csv --genome=GRCh38 -bg
#+end_src
**** TODO Données simuscop 200x
SCHEDULED: <2023-08-22 Tue>
* Résultats
** TODO Speed-up BWA-mem
SCHEDULED: <2023-08-26 Sat>
** TODO Speed-up Hapotypecaller
SCHEDULED: <2023-08-26 Sat>