*** Best practices for the analytical validation of clinical whole-genome sequencing intended for the diagnosis of germline disease

*** NA12878
**** KILL [[https://precision.fda.gov/challenges/truth/results][fdaPrecision challenge]]
Attention, génome et en hg19 donc comparaison non adaptée ...
**** TODO Best practices for the analytical validation of clinical whole-genome sequencing intended for the diagnosis of germline disease
SCHEDULED: <2023-03-04 Sat>

Peu d’info sur la bionfo, sauf Supplementary :

1. variant calling seul
2. NA12878 + données simulées
3. exome
4. évalué via F-score
Résultat: BWA/Novoalign_DeepVariant
Aligneurs
- BWA-MEM 0.7.16
- Bowtie2 2.2.6
- Novoalign 3.08.02
- SOAP 2.21
- MOSAIK 2.2.3
Variantcalling
- GATK HaplotypeCaller 4
- FreeBayes 1.1.0
- SAMtools mpileup 1.7
- DeepVariant r0.4
  SNV
| Exome | Pipeline |    TP |   FP |  FN | Sensitivity | Precision | F-Score |   FDR |
|     1 | BWA_GATK | 23689 | 1397 | 613 |       0.975 |     0.944 |   0.959 | 0.057 |
|     2 | BWA_GATK | 23946 |  865 | 356 |       0.985 |     0.965 |   0.975 | 0.036 |
indel
 |   TP | FP | FN | Sensitivity | Precision | F-Score |   FDR |   |
 | 1254 | 72 | 75 |       0.944 |     0.946 |   0.945 | 0.054 |   |
 | 1309 | 10 | 20 |       0.985 |     0.992 |   0.989 | 0.008 |   |
Valeur brutes :
https://static-content.springer.com/esm/art%3A10.1186%2Fs12859-019-2928-9/MediaObjects/12859_2019_2928_MOESM8_ESM.pdf
Autres articles avec même comparaison en exome sur NA12878
- Hwang et al., 2015 studyi
- Highnam et al, 2015
-  Cornish and Guda, 2015

:PROPERTIES:
:CATEGORY: tests
:END:

*** TODO 63003856_S135

*** HOLD 63003856_S135

:PROPERTIES:
:CATEGORY: NA18278
:END:

:PROPERTIES:
:CATEGORY: GIAB
:END:

**** Comparaison

**** DONE Comparaison
CLOSED: [2023-03-04 Sat 11:14]

***** TODO Restreindre genome de référence
****** Discussion Alexis

**** DONE Restreindre genome de référence
CLOSED: [2023-03-04 Sat 11:15]
***** Discussion Alexis

****** DONE -f, -R ou -T ?

***** DONE -f, -R ou -T ?

******* rtg tools

****** rtg tools

****** TODO Utiliser seulement exons
******* DONE BestRefSeq 

**** TODO Exons seuls
***** DONE BestRefSeq 

******** DONE Tests après correction bug dans noms de chromosome : precision ~ ok, recall très mauvais -> trop de FN ?

****** DONE Tests après correction bug dans noms de chromosome : precision ~ ok, recall très mauvais -> trop de FN ?

******** DONE Vérifier exons: on a l'union des exons de tous les transcripts...

****** DONE Vérifier exons: on a l'union des exons de tous les transcripts...

******* KILL Bed illumina

***** KILL Bed illumina

******** KILL Sans filtre vep: Inversion truth et query... on recommence

****** KILL Sans filtre vep: Inversion truth et query... on recommence

******** KILL Sans filtre vep: mieux mais pas exceptionnel

****** KILL Sans filtre vep: mieux mais pas exceptionnel

******** KILL Comprendre pour les FN explosent avec vep annot

****** KILL Comprendre pour les FN explosent avec vep annot

******* KILL Obtenir les zones couvertes depuis le bam directement

***** KILL Obtenir les zones couvertes depuis le bam directement

******* KILL Vérifier un bam de référence de NA12878

***** KILL Vérifier un bam de référence de NA12878

******** Étude

****** Étude

******* KILL Filtrer variants introniques de référence avec vep

***** KILL Filtrer variants introniques de référence avec vep

******** KILL variant calling seulf + seulement -f: nombreux FP

****** KILL variant calling seulf + seulement -f: nombreux FP

******* TODO Avec les exons donnése par GIAB (en hg19) : bon résultats !!

***** TODO Bed donnés par GIAB (en hg19)
****** DONE résultats corrects
SCHEDULED: <2023-03-02 Thu>

***** TODO Comprendre pourquoi on est inférieur à Kumaran et al 2019
SCHEDULED: <2023-03-05 Sun>
****** TODO Comprende/améliorer Recall SNP 0.855
SCHEDULED: <2023-03-04 Sat>
******* TODO Regarder les FN (SNP)
SCHEDULED: <2023-03-04 Sat>
******** Manuel:
NC_000001.11:1385919 pas de read 1/1:FN:.:i1_5:INDEL:homalt:.
NC_000001.11:1623412 1 read     1/1:FN:.:ti:SNP:homalt:.
NC_000001.11:1668449 33 read sur 160 voient l'allèle alternative 1/1:FN:am:ti:SNP:homalt:.
NC_000001.11:1676135 67 reads, non vu 0/1:FN:.:ti:SNP:het:.
NC_000001.11:1734812 1/1:FN:.:ti:SNP:homalt:.
NC_000001.11:1745808 1/1:FN:.:ti:SNP:homalt:.
NC_000001.11:1745814 1/1:FN:.:ti:SNP:homalt:.
NC_000001.11:1953616 1/1:FN:.:ti:SNP:homalt:.
NC_000001.11:2512975 0/1:FN:.:ti:SNP:het:

******* TODO Statistiques
SCHEDULED: <2023-03-04 Sat>
NB: awk car bcftools query ne semble pas prendre en compte les nouvelles colonnes

****** TODO Comparer avec stats de NA12878 dans example/happy sur chr21 (exons fournis)
******* DONE DP_over_30_not_SNP_consensual_sequence.vcf: horrible

On confirnme le nombre de SNP:
❯ awk '!/^#/ && $10~/:FN:/ && $10~/SNP/' test-allchr.vcf  | wc -l
6665
Une minorité concerne des problème d'haploides
❯ awk '!/^#/ && $10~/:FN:/ && $11!~/NOCALL/ && $10~/SNP/' test-allchr.vcf
304
avec 1/3 où l'exome manque une allèles
❯ awk '!/^#/ && $10~/:FN:/ && $11!~/NOCALL/ && $10~/SNP/ && $10~/homalt/' test-allchr.vcf  | wc -l
101
et 2/3 où il y a une allèle "en trop"
La majorité ne sont pas vu
❯ awk '!/^#/ && $10~/:FN:/ && $11~/NOCALL/ && $10~/SNP/' test-allchr.vcf  | wc -l
6361
Nombre de reads pour chaque position en bash (!)
#+begin_src bash
 awk '!/^#/ && $10~/:FN:/ && $11~/NOCALL/ && $10~/SNP/ {print $1":"$2"-"$2}' test-allchr.vcf  |  xargs -I {} sh -c 'echo -n {}";"; samtools view ../NA12878_NIST.b
am {} | wc -l' > count.csv
#+end_src
On examine avec Julia
#+begin_src julia
using CSV, DataFrames
f = DataFrame(CSV.File("count.csv", header=false))
size(f[f.Column2 .< 10,:])
#+end_src
(4430, 2)
Est-on dans des exons ?
Ça a l'air sur des exemples simples ...
****** TODO Vérifier qu'il ne manque pas des exons (avec bam ?)
SCHEDULED: <2023-03-05 Sun>
***** KILL Comparer avec stats de NA12878 dans example/happy sur chr21 (exons fournis)
CLOSED: [2023-03-04 Sat 11:01]
****** DONE DP_over_30_not_SNP_consensual_sequence.vcf: horrible

******* DP _over_30 : mieux

****** DONE DP _over_30 : mieux
CLOSED: [2023-03-04 Sat 11:01]

******** STRT NC_000021.9:45515163 TA -> T : référence non normalisée ??

******* DONE NC_000021.9:45515163 TA -> T : référence non normalisée ??
CLOSED: [2023-03-04 Sat 11:01]

******** NC_000021.9     14108836  T -> C

******* DONE NC_000021.9     14108836  T -> C
CLOSED: [2023-03-04 Sat 11:01]

******* variant calling seul : meilleur score pour l'instant (77% recall, 95% precision)

****** DONE variant calling seul : meilleur score pour l'instant (77% recall, 95% precision)
CLOSED: [2023-03-04 Sat 11:01]

******** NC_000021.9:14144627 : FN, dans le bam mais pas dans le vcf (2 reads/9)

******* NC_000021.9:14144627 : FN, dans le bam mais pas dans le vcf (2 reads/9)

********* --debug

******** --debug

********* DONE --linked-de-bruijn-graph : idem

******** DONE --linked-de-bruijn-graph : idem

********* DONE examine sortie --bamout : non présent

******** DONE examine sortie --bamout : non présent

********* 14582339: FN mais pas de reads...
********* 14583327 idem
********* 17512551 idem
********* 17567111: difference d'haplotype
********* 17567621 pas de reads

******** 14582339: FN mais pas de reads...
******** 14583327 idem
******** 17512551 idem
******** 17567111: difference d'haplotype
******** 17567621 pas de reads