ence est TAAAA donc notre TAA -> TA aurait du être reconnu comme identique. Faut-il utiliser un autre fasta ?
La décomposition de xcmp est
45515163 TAA ->
T
45515164
AA -> A
NB: testée avec pre.py
#+begin_src
grep '^#' NA12878_NIST7035_DP_over_30.vcf > NA12878_NIST7035_DP_over_30_chr21.vcf
grep '^NC_000021' NA12878_NIST7035_DP_over_30.vcf >> NA12878_NIST7035_DP_over_30_chr21.vcf
pre.py NA12878_NIST7035_DP_over_30_chr21.vcf out.vcf.gz -r GCA_000001405.15_GRCh38_no_alt_analysis_set.fasta
#+end_src
Selon l'algorithm, on devrait avoir TA -> T et non AA -> A (non left-align)
Si on débug:
#+begin_src
grep "^#" NA12878_NIST7035_DP_over_30_chr21.vcf > test_align.vcf
grep 45515163 NA12878_NIST7035_DP_over_30_chr21.vcf >> test_align.vcf.gz
bgzip test_align.vcf.gz
bcftools index test_align.vcf.gz
multimerge test_align.vcf.gz -o out.vcf.gz -r GCA_000001405.15_GRCh38_no_alt_analysis_set.fasta
#+end_src
Idem... Si on utile un genome de référence plus récent ?
#+begin_src
multimerge test_align.vcf.gz -o out.vcf.gz -r /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.fna --process-full=1
#+end_src
Idem. Et vcfeval ?
#+begin_src
tabix HG001_GRCh38_1_22_v4.2.1_benchmark_chr21.vcf.gz
tabix test_align.vcf.gz
rtg vcfeval -b HG001_GRCh38_1_22_v4.2.1_benchmark_chr21.vcf.gz -c test_align.vcf.gz -t /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.sdf -o vcfeval-test
#+end_src
Selected score threshold using: maximized F-measure
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
99.000 0 0 1 54828 0.0000 0.0000 0.0000
None 0 0 1 54828 0.0000 0.0000 0.0000
On essaie les différentes étapes
#+begin_src
multimerge test_align.vcf.gz -o out.vcf.gz -r /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.fna --homref-split 1 --homref-vcf-out 1 --trimalleles 1 --splitalleles 1
#+end_src
1er changement
TAA T
TAA TA
Après réflexion, c'est la référence qui n'est pas normalisée ! (left-trimmed)
******** DONE NC_000021.9 14108836 T -> C
CLOSED: [2023-03-04 Sat 11:01]
Dans le bam mais filtré par DP
******* DONE variant calling seul : meilleur score pour l'instant (77% recall, 95% precision)
CLOSED: [2023-03-04 Sat 11:01]
tests/chr21-alexis
Type Filter TRUTH.TOTAL TRUTH.TP TRUTH.FN QUERY.TOTAL QUERY.FP QUERY.UNK FP.gt FP.al METRIC.Recall METRIC.Precision
INDEL ALL 76 43 33 82 18 20 3 5 0.565789 0.709677
SNP ALL 582 448 134 530 25 57 6 1 0.769759 0.947146
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.243902 0.629617 NaN NaN 1.181818 1.612903
0.107547 0.849289 3.098592 2.925926 1.530435 1.774869
******** NC_000021.9:14144627 : FN, dans le bam mais pas dans le vcf (2 reads/9)
On récupére tout le vcf: pas dedans
Dans le bam : 9/2
Idem pour le bam dans notre pipeline
- base qualité 33 sur les 2 reads
https://gatk.broadinstitute.org/hc/en-us/articles/360043491652-When-HaplotypeCaller-and-Mutect2-do-not-call-an-expected-variant
https://gatk.broadinstitute.org/hc/en-us/articles/360035891111-Expected-variant-at-a-specific-site-was-not-called
On debug
#+begin_src
cd /Work/Users/apraga/bisonex/out/NA12878_NIST7035/preprocessing/applybqsr
samtools view -b NA12878_NIST7035.bam NC_000021.9 -o NA12878_NIST7035_chr21.bam
samtools index NA12878_NIST7035_chr21.bam
#+end_src
********* --debug
#+begin_src
gatk --java-options "-Xmx3g" HaplotypeCaller --input NA12878_NIST7035_chr21.bam \
--output debug_chr1.vcf.gz \
--reference /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.fna \
--dbsnp /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP.gz \
--tmp-dir . \
--max-mnp-distance 2 --debug &> lol.txt
#+end_src
Les allèles sont bien retrouvées
#+begin_quote
14:41:05.530 INFO EventMap - === Best Haplotypes ===
14:41:05.530 INFO EventMap - AATTTAATTTTCTTACCTTTCTGGGTATGTAAGTGATTTTA
14:41:05.530 INFO EventMap - > Cigar = 41M
14:41:05.530 INFO EventMap - >> Events = EventMap{}
14:41:05.530 INFO EventMap - AATTTAATTTTCTTACCTTTTTGGGTATGTAAGTGATTTTA
14:41:05.530 INFO EventMap - > Cigar = 41M
14:41:05.530 INFO EventMap - >> Events = EventMap{NC_000021.9:14144627-14144627 [C*, T],}
14:41:05.530 INFO HaplotypeCallerGenotypingEngine - Genotyping event at 14144627 with alleles = [C*, T]
#+end_quote
NB: même en filtranrt sur le chromosome 21 avec julia, haplotypecaller parcourt tous les chromosomes
********* DONE --linked-de-bruijn-graph : idem
CLOSED: [2023-02-26 Sun 17:26]
********* DONE examine sortie --bamout : non présent
CLOSED: [2023-02-26 Sun 19:53]
#+begin_src
cd test/chr21-alexis
gatk --java-options "-Xmx3g" HaplotypeCaller \
--input /Work/Users/apraga/bisonex/script/files/bam/NA12878_chr21.bam \
--output debug_chr1.vcf.gz \
--reference /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.fna \
--dbsnp /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP.gz \
--tmp-dir . \
--max-mnp-distance 2 -bamout debug.bam
#+end_src
Pas de reads
#+begin_quote
If you see nothing overlapping your region, then it might not have been flagged as active, or could have failed to assemble.
#+end_quote
********* 14582339: FN mais pas de reads...
********* 14583327 idem
********* 17512551 idem
********* 17567111: difference d'haplotype
********* 17567621 pas de reads
****** DONE Comparer avec sortie du variant calling vcf donné par GIAB
CLOSED: [2023-04-02 Sun 17:11]
******* DONE vcfeval
CLOSED: [2023-04-01 Sat 11:59] SCHEDULED: <2023-04-01 Sat>
#+begin_src sh
nextflow run workflows/test.nf -profile standard,helios -resume --test.vcfeval --test.giabVCF --outdir=test-giabVCF
cat test-giabVCF/vcfeval/output/summary.txt
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
1.000 44818 44818 2892 6087 0.9394 0.8804 0.9089
None 44819 44819 2896 6086 0.9393 0.8804 0.9089
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
1.000 44818 44818 2892 6087 0.9394 0.8804 0.9089
None 44819 44819 2896 6086 0.9393 0.8804 0.9089
******* DONE happy
CLOSED: [2023-04-01 Sat 11:56]
Type Filter TRUTH.TOTAL TRUTH.TP TRUTH.FN QUERY.TOTAL QUERY.FP QUERY.UNK FP.gt FP.al METRIC.Recall METRIC.PrecisioN
INDEL PASS 4871 3678 1193 7036 1299 2011 208 217 0.755081 0.741493
SNP PASS 46032 41138 4894 47694 1622 4930 362 31 0.893683 0.962071
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.285816 0.748225 NaN NaN 1.617499 2.524051
0.103367 0.926617 2.529552 2.412446 1.620686 1.688868
****** DONE Statistiques avec vcfeval
CLOSED: [2023-04-02 Sun 17:10] SCHEDULED: <2023-04-01 Sat>
**** TODO HG002 :hg002:
SCHEDULED: <2023-04-10 Mon>
#+begin_src
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/compareVCF.nf -profile standard,helios -resume --outdir=compareHG002 --test.compare=vcfeval --test.query=out/HG002/variantCalling/haplotypecaller/HG002.vcf.gz --test.id=HG002
#+end_src
***** Mauvais résultats
avec vcfeval
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
0.000 24585 24390 10060 39415 0.7080 0.3841 0.4980
None 24585 24390 10060 39415 0.7080 0.3841 0.4980
La sortie du variantCalling est celle d'happy ???
On relance...
***** DONE Vérifier vcf en hg38
CLOSED: [2023-04-12 Wed 10:33] SCHEDULED: <2023-04-12 Wed>
***** TODO Capture en hg19 ?
SCHEDULED: <2023-04-12 Wed>
***** TODO Vraiment fichier de capture ou zone d'intérêt ?
SCHEDULED: <2023-04-12 Wed>
"target region" +/- 50bp
[[https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/README.txt][README]]
list file describing the variant calling regions (target regions extended with 50 bp on each end)
***** TODO Réessayer avec .bed fourni par AGilent
SCHEDULED: <2023-04-13 Thu>
Agilent SureSelect Human All Exon V5 kit
Disponible en hg38
***** TODO Bam -> fastq : impact du génome d'alignement ?
***** TODO Impact du joint calling ?
**** DONE Résumer résultats pour Paul + article :resultats:
CLOSED: [2023-04-06 Thu 21:41] SCHEDULED: <2023-04-02 Sun>
***** DONE HG001 :
CLOSED: [2023-04-06 Thu 21:41] SCHEDULED: <2023-04-02 Sun>
****** Données brutes
Version GIAB avec hap.py + vcfeval:
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/compareVCF.nf -profile standard,helios -resume --outdir=compareNA12878-giab --test.compare=happy,vcfeval --test.query=giab --test.id=HG001
#+end_src
Notre version avec hap.py + vcfeval
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/compareVCF.nf -profile standard,helios -resume --outdir=compareNA12878 --test.vcfeval --test.query="out/NA12878_NIST/variantCalling/haplotypecaller/NA12878_NIST.vcf.gz" --test.happy
#+end_src
On concatene les csv avec une colonne indicant le type
# awk '{if (NR==1) {print "Data,Algorithm" $0} else {print "bisonx,happy,"$0}}' compareNA12878/happy/NA12878.summary.csv
compareNA12878/happy/NA12878.summary.csv
| 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 | 4871 | 3461 | 1410 | 7048 | 1554 | 1987 | 193 | 346 | 0.710532 | 0.692946 | 0.281924 | 0.701629 | | | 1.6174985978687606 | 3.0674091441969518 |
| INDEL | PASS | 4871 | 3461 | 1410 | 7048 | 1554 | 1987 | 193 | 346 | 0.710532 | 0.692946 | 0.281924 | 0.701629 | | | 1.6174985978687606 | 3.0674091441969518 |
| SNP | ALL | 46032 | 39367 | 6665 | 44599 | 1186 | 4042 | 304 | 30 | 0.855209 | 0.970757 | 0.09063 | 0.909327 | 2.529551552318896 | 2.402150701647346 | 1.6206857273037931 | 1.6273423688862698 |
| SNP | PASS | 46032 | 39367 | 6665 | 44599 | 1186 | 4042 | 304 | 30 | 0.855209 | 0.970757 | 0.09063 | 0.909327 | 2.529551552318896 | 2.402150701647346 | 1.6206857273037931 | 1.6273423688862698 |
compareNA12878/vcfeval/NA12878.summary.txt
| Threshold | True-pos-baseline | True-pos-call | False-pos | False-neg | Precision | Sensitivity | F-measure |
|-----------+-------------------+---------------+-----------+-----------+-----------+-------------+-----------|
| 3.000 | 42789 | 42416 | 2598 | 8080 | 0.9423 | 0.8412 | 0.8889 |
| None | 42798 | 42425 | 2616 | 8071 | 0.9419 | 0.8413 | 0.8888 |
Indel avec le plus petit seuil : zcat NA12878.non_snp_roc.tsv.gz
Attention à inverser precision et recall !
zcat NA12878.non_snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
0.71390.7136
SNP avec le plus petit seuil : zcat NA12878.non_snp_roc.tsv.gz
Attention à inverser precision et recall !
$ zcat NA12878.snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
0.85470.9727
compareNA12878-giab/vcfeval/NA12878.summary.txt
| Threshold | True-pos-baseline | True-pos-call | False-pos | False-neg | Precision | Sensitivity | F-measure |
| 1.000 | 44812 | 44812 | 2878 | 6057 | 0.9397 | 0.8809 | 0.9093 |
| None | 44813 | 44813 | 2882 | 6056 | 0.9396 | 0.8809 | 0.9093 |
SNP:
$ zcat NA12878.snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
0.89370.9621
indel
$ zcat NA12878.non_snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
0.75980.7445
compareNA12878-giab/happy/NA12878.summary.csv
| 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 | 4871 | 3678 | 1193 | 7036 | 1299 | 2011 | 208 | 217 | 0.755081 | 0.741493 | 0.285816 | 0.748225 | | | 1.6174985978687606 | 2.5240506329113925 |
| INDEL | PASS | 4871 | 3678 | 1193 | 7036 | 1299 | 2011 | 208 | 217 | 0.755081 | 0.741493 | 0.285816 | 0.748225 | | | 1.6174985978687606 | 2.5240506329113925 |
| SNP | ALL | 46032 | 41138 | 4894 | 47694 | 1622 | 4930 | 362 | 31 | 0.893683 | 0.962071 | 0.103367 | 0.926617 | 2.529551552318896 | 2.4124463519313304 | 1.6206857273037931 | 1.6888675840288743 |
| SNP | PASS | 46032 | 41138 | 4894 | 47694 | 1622 | 4930 | 362 | 31 | 0.893683 | 0.962071 | 0.103367 | 0.926617 | 2.529551552318896 | 2.4124463519313304 | 1.6206857273037931 | 1.688867584028874 |
****** Résumé
| Données | Algorithm | Type | Recall | Precision |
|---------+-----------+---------+--------+-----------|
| Bisonex | Happy | SNP | 0.8552 | 0.9708 |
| Bisonex | vcfeval | SNP | 0.8547 | 0.9727 |
| Bisonex | Happy | INDEL | 0.7105 | 0.6929 |
| Bisonex | vcfeval | Non-SNP | 0.7139 | 0.7136 |
|---------+-----------+---------+--------+-----------|
| GIAB | happy | INDEL | 0.7551 | 0.7415 |
| GIAB | vcfeval | INDEL | 0.7598 | 0.7445 |
| GIAB | happy | SNP | 0.8937 | 0.9621 |
| giab | vcfeval | SNP | 0.8937 | 0.9621 |
*** TODO Platinum genome
https://emea.illumina.com/platinumgenomes.html
*** TODO Séquencer NA12878
Discussion avec Paul : sous-traitant ne nous donnera pas les données, il faut commander l'ADN
** TODO Fastq avec tous les variants centogène
*** TODO Extraire liste des variants
*** TODO Générer fastq
*** TODO Vérifier qu
ence est TAAAA donc notre TAA -> TA aurait du être reconnu comme identique. Faut-il utiliser un autre fasta ?
La décomposition de xcmp est
45515163 TAA -> T
45515164 AA -> A
NB: testée avec pre.py
#+begin_src
grep '^#' NA12878_NIST7035_DP_over_30.vcf > NA12878_NIST7035_DP_over_30_chr21.vcf
grep '^NC_000021' NA12878_NIST7035_DP_over_30.vcf >> NA12878_NIST7035_DP_over_30_chr21.vcf
pre.py NA12878_NIST7035_DP_over_30_chr21.vcf out.vcf.gz -r GCA_000001405.15_GRCh38_no_alt_analysis_set.fasta
#+end_src
Selon l'algorithm, on devrait avoir TA -> T et non AA -> A (non left-align)
Si on débug:
#+begin_src
grep "^#" NA12878_NIST7035_DP_over_30_chr21.vcf > test_align.vcf
grep 45515163 NA12878_NIST7035_DP_over_30_chr21.vcf >> test_align.vcf.gz
bgzip test_align.vcf.gz
bcftools index test_align.vcf.gz
multimerge test_align.vcf.gz -o out.vcf.gz -r GCA_000001405.15_GRCh38_no_alt_analysis_set.fasta
#+end_src
Idem... Si on utile un genome de référence plus récent ?
#+begin_src
multimerge test_align.vcf.gz -o out.vcf.gz -r /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.fna --process-full=1
#+end_src
Idem. Et vcfeval ?
#+begin_src
tabix HG001_GRCh38_1_22_v4.2.1_benchmark_chr21.vcf.gz
tabix test_align.vcf.gz
rtg vcfeval -b HG001_GRCh38_1_22_v4.2.1_benchmark_chr21.vcf.gz -c test_align.vcf.gz -t /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.sdf -o vcfeval-test
#+end_src
Selected score threshold using: maximized F-measure
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
99.000 0 0 1 54828 0.0000 0.0000 0.0000
None 0 0 1 54828 0.0000 0.0000 0.0000
On essaie les différentes étapes
#+begin_src
multimerge test_align.vcf.gz -o out.vcf.gz -r /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.fna --homref-split 1 --homref-vcf-out 1 --trimalleles 1 --splitalleles 1
#+end_src
1er changement
TAA T
TAA TA
Après réflexion, c'est la référence qui n'est pas normalisée ! (left-trimmed)
******** DONE NC_000021.9 14108836 T -> C
CLOSED: [2023-03-04 Sat 11:01]
Dans le bam mais filtré par DP
******* DONE variant calling seul : meilleur score pour l'instant (77% recall, 95% precision)
CLOSED: [2023-03-04 Sat 11:01]
tests/chr21-alexis
Type Filter TRUTH.TOTAL TRUTH.TP TRUTH.FN QUERY.TOTAL QUERY.FP QUERY.UNK FP.gt FP.al METRIC.Recall METRIC.Precision
INDEL ALL 76 43 33 82 18 20 3 5 0.565789 0.709677
SNP ALL 582 448 134 530 25 57 6 1 0.769759 0.947146
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.243902 0.629617 NaN NaN 1.181818 1.612903
0.107547 0.849289 3.098592 2.925926 1.530435 1.774869
******** NC_000021.9:14144627 : FN, dans le bam mais pas dans le vcf (2 reads/9)
On récupére tout le vcf: pas dedans
Dans le bam : 9/2
Idem pour le bam dans notre pipeline
- base qualité 33 sur les 2 reads
https://gatk.broadinstitute.org/hc/en-us/articles/360043491652-When-HaplotypeCaller-and-Mutect2-do-not-call-an-expected-variant
https://gatk.broadinstitute.org/hc/en-us/articles/360035891111-Expected-variant-at-a-specific-site-was-not-called
On debug
#+begin_src
cd /Work/Users/apraga/bisonex/out/NA12878_NIST7035/preprocessing/applybqsr
samtools view -b NA12878_NIST7035.bam NC_000021.9 -o NA12878_NIST7035_chr21.bam
samtools index NA12878_NIST7035_chr21.bam
#+end_src
********* --debug
#+begin_src
gatk --java-options "-Xmx3g" HaplotypeCaller --input NA12878_NIST7035_chr21.bam \
--output debug_chr1.vcf.gz \
--reference /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.fna \
--dbsnp /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP.gz \
--tmp-dir . \
--max-mnp-distance 2 --debug &> lol.txt
#+end_src
Les allèles sont bien retrouvées
#+begin_quote
14:41:05.530 INFO EventMap - === Best Haplotypes ===
14:41:05.530 INFO EventMap - AATTTAATTTTCTTACCTTTCTGGGTATGTAAGTGATTTTA
14:41:05.530 INFO EventMap - > Cigar = 41M
14:41:05.530 INFO EventMap - >> Events = EventMap{}
14:41:05.530 INFO EventMap - AATTTAATTTTCTTACCTTTTTGGGTATGTAAGTGATTTTA
14:41:05.530 INFO EventMap - > Cigar = 41M
14:41:05.530 INFO EventMap - >> Events = EventMap{NC_000021.9:14144627-14144627 [C*, T],}
14:41:05.530 INFO HaplotypeCallerGenotypingEngine - Genotyping event at 14144627 with alleles = [C*, T]
#+end_quote
NB: même en filtranrt sur le chromosome 21 avec julia, haplotypecaller parcourt tous les chromosomes
********* DONE --linked-de-bruijn-graph : idem
CLOSED: [2023-02-26 Sun 17:26]
********* DONE examine sortie --bamout : non présent
CLOSED: [2023-02-26 Sun 19:53]
#+begin_src
cd test/chr21-alexis
gatk --java-options "-Xmx3g" HaplotypeCaller \
--input /Work/Users/apraga/bisonex/script/files/bam/NA12878_chr21.bam \
--output debug_chr1.vcf.gz \
--reference /Work/Groups/bisonex/data/genome/GRCh38.p13/genomeRef.fna \
--dbsnp /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP.gz \
--tmp-dir . \
--max-mnp-distance 2 -bamout debug.bam
#+end_src
Pas de reads
#+begin_quote
If you see nothing overlapping your region, then it might not have been flagged as active, or could have failed to assemble.
#+end_quote
********* 14582339: FN mais pas de reads...
********* 14583327 idem
********* 17512551 idem
********* 17567111: difference d'haplotype
********* 17567621 pas de reads
****** DONE Comparer avec sortie du variant calling vcf donné par GIAB
CLOSED: [2023-04-02 Sun 17:11]
******* DONE vcfeval
CLOSED: [2023-04-01 Sat 11:59] SCHEDULED: <2023-04-01 Sat>
#+begin_src sh
nextflow run workflows/test.nf -profile standard,helios -resume --test.vcfeval --test.giabVCF --outdir=test-giabVCF
cat test-giabVCF/vcfeval/output/summary.txt
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
1.000 44818 44818 2892 6087 0.9394 0.8804 0.9089
None 44819 44819 2896 6086 0.9393 0.8804 0.9089
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
1.000 44818 44818 2892 6087 0.9394 0.8804 0.9089
None 44819 44819 2896 6086 0.9393 0.8804 0.9089
******* DONE happy
CLOSED: [2023-04-01 Sat 11:56]
Type Filter TRUTH.TOTAL TRUTH.TP TRUTH.FN QUERY.TOTAL QUERY.FP QUERY.UNK FP.gt FP.al METRIC.Recall METRIC.PrecisioN
INDEL PASS 4871 3678 1193 7036 1299 2011 208 217 0.755081 0.741493
SNP PASS 46032 41138 4894 47694 1622 4930 362 31 0.893683 0.962071
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.285816 0.748225 NaN NaN 1.617499 2.524051
0.103367 0.926617 2.529552 2.412446 1.620686 1.688868
****** DONE Statistiques avec vcfeval
CLOSED: [2023-04-02 Sun 17:10] SCHEDULED: <2023-04-01 Sat>
***** DONE Résultats finaux
CLOSED: [2023-04-14 Fri 09:53]
Version GIAB avec hap.py + vcfeval:
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/compareVCF.nf -profile standard,helios -resume --outdir=compareNA12878-giab --test.compare=happy,vcfeval --test.query=giab --test.id=HG001
#+end_src
Notre version avec hap.py + vcfeval
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/compareVCF.nf -profile standard,helios -resume --outdir=compareNA12878 --test.vcfeval --test.query="out/NA12878_NIST/variantCalling/haplotypecaller/NA12878_NIST.vcf.gz" --test.happy
#+end_src
On concatene les csv avec une colonne indicant le type
# awk '{if (NR==1) {print "Data,Algorithm" $0} else {print "bisonx,happy,"$0}}' compareNA12878/happy/NA12878.summary.csv
compareNA12878/happy/NA12878.summary.csv
| 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 | 4871 | 3461 | 1410 | 7048 | 1554 | 1987 | 193 | 346 | 0.710532 | 0.692946 | 0.281924 | 0.701629 | | | 1.6174985978687606 | 3.0674091441969518 |
| INDEL | PASS | 4871 | 3461 | 1410 | 7048 | 1554 | 1987 | 193 | 346 | 0.710532 | 0.692946 | 0.281924 | 0.701629 | | | 1.6174985978687606 | 3.0674091441969518 |
| SNP | ALL | 46032 | 39367 | 6665 | 44599 | 1186 | 4042 | 304 | 30 | 0.855209 | 0.970757 | 0.09063 | 0.909327 | 2.529551552318896 | 2.402150701647346 | 1.6206857273037931 | 1.6273423688862698 |
| SNP | PASS | 46032 | 39367 | 6665 | 44599 | 1186 | 4042 | 304 | 30 | 0.855209 | 0.970757 | 0.09063 | 0.909327 | 2.529551552318896 | 2.402150701647346 | 1.6206857273037931 | 1.6273423688862698 |
compareNA12878/vcfeval/NA12878.summary.txt
| Threshold | True-pos-baseline | True-pos-call | False-pos | False-neg | Precision | Sensitivity | F-measure |
|-----------+-------------------+---------------+-----------+-----------+-----------+-------------+-----------|
| 3.000 | 42789 | 42416 | 2598 | 8080 | 0.9423 | 0.8412 | 0.8889 |
| None | 42798 | 42425 | 2616 | 8071 | 0.9419 | 0.8413 | 0.8888 |
Indel avec le plus petit seuil : zcat NA12878.non_snp_roc.tsv.gz
Attention à inverser precision et recall !
zcat NA12878.non_snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
0.71390.7136
SNP avec le plus petit seuil : zcat NA12878.non_snp_roc.tsv.gz
Attention à inverser precision et recall !
$ zcat NA12878.snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
0.85470.9727
compareNA12878-giab/vcfeval/NA12878.summary.txt
| Threshold | True-pos-baseline | True-pos-call | False-pos | False-neg | Precision | Sensitivity | F-measure |
| 1.000 | 44812 | 44812 | 2878 | 6057 | 0.9397 | 0.8809 | 0.9093 |
| None | 44813 | 44813 | 2882 | 6056 | 0.9396 | 0.8809 | 0.9093 |
SNP:
$ zcat NA12878.snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
0.89370.9621
indel
$ zcat NA12878.non_snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
0.75980.7445
compareNA12878-giab/happy/NA12878.summary.csv
| 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 | 4871 | 3678 | 1193 | 7036 | 1299 | 2011 | 208 | 217 | 0.755081 | 0.741493 | 0.285816 | 0.748225 | | | 1.6174985978687606 | 2.5240506329113925 |
| INDEL | PASS | 4871 | 3678 | 1193 | 7036 | 1299 | 2011 | 208 | 217 | 0.755081 | 0.741493 | 0.285816 | 0.748225 | | | 1.6174985978687606 | 2.5240506329113925 |
| SNP | ALL | 46032 | 41138 | 4894 | 47694 | 1622 | 4930 | 362 | 31 | 0.893683 | 0.962071 | 0.103367 | 0.926617 | 2.529551552318896 | 2.4124463519313304 | 1.6206857273037931 | 1.6888675840288743 |
| SNP | PASS | 46032 | 41138 | 4894 | 47694 | 1622 | 4930 | 362 | 31 | 0.893683 | 0.962071 | 0.103367 | 0.926617 | 2.529551552318896 | 2.4124463519313304 | 1.6206857273037931 | 1.688867584028874 |
**** DONE HG002 :hg002:
CLOSED: [2023-04-14 Fri 09:54] SCHEDULED: <2023-04-10 Mon>
#+begin_src
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/giabFastq.nf -profile standard,helios
NXF_OPTS=-D"user.name=${USER}" nextflow run main.nf -profile standard,helios -resume --input="/Work/Groups/bisonex/data/giab/GRCh38/HG002_{1,2}.fq.gz --test.id=HG002
Only the capture file differs. Results are better using the capture file given by Agilent, stored in data/
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/compareVCF.nf -profile standard,helios -resume --outdir=compareHG002 --test.id=HG002 --test.query=out/HG002_1/variantCalling/haplotypecaller/HG002_1.vcf.gz --test.compare=vcfeval,happy --test.capture=data/AgilentSureSelectv05_hg38.bed
#
#+end_src
***** DONE Mauvais résultats
CLOSED: [2023-04-14 Fri 09:42]
avec vcfeval
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
0.000 24585 24390 10060 39415 0.7080 0.3841 0.4980
None 24585 24390 10060 39415 0.7080 0.3841 0.4980
La sortie du variantCalling est celle d'happy ???
On relance...
***** DONE Vérifier vcf en hg38
CLOSED: [2023-04-12 Wed 10:33] SCHEDULED: <2023-04-12 Wed>
***** KILL Capture en hg19 ?
CLOSED: [2023-04-13 Thu 09:46] SCHEDULED: <2023-04-12 Wed>
***** KILL Vraiment fichier de capture ou zone d'intérêt ?
CLOSED: [2023-04-13 Thu 09:45] SCHEDULED: <2023-04-12 Wed>
"target region" +/- 50bp
[[https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/README.txt][README]]
list file describing the variant calling regions (target regions extended with 50 bp on each end)
***** DONE .bed fourni par AGilent: sensbilité très mauvaise
CLOSED: [2023-04-13 Thu 09:46] SCHEDULED: <2023-04-13 Thu>
Agilent SureSelect Human All Exon V5 kit
Disponible en hg38
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
0.000 19653 19501 6410 21657 0.7526 0.4757 0.5830
None 19653 19501 6410 21657 0.7526 0.4757 0.5830
***** DONE Trier par nom avec samtools sort : bons résultats
CLOSED: [2023-04-14 Fri 09:25] SCHEDULED: <2023-04-13 Thu>
Avec capture fourni par GIAB
vcf eval
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
5.000 57443 57032 984 6557 0.9830 0.8975 0.9383
None 57457 57046 1009 6543 0.9826 0.8978 0.9383
Happy
| 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 | 6150 | 5007 | 1143 | 6978 | 556 | 1346 | 151 | 168 | 0.814146 | 0.901278 | 0.192892 | 0.8555 | | | 1.5434221840068787 | 1.9467178175618074 |
| INDEL | PASS | 6150 | 5007 | 1143 | 6978 | 556 | 1346 | 151 | 168 | 0.814146 | 0.901278 | 0.192892 | 0.8555 | | | 1.5434221840068787 | 1.9467178175618074 |
| SNP | ALL | 57818 | 52464 | 5354 | 56016 | 500 | 3046 | 90 | 30 | 0.907399 | 0.990561 | 0.054377 | 0.947158 | 2.4892012548262548 | 2.426824047458871 | 1.5904527117884357 | 1.6107795598657217 |
| SNP | PASS | 57818 | 52464 | 5354 | 56016 | 500 | 3046 | 90 | 30 | 0.907399 | 0.990561 | 0.054377 | 0.947158 | 2.4892012548262548 | 2.426824047458871 | 1.5904527117884357 | 1.6107795598657217 |
***** DONE Capture agilent légment meilleur que celui fourni par GIAB (padding ?)
CLOSED: [2023-04-14 Fri 09:48]
GIAB:
vcf eval
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
5.000 57443 57032 984 6557 0.9830 0.8975 0.9383
None 57457 57046 1009 6543 0.9826 0.8978 0.9383
Happy
| 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 | 6150 | 5007 | 1143 | 6978 | 556 | 1346 | 151 | 168 | 0.814146 | 0.901278 | 0.192892 | 0.8555 | | | 1.5434221840068787 | 1.9467178175618074 |
| INDEL | PASS | 6150 | 5007 | 1143 | 6978 | 556 | 1346 | 151 | 168 | 0.814146 | 0.901278 | 0.192892 | 0.8555 | | | 1.5434221840068787 | 1.9467178175618074 |
| SNP | ALL | 57818 | 52464 | 5354 | 56016 | 500 | 3046 | 90 | 30 | 0.907399 | 0.990561 | 0.054377 | 0.947158 | 2.4892012548262548 | 2.426824047458871 | 1.5904527117884357 | 1.6107795598657217 |
| SNP | PASS | 57818 | 52464 | 5354 | 56016 | 500 | 3046 | 90 | 30 | 0.907399 | 0.990561 | 0.054377 | 0.947158 | 2.4892012548262548 | 2.426824047458871 | 1.5904527117884357 | 1.6107795598657217 |
Agilent
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
6.000 37241 36965 449 4069 0.9880 0.9015 0.9428
None 37248 36972 461 4062 0.9877 0.9017 0.9427
| 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 | 2909 | 2477 | 432 | 3229 | 207 | 519 | 52 | 50 | 0.851495 | 0.923616 | 0.160731 | 0.886091 | | | 1.4964850615114236 | 1.8339222614840989 |
| INDEL | PASS | 2909 | 2477 | 432 | 3229 | 207 | 519 | 52 | 50 | 0.851495 | 0.923616 | 0.160731 | 0.886091 | | | 1.4964850615114236 | 1.8339222614840989 |
| SNP | ALL | 38406 | 34793 | 3613 | 36935 | 275 | 1868 | 37 | 15 | 0.905926 | 0.992158 | 0.050575 | 0.947083 | 2.6247759222568168 | 2.5752854654538417 | 1.588953331534934 | 1.6192536889897844 |
| SNP | PASS | 38406 | 34793 | 3613 | 36935 | 275 | 1868 | 37 | 15 | 0.905926 | 0.992158 | 0.050575 | 0.947083 | 2.6247759222568168 | 2.5752854654538417 | 1.588953331534934 | 1.6192536889897844 |
**** DONE Résumer résultats pour Paul + article :resultats:
CLOSED: [2023-04-06 Thu 21:41] SCHEDULED: <2023-04-02 Sun>
*** TODO Platinum genome
https://emea.illumina.com/platinumgenomes.html
*** TODO Séquencer NA12878
Discussion avec Paul : sous-traitant ne nous donnera pas les données, il faut commander l'ADN
** TODO Fastq avec tous les variants centogène
*** TODO Extraire liste des variants
*** TODO Générer fastq
*** TODO Vérifier qu