apraga/org - Change 6CRMCEKX7FPZ6JDUEJ3BEAHO5EX5HUNI7FALEHA2FFSPRYHS63BQC

Merge branch 'master' of git.sr.ht:~scut/org

Created by Alexis Praga on July 18, 2023

6CRMCEKX7FPZ6JDUEJ3BEAHO5EX5HUNI7FALEHA2FFSPRYHS63BQC

Dependencies

In channels

main

Change contents

Replacement in projects/bisonex.org at line 34 [5.35]

B:BD[4.8003] → [4.8003:8134]

B:BD[4.8134] → [2.57506:65698]

BD:BD[4.8134] → [2.57506:65698]

B:BD[4.8134] → [3.29:7940]

lter  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]
******* HOLD Examiner les FP
******* HOLD 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
******* 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
******* 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.001231                  2.0785           
_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
******* 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
******* 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 !
******** Algorithmes
 - 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
******* 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

[4.8003]

[3.7940]

lter  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]
******* HOLD Examiner les FP
******* HOLD 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
******* 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
******* 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 !
******** Algorithmes
 - 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
******* 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

Replacement in projects/bisonex.org at line 36 [5.35]

B:BD[3.16132] → [3.16132:16335]

∅:D[3.16335] → [6.43625:51603]

B:BD[6.43625] → [6.43625:51603]

  922       490       432          942       439          0     32     56       0.531453          0.533970
INDEL   PASS          922       490       432          942       439          0     32     56  
     0.531453          0.533970
  SNP    ALL        24618     18689      5929        21257      2571          0    239     17       0.759160          0.879052
  SNP   PASS        24618     18689      5929        21257      2571          0    239     17       0.759160          0.879052
 METRIC.Frac_NA  METRIC.F1_Score  TRUTH.TOTAL.TiTv_ratio  QUERY.TOTAL.TiTv_ratio  TRUTH.TOTAL.het_hom_ratio  QUERY.TOTAL.het_hom_ratio
            0.0         0.532709                     NaN                     NaN                   1.628743                   2.799180
            0.0         0.532709                     NaN                     NaN                   1.628743                   2.799180
            0.0         0.814719                2.899604                2.881526                   1.598479                   1.564378
            0.0         0.814719                2.899604                2.881526                   1.598479                   1.564378
******* KILL Comprendre pour les FN explosent avec vep annot
CLOSED: [2023-02-24 Fri 23:44]
NC_000001.11	924024	.	C	G
non couverte dans bam -> 1er exons de SAMD11 mais non couverte
Il est dans NM_001385641.1 mais pas dans NM_152486.4
****** KILL Obtenir les zones couvertes depuis le bam directement
CLOSED: [2023-02-24 Fri 23:44]
#+begin_src  sh
samtools sort NA12878_chr1.bam -o NA12878_chr1_sorted.bam
bedtools genomecov -ibam NA12878_chr1_sorted.bam  -bg > chr1.bedgraph
head chr1.bedgraph
#+end_src
#+RESULTS:
: NC_000001.11	10001	10021	1
: NC_000001.11	10021	10059	2
: NC_000001.11	10059	10063	1
: NC_000001.11	10063	10087	2
: NC_000001.11	10087	10110	1
: NC_000001.11	10110	10113	3
: NC_000001.11	10113	10121	2
On prend les régions avec 20 reads
awk '$4 > 20 {print $1"\t"$2"\t"$3}' chr1.bedgraph  > tomerge.bed
On fusionne les régions
bedtools merge -i tomerge.bed > exons_from_bam.bed
Problème : on manque parfois le bord des exons...
****** KILL Vérifier un bam de référence de NA12878
CLOSED: [2023-02-24 Fri 23:44]
Notamment pour la définitions des exons, on a des reads qui ne semblent pas alignés sur les bons exons selon IGV
On donne directemet à IGV le bam d'illuma WES ( https://github.com/genome-in-a-bottle/giab_data_indexes )
IGV n'aime pas le ftp apparement...
On récupére 2 échantillons pour ce patient sur HiSeq Exome
#+begin_src sh :dir /ssh:meso:/Work/Groups/bisonex/data/giab/GRCh38
wget https://raw.githubusercontent.com/genome-in-a-bottle/giab_data_indexes/master/NA12878/alignment.index.NA12878_HiSeq_Exome_Garvan_GRCh37_09252015 -O todl.txt
awk '{print $1}' todl.txt | wget -i -
#+end_src
On récupére le premier échantillons en prenant just le chromosome 1
#+begin_src
samtools merge NA12878-NIST7035-HiSeq_Exome_Garan_GRCh37.bam project
.NIST_NIST7035_H7AP8ADXX_TAAGGCGA_1_NA12878.bwa.markDuplicates.bam project.NIST_NIST7035_H7AP8ADXX_TA
AGGCGA_2_NA12878.bwa.markDuplicates.bam
samtools index NA12878-NIST7035-HiSeq_Exome_Garan_GRCh37.bam
#+end_src
******* Étude
NC_000001.11:924024 :
- bed d’illumina : pas correct
- bed refseq : meux
- bam :  exon non ciblé, probablement parce que Mane select n’a pas été utilisé mais refseq (https://www.genenames.org/data/gene-symbol-report/#!/hgnc_id/HGNC:28706)
NC_000001.11:924310: idem
NC_000001.11:924321: idem
NC_000001.11:924533: idem
NC_000001.11:966227: le bed ne couvre pas le bon exon !
https://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr1%3A965863%2D966591&hgsid=295298291_vhDbIv5YIckCKLpGB7UUaZi2cAkr
NC_000001.11:966272
****** KILL Filtrer variants introniques de référence avec vep
CLOSED: [2023-02-24 Fri 23:44]
******* KILL variant calling seulf + seulement -f: nombreux FP
CLOSED: [2023-02-24 Fri 23:44]
/Work/Users/apraga/bisonex/work/68/f1cf72a5a4078fdf743fb3844b369a
 Type Filter  TRUTH.TOTAL  TRUTH.TP  TRUTH.FN  QUERY.TOTAL  QUERY.FP  QUERY.UNK  FP.gt  FP.al  METRIC.Recall  METRIC.Precision
INDEL    ALL          519       284       235        52169     37789      14094     12     64       0.547206          0.007511
INDEL   PASS          519       284       235        52169     37789      14094     12     64       0.547206          0.007511
  SNP    ALL        22131     17434      4697       357652    305313      34904    189     32       0.787764          0.054020
  SNP   PASS        22131     17434      4697       357652    305313      34904    189     32       0.787764          0.054020
 METRIC.Frac_NA  METRIC.F1_Score  TRUTH.TOTAL.TiTv_ratio  QUERY.TOTAL.TiTv_ratio  TRUTH.TOTAL.het_hom_ratio  QUERY.TOTAL.het_hom_ratio
       0.270160         0.014820                     NaN                     NaN                   1.775956                   0.509510
       0.270160         0.014820                     NaN                     NaN                   1.775956                   0.509510
       0.097592         0.101108                2.971834                1.429533                   1.579776                   0.324923
       0.097592         0.101108                2.971834                1.429533                   1.579776                   0.324923
****** DONE Bed donnés par GIAB (en hg19) : résultats corrects :resultats:
CLOSED: [2023-03-09 Thu 22:42] SCHEDULED: <2023-03-02 Thu>
Téléchargé à la main et converti via UCSCS. À automatiser, voir : *hg38
Sans oublier de renommer les chromosomes !
| Type  | TRUTH.TOTAL | TRUTH.TP | TRUTH.FN | QUERY.TOTAL | QUERY.FP | QUERY.UNK | FP.gt | FP.al | METRIC.Recall | METRIC.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 |
 Type   Ssensibilité    VPP
 INDEL       0.710532   0.692946 
 SNP         0.855253   0.970759
 | METRIC.Frac_NA | METRIC.F1_Score | TRUTH.TOTAL.TiTv_ratio | QUERY.TOTAL.TiTv_ratio | TRUTH.TOTAL.het_hom_ratio | QUERY.TOTAL.het_hom_ratio |
 |----------------+-----------------+------------------------+------------------------+---------------------------+---------------------------|
 |       0.281924 |        0.701629 |                        |                        |        1.6174985978687606 |        3.0674091441969518 |
 |       0.090605 |        0.909353 |      2.529551552318896 |     2.4019675131548843 |        1.6206857273037931 |        1.6274012964054214 |
****** DONE Re-tester avec exons Refseq et option -T: résultats un peu moins bons
CLOSED: [2023-03-09 Thu 22:42] SCHEDULED: <2023-03-08 Wed>
On utilise directement les coordonées données par refseq
 Type Filter  TRUTH.TOTAL  TRUTH.TP  TRUTH.FN  QUERY.TOTAL  QUERY.FP  QUERY.UNK  FP.gt  FP.al  METRIC.Recall  METRIC.Precision
INDEL    ALL         7226      3417      3809         6978      1599       1918    228    353       0.472876          0.683992
  SNP    ALL        59052     37825     21227        43480      1913       3740    675     35       0.640537          0.951862
 METRIC.Frac_NA  METRIC.F1_Score  TRUTH.TOTAL.TiTv_ratio  QUERY.TOTAL.TiTv_ratio  TRUTH.TOTAL.het_hom_ratio  QUERY.TOTAL.het_hom_ratio
       0.274864         0.559171                     NaN                     NaN                   1.547733                   2.756151
       0.086017         0.765767                2.433271                2.350281                   1.575230                   1.492346
****** DONE Vérifier que la zone de capture a vraiment besoin d'un liftover...
CLOSED: [2023-04-02 Sun 14:13]
D'après la documentation, oui
https://support.illumina.com/sequencing/sequencing_kits/nextera-rapid-captur

[3.16132]

[6.51603]

  922       490       432          942       439          0     32     56       0.531453          0.533970
INDEL   PASS          922       490       432          942       439          0     32     56       0.531453          0.533970
  SNP    ALL        24618     18689      5929        21257      2571          0    239     17       0.759160          0.879052
  SNP   PASS        24618     18689      5929        21257      2571          0    239     17       0.759160          0.879052
 METRIC.Frac_NA  METRIC.F1_Score  TRUTH.TOTAL.TiTv_ratio  QUERY.TOTAL.TiTv_ratio  TRUTH.TOTAL.het_hom_ratio  QUERY.TOTAL.het_hom_ratio
            0.0         0.532709                     NaN                     NaN                   1.628743                   2.799180
            0.0         0.532709                     NaN                     NaN                   1.628743                   2.799180
            0.0         0.814719                2.899604                2.881526                   1.598479                   1.564378
            0.0         0.814719                2.899604                2.881526                   1.598479                   1.564378
******* KILL Comprendre pour les FN explosent avec vep annot
CLOSED: [2023-02-24 Fri 23:44]
NC_000001.11	924024	.	C	G
non couverte dans bam -> 1er exons de SAMD11 mais non couverte
Il est dans NM_001385641.1 mais pas dans NM_152486.4
****** KILL Obtenir les zones couvertes depuis le bam directement
CLOSED: [2023-02-24 Fri 23:44]
#+begin_src  sh
samtools sort NA12878_chr1.bam -o NA12878_chr1_sorted.bam
bedtools genomecov -ibam NA12878_chr1_sorted.bam  -bg > chr1.bedgraph
head chr1.bedgraph
#+end_src
#+RESULTS:
: NC_000001.11	10001	10021	1
: NC_000001.11	10021	10059	2
: NC_000001.11	10059	10063	1
: NC_000001.11	10063	10087	2
: NC_000001.11	10087	10110	1
: NC_000001.11	10110	10113	3
: NC_000001.11	10113	10121	2
On prend les régions avec 20 reads
awk '$4 > 20 {print $1"\t"$2"\t"$3}' chr1.bedgraph  > tomerge.bed
On fusionne les régions
bedtools merge -i tomerge.bed > exons_from_bam.bed
Problème : on manque parfois le bord des exons...
****** KILL Vérifier un bam de référence de NA12878
CLOSED: [2023-02-24 Fri 23:44]
Notamment pour la définitions des exons, on a des reads qui ne semblent pas alignés sur les bons exons selon IGV
On donne directemet à IGV le bam d'illuma WES ( https://github.com/genome-in-a-bottle/giab_data_indexes )
IGV n'aime pas le ftp apparement...
On récupére 2 échantillons pour ce patient sur HiSeq Exome
#+begin_src sh :dir /ssh:meso:/Work/Groups/bisonex/data/giab/GRCh38
wget https://raw.githubusercontent.com/genome-in-a-bottle/giab_data_indexes/master/NA12878/alignment.index.NA12878_HiSeq_Exome_Garvan_GRCh37_09252015 -O todl.txt
awk '{print $1}' todl.txt | wget -i -
#+end_src
On récupére le premier échantillons en prenant just le chromosome 1
#+begin_src
samtools merge NA12878-NIST7035-HiSeq_Exome_Garan_GRCh37.bam project
.NIST_NIST7035_H7AP8ADXX_TAAGGCGA_1_NA12878.bwa.markDuplicates.bam project.NIST_NIST7035_H7AP8ADXX_TA
AGGCGA_2_NA12878.bwa.markDuplicates.bam
samtools index NA12878-NIST7035-HiSeq_Exome_Garan_GRCh37.bam
#+end_src
******* Étude
NC_000001.11:924024 :
- bed d’illumina : pas correct
- bed refseq : meux
- bam :  exon non ciblé, probablement parce que Mane select n’a pas été utilisé mais refseq (https://www.genenames.org/data/gene-symbol-report/#!/hgnc_id/HGNC:28706)
NC_000001.11:924310: idem
NC_000001.11:924321: idem
NC_000001.11:924533: idem
NC_000001.11:966227: le bed ne couvre pas le bon exon !
https://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr1%3A965863%2D966591&hgsid=295298291_vhDbIv5YIckCKLpGB7UUaZi2cAkr
NC_000001.11:966272
****** KILL Filtrer variants introniques de référence avec vep
CLOSED: [2023-02-24 Fri 23:44]
******* KILL variant calling seulf + seulement -f: nombreux FP
CLOSED: [2023-02-24 Fri 23:44]
/Work/Users/apraga/bisonex/work/68/f1cf72a5a4078fdf743fb3844b369a
 Type Filter  TRUTH.TOTAL  TRUTH.TP  TRUTH.FN  QUERY.TOTAL  QUERY.FP  QUERY.UNK  FP.gt  FP.al  METRIC.Recall  METRIC.Precision
INDEL    ALL          519       284       235        52169     37789      14094     12     64       0.547206          0.007511
INDEL   PASS          519       284       235        52169     37789      14094     12     64       0.547206          0.007511
  SNP    ALL        22131     17434      4697       357652    305313      34904    189     32       0.787764          0.054020
  SNP   PASS        22131     17434      4697       357652    305313      34904    189     32       0.787764          0.054020
 METRIC.Frac_NA  METRIC.F1_Score  TRUTH.TOTAL.TiTv_ratio  QUERY.TOTAL.TiTv_ratio  TRUTH.TOTAL.het_hom_ratio  QUERY.TOTAL.het_hom_ratio
       0.270160         0.014820                     NaN                     NaN                   1.775956                   0.509510
       0.270160         0.014820                     NaN                     NaN                   1.775956                   0.509510
       0.097592         0.101108                2.971834                1.429533                   1.579776                   0.324923
       0.097592         0.101108                2.971834                1.429533                   1.579776                   0.324923
****** DONE Bed donnés par GIAB (en hg19) : résultats corrects :resultats:
CLOSED: [2023-03-09 Thu 22:42] SCHEDULED: <2023-03-02 Thu>
Téléchargé à la main et converti via UCSCS. À automatiser, voir : *hg38
Sans oublier de renommer les chromosomes !
| Type  | TRUTH.TOTAL | TRUTH.TP | TRUTH.FN | QUERY.TOTAL | QUERY.FP | QUERY.UNK | FP.gt | FP.al | METRIC.Recall | METRIC.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 |
 Type   Ssensibilité    VPP
 INDEL       0.710532   0.692946 
 SNP         0.855253   0.970759
 | METRIC.Frac_NA | METRIC.F1_Score | TRUTH.TOTAL.TiTv_ratio | QUERY.TOTAL.TiTv_ratio | TRUTH.TOTAL.het_hom_ratio | QUERY.TOTAL.het_hom_ratio |
 |----------------+-----------------+------------------------+------------------------+---------------------------+---------------------------|
 |       0.281924 |        0.701629 |                        |                        |        1.6174985978687606 |        3.0674091441969518 |
 |       0.090605 |        0.909353 |      2.529551552318896 |     2.4019675131548843 |        1.6206857273037931 |        1.6274012964054214 |
****** DONE Re-tester avec exons Refseq et option -T: résultats un peu moins bons
CLOSED: [2023-03-09 Thu 22:42] SCHEDULED: <2023-03-08 Wed>
On utilise directement les coordonées données par refseq
 Type Filter  TRUTH.TOTAL  TRUTH.TP  TRUTH.FN  QUERY.TOTAL  QUERY.FP  QUERY.UNK  FP.gt  FP.al  METRIC.Recall  METRIC.Precision
INDEL    ALL         7226      3417      3809         6978      1599       1918    228    353       0.472876          0.683992
  SNP    ALL        59052     37825     21227        43480      1913       3740    675     35       0.640537          0.951862
 METRIC.Frac_NA  METRIC.F1_Score  TRUTH.TOTAL.TiTv_ratio  QUERY.TOTAL.TiTv_ratio  TRUTH.TOTAL.het_hom_ratio  QUERY.TOTAL.het_hom_ratio
       0.274864         0.559171                     NaN                     NaN                   1.547733                   2.756151
       0.086017         0.765767                2.433271                2.350281                   1.575230                   1.492346
****** DONE Vérifier que la zone de capture a vraiment besoin d'un liftover...
CLOSED: [2023-04-02 Sun 14:13]
D'après la documentation, oui
https://support.illumina.com/sequencing/sequencing_kits/nextera-rapid-captur