B:BD[
6.24367] → [
6.24367:24496]
B:BD[
6.24496] → [
4.24733:33072]
318565 g.89318565A>G snv heterozygous A G 303 98 404
24 │ NC_000015.10 89
323426 g.89323426C>G snv heterozygous C G 93 109 202
25 │ NC_000015.10 89318595 g.89318595T>C snv heterozygous T C 321 128 453
26 │ NC_000015.10 48488437 g.48488437T>C snv heterozygous T C 356 132 488
CLOSED: [2023-05-01 Mon 17:18]
***** KILL Chromosome1 15 :Test haplotype caller : échec car CIGARE non mis à jour
CLOSED: [2023-05-13 Sat 18:29] SCHEDULED: <2023-05-01 Mon>
#+begin_src
julia -Jbisonex.so --project=. insertVariants.jl `63003856_S135_chr15.bam` 63003856_S135_chr15_inserted.bam
scp 63003856_S135_chr15_inserted.bam* meso:/Work/Users/apraga/bisonex/tests/synthetic/
#+end_src
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/synthetic :results silent
ln -s /Work/Projects/bisonex/data/dbSNP/GRCh38.p13/dbSNP.gz .
ln -s /Work/Projects/bisonex/data/dbSNP/GRCh38.p13/dbSNP.gz.tbi
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef.dict .
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef.fna .
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef.fna.fai .
#+end_src
puis
#+begin_src
gatk --java-options "-Xmx3072M" HaplotypeCaller --input 63003856_S135_chr15_inserted.bam --output testchr15.vcf.gz --reference genomeRef.fna --tmp-dir . -L NC_000015.10
#+end_src
scp meso:/Work/Users/apraga/bisonex/tests/synthetic/testchr15.vcf.gz haplotypecaller-chr15.vcf.gz
Aucun variant inséré
- base quality ok
-
****** DONE bam out : non appelé
CLOSED: [2023-05-01 Mon 21:57]
gatk --java-options "-Xmx3072M" HaplotypeCaller --input 63003856_S135_chr15_inserted.bam --output haplotypecaller-chr15.vcf.gz --reference genomeRef.f
na --tmp-dir . -L NC_000015.10 --bam-output debug.bam
****** DONE --linked-de-bruijn-graph : idem
CLOSED: [2023-05-01 Mon 21:57]
readlink testchr15.vcf.gz -f^C
[apraga@mesointeractive synthetic]$ gatk --java-options "-Xmx3072M" HaplotypeCaller --input 63003856_S135_chr15_inserted.bam --output haplotypecaller-chr15.vcf.gz --reference genomeRef.fna --tmp-dir . -L NC_000015.10 --linked-de-bruijn-graph
****** KILL regénérer fastq
CLOSED: [2023-05-13 Sat 18:29]
Non
***** KILL Générer bam données pour tous les chromosomes
CLOSED: [2023-05-13 Sat 18:29]
timeit julia -Jbisonex.so --project=. insertVariants.jl ~/code/bisonex/out/63003856/preprocessing/63003856_S135.bam 63003856_S135_inserted.bam
40min 516ms 835µs 405ns
Avertissement:
[W::bam_hdr_read] EOF marker is absent. The input is probably truncated
Inserted.bam et excluded.bam (fichier avant le merge) ont l'air ok...
On réessaie à la main : ça passe
#+begin_src
samtools merge test-all.bam inserted.bam excluded.bam
❯ mv test-all.bam `63003856_S135_inserted.bam` -f
❯ mv test-all.bam.bai `63003856_S135_chr15_inserted.bam.bai` -f
#+end_src
***** DONE BAm2fastq pour avoir CIGAR à jour : échec (variants "cachés")
CLOSED: [2023-05-04 Thu 20:30] SCHEDULED: <2023-05-01 Mon>
On lance la génération de bam depuis le mesocentro (la copie plante via le VPN)
#+begin_src sh
cd /Work/Users/apraga/recherche/bisonex/generate
julia --project=. insertVariants.jl ../../../bisonex/out/63003856_S135/preprocessing/applybqsr/63003856_S135.bam 63003856_S135_inserted.bam
#+end_src
Workflow après avec désactivé storeDir pour SAMTOOLS_BAM2FQ dans nextflow.config (pourquoi ??)
#+begin_src nextflow
include { SAMTOOLS_BAM2FQ } from "${params.modulesDir}/samtools/bam2fq/main"
include { SAMTOOLS_SORT as sortBamByName } from "${params.modulesDir}/samtools/sort/main"
workflow {
f = Channel.fromPath("${params.dataDir}/synthetic/63003856_S135_inserted.bam",
checkIfExists: true).map{it -> [["id": "synthetic_63003856"], it]}
// Important: use "-n" option !!
sortBamByName(f)
SAMTOOLS_BAM2FQ(sortBamByName.out.bam, true)
}
#+end_src
Puis
#+begin_src
cp work/34/fb2fc136f6f6d7f42d0960512f06de/*.fq.gz /Work/Groups/bisonex/data/synthetic/
#+end_src
***** KILL Lancer pipeline
CLOSED: [2023-05-04 Thu 20:30] SCHEDULED: <2023-05-01 Mon>
NXF_OPTS=-D"user.name=apraga" nextflow run main.nf -c nextflow.config -profile standard,helios -bg --input="/Work/Groups/bisonex/data/synthetic/synthetic_63003856_{1,2}.fq.gz" --outdir out/synthetic_63003856
*** TODO Bamsurgeon :bamsurgeon:
**** TODO Package nix
1. Patcher la recherche du génome de référence pour bien trouver les index (en utilisant une regexp comme nf-core)
2. Rajouter le chemin de picard dans les arguments
3. Option -O3 pour performance
***** DONE Erreur ValueError: quality and sequence mismatch
CLOSED: [2023-05-19 Fri 18:44]
****** DONE Idem avec dernière version sur github
CLOSED: [2023-05-18 Thu 14:36]
****** DONE Version 1.3: ok mais
CLOSED: [2023-05-19 Fri 18:44]
Test sur chr22: variants ok mais VAF=1...
S'exécute "normalement" (échec selon nextflow) mais le bam de sorstie quasiment vide
La fusion du bam avec les variants et du fichier de référence n'a fonctionné correctement.
******* DONE Lancer replacereads.py à la main
CLOSED: [2023-05-18 Thu 21:41]
Dans /Work/Users/apraga/bisonex/work/f3/ce044f80ca91016d68d1bc4f4f5301
#+begin_src sh
/nix/store/xw277la6w4sjqlsvw9h32cvrlacrfkgm-python3-3.10.9-env/bin/python3.10 /nix/store/abzangf0q8k37053p776cfkw181dzjn3-bamsurgeon-1.3/bin/bamsurgeon/replacereads.py -b cento.bam -r addsnv.e14561be-4fdd-45cc-9989-048ab6da6cc6.muts.bam -o snv-manual.bam
#+end_src
Puis
#+begin_src
samtools sort snv-manual.bam -o snv-manual-sorted.bam
#+end_src
A l'air de fonctionne
***** TODO Corriger run nextflow pour éviter les erreurs
SCHEDULED: <2023-05-21 Sun>
Trop de message d'erreur en sortie ?
**** DONE Test sur mini-bam: échec
CLOSED: [2023-05-14 Sun 21:12]
❯ samtools view -h ~/code/bisonex/simuscop-centogene-200x/cento/preprocessing/mapped/cento.bam | head -n1000 | samtools view -Sb - > mini.bam
❯ samtools index mini.bam
Sans spécfier le variant:
#+begin_quote
NC_000001.11 17651 17651
#+end_quote
./result/bin/addsnv -v snv.txt -f mini.bam -r ../data/genomeRef.fna -o test.bam
**** DONE Test chr22
CLOSED: [2023-05-15 Mon 23:24]
Pas assez de reads, on prend le chromosome 22
#+begin_src sh
samtools view ../simuscop-centogene-200x/cento/preprocessing/mapped/cento.bam NC_000022.11 -b -o chr22.bam
samtools index chr22.bam
#+end_src
Mésocentre
dans tests/bamsurgeno
#+begin_src
addsnv -v snv.txt -f chr22.bam -r ../genomeRef.fna -o test.bam --aligner mem
#+end_src
***** DONE SNV aléatoire:
CLOSED: [2023-05-15 Mon 23:13]
NC_000022.11 17499704 17499704 0.2
On retrouve bien un variant à cette position A > T
***** DONE SNV avec ALT prédéfini : retrouvée dans IGV (mais pas dans pileup)
CLOSED: [2023-05-15 Mon 23:13]
NC_000022.11 17499704 17499704 0.2 G
***** DONE Variants patients chr22: ok IGV
CLOSED: [2023-05-15 Mon 23:23]
Fichier non trié donc
samtools sort test.bam -o test-sorted.bam
samtools index test-sorted.bam
***** DONE Vérifier qu'il faut POS et POS+1: non
CLOSED: [2023-05-14 Sun 21:21]
**** TODO Variants cento
***** STRT SNV
SCHEDULED: <2023-05-15 Mon>
Attention à la mémoire: 32G ne semble pas suffire avec 12 threads
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/bamsurgeon.nf -profile standard,helios --input=tests/bamsurgeon/snv-cento.tsv -bg
#+end_src
ET
#+begin_src nextflow
workflow {
f = Channel.fromPath(params.input, checkIfExists: true)
bam = Channel.fromPath("simuscop-centogene-200x/cento/preprocessing/mapped/cento.bam",
checkIfExists: true)
bamIndex = bam.map { it -> it + ".bai" }
downloadGenome | indexGenome
indexGenome.out.index | view
addSNV(f, bam, bamIndex, downloadGenome.out, indexGenome.out.index, indexGenome.out.dict, indexGenome.out.fai)
}
#+end_src
****** DONE v1.3: Lancer le pipeline pour vérifier qu'on retrouve les variants
CLOSED: [2023-05-19 Fri 18:41] SCHEDULED: <2023-05-18 Thu>
****** TODO Forcer la VAF
SCHEDULED: <2023-05-19 Fri>
Est à 1
***** TODO del
***** TODO ins
*** Divers
**** DONE Vérifier nombre de reads fastq - bam
CLOSED: [2022-10-09 Sun 22:31]
318565 g.89318565A>G snv heterozygous A G 303 98 404
24 │ NC_000015.10 89323426 g.89323426C>G snv heterozygous C G 93 109 202
25 │ NC_000015.10 89318595 g.89318595T>C snv heterozygous T C 321 128 453
26 │ NC_000015.10 48488437 g.48488437T>C snv heterozygous T C 356 132 488
CLOSED: [2023-05-01 Mon 17:18]
***** KILL Chromosome1 15 :Test haplotype caller : échec car CIGARE non mis à jour
CLOSED: [2023-05-13 Sat 18:29] SCHEDULED: <2023-05-01 Mon>
#+begin_src
julia -Jbisonex.so --project=. insertVariants.jl `63003856_S135_chr15.bam` 63003856_S135_chr15_inserted.bam
scp 63003856_S135_chr15_inserted.bam* meso:/Work/Users/apraga/bisonex/tests/synthetic/
#+end_src
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/synthetic :results silent
ln -s /Work/Projects/bisonex/data/dbSNP/GRCh38.p13/dbSNP.gz .
ln -s /Work/Projects/bisonex/data/dbSNP/GRCh38.p13/dbSNP.gz.tbi
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef.dict .
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef.fna .
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef.fna.fai .
#+end_src
puis
#+begin_src
gatk --java-options "-Xmx3072M" HaplotypeCaller --input 63003856_S135_chr15_inserted.bam --output testchr15.vcf.gz --reference genomeRef.fna --tmp-dir . -L NC_000015.10
#+end_src
scp meso:/Work/Users/apraga/bisonex/tests/synthetic/testchr15.vcf.gz haplotypecaller-chr15.vcf.gz
Aucun variant inséré
- base quality ok
-
****** DONE bam out : non appelé
CLOSED: [2023-05-01 Mon 21:57]
gatk --java-options "-Xmx3072M" HaplotypeCaller --input 63003856_S135_chr15_inserted.bam --output haplotypecaller-chr15.vcf.gz --reference genomeRef.f
na --tmp-dir . -L NC_000015.10 --bam-output debug.bam
****** DONE --linked-de-bruijn-graph : idem
CLOSED: [2023-05-01 Mon 21:57]
readlink testchr15.vcf.gz -f^C
[apraga@mesointeractive synthetic]$ gatk --java-options "-Xmx3072M" HaplotypeCaller --input 63003856_S135_chr15_inserted.bam --output haplotypecaller-chr15.vcf.gz --reference genomeRef.fna --tmp-dir . -L NC_000015.10 --linked-de-bruijn-graph
****** KILL regénérer fastq
CLOSED: [2023-05-13 Sat 18:29]
Non
***** KILL Générer bam données pour tous les chromosomes
CLOSED: [2023-05-13 Sat 18:29]
timeit julia -Jbisonex.so --project=. insertVariants.jl ~/code/bisonex/out/63003856/preprocessing/63003856_S135.bam 63003856_S135_inserted.bam
40min 516ms 835µs 405ns
Avertissement:
[W::bam_hdr_read] EOF marker is absent. The input is probably truncated
Inserted.bam et excluded.bam (fichier avant le merge) ont l'air ok...
On réessaie à la main : ça passe
#+begin_src
samtools merge test-all.bam inserted.bam excluded.bam
❯ mv test-all.bam `63003856_S135_inserted.bam` -f
❯ mv test-all.bam.bai `63003856_S135_chr15_inserted.bam.bai` -f
#+end_src
***** DONE BAm2fastq pour avoir CIGAR à jour : échec (variants "cachés")
CLOSED: [2023-05-04 Thu 20:30] SCHEDULED: <2023-05-01 Mon>
On lance la génération de bam depuis le mesocentro (la copie plante via le VPN)
#+begin_src sh
cd /Work/Users/apraga/recherche/bisonex/generate
julia --project=. insertVariants.jl ../../../bisonex/out/63003856_S135/preprocessing/applybqsr/63003856_S135.bam 63003856_S135_inserted.bam
#+end_src
Workflow après avec désactivé storeDir pour SAMTOOLS_BAM2FQ dans nextflow.config (pourquoi ??)
#+begin_src nextflow
include { SAMTOOLS_BAM2FQ } from "${params.modulesDir}/samtools/bam2fq/main"
include { SAMTOOLS_SORT as sortBamByName } from "${params.modulesDir}/samtools/sort/main"
workflow {
f = Channel.fromPath("${params.dataDir}/synthetic/63003856_S135_inserted.bam",
checkIfExists: true).map{it -> [["id": "synthetic_63003856"], it]}
// Important: use "-n" option !!
sortBamByName(f)
SAMTOOLS_BAM2FQ(sortBamByName.out.bam, true)
}
#+end_src
Puis
#+begin_src
cp work/34/fb2fc136f6f6d7f42d0960512f06de/*.fq.gz /Work/Groups/bisonex/data/synthetic/
#+end_src
***** KILL Lancer pipeline
CLOSED: [2023-05-04 Thu 20:30] SCHEDULED: <2023-05-01 Mon>
NXF_OPTS=-D"user.name=apraga" nextflow run main.nf -c nextflow.config -profile standard,helios -bg --input="/Work/Groups/bisonex/data/synthetic/synthetic_63003856_{1,2}.fq.gz" --outdir out/synthetic_63003856
*** TODO Bamsurgeon :bamsurgeon:
**** TODO Package nix
1. Patcher la recherche du génome de référence pour bien trouver les index (en utilisant une regexp comme nf-core)
2. Rajouter le chemin de picard dans les arguments
3. Option -O3 pour performance
***** DONE Erreur ValueError: quality and sequence mismatch
CLOSED: [2023-05-19 Fri 18:44]
****** DONE Idem avec dernière version sur github
CLOSED: [2023-05-18 Thu 14:36]
****** DONE Version 1.3: ok mais
CLOSED: [2023-05-19 Fri 18:44]
Test sur chr22: variants ok mais VAF=1...
S'exécute "normalement" (échec selon nextflow) mais le bam de sorstie quasiment vide
La fusion du bam avec les variants et du fichier de référence n'a fonctionné correctement.
******* DONE Lancer replacereads.py à la main
CLOSED: [2023-05-18 Thu 21:41]
Dans /Work/Users/apraga/bisonex/work/f3/ce044f80ca91016d68d1bc4f4f5301
#+begin_src sh
/nix/store/xw277la6w4sjqlsvw9h32cvrlacrfkgm-python3-3.10.9-env/bin/python3.10 /nix/store/abzangf0q8k37053p776cfkw181dzjn3-bamsurgeon-1.3/bin/bamsurgeon/replacereads.py -b cento.bam -r addsnv.e14561be-4fdd-45cc-9989-048ab6da6cc6.muts.bam -o snv-manual.bam
#+end_src
Puis
#+begin_src
samtools sort snv-manual.bam -o snv-manual-sorted.bam
#+end_src
A l'air de fonctionne
***** TODO Corriger run nextflow pour éviter les erreurs
SCHEDULED: <2023-05-21 Sun>
Trop de message d'erreur en sortie ?
***** DONE Test sur mini-bam: échec
CLOSED: [2023-05-14 Sun 21:12]
❯ samtools view -h ~/code/bisonex/simuscop-centogene-200x/cento/preprocessing/mapped/cento.bam | head -n1000 | samtools view -Sb - > mini.bam
❯ samtools index mini.bam
Sans spécfier le variant:
#+begin_quote
NC_000001.11 17651 17651
#+end_quote
./result/bin/addsnv -v snv.txt -f mini.bam -r ../data/genomeRef.fna -o test.bam
***** DONE Test chr22
CLOSED: [2023-05-15 Mon 23:24]
Pas assez de reads, on prend le chromosome 22
#+begin_src sh
samtools view ../simuscop-centogene-200x/cento/preprocessing/mapped/cento.bam NC_000022.11 -b -o chr22.bam
samtools index chr22.bam
#+end_src
Mésocentre
dans tests/bamsurgeno
#+begin_src
addsnv -v snv.txt -f chr22.bam -r ../genomeRef.fna -o test.bam --aligner mem
#+end_src
****** DONE SNV aléatoire:
CLOSED: [2023-05-15 Mon 23:13]
NC_000022.11 17499704 17499704 0.2
On retrouve bien un variant à cette position A > T
****** DONE SNV avec ALT prédéfini : retrouvée dans IGV (mais pas dans pileup)
CLOSED: [2023-05-15 Mon 23:13]
NC_000022.11 17499704 17499704 0.2 G
****** DONE Variants patients chr22: ok IGV
CLOSED: [2023-05-15 Mon 23:23]
Fichier non trié donc
samtools sort test.bam -o test-sorted.bam
samtools index test-sorted.bam
****** DONE Vérifier qu'il faut POS et POS+1: non
CLOSED: [2023-05-14 Sun 21:21]
**** TODO Variants cento
***** STRT SNV
SCHEDULED: <2023-05-15 Mon>
Attention à la mémoire: 32G ne semble pas suffire avec 12 threads
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/bamsurgeon.nf -profile standard,helios --input=tests/bamsurgeon/snv-cento.tsv -bg
#+end_src
ET
#+begin_src nextflow
workflow {
f = Channel.fromPath(params.input, checkIfExists: true)
bam = Channel.fromPath("simuscop-centogene-200x/cento/preprocessing/mapped/cento.bam",
checkIfExists: true)
bamIndex = bam.map { it -> it + ".bai" }
downloadGenome | indexGenome
indexGenome.out.index | view
addSNV(f, bam, bamIndex, downloadGenome.out, indexGenome.out.index, indexGenome.out.dict, indexGenome.out.fai)
}
#+end_src
****** DONE v1.3: Lancer le pipeline pour vérifier qu'on retrouve les variants
CLOSED: [2023-05-19 Fri 18:41] SCHEDULED: <2023-05-18 Thu>
****** TODO Corrigier position pour avoir une bonne VAF
SCHEDULED: <2023-05-19 Fri>
POS POS+1 VAF ALT
Attention, la base corrigée est à POS+1...
****** DONE Comparaison manuelle avec julia (VAF = 1...)
CLOSED: [2023-05-19 Fri 21:58]
552 found over 585
***** TODO del
***** TODO ins
*** TODO Julia : Bamscissors :bamscissors:
Modification de la séquence dans BAM.
*Pas de mise à jour de CIGAR*
On convertit en fastq et on lance le pipeline pour "corriger"
#+begin_src sh
cd /home/alex/code/bisonex/out/63003856/preprocessing/mapped
samtools view 63003856_S135.bam NC_000022.11 -o 63003856_S135_chr22.bam
cd /home/alex/recherche/bisonex/code/BamScissors.jl
cp ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135_chr22.bam .
samtools index 63003856_chr22.bam
#+end_src
Le script va modifier le bam, le trier et générer le fastq. !!!
Attention: ne pas oublier l'option -n !!!
#+begin_src sh
time julia --project=.. insertVariant.jl
scp 63003856_S135_chr22_{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
#+end_src
**** TODO Phase 1 : chr22, VAF=1
***** DONE 1 SNV : ok !
CLOSED: [2023-05-20 Sat 19:35] SCHEDULED: <2023-05-20 Sat>
#+begin_src
make run READS="tests/bamscissors/corrected_{1,2}.fq.gz"
Puis on lance le pipeline sur correct_1
- [X] Variant visible dans IGV
- [X] Variant visible après alignement
- [X] Variant visible après appel de variant
***** TODO Tous les SNV du chromosome
SCHEDULED: <2023-05-20 Sat>
**** TODO PHase 2 : chr22, VAF variable: de nombreux reads sont perdus
Après l'alignement donc soit à cause de bam2fq, soit à cause de l'alignement.
Plutôt bam2fq ?
ex: chr22:42,202,582-42,223,574
SCHEDULED: <2023-05-22 Mon>
*** Divers
**** DONE Vérifier nombre de reads fastq - bam
CLOSED: [2022-10-09 Sun 22:31]