apraga/org - Change GGUC6BRBVBMIWZAYAYHNLL37QXWKTUUANB4PBMHIJIBJWBWP7JPAC

Bisonex : finally some results

Created by Alexis Praga on February 26, 2023

GGUC6BRBVBMIWZAYAYHNLL37QXWKTUUANB4PBMHIJIBJWBWP7JPAC

Dependencies

In channels

main

Change contents

Insertion in workout.org at line 5898 [6.1]

[5.34]

* <2023-02-26 Sun> Workout
** RTO:
20-10-10 (peu de motivation)
** L-sit
- 5
- 5
- 5
** Muscle-up
- 1-1-1 - 5 négatif
- 2-1-1 - 4 négatif
- 1-1-1 - 4 négatif
** Extension:
- 22
- 22
- 22
** FL tucked row :
- 3-2
- 3-2
- 3-2
** Pistols :
- 4
- 4
- 4
** Planche tucked push-up:
- 5-5
- 5-5
- 5-5
** Compression:
Jambes jointes en pensant à rapproche les genoux = difficile !
- 4-4
- 4-4
- 4-4
** Norwegian roll
- 4
- 4
- 4

Insertion in projects.org at line 136 [6.123895]

[7.431]

[8.8973]

* Julia
Mettre à jour juliastats
https://discourse.julialang.org/t/suggestion-for-homepage-of-https-juliastats-org/94948/2

Insertion in projects/bisonex.org at line 10 [10.35]

[3.140]

[9.176]

** Étapes du pipeline
*** Variant calling: Haplotype caller
https://gatk.broadinstitute.org/hc/en-us/articles/360035531412
Définis l'algorithme + image

Replacement in projects/bisonex.org at line 49 [10.35]

B:BD[3.1270] → [3.1270:1271]

[3.1270]

[11.14]

*** Débugger variant calling (haplotypecaller)
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

Insertion in projects/bisonex.org at line 985 [10.35]

[12.32132]

*** TODO Utilise AVX pour accélerer l'exécution
Sans cela, on a l'avertissement
#+begin_quote
17:28:00.720 INFO  PairHMM - OpenMP multi-threaded AVX-accelerated native PairHMM implementation is not supported
17:28:00.721 INFO  NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/nix/store/cy9ckxqwrkifx7wf02hm4ww1p6lnbxg9-gatk-4.2.4.1/bin/gatk-package-4.2.4.1-local.jar!/com/intel/gkl/native/libgkl_utils.so
17:28:00.733 WARN  NativeLibraryLoader - Unable to load libgkl_utils.so from native/libgkl_utils.so (/Work/Users/apraga/bisonex/out/NA12878_NIST7035/preprocessing/applybqsr/libgkl_utils821485189051585397.so: libgomp.so.1: cannot open shared object file: No such file or directory)
17:28:00.733 WARN  IntelPairHmm - Intel GKL Utils not loaded
17:28:00.733 WARN  PairHMM - ***WARNING: Machine does not have the AVX instruction set support needed for the accelerated AVX PairHmm. Falling back to the MUCH slower LOGLESS_CACHING implementation!
17:28:00.763 INFO  ProgressMeter - Starting traversal
#+end_quote
libgomp.so est fourni par gcc donc il faut charger le module
 module load gcc@11.3.0/gcc-12.1.0

Replacement in projects/bisonex.org at line 2825 [10.35]

B:BD[13.28] → [14.6268:6282]

B:BD[14.6282] → [4.15:80]

**** Nextflow
***** TODO Télécharger NA12878 (HG001)
****** TODO Fastq HiSeq

[13.28]

[4.80]

**** DONE télécharger données avec Nextflow
CLOSED: [2023-02-25 Sat 19:47]
***** DONE Télécharger NA12878 (HG001)
CLOSED: [2023-02-25 Sat 19:46]
****** DONE Fastq HiSeq
CLOSED: [2023-02-25 Sat 19:46]

Replacement in projects/bisonex.org at line 2834 [10.35]

B:BD[4.486] → [4.486:712]


#+begin_src
for j in 1 2; do
for i in 1 2; do wget "ftp://ftp-trace.ncbi.nih.gov/ReferenceSamples/giab/data/NA12878/Garvan_NA12878_HG001_HiSeq_Exome/NIST7035_TAAGGCGA_L00${j}_R${i}_001_trimmed.fastq.gz"; done; done
#+end_src

[4.486]

[4.712]

Informations:
- https://ftp-trace.ncbi.nih.gov/ReferenceSamples/giab/data/NA12878/Garvan_NA12878_HG001_HiSeq_Exome/Garvan_NA12878_HG001_HiSeq_Exome.README
- Sequencer: HiSeq2500
- kit: Nextera Rapid Capture Exome and Expanded Exome

Replacement in projects/bisonex.org at line 2839 [10.35]

B:BD[4.713] → [4.713:863]


Il faut concatener les 4 fastq en 2
#+begin_src
 zcat NIST7035_TAAGGCGA_L00*_R1_* | gzip -c > NIST7035_TAAGGCGA_R1_001_trimm
#+end_src
ed.fastq.gz

[4.713]

[4.863]

Il y a 2 samples (NIST7035 et NIST7086), chacun sur 2 lanes -> à concaténer

Replacement in projects/bisonex.org at line 2843 [10.35]
B:BD[4.980] → [4.980:1004]
```
****** TODO Exons (bed)
```
[4.980]
[4.1004]
```
****** DONE Exons (bed)
CLOSED: [2023-02-25 Sat 19:46]
```

Replacement in projects/bisonex.org at line 2851 [10.35]

B:BD[14.6498] → [14.6498:6587]

***** TODO Genome de reference NCBI
Téléchargé mais index non stocké dars /Work TODO

[14.6498]

[14.6587]

***** DONE Genome de reference NCBI
CLOSED: [2023-02-25 Sat 19:46]
***** TODO Bed avec les exons
****** DONE hg19
CLOSED: [2023-02-26 Sun 22:37]
****** TODO hg38
- [ ] Télécharger hg19 : ok
- [ ] convertir bed en interval list
picard BedToIntervalList -I exons_illumina.bed  -O exons_illumina.list -SD  ../../genome/GRCh19/genomeRef.dict
- [ ] puis en hg38
picard LiftOverIntervalList -I exons_illumina.list  -O exons_illumina_hg38.list --CHAIN hg19ToHg38.over.chain -SD  ../../genome/GRCh38.p13/genomeRef.dict
- [ ] puis en bed

Replacement in projects/bisonex.org at line 2938 [10.35]
B:BD[15.834] → [16.15:42]
```
****** TODO -f, -R ou -T ?
```
[15.834]
[16.42]
```
****** DONE -f, -R ou -T ?
CLOSED: [2023-02-25 Sat 19:47]
```

Replacement in projects/bisonex.org at line 3087 [10.35]

B:BD[4.1681] → [4.1681:1726]

∅:D[4.1726] → [5.180:995]

B:BD[17.2203] → [5.180:995]

******* TODO GIAB donne les exons en hg19 !!
****** TODO Comparer avec stats de NA12878 dans example/happy sur chr21 (exons fournis)
Résultats
| Type  | Filter | TRUTH.TOTAL | TRUTH.TP | TRUTH.FN | QUERY.TOTAL | QUERY.FP | QUERY.UNK | FP.gt | FP.al | METRIC.Recall | METRIC.Precision |
| INDEL | ALL    |        8530 |     8247 |      283 |       14509 |      231 |      5835 |    18 |   141 |      0.966823 |         0.973369 |
| INDEL | PASS   |        8530 |     8247 |      283 |       14483 |      231 |      5809 |    18 |   141 |      0.966823 |         0.973369 |
| SNP   | ALL    |       49913 |    49466 |      447 |       70641 |       69 |     21118 |    40 |     9 |      0.991044 |         0.998607 |
| SNP   | PASS   |       49913 |    49465 |      448 |       70551 |       69 |     21029 |    40 |     9 |      0.991024 |         0.998607 |

[4.1681]

[5.995]

******* TODO Avec les exons donnése par GIAB (en hg19) : bon résultats !!
Téléchargé à la main et converti via UCSCS. À automatiser, voir : *hg38
Sans oublier de renommer les chromosomes !

Replacement in projects/bisonex.org at line 3091 [10.35]

B:BD[5.996] → [5.996:1159]

Commande
#+begin_src
hap.py pg-hg38.bcf NA12878-GATK3-chr21.vcf.gz -f pg-hg38-conf.bed.gz -r hg38.chr21.fa -o GATK3-xcmp --roc QUAL --roc-filter LowQual
#+end_src

[5.996]

[18.666]

| 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 |

Replacement in projects/bisonex.org at line 3096 [10.35]

B:BD[16.114] → [5.1160:1338]

En filtrant les exons avec le bed fourni
#+begin_src
 hap.py pg-hg38.bcf NA12878-GATK3-chr21.vcf.gz -f pg-hg38-conf.bed.gz -r hg38.chr21.fa -T exons_hg38.bed.gz -o lol
#+end_src

[16.114]

[5.1338]

 | 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 |

Deletion in projects/bisonex.org at line 3101 [10.35]

B:BD[5.1339] → [5.1339:2110]

  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           97        96         1          138         2         37      0      0       0.989691          0.980198
INDEL   PASS           97        96         1          138         2         37      0      0       0.989691          0.980198
  SNP    ALL          638       636         2          674         0         38      0      0       0.996865          1.000000
  SNP   PASS          638       636         2          674         0         38      0      0       0.996865          1.000000

Replacement in projects/bisonex.org at line 3102 [10.35]

B:BD[16.378] → [19.15:103]

B:BD[19.103] → [20.15:76]

****** TODO Utiliser les filtre pour éliminer FP
On a déjà des DP à éliminer
reads
******* DP_over_30_not_SNP_consensual_sequence.vcf: horrible

[16.378]

[19.154]

****** TODO Comparer avec stats de NA12878 dans example/happy sur chr21 (exons fournis)
******* DONE DP_over_30_not_SNP_consensual_sequence.vcf: horrible
CLOSED: [2023-02-25 Sat 19:47]

Deletion in projects/bisonex.org at line 3111 [10.35]

B:BD[19.1476] → [4.1727:1826]

******* NC_000021.9     14144627
Dans le bam mais pas dans le vcf...
Cela vient du variant calling

Replacement in projects/bisonex.org at line 3176 [10.35]

B:BD[4.1944] → [4.1944:2728]

cd /Work/Users/apraga/bisonex/work/3a/ebf0249db81166c5b12f03cc8167b2
| 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 |
| INDEL | PASS   |          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 |
| SNP   | PASS   |         582 |      448 |      134 |         530 |       25 |        57 |     6 |     1 |      0.769759 |         0.947146 |

[4.1944]

[21.2042]

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
INDEL   PASS           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
  SNP   PASS          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.243902         0.629617                     NaN                     NaN                   1.181818                   1.612903
       0.107547         0.849289                3.098592                2.925926                   1.530435                   1.774869
       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

Replacement in projects/bisonex.org at line 3211 [10.35]

B:BD[21.2043] → [5.5118:5838]

 | 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 |                        |                        |        1.1818181818181819 |        1.6129032258064515 |
 |       0.243902 |        0.629617 |                        |                        |        1.1818181818181819 |        1.6129032258064515 |
 |       0.107547 |        0.849289 |     3.0985915492957745 |      2.925925925925926 |        1.5304347826086957 |         1.774869109947644 |
 |       0.107547 |        0.849289 |     3.0985915492957745 |      2.925925925925926 |        1.5304347826086957 |         1.774869109947644 |

[21.2043]

[4.2729]

#+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]

Replacement in projects/bisonex.org at line 3229 [10.35]

B:BD[4.2730] → [4.2730:2975]

******** NC_000021.9     14109044:
Dans bam filtré sur le chromosome 21 mais non dans vcf...
On appelle haplotype caller sans --max-mnp-distance
#+begin_src slurm
!/bin/bash
#SBATCH -c 4
#SBATCH -p smp
#SBATCH --time=08:00:00
#SBATCH --mem=32G

[4.2730]

[4.2975]

#+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

Deletion in projects/bisonex.org at line 3239 [10.35]

B:BD[4.2976] → [4.2976:3473]

module load nix/2.11.0
dir=/Work/Groups/bisonex/data/giab/
genomeRef=/Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef.fna
dbsnpDir=/Work/Projects/bisonex/data-alexis-reference/dbSNP
bam=/Work/Users/apraga/bisonex/script/files/bam/NA12878_NIST7035_recalibrated_hg38.bam
vcf=/Work/Users/apraga/bisonex/script/files/vcf/NA12878_NIST7035_nodistance.vcf
gatk --java-options "-Xmx32g" HaplotypeCaller \
     -R $genomeRef \
     -I $bam \
     -O $vcf \
     -D "$dbsnpDir"/GCF_000001405.39.gz \

Replacement in projects/bisonex.org at line 3240 [10.35]

B:BD[4.3483] → [4.3483:3539]

Idem
Sur le bam en entier, on le retrouve bien aussi...

[4.3483]

[13.818]

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

Insertion in notes/japanese.org at line 486 [6.353504]

[22.1363]


* <2023-02-26 Sun> - Aya-sensei : discussion
Grammaire
-ta koto ga arimasu : déjà fait
 ちゃんと : suffisament
 飼う(ka): avoir un animal
  いがくせいぶつがく 医学生物学 = biologie médicale
  sample = sample
  かんじゃ 患者 patient
  けんきゅう 研究 recherche
  けっか 結果 résultat
  はんしょくき 繁殖期 saison de reproduction

Insertion in notes/emacs_tips.org at line 182 [6.577073]
[2.113]
[2.113]
```
org-insert-structure-template
```
Replacement in notes/books.org at line 14 [6.675275]
B:BD[23.596] → [14.9355:9385]
```
*** TODO Benefits of Old Laws
```
[23.596]
[14.9385]
```
*** STRT Benefits of Old Laws
```
Replacement in notes/books.org at line 30 [6.675275]
B:BD[23.962] → [24.1328:1336]
```
** Pern
```
[23.962]
[23.977]
```
** HOLD Pern
```
Replacement in notes/books.org at line 39 [6.675275]
∅:D[25.268] → [23.1506:1540]
B:BD[23.1506] → [23.1506:1540]
```
*** TODO The Skies of Pern (2001)
```
[25.268]
[23.1870]
```
*** DONE The Skies of Pern (2001)
CLOSED: [2023-02-25 Sat 18:57]
```