B:BD[
5.58110] → [
6.29:8221]
B:BD[
6.8221] → [
3.16415:24607]
uvais résultat
SCHEDULED: <2023-07-03 Mon>
****** DONE 2x moins de variants, beaucoup trop de FP : erreur de FILTER
CLOSED: [2023-07-08 Sat 11:17] SCHEDULED: <2023-07-07 Fri>
******* DONE [#A] Certains FP sont des vrais FP: erreur de flag et non d'happy
CLOSED: [2023-07-08 Sat 10:43] SCHEDULED: <2023-07-07 Fri>
******** Ex: chr1:408793
G C 600.64 . BS=408793;Regions=CONF,TS_contained GT:BD:BK:BI:BVT:BLT:QQ ./.:.:.:.:NOCALL:nocall:. 0/1:FP:.:tv:SNP:het:600.64
CLOSED: [2023-07-07 Fri 20:03]
********* DONE Vérifier que l'on est bien dans la zone "high confidence" définie dans le bed
oui:
chr1 408514 414397 . 500 +
********* DONE Vérifier que le variant est bien dans le truth et query
CLOSED: [2023-07-07 Fri 20:14]
truth
chr1 408793 . G C 50 MismatchedRefAllele AttemptedAlleles=C*->G;AttemptedLocus=chr1:408793-408793;SwappedAlleles;callable=CS_HiSeqPE300xGATK_callable,CS_10XLRGATK_callable,CS_CCS15kb_20kbGATK4_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbGATK4,CGnormal,HiSeqPE300xfreebayes,CCS15kb_20kbDV,10XLRGATK;callsets=6;datasetnames=HiSeqPE300x,CCS15kb_20kb,CGnormal,10XChromiumLR;datasets=4;datasetsmissingcall=IonExome,SolidSE75bp;filt=CS_CCS15kb_20kbDV_filt,CS_CCS15kb_20kbGATK4_filt;platformnames=Illumina,PacBio,CG,10X;platforms=4 GT:AD:ADALL:DP:GQ 1/0:153,159:141,132:689:99
query
chr1 408793 . G C 600.64 . AC=1;AF=0.500;AN=2;BaseQRankSum=-1.390;DP=64;ExcessHet=0.0000;FS=0.961;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=9.53;ReadPosRankSum=-3.146;SOR=0.906 GT:AD:DP:GQ:PL 0/1:33,30:63:99:608,0,798
********* DONE Tester happy sur ce variant: non retrouvé
CLOSED: [2023-07-08 Sat 10:42]
Il faut au moins 2 variant + les index (sinon "no chromosome in common")
Donc on augmente la taille de la région
#+begin_src sh :dir ~/code/bisonex/out
bcftools view HG001_GRCh38_1_22_v4_lifted_merged.vcf.gz -o lifted_test.vcf.gz chr1:408793-408900
bcftools view HG001-SRX11061486_SRR14724513-T2T.vcf.gz -o run_test.vcf.gz chr1:408793-408900
bcftools index lifted_test.vcf.gz
bcftools index run_test.vcf.gz
#+end_src
#+begin_src sh
hap.py lifted_test.vcf.gz run_test.vcf.gz --reference ../data/chm13v2.0.fa -o test
#+end_src
On trouve bien les 2 variants. L'indel est retrouvé mais pas le SNP
chr1 408793 . G C 600.64 . BS=408793 GT:BD:BK:BI:BVT:BLT:QQ ./.:.:.:.:NOCALL:nocall:. 0/1:FP:.:tv:SNP:het:600.64
chr1 408874 . CGAA C 1476.6 . BS=408874 GT:BD:BK:BI:BVT:BLT:QQ 0/1:TP:gm:d1_5:INDEL:het:1476.6 0/1:TP:gm:d1_5:INDEL:het:1476.6
********* DONE Changer flag : mismatch -> PASS: OUI
CLOSED: [2023-07-08 Sat 10:42]
********* DONE Changer haplotype 1/0 -> 0/1 truth: idem
CLOSED: [2023-07-08 Sat 10:35]
********* DONE Tester avec vcfeval : idemnt
CLOSED: [2023-07-08 Sat 10:43]
******* DONE Comparer hg38 et T2T: 2x moinsr de variants, trop de FP et FN
CLOSED: [2023-07-07 Fri 18:38]
T2T
| 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 | 93 | 0.595642 | 0.460821 | 0.286285 | 0.519629 | NaN | NaN | 2.428571 | 2.465116 |
| INDEL | PASS | 413 | 246 | 167 | 751 | 289 | 215 | 2 | 93 | 0.595642 | 0.460821 | 0.286285 | 0.519629 | NaN | NaN | 2.428571 | 2.465116 |
| SNP | ALL | 11236 | 10985 | 251 | 23597 | 9771 | 2841 | 26 | 58 | 0.977661 | 0.529245 | 0.120397 | 0.686734 | 3.11461 | 2.85705 | 3.640644 | 2.114633 |
| SNP | PASS | 11236 | 10985 | 251 | 23597 | 9771 | 2841 | 26 | 58 | 0.977661 | 0.529245 | 0.120397 | 0.686734 | 3.11461 | 2.85705 | 3.640644 | 2.114633 |
Hg38
| 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 | 549 | 489 | 60 | 899 | 64 | 340 | 8 | 17 | 0.890710 | 0.885510 | 0.378198 | 0.888102 | NaN | NaN | 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 ?
SCHEDULED: <2023-07-08 Sat>
******* DONE Regarder avec Julia si ce sont vraiment des FP: 61/5277 qui ne le sont pas
CLOSED: [2023-07-09 Sun 12:09]
******* TODO Examiner les FP
******* TODO Tester un FP
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
******* TODO La plupart des FP (4705/5566) sont homozygotes: erreur de référence ?
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
******* TODO Enlever les FP qui correspondent à un changement dans le génome
SCHEDULED: <2023-07-09 Sun>
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
Algorithm
1. Définir un ID sur le VCF en hg38 :
bcftools annotate --set-id +'hg38\_%CHROM\_%POS\_%REF\_%FIRST_ALT' file.vcf
2. Pour chaque FP, récupérer depuis le VCF généré par happy
- l'identifiant en hg38
- les coordonnées en T2T
- séquence alternative
3. Pour chaque identifiant, récupérer depuis le VCF en hg38
- séquence de réference
- les coordonnées
4. Comparer ALT en T2T à la REF en hg38
******* 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]
******* TODO 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]
**** DONE NA12878 :na12878:hg38:
CLOSED: [2023-06-30 Fri 22:30]
***** DONE Discussion alexis : Mail
CLOSED: [2023-03-29 Wed 22:40]
Avec le patient NA12878 et comparaison avec hap.py du VCF de Genome In A Bottle ("gold" standard), on avait pour rappel
- sensibilité (=recall) 71% pour indel, 85% SNP
- précision (= VPP) 69 et 97% respectivement
| Type | TRUTH | TP | FN | QUERY | FP | UNK | FP.gt | FP.al | Recall | Precision |
| INDEL | 4871 | 3461 | 1410 | 7048 | 1554 | 1987 | 193 | 346 | 0.710532 | 0.692946 |
| SNP | 46032 | 39369 | 6663 | 44600 | 1186 | 4041 | 304 | 30 | 0.855253 | 0.970759 |
Les statistiques sur les génomes sont bien meilleurs (cf precisionFDA challenge).
Pour les exome, un article [1] a fait a des meilleures stats sur ce patient avec BWA et GATK mais ils ont moins de variant (on a presque un facteur 2 !).
Je soupçonne qu'on ne travaille pas sur les mêmes zones de capture (pas réussi à récupérer leur .bed)
| Exome | Type | TP | FP | FN | Sensitivity | Precision | F-Score | FDR |
| 1 | SNV | 23689 | 1397 | 613 | 0.975 | 0.944 | 0.959 | 0.057 |
| 2 | SNV | 23946 | 865 | 356 | 0.985 | 0.965 | 0.975 | 0.036 |
| 1 | indel | 1254 | 72 | 75 | 0.944 | 0.946 | 0.945 | 0.054 |
| 2 | indel | 1309 | 10 | 20 | 0.985 | 0.992 | 0.989 | 0.008 |
Pour essayer d'améliorer les statistiques :
- La version du génome GRC38 vs GRCh38.p13 ne change quasiment rien
- Désactiver dbSNP ne change strictement rien pour le variant calling
J'ai exploré les faux négatifs :
- la grande majorité n'est juste pas vue (ce n'est pas un problème d'haploïde/génotype)
- la répartition par chromosome est relativement homogène, sauf sur le 6 ()
- la majorité est en 5' et 3'UTR (selon Best refseq)
Conclusion: je pense m'arrêter là pour la validation du variant calling par manque de temps. Il faudrait creuser pour savoir pourquoi certains variants ne sont pas vus par GATK mais ce n'est pas la majorité. En tout cas, je peux justifier d'une première analyse pour la thèse.
Ça te va ?
[1]
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-019-2928-9
Résultats ici https://static-content.springer.com/esm/art%3A10.1186%2Fs12859-019-2928-9/MediaObjects/12859_2019_2928_MOESM8_ESM.pdf
***** DONE Comparaison
CLOSED: [2023-03-04 Sat 11:14]
HGREF=/Work/Groups/bisonex/data-alexis-reference/genome/GRCh38_latest_genomic.fna ./result/bin/hap.py /Work/Groups/bisonex/NA12878/HG001_GRCh38_1_22_v4.2.1
_benchmark_renamed.vcf.gz script/files/vcf/NA12878_NIST7035_vep_annot.vcf -f /Work/Groups/bison
ex/NA12878/HG001_GRCh38_1_22_v4.2.1_benchmark.bed -o test
na1878.slurm
#+begin_src slurm
#!/bin/bash
#SBATCH -c 4
#SBATCH -p smp
#SBATCH --time=01:00:00
#SBATCH --mem=32G
module load nix/2.11.0
export HGREF=/Work/Groups/bisonex/data-alexis-reference/genome/GRCh38_latest_genomic.fna
dir=/Work/Groups/bisonex/data/NA12878/GRCh38
hap.py ${dir}/HG001_GRCh38_1_22_v4.2.1_benchmark.vcf.gz script/files/vcf/NA12878_NIST7035.vcf -f ${dir}/HG001_GRCh38_1_22_v4.2.1_benchmark.bed -o test
#+end_src
****** KILL beaucoup trop de faux négatifs
CLOSED: [2023-02-17 Fri 19:37]
******* DONE Test 1 : vep annot : beaucoup trop de faux négatif
CLOSED: [2023-02-06 lun. 13:40]
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 276768 274 276494 1500 257 968 26 15 0.000990 0.516917 0.645333 0.001976 NaN NaN 1.483361 6.129187
INDEL PASS 276768 274 276494 1500 257 968 26 15 0.000990 0.516917 0.645333 0.001976 NaN NaN 1.483361 6.129187
SNP ALL 1937706 1193 1936513 3338 106 2037 11 2 0.000616 0.918524 0.610246 0.001231 2.0785 1.861183 1.539064 2.703663
SNP PASS 1937706 1193 1936513 3338 106 2037 11 2 0.000616 0.918524 0.610246 0.00123
uvais résultat
SCHEDULED: <2023-07-03 Mon>
****** DONE 2x moins de variants, beaucoup trop de FP : erreur de FILTER
CLOSED: [2023-07-08 Sat 11:17] SCHEDULED: <2023-07-07 Fri>
******* DONE [#A] Certains FP sont des vrais FP: erreur de flag et non d'happy
CLOSED: [2023-07-08 Sat 10:43] SCHEDULED: <2023-07-07 Fri>
******** Ex: chr1:408793
G C 600.64 . BS=408793;Regions=CONF,TS_contained GT:BD:BK:BI:BVT:BLT:QQ ./.:.:.:.:NOCALL:nocall:. 0/1:FP:.:tv:SNP:het:600.64
CLOSED: [2023-07-07 Fri 20:03]
********* DONE Vérifier que l'on est bien dans la zone "high confidence" définie dans le bed
oui:
chr1 408514 414397 . 500 +
********* DONE Vérifier que le variant est bien dans le truth et query
CLOSED: [2023-07-07 Fri 20:14]
truth
chr1 408793 . G C 50 MismatchedRefAllele AttemptedAlleles=C*->G;AttemptedLocus=chr1:408793-408793;SwappedAlleles;callable=CS_HiSeqPE300xGATK_callable,CS_10XLRGATK_callable,CS_CCS15kb_20kbGATK4_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbGATK4,CGnormal,HiSeqPE300xfreebayes,CCS15kb_20kbDV,10XLRGATK;callsets=6;datasetnames=HiSeqPE300x,CCS15kb_20kb,CGnormal,10XChromiumLR;datasets=4;datasetsmissingcall=IonExome,SolidSE75bp;filt=CS_CCS15kb_20kbDV_filt,CS_CCS15kb_20kbGATK4_filt;platformnames=Illumina,PacBio,CG,10X;platforms=4 GT:AD:ADALL:DP:GQ 1/0:153,159:141,132:689:99
query
chr1 408793 . G C 600.64 . AC=1;AF=0.500;AN=2;BaseQRankSum=-1.390;DP=64;ExcessHet=0.0000;FS=0.961;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=9.53;ReadPosRankSum=-3.146;SOR=0.906 GT:AD:DP:GQ:PL 0/1:33,30:63:99:608,0,798
********* DONE Tester happy sur ce variant: non retrouvé
CLOSED: [2023-07-08 Sat 10:42]
Il faut au moins 2 variant + les index (sinon "no chromosome in common")
Donc on augmente la taille de la région
#+begin_src sh :dir ~/code/bisonex/out
bcftools view HG001_GRCh38_1_22_v4_lifted_merged.vcf.gz -o lifted_test.vcf.gz chr1:408793-408900
bcftools view HG001-SRX11061486_SRR14724513-T2T.vcf.gz -o run_test.vcf.gz chr1:408793-408900
bcftools index lifted_test.vcf.gz
bcftools index run_test.vcf.gz
#+end_src
#+begin_src sh
hap.py lifted_test.vcf.gz run_test.vcf.gz --reference ../data/chm13v2.0.fa -o test
#+end_src
On trouve bien les 2 variants. L'indel est retrouvé mais pas le SNP
chr1 408793 . G C 600.64 . BS=408793 GT:BD:BK:BI:BVT:BLT:QQ ./.:.:.:.:NOCALL:nocall:. 0/1:FP:.:tv:SNP:het:600.64
chr1 408874 . CGAA C 1476.6 . BS=408874 GT:BD:BK:BI:BVT:BLT:QQ 0/1:TP:gm:d1_5:INDEL:het:1476.6 0/1:TP:gm:d1_5:INDEL:het:1476.6
********* DONE Changer flag : mismatch -> PASS: OUI
CLOSED: [2023-07-08 Sat 10:42]
********* DONE Changer haplotype 1/0 -> 0/1 truth: idem
CLOSED: [2023-07-08 Sat 10:35]
********* DONE Tester avec vcfeval : idemnt
CLOSED: [2023-07-08 Sat 10:43]
******* DONE Comparer hg38 et T2T: 2x moinsr de variants, trop de FP et FN
CLOSED: [2023-07-07 Fri 18:38]
T2T
| 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 | 93 | 0.595642 | 0.460821 | 0.286285 | 0.519629 | NaN | NaN | 2.428571 | 2.465116 |
| INDEL | PASS | 413 | 246 | 167 | 751 | 289 | 215 | 2 | 93 | 0.595642 | 0.460821 | 0.286285 | 0.519629 | NaN | NaN | 2.428571 | 2.465116 |
| SNP | ALL | 11236 | 10985 | 251 | 23597 | 9771 | 2841 | 26 | 58 | 0.977661 | 0.529245 | 0.120397 | 0.686734 | 3.11461 | 2.85705 | 3.640644 | 2.114633 |
| SNP | PASS | 11236 | 10985 | 251 | 23597 | 9771 | 2841 | 26 | 58 | 0.977661 | 0.529245 | 0.120397 | 0.686734 | 3.11461 | 2.85705 | 3.640644 | 2.114633 |
Hg38
| 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 | 549 | 489 | 60 | 899 | 64 | 340 | 8 | 17 | 0.890710 | 0.885510 | 0.378198 | 0.888102 | NaN | NaN | 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 ?
SCHEDULED: <2023-07-08 Sat>
******* DONE Regarder avec Julia si ce sont vraiment des FP: 61/5277 qui ne le sont pas
CLOSED: [2023-07-09 Sun 12:09]
******* TODO Examiner les FP
******* TODO Tester un FP
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
******* TODO La plupart des FP (4705/5566) sont homozygotes: erreur de référence ?
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
******* TODO Enlever les FP qui correspondent à un changement dans le génome
SCHEDULED: <2023-07-09 Sun>
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 !
Algorithme
1. 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
******* 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]
******* TODO 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]
**** DONE NA12878 :na12878:hg38:
CLOSED: [2023-06-30 Fri 22:30]
***** DONE Discussion alexis : Mail
CLOSED: [2023-03-29 Wed 22:40]
Avec le patient NA12878 et comparaison avec hap.py du VCF de Genome In A Bottle ("gold" standard), on avait pour rappel
- sensibilité (=recall) 71% pour indel, 85% SNP
- précision (= VPP) 69 et 97% respectivement
| Type | TRUTH | TP | FN | QUERY | FP | UNK | FP.gt | FP.al | Recall | Precision |
| INDEL | 4871 | 3461 | 1410 | 7048 | 1554 | 1987 | 193 | 346 | 0.710532 | 0.692946 |
| SNP | 46032 | 39369 | 6663 | 44600 | 1186 | 4041 | 304 | 30 | 0.855253 | 0.970759 |
Les statistiques sur les génomes sont bien meilleurs (cf precisionFDA challenge).
Pour les exome, un article [1] a fait a des meilleures stats sur ce patient avec BWA et GATK mais ils ont moins de variant (on a presque un facteur 2 !).
Je soupçonne qu'on ne travaille pas sur les mêmes zones de capture (pas réussi à récupérer leur .bed)
| Exome | Type | TP | FP | FN | Sensitivity | Precision | F-Score | FDR |
| 1 | SNV | 23689 | 1397 | 613 | 0.975 | 0.944 | 0.959 | 0.057 |
| 2 | SNV | 23946 | 865 | 356 | 0.985 | 0.965 | 0.975 | 0.036 |
| 1 | indel | 1254 | 72 | 75 | 0.944 | 0.946 | 0.945 | 0.054 |
| 2 | indel | 1309 | 10 | 20 | 0.985 | 0.992 | 0.989 | 0.008 |
Pour essayer d'améliorer les statistiques :
- La version du génome GRC38 vs GRCh38.p13 ne change quasiment rien
- Désactiver dbSNP ne change strictement rien pour le variant calling
J'ai exploré les faux négatifs :
- la grande majorité n'est juste pas vue (ce n'est pas un problème d'haploïde/génotype)
- la répartition par chromosome est relativement homogène, sauf sur le 6 ()
- la majorité est en 5' et 3'UTR (selon Best refseq)
Conclusion: je pense m'arrêter là pour la validation du variant calling par manque de temps. Il faudrait creuser pour savoir pourquoi certains variants ne sont pas vus par GATK mais ce n'est pas la majorité. En tout cas, je peux justifier d'une première analyse pour la thèse.
Ça te va ?
[1]
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-019-2928-9
Résultats ici https://static-content.springer.com/esm/art%3A10.1186%2Fs12859-019-2928-9/MediaObjects/12859_2019_2928_MOESM8_ESM.pdf
***** DONE Comparaison
CLOSED: [2023-03-04 Sat 11:14]
HGREF=/Work/Groups/bisonex/data-alexis-reference/genome/GRCh38_latest_genomic.fna ./result/bin/hap.py /Work/Groups/bisonex/NA12878/HG001_GRCh38_1_22_v4.2.1
_benchmark_renamed.vcf.gz script/files/vcf/NA12878_NIST7035_vep_annot.vcf -f /Work/Groups/bison
ex/NA12878/HG001_GRCh38_1_22_v4.2.1_benchmark.bed -o test
na1878.slurm
#+begin_src slurm
#!/bin/bash
#SBATCH -c 4
#SBATCH -p smp
#SBATCH --time=01:00:00
#SBATCH --mem=32G
module load nix/2.11.0
export HGREF=/Work/Groups/bisonex/data-alexis-reference/genome/GRCh38_latest_genomic.fna
dir=/Work/Groups/bisonex/data/NA12878/GRCh38
hap.py ${dir}/HG001_GRCh38_1_22_v4.2.1_benchmark.vcf.gz script/files/vcf/NA12878_NIST7035.vcf -f ${dir}/HG001_GRCh38_1_22_v4.2.1_benchmark.bed -o test
#+end_src
****** KILL beaucoup trop de faux négatifs
CLOSED: [2023-02-17 Fri 19:37]
******* DONE Test 1 : vep annot : beaucoup trop de faux négatif
CLOSED: [2023-02-06 lun. 13:40]
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 276768 274 276494 1500 257 968 26 15 0.000990 0.516917 0.645333 0.001976 NaN NaN 1.483361 6.129187
INDEL PASS 276768 274 276494 1500 257 968 26 15 0.000990 0.516917 0.645333 0.001976 NaN NaN 1.483361 6.129187
SNP ALL 1937706 1193 1936513 3338 106 2037 11 2 0.000616 0.918524 0.610246 0.001231 2.0785 1.861183 1.539064 2.703663
SNP PASS 1937706 1193 1936513 3338 106 2037 11 2 0.000616 0.918524 0.610246 0.00123