B:BD[
5.63935] → [
5.63935:64841]
B:BD[
5.64841] → [
3.57654:65964]
∅:D[
3.65964] → [
5.73033:73761]
B:BD[
5.73033] → [
5.73033:73761]
B:BD[
5.73761] → [
6.32407:40637]
∅:D[
6.40637] → [
5.81953:89379]
B:BD[
5.81953] → [
5.81953:89379]
B:BD[
5.89379] → [
3.65965:82349]
∅:D[
3.82349] → [
6.56256:63005]
B:BD[
6.56256] → [
6.56256:63005]
∅:D[
6.63005] → [
5.98337:105848]
B:BD[
5.98337] → [
5.98337:105848]
∅:D[
5.105848] → [
7.21561:24047]
B:BD[
7.21561] → [
7.21561:24047]
B:BD[
7.24047] → [
8.2320:8026]
B:BD[
8.8026] → [
5.105849:109094]
B:BD[
5.109094] → [
9.8932:17155]
∅:D[
9.17155] → [
8.19463:26896]
B:BD[
8.19463] → [
8.19463:26896]
B:BD[
8.26896] → [
10.27:89]
B:BD[
10.89] → [
11.16869:25064]
∅:D[
11.25064] → [
9.24651:32385]
B:BD[
9.24651] → [
9.24651:32385]
B:BD[
9.32385] → [
3.82350:86775]
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 HG003 :hg003:hg38:
CLOSED: [2023-04-16 Sun 00:20]
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run main.nf -profile standard,helios --input /Work/Groups/bisonex/data/giab/GRCh38/HG003_{1,2}.fq.gz -bg
#+end_src
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/compareVCF.nf -profile standard,helios -resume --outdir=compareHG003 --test.id=HG003 --test.query=out/HG003_1/variantCalling/haplotypecaller/HG003_1.vcf.gz --test.compare=vcfeval,happy --test.capture=data/AgilentSureSelectv05_hg38.bed
#+end_src
vcfeval
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
5.000 36745 36473 486 3988 0.9869 0.9021 0.9426
None 36748 36476 495 3985 0.9866 0.9022 0.9425
$ zcat NA12878.snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
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 2731 2290 441 3092 208 577 62 53 0.838521 0.917296 0.186611 0.876141 NaN NaN 1.505145 1.888993
INDEL PASS 2731 2290 441 3092 208 577 62 53 0.838521 0.917296 0.186611 0.876141 NaN NaN 1.505145 1.888993
SNP ALL 37997 34481 3516 36861 306 2074 33 13 0.907466 0.991204 0.056265 0.947488 2.611269 2.565915 1.555780 1.621727
SNP PASS 37997 34481 3516 36861 306 2074 33 13 0.907466 0.991204 0.056265 0.947488 2.611269 2.5659
**** DONE HG004 :hg38:
CLOSED: [2023-04-16 Sun 00:20]
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run main.nf -profile standard,helios --input /Work/Groups/bisonex/data/giab/GRCh38/HG004_{1,2}.fq.gz -bg
#+end_src
vcfeval
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
6.000 36938 36678 421 4040 0.9887 0.9014 0.9430
None 36942 36682 432 4036 0.9884 0.9015 0.9429
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 2787 2388 399 3183 195 580 53 38 0.856835 0.925086 0.182218 0.889654 NaN NaN 1.507834 1.848649
INDEL PASS 2787 2388 399 3183 195 580 53 38 0.856835 0.925086 0.182218 0.889654 NaN NaN 1.507834 1.848649
SNP ALL 38185 34560 3625 36921 254 2107 46 7 0.905067 0.992704 0.057068 0.946862 2.589175 2.553546 1.632595 1.653534
SNP PASS 38185 34560 3625 36921 254 2107 46 7 0.905067 0.992704 0.057068 0.946862 2.589175 2.553546 1.632595 1.653534
**** DONE Résumer résultats pour Paul + article :resultats:hg38:
CLOSED: [2023-04-06 Thu 21:41] SCHEDULED: <2023-04-02 Sun>
**** DONE Plot : ashkenazim trio :hg38:
CLOSED: [2023-04-18 Tue 21:27] SCHEDULED: <2023-04-16 Sun>
/Entered on/ [2023-04-16 Sun 17:29]
**** TODO NA12878 :na12878:T2T:
SCHEDULED: <2023-06-15 Thu>
*** KILL Platinum genome
CLOSED: [2023-06-14 Wed 22:37]
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 :centogene:
*** DONE Extraire liste des SNVs
CLOSED: [2023-04-22 Sat 17:32] SCHEDULED: <2023-04-17 Mon>
**** DONE Corriger manquant à la main
CLOSED: [2023-04-22 Sat 17:31]
La sortie est sauvegardé dans git-annex : variants_success.csv
**** DONE Automatique
CLOSED: [2023-04-22 Sat 17:31]
*** DONE Convert SNVs : transcript -> génomique
CLOSED: [2023-06-03 Sat 17:16]
**** DONE Variant_recoder
CLOSED: [2023-04-26 Wed 21:21] SCHEDULED: <2023-04-22 Sat>
***** KILL Haskell: 160 manquant : recoded-success.csv
CLOSED: [2023-04-25 Tue 18:32]
La liste des variants a été générée en Haskel l et nettoyée à la main.
On générer une liste de variant pour variant_rec oder et on soumet tout d'un coup.
[[file:~/recherche/bisonex/parsevariants/app/Main.hs][parsevariant]]
#+begin_src haskell
recodeVariant = do
prepareVariantRecod er "variant_success.csv" "renamed.csv"
runVariantRecoder "renamed.csv" "recoded.json"
#+end_src
#+RESULTS:
: <interactive>:4:3-19: error:
: Variable not in scope: runVariantRecoder :: String -> String -> t
: gh
Problème : 160 n'ont pas pu être lu sur 820, probablement à cause du numéro mineur de transcrit
La sortie est sauvegardé dans git-annex : variants-recoded-raw.json.
***** KILL Julia
CLOSED: [2023-04-25 Tue 18:32]
On regénère la liste de variant et on passe à Julia pour préparer l'appel en parallèle à variant recoder
[[file:~/recherche/bisonex/parsevariants/variantRecoder.jl][variantRecoder.jl]]
#+begin_src julia
setupVariantRecoder(unique(init), n)
#+end_src
Puis
#+begin_src sh
parallel -a parallel-recoder.sh --jobs 10
#+end_src
On récupère les résultats
#+begin_src julia
(fails, success) = mergeVariantRecoder(n)
CSV.write(fSuccess, success)
CSV.write(fFailures, fails)
#+end_src
Certains variants ne sont pas trouvé, donc on prépare un nouveau job en enlevant les versionrs mineures des transcrits
#+begin_src julia
# Cleanup json and txt
if isfile(fSuccess) && isfile(fFailures)
foreach(rm, variantRecoderInput())
foreach(rm, variantRecoderOutput())
end
redoFails(fFailures)
#+end_src
Puis
#+begin_src sh
parallel -a parallel-recoder.sh --jobs 3
#+end_src
Il manque encore 70 transcrits
**** DONE Julia avec mobidetails: recode-failures-mobidetails.csv
CLOSED: [2023-04-25 Tue 18:58]
Nouvelle stratégie : on essaie une fois variant recoder.
Pour tous les échecs, on utilise mobidetails (~170).
Si l'ID n'est pas trouvé, on incrémente le numéro de version 2 fois
**** DONE Reste une dizaine à corriger à la main
CLOSED: [2023-04-26 Wed 21:21]
- [X] certains transcrits ont juste été supprimé
- [X] Erreur de parsing, manque souvent un -
#+begin_src julia
lastTryMobidetails("recoded-failures-mobidetails.csv")
#+end_src
**** DONE Fusionner données
CLOSED: [2023-04-26 Wed 22:35]
#+begin_src julia
function merge
AllGenomic()
dNew = mergeAll("recoded-success.csv",
"recoded-failures-mobidetails.csv",
"recoded-failures-mobidetails-redo.csv")
dInit = @chain DataFrame(CSV.File("variant_success.csv")) begin
@transform :transcript = :transcript .* ":" .* :coding .* :codingPos .* :codingChange
@select :file :transcript :classification :zygosity
@rename :classificationCentogene = :classification
end
dTmp = outerjoin(dInit, dNew, on = :transcript)
CSV.write("variant_genomic.csv", dTmp)
end
fSuccess = "recoded-success.csv"
fFailures = "recoded-failures.csv"
# variantRecoder(fSuccess, fFailures)
# mobidetailsOnFailures(fFailures)
# la
es filtrées ? apparement non
CLOSED: [2023-04-30 Sun 16:41]
183122 read(s) filtered by: MappingQualityReadFilter
3674 read(s) filtered by: NotDuplicateReadFilter
********* DONE --disable-read-filter MappingQualityReadFilter: idem
CLOSED: [2023-04-30 Sun 16:34]
On a bien - 0 read(s) filtered by: MappingQualityAvailableReadFilter
********* DONE --disable-read-filter NotDuplicateReadFilter: idem
CLOSED: [2023-04-30 Sun 16:40]
******** DONE Essayer freebayes : idem
CLOSED: [2023-04-30 Sun 16:22]
freebayes -f genomeRef.fna -r NC_000015.10 cento_chr15.bam > freebayes-test-chr15.vcf
******** DONE Avec toutes les options : idem
--linked-de-bruijn-graph --recover-all-dangling-branches --min-pruning 0 --bam-output debug.bam
CLOSED: [2023-04-30 Sun 16:50]
******** DONE Vérifier qu'on regarde le même bam : oui
CLOSED: [2023-04-30 Sun 16:50]
******** DONE Désactiver dbSNP : idem
CLOSED: [2023-04-30 Sun 16:52]
******** DONE Changer kmer size : idem
CLOSED: [2023-04-30 Sun 16:56]
par exemple[[https://gatk.broadinstitute.org/hc/en-us/community/posts/360075653152-REAL-Variant-not-called-by-HaplotypeCaller][forum gatk]] --kmer-size 18 --kmer-size 22
******** DONE --adaptive-pruning true
CLOSED: [2023-05-01 Mon 19:57]
******* DONE Mapping quality : est à 0 !!!!
CLOSED: [2023-05-01 Mon 19:58]
****** KILL Comparer VCF avec vcfeval :haplotypecaller:
CLOSED: [2023-06-12 Mon 23:24]
On prépare les données en julia
#+begin_src ~/recherche/bisonex/simuscop
julia --project=. toVCF.jl
#+end_src
Puis on export sur le mésocentre
#+begin_src
scp variants_for_vcfeval.tsv.gz* meso:centogene_variants/
#+end_src
#+begin_src
z bis
cd simuscop-200x
rtg vcfeval -b ~/centogene_variants/variants_for_vcfeval.tsv.gz -c cento/variantCalling/haplotypecaller/cento.vcf.gz -o compare-haplotypecaller -t /Work/Groups/bisonex/data/giab/GRCh38/genomeRef.sdf
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
82.000 540 540 60 45 0.9000 0.9231 0.9114
None 546 546 329 39 0.6240 0.9333 0.7479
****** KILL Comparer avec hap.py :haplotypecaller:
CLOSED: [2023-06-12 Mon 23:24]
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/checkInserted.nf -profile standard,helios --outdir=compare-simuscop-200x --query=out/simuscop-centogene-200x/cento/callVariant/haplotypecaller/cento.vcf.gz --truth=centogene_variants/variants_for_vcfeval.tsv.gz --id=simuscop-200x-check
****** DONE Méthode naïve 549/585
CLOSED: [2023-05-04 Thu 21:57]
Haplotypecaller: Nb reference SNV 692 vs found 585
Variant calling, filter technical: reference SNV 692 vs found 521
***** KILL Avant annotation
CLOSED: [2023-06-12 Mon 23:25] SCHEDULED: <2023-04-28 Fri>
#+begin_src
cd cento/variantCalling
bgzip filter-technical.vcf
tabix -p vcf filter-technical.vcf.gz -f
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
12.000 519 519 55 66 0.9042 0.8872 0.8956
None 519 519 55 66 0.9042 0.8872 0.8956
****** DONE Méthode naïve 521/585
CLOSED: [2023-05-04 Thu 21:57]
Haplotypecaller: Nb reference SNV 692 vs found 585
Variant calling, filter technical: reference SNV 692 vs found 521
****** KILL Comparer avec hap.py
CLOSED: [2023-06-12 Mon 23:24]
***** KILL Après filtre annotation
CLOSED: [2023-06-12 Mon 23:25]
****** DONE Méthode naïve : 493/585
CLOSED: [2023-05-04 Thu 22:09]
****** KILL Comparer avec hap.py
CLOSED: [2023-06-12 Mon 23:25]
****** KILL VCf eval
CLOSED: [2023-06-12 Mon 23:25]
cd cento/annotation/
bgzip postvep-filter.vcf
tabix postvep-filter.vcf.gz
cd ../..
rtg vcfeval -b ~/centogene_variants/variants_for_vcfeval.tsv.gz -c cento/annotation/postvep-filter.vcf.gz -o compare-vepfilter -t /Work/Groups/bisonex/data/giab/GRCh38/genomeRef.sdf
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
12.000 491 491 50 94 0.9076 0.8393 0.8721
None 491 491 50 94 0.9076 0.8393 0.8721
**** TODO Vérifier tous les variants sont retrouvés en 200x: hg38 :T2T:
***** DONE Après alignement
CLOSED: [2023-04-29 Sat 18:27] SCHEDULED: <2023-04-28 Fri>
****** DONE SNV: avec doublons
CLOSED: [2023-04-28 Fri 18:12]
On utilise [[file:~/recherche/bisonex/simuscop/checkBam.jl][checkBam.jl]]
#+begin_src julia
d = prepareVariant("../parsevariants/variant_genomic.csv")
root = "/home/alex/code/bisonex/simuscop-centogene/cento"
bam = root * "/preprocessing/applybqsr/cento.bam"
bai = root * "/preprocessing/recalibrated/cento.bam.bai"
snv = getSNV(d, bam, bai)
#+end_src
Nombreux faux homozygouteS
Vérification avec checkFalseHemizygous(snv) : nombreux doublons dans le fichier pour simuscop...
****** DONE SNV sans doublons
CLOSED: [2023-04-29 Sat 18:27]
******* DONE 18 faux homozygote mais avec peu de reads
CLOSED: [2023-04-29 Sat 18:27]
julia> @subset snv :refCount .== 0 :altCount .> 0 :zygosity .== "heterozygous"
18×10 DataFrame
Row │ chrom pos variant variantType zygosity ref alt refCount altCount readsCount
│ SubStrin…? Int64 SubStrin…? String? String15 SubStrin… SubStrin… Int64 Int64 Int64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ NC_000022.11 42213078 g.42213078T>G snv heterozygous T G 0 1 1
2 │ NC_000012.12 101680427 g.101680427C>A snv heterozygous C A 0 3 3
3 │ NC_000014.9 105385684 g.105385684G>C snv heterozygous G C 0 4 4
4 │ NC_000011.10 125978299 g.125978299C>T snv heterozygous C T 0 3 3
5 │ NC_000023.11 77998618 g.77998618C>T snv heterozygous C T 0 2 2
6 │ NC_000015.10 66703292 g.66703292C>T snv heterozygous C T 0 3 3
7 │ NC_000010.11 87961118 g.87961118G>A snv heterozygous G A 0 3 3
8 │ NC_000012.12 112477719 g.112477719A>G snv heterozygous A G 0 2 2
9 │ NC_000020.11 6778406 g.6778406C>T snv heterozygous C T 0 3 3
10 │ NC_000023.11 68192943 g.68192943G>A snv heterozygous G A 0 2 2
11 │ NC_000004.12 987858 g.987858C>T snv heterozygous C T 0 3 4
12 │ NC_000015.10 66435145 g.66435145G>A snv heterozygous G A 0 1 2
13 │ NC_000002.12 47809595 g.47809595C>T snv heterozygous C T 0 2 2
14 │ NC_000003.12 136477305 g.136477305C>G snv heterozygous C G 0 4 4
15 │ NC_000005.10 157285458 g.157285458C>T snv heterozygous C T 0 3 3
16 │ NC_000012.12 23604413 g.23604413T>G snv heterozygous T G 0 5 5
17 │ NC_000019.10 52219703 g.52219703C>T snv heterozygous C T 0 1 1
18 │ NC_000016.10 88856757 g.88856757C>T snv heterozygous C T 0 8 8
******* DONE 8 non retrouvé => probablement hors de la zjone de capture
CLOSED: [2023-04-28 Fri 19:49]
julia> @subset snv :refCount .== 0 :altCount .== 0
8×10 DataFrame
Row │ chrom pos variant variantType zygosity ref alt refCount altCount readsCount
│ SubStrin…? Int64 SubStrin…? String? String15 SubStrin… SubStrin… Int64 Int64 Int64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ NC_000015.10 74343027 g.74343027C>T snv heterozygous C T 0 0 0
2 │ NC_000011.10 20638345 g.20638345A>G snv heterozygous A G 0 0 0
3 │ NC_000004.12 139370252 g.139370252C>T snv heterozygous C T 0 0 2
4 │ NC_000017.11 61966475 g.61966475G>T snv heterozygous G T 0 0 0
5 │ NC_000019.10 54144058 g.54144058G>A snv heterozygous G A 0 0 0
6 │ NC_000023.11 77635947 g.77635947A>G snv hemizygous A G 0 0 0
7 │ NC_000005.10 1258495 g.1258495G>A snv heterozygous G A 0 0 0
8 │ NC_000012.12 2449086 g.2449086C>G snv heterozygous C G 0 0 0
***** KILL Après haplotypecaller
CLOSED: [2023-06-12 Mon 23:24]
****** KILL 20x
CLOSED: [2023-04-29 Sat 15:39]
Manque 183 sur 766
[[file:~/recherche/bisonex/simuscop/checkVCF.jl][checkVCF.jl]]
#+begin_src julia
@subset leftjoin(d2, dHaplo2, on=:genomic) ismissing.(:Column1)
#+end_src
Problème de profondeur ?
Ex: chr13 nombre de 101081606
NC_000011.10 16014966 g.16014966G>A
1 read sur 11 pour allèle alternative
Sur le patient de référence, 202 reads!
Celui-ci n'est pas le fichier de capture (ni dans le bam !)
ex: NC_000015.10 74343027 g.74343027C>T
Pour les autres, on devrait les retrouver...
Vérifier le nombre de reads sur 63003856
Vérifier la paramétrisation du modèle également
****** DONE [#B] 200x
CLOSED: [2023-05-18 Thu 11:04] SCHEDULED: <2023-04-30 Sun>
120 manquants (99 sans doublon)!
On vérifie dans IGV (vcf + bam après alignement) :
******* snv NC_000015.10 74343027
- rien d'appelé
- pas une région répétée
- base quality (voir [[*Phred score][Phred score]] ) à 37 donc ok
- variant retrouvé à 26/42
- Bam après aplybqsr: base qualità 35 donc ok
chr15 également à 89318565, variant retrouvé à 25/33 avec basequal de 37
Sans oublier de charger les instructions avx
#+begin_src sh
module load gcc@11.3.0/gcc-12.1.0
#+end_src
On coupe le .bam par chromosome pour débugger (sur le mesocentre)
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/simuscop-centogene-200x/cento/testing :results silent
ln -s ../preprocessing/applybqsr/cento.bam .
ln -s ../preprocessing/recalibrated/cento.bam.bai .
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
On doit lancer à la main (org-mode ne connait pas le chemin de samtools)
samtools view -b cento.bam NC_000015.10 > cento_chr15.bam
samtools index cento_chr15.bam
Puis on se restreint au chronmosome 15
samtools faidx genomeRef.fna NC_000015.10 > genomeRef_chr15.fa
samtools faidx genomeRef_chr15.fa
gatk CreateSequenceDictionary -R genomeRef_chr15.fa -O genomeRef_chr15.dict
On restreint au chromosome 15 avec l'option -L (dure = 1min)
gatk --java-options "-Xmx3072M" HaplotypeCaller --input cento_chr15.bam \
--output test.vcf.gz --reference genomeRef.fna --dbsnp dbSNP.gz --tmp-dir . --max-mnp-distance 2 -L NC_000015.10
******* DONE Tutorial haplotycaller
CLOSED: [2023-05-01 Mon 19:58]
Procédure : https://gatk.broadinstitute.org/hc/en-us/articles/360043491652-When-HaplotypeCaller-and-Mutect2-do-not-call-an-expected-variant
******** DONE Supprimer --max-mnp-distance = 2: idem
CLOSED: [2023-04-30 Sun 15:42]
******** DONE --debug &> run.log : Non appelé...
CLOSED: [2023-04-30 Sun 15:52]
******** DONE --linked-de-bruijn-graph: idem
CLOSED: [2023-04-30 Sun 15:55]
******** DONE --recover-all-dangling-branches
CLOSED: [2023-04-30 Sun 16:01]
******** DONE --min-pruning 0 : plus mais pas celui là
CLOSED: [2023-04-30 Sun 15:59]
******** DONE --bam-output
CLOSED: [2023-04-30 Sun 16:50]
********* DONE : rien !
CLOSED: [2023-04-30 Sun 16:08]
********* DONE + --recover-all-dangling-branches : rien !
CLOSED: [2023-04-30 Sun 16:08]
******** DONE Données filtrées ? apparement non
CLOSED: [2023-04-30 Sun 16:41]
183122 read(s) filtered by: MappingQualityReadFilter
3674 read(s) filtered by: NotDuplicateReadFilter
********* DONE --disable-read-filter MappingQualityReadFilter: idem
CLOSED: [2023-04-30 Sun 16:34]
On a bien - 0 read(s) filtered by: MappingQualityAvailableReadFilter
********* DONE --disable-read-filter NotDuplicateReadFilter: idem
CLOSED: [2023-04-30 Sun 16:40]
******** DONE Essayer freebayes : idem
CLOSED: [2023-04-30 Sun 16:22]
freebayes -f genomeRef.fna -r NC_000015.10 cento_chr15.bam > freebayes-test-chr15.vcf
******** DONE Avec toutes les options : idem
--linked-de-bruijn-graph --recover-all-dangling-branches --min-pruning 0 --bam-output debug.bam
CLOSED: [2023-04-30 Sun 16:50]
******** DONE Vérifier qu'on regarde le même bam : oui
CLOSED: [2023-04-30 Sun 16:50]
******** DONE Désactiver dbSNP : idem
CLOSED: [2023-04-30 Sun 16:52]
******** DONE Changer kmer size : idem
CLOSED: [2023-04-30 Sun 16:56]
par exemple[[https://gatk.broadinstitute.org/hc/en-us/community/posts/360075653152-REAL-Variant-not-called-by-HaplotypeCaller][forum gatk]] --kmer-size 18 --kmer-size 22
******** DONE --adaptive-pruning true
CLOSED: [2023-05-01 Mon 19:57]
******* DONE Mapping quality : est à 0 !!!!
CLOSED: [2023-05-01 Mon 19:58]
****** KILL Comparer VCF avec vcfeval :haplotypecaller:
CLOSED: [2023-06-12 Mon 23:24]
On prépare les données en julia
#+begin_src ~/recherche/bisonex/simuscop
julia --project=. toVCF.jl
#+end_src
Puis on export sur le mésocentre
#+begin_src
scp variants_for_vcfeval.tsv.gz* meso:centogene_variants/
#+end_src
#+begin_src
z bis
cd simuscop-200x
rtg vcfeval -b ~/centogene_variants/variants_for_vcfeval.tsv.gz -c cento/variantCalling/haplotypecaller/cento.vcf.gz -o compare-haplotypecaller -t /Work/Groups/bisonex/data/giab/GRCh38/genomeRef.sdf
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
82.000 540 540 60 45 0.9000 0.9231 0.9114
None 546 546 329 39 0.6240 0.9333 0.7479
****** KILL Comparer avec hap.py :haplotypecaller:
CL
OSED: [2023-06-12 Mon 23:24]
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/checkInserted.nf -profile standard,helios --outdir=compare-simuscop-200x --query=out/simuscop-centogene-200x/cento/callVariant/haplotypecaller/cento.vcf.gz --truth=centogene_variants/variants_for_vcfeval.tsv.gz --id=simuscop-200x-check
****** DONE Méthode naïve 549/585
CLOSED: [2023-05-04 Thu 21:57]
Haplotypecaller: Nb reference SNV 692 vs found 585
Variant calling, filter technical: reference SNV 692 vs found 521
***** KILL Avant annotation
CLOSED: [2023-06-12 Mon 23:25] SCHEDULED: <2023-04-28 Fri>
#+begin_src
cd cento/variantCalling
bgzip filter-technical.vcf
tabix -p vcf filter-technical.vcf.gz -f
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
12.000 519 519 55 66 0.9042 0.8872 0.8956
None 519 519 55 66 0.9042 0.8872 0.8956
****** DONE Méthode naïve 521/585
CLOSED: [2023-05-04 Thu 21:57]
Haplotypecaller: Nb reference SNV 692 vs found 585
Variant calling, filter technical: reference SNV 692 vs found 521
****** KILL Comparer avec hap.py
CLOSED: [2023-06-12 Mon 23:24]
***** KILL Après filtre annotation
CLOSED: [2023-06-12 Mon 23:25]
****** DONE Méthode naïve : 493/585
CLOSED: [2023-05-04 Thu 22:09]
****** KILL Comparer avec hap.py
CLOSED: [2023-06-12 Mon 23:25]
****** KILL VCf eval
CLOSED: [2023-06-12 Mon 23:25]
cd cento/annotation/
bgzip postvep-filter.vcf
tabix postvep-filter.vcf.gz
cd ../..
rtg vcfeval -b ~/centogene_variants/variants_for_vcfeval.tsv.gz -c cento/annotation/postvep-filter.vcf.gz -o compare-vepfilter -t /Work/Groups/bisonex/data/giab/GRCh38/genomeRef.sdf
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
12.000 491 491 50 94 0.9076 0.8393 0.8721
None 491 491 50 94 0.9076 0.8393 0.8721
*** KILL NEAT : trop lent :neat:
CLOSED: [2023-04-29 Sat 22:06]
**** KILL Génération fastq sur exno 5 GATAD2B
CLOSED: [2023-04-29 Sat 22:06]
Trop lent : pour 1 exon : 1500 secondes !
#+begin_src sh
samtools faidx genomeRef.fna NC_000001.11 | save -f genomeRef_chr1.fna
python gen_reads.py -r ../test-simuscop/genomeRef_chr1.fna -o lol -tr ../test-simuscop/gatad2b-exon6.bed -R 147 --pe 150 10
#+end_src
*** KILL ReSeq : exome avec exons comme fasta mais ne gère pas des exons trop petits :reseq:
CLOSED: [2023-04-30 Sun 19:44] SCHEDULED: <2023-04-29 Sat>
#+begin_quote
Can I simulate exome sequencing? Yes. You need to use a reference that only contains the exons as individual scaffolds. Using --refBiasFile you can specify the coverage of individual exons. To simulate intron contamination you can add the whole reference to the reference containing the exons and strongly reduce the coverage for these scaffolds using --refBiasFile.
#+end_quote
Par contre, rapide
**** DONE Fasta pour exons seuls
CLOSED: [2023-04-30 Sun 19:25]
Depuis le GFF
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_genomic.gff.gz
#+end_src
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
gunzip -c GCF_000001405.39_GRCh38.p13_genomic.gff.gz | grep -w "exon" > exons.gff
#+end_src
On génère les exons
#+begin_src sh :dir ~/code/bisonex/test-reseq
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed exons.gff -fo exons.fna
#+end_src
A tester avec un profile déjà fait :
https://github.com/schmeing/ReSeq-profiles/tree/master/profiles
On cherche l'exons qui nous intéresse
NC_000001.11 g.153817496 A>T
N'y est pas ??
***** DONE On test sur les 2 premiers : exec
CLOSED: [2023-04-30 Sun 18:39]
#+begin_src
head exons.fa -n 2 > 2exons.fna
#+end_src
#+begin_src sh
../ReSeq/bin/reseq illuminaPE -j 32 -R exons.fa -s Ec-Hi2000-TruSeq.reseq --ipfIterations 0 -1 reseq-sim_1.fq reseq_sim_2.fq
#+end_src
#+begin_quote
error: All reference sequences are too short for simulating. They should have at least 1991 bases
#+end_quote
#+begin_src sh
grep '^>NC_000001.10' exons.fa | sed 's/:/,/;s/-/,/;s/^>//' > exons.csv
#+end_src
***** DONE Sur 200 premiers exons du chr1
CLOSED: [2023-04-30 Sun 19:17]
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
head -n200 exons.fna > exons-200.fna
bwa index exons-200.fna
#+end_src
Simulation avec 30x
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
../ReSeq/bin/reseq illuminaPE -R exons-200.fna -s Ec-Hi2000-TruSeq.reseq --ipfIterations 0 -1 reseq1.fq -2 reseq2.fq -c 30
#+end_src
Attention, pour l'alignement, il faut le nfa complet ! Sinon erreur du type
Erreurs:::sam_hdr_create] Duplicated sequence "NC_000001.10:762970-763155" in file "-"
Et pas de bam avec
samtools sort: failed to change sort order header to 'coordinate'
#+begin_src
bwa mem ../test-simuscop/bwa/genomeRef.fna reseq1.fq reseq2.fq | samtools sort -o reseq.bam
#+end_src
Manque des exons et l'allure ne correspond pas...
***** DONE Utiliser le fichier de capture : exons trop petits
CLOSED: [2023-04-30 Sun 19:25]
Comme pour ART
Trop court avec
echo -e "NC_000001.11\t153817371\t153817542" > gatad2b-exon6.bed
Donc on ajoute 1000 de chaque côté
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
echo -e "NC_000001.11\t153816371\t153818542" > gatad2b-exon6.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6.bed -fo gatad2b-exon6.fna
bwa index gatad2b-exon6.bed
../ReSeq/bin/reseq illuminaPE -R gatad2b-exon6.fna -s Ec-Hi2000-TruSeq.reseq --ipfIterations 0 -1 reseq1.fq -2 reseq2.fq -c 30
bwa mem ../test-simuscop/bwa/genomeRef.fna reseq1.fq reseq2.fq | samtools sort -o reseq.bam
samtools index reseq.bam
#+end_src
**** KILL Sur le chromosome 15 puis trier à la main sur les zones de capture ?
CLOSED: [2023-04-30 Sun 19:44]
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
samtools faidx ../test-simuscop/genomeRef.fna NC_000015.10 > chr15.fna
../ReSeq/bin/reseq illuminaPE -R chr15.fna -s Ec-Hi2000-TruSeq.reseq --ipfIterations 0 -1 reseq1.fq -2 reseq2.fq -c 30
#+end_src
*** DONE ART : fonctionne très mal en targeted
CLOSED: [2023-04-30 Sun 11:49]
**** DONE
Génération de reads
CLOSED: [2023-04-30 Sun 11:49]
***** DONE Avec seulement les exons en séquence
CLOSED: [2023-04-30 Sun 10:24]
head -n6 exons.fa | save three-exons.fna
../art_bin_MountRainier/art_illumina -ss HS25 -i three-exons.fna -o ./paired_end_com -l 150 -f 10 -p -m 500 -s 10 -sam
Le sam n'est pas visible sur igv mais si on aligne avec bwa mem, on a quelques reads
***** DONE Extraire une zone de capture dans le fasta
CLOSED: [2023-04-30 Sun 11:49]
NC_000001.11 g.153817496 A>T
****** DONE Essai 1: ne dépasse pas la zone
CLOSED: [2023-04-30 Sun 10:49]
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153817371\t153817542" > gatad2b-exon6.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6.bed -fo gatad2b-exon6.fa
#+end_src
-ss HS25 : nom du profile illumina
-l 150 : reads de 150
-f 10 : coverage de 10
-p : paired end
-m 500 : longueur moyenne des fragment d'ADN
-s 10 : déviation standard
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6.fa -o ./gatad2b -l 150 -f 100 -p -m 500 -s 10
#+end_src
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.b
am
samtools index gatad2b.bam
#+end_src
#+attr_html: :width 800px
[[./art-capture-1.png]]
****** Avec offset
50bp idem
NC_000001.11 153817371 153817542 -> NC_000001.11 153817321 153817592
On essaie 1000
NC_000001.11 153817371 153817542 -> NC_000001.11 153816371 153818542 ->
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153816371\t153818542" > gatad2b-exon6-offset.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6-offset.bed -fo gatad2b-exon6-offset.fa
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6-offset.fa -o ./gatad2b -l 150 -f 100 -p -m 500 -s 10
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.bam
samtools index gatad2b.bam
#+end_src
mieux mais trop large
#+attr_html: :width 800px
[[./art-exon6-offset1000.png]]
Sur les vraies données, on a une large de 500bp environ
#+attr_html: :width 800px
[[./illumina-ref-exon6.png]]
Ici l'exon fait 250bp donc on rajouter 125bp de chaque côté
NC_000001.11 153817371 153817542
NC_000001.11 153817246 153817667
Résulats incohérents :on a 2 colonnes séparées !
Essai avec +200bp de chaque côt
NC_000001.11 153817371 153817542
NC_000001.11 153817171 153817742
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153817171\t153817742"> gatad2b-exon6-offset.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6-offset.bed -fo gatad2b-exon6-offset.fa
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6-offset.fa -o ./gatad2b -l 150 -f 100 -p -m 500 -s 50
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.bam
samtools index gatad2b.bam
#+end_src
En changeant la longueur des reads
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153817171\t153817742"> gatad2b-exon6-offset.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6-offset.bed -fo gatad2b-exon6-offset.fa
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6-offset.fa -o ./gatad2b -l 126 -f 100 -p -m 500 -s 50
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.bam
samtools index gatad2b.bam
#+end_src
Idem, fait 2 "tas" séparés de chaque côté de l'exon
+/- 100
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153817271\t153817642" > gatad2b-exon6-offset.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6-offset.bed -fo gatad2b-exon6-offset.fa
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6-offset.fa -o ./gatad2b -l 125 -f 100 -p -m 500 -s 30
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.bam
samtools index gatad2b.bam
#+end_src
*** KILL NGSNG : plante
CLOSED: [2023-05-01 Mon 09:27]
https://github.com/RAHenriksen/NGSNGS
**** KILL Essai avec fasta
CLOSED: [2023-05-01 Mon 09:27]
#+begin_src sh :dir ~/code/bisonex/test-ngsngs :results silent
echo -e "NC_000001.11\t153817371\t153817542" > gatad2b-exon6.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6.bed -fo gatad2b-exon6.fa
#+end_src
-ss HS25 : nom du profile illumina
-l 150 : reads de 150
-f 10 : coverage de 10
-p : paired end
-m 500 : longueur moyenne des fragment d'ADN
-s 10 : déviation standard
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6.fa -o ./gatad2b -l 150 -f 100 -p -m 500 -s 10
#+end_src
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.b
am
samtools index gatad2b.bam
#+end_src
*** KILL Script maison :generate:
CLOSED: [2023-05-13 Sat 18:29] SCHEDULED: <2023-05-01 Mon>
**** KILL SNV
CLOSED: [2023-05-13 Sat 18:29] SCHEDULED: <2023-05-01 Mon>
***** DONE Script python: ok seulement pour reads corrigé, trop long sinon
CLOSED: [2023-05-01 Mon 13:32]
Même sur un seul chromosome (15)...
***** DONE Couper les données avec bedtools intersect ?
CLOSED: [2023-05-01 Mon 20:33]
-wa pour avoir les reads qui corresponds
-v pour ceux qui n'intersecte pass
ok sur un exemple
****** DONE Test simple : ok
CLOSED: [2023-05-01 Mon 13:43]
#+attr_html: :width 500px
[[./test-intersect-chr15.png]]
#+attr_html: :width 500px
[[./test-nointersect-chr15.png]]
****** DONE Chromosome 15 : vérifier BAM
CLOSED: [2023-05-01 Mon 20:33] SCHEDULED: <2023-05-01 Mon>
#+begin_src
samtools view -b `63003856_S135.bam` NC_000015.10 > `63003856_S135_chr15.bam`
#+end_src
On génère un python avec les dépendances
#+begin_src
nix-build
#+end_src
Puis on lance un script julia qui va couper le bam en 2, lancer le script python sur l'intersection et fusionner le résultat
#+begin_src
julia -Jbisonex.so --project=. insertVariants.jl `63003856_S135_chr15.bam` test_new.bam
#+end_src
NB: pour accélerer l'exécution, générer une sysimage :
#+begin_src julia
(1.7) pkg> activate .
julia> create_sysimage(; sysimage_path="bisonex.so")
#+end_src
3 variants: Ok sur le nombre de reads et les variants
Ok pour homozygote
#+attr_html: :width 500px
[[./check-snv-chr15.png]]
****** DONE Chromosome 15 Vérifier VAF avec checkBam.jl: ok
julia> @subset snv :chrom .== "NC_000015.10"
26×10 DataFrame
Row │ chrom pos variant variantType zygosity ref alt refCount altCount readsCount
│ SubStrin…? Int64 SubStrin…? String? String15 SubStrin… SubStrin… Int64 Int64 Int64
─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ NC_000015.10 74343027 g.74343027C>T snv
heterozygous C T 61 58 120
2 │ NC_000015.10 75400778 g.75400778C>G snv heterozygous C G 108 79 187
3 │ NC_000015.10 89327201 g.89327201C>T snv heterozygous C T 243 241 486
4 │ NC_000015.10 48767448 g.48767448A>C snv heterozygous A C 72 70 142
5 │ NC_000015.10 75411685 g.75411685T>C snv heterozygous T C 79 81 160
6 │ NC_000015.10 66703292 g.66703292C>T snv heterozygous C T 70 60 130
7 │ NC_000015.10 89325639 g.89325639G>A snv heterozygous G A 257 267 524
8 │ NC_000015.10 89330184 g.89330184G>A snv heterozygous G A 258 287 548
9 │ NC_000015.10 89330184 g.89330184G>A snv heterozygous G A 258 287 548
10 │ NC_000015.10 89325639 g.89325639G>A snv heterozygous G A 257 267 524
11 │ NC_000015.10 42401752 g.42401752G>A snv homozygous G A 61 212 273
12 │ NC_000015.10 89327201 g.89327201C>T snv heterozygous C T 243 241 486
13 │ NC_000015.10 38339896 g.38339896G>A snv heterozygous G A 56 86 144
14 │ NC_000015.10 26869324 g.26869324A>T snv heterozygous A T 62 49 113
15 │ NC_000015.10 66435145 g.66435145G>A snv heterozygous G A 98 95 193
16 │ NC_000015.10 60514655 g.60514655G>A snv heterozygous G A 94 99 194
17 │ NC_000015.10 42410947 g.42410947A>G snv heterozygous A G 153 123 276
18 │ NC_000015.10 75430368 g.75430368C>T snv heterozygous C T 80 62 142
19 │ NC_000015.10 25375494 g.25375494T>C snv heterozygous T C 103 104 207
20 │ NC_000015.10 60497497 g.60497497C>A snv heterozygous C A
61 65 126
21 │ NC_000015.10 74891539 g.74891539C>T snv heterozygous C T 118 124 242
22 │ NC_000015.10 48488433 g.48488433A>G snv heterozygous A G 367 122 492
23 │ NC_000015.10 89318565 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 /Wor
k/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
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]
**** HOLD Variants cento
***** STRT SNV
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>
****** HOLD Corrigier position pour avoir une bonne VAF
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 [[id:966a298c-948a-4694-a6f5-c326b1046a05][XAMscissors.jl]] :xamscissors:
**** TODO Test SNV
***** DONE Phase 1 : chr22, VAF=1
CLO
SED: [2023-05-29 Mon 15:36]
****** 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
****** KILL Tester SNV chromosome 22
CLOSED: [2023-05-29 Mon 15:36] SCHEDULED: <2023-05-20 Sat>
***** KILL PHase 2 : chr22, VAF variable
CLOSED: [2023-06-12 Mon 23:27] SCHEDULED: <2023-06-03 Sat>
****** DONE de nombreux reads sont perdus -> ok sur un SNV en alignant sur chromosome 22
CLOSED: [2023-05-29 Mon 15:38]
Problème dansr le BAM car sans les insertion
******* DONE Filtrer les reads sans pair : idem
CLOSED: [2023-05-23 Tue 00:01]
FOUND:Found 16662 unpaired mates
at htsjdk.samtools.SAMUtils.processValidationError(SAMUtils.java:470)
at picard.sam.SamToFastq.doWork(SamToFastq.java:224)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:308)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)
ex: chr22:42213078
On essaie
samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -hf 0x2 -o 63003856_chr22.bam
Ne rale pas...
SCHEDULED: <2023-05-22 Mon>
******* DONE Problème de la conversion BAM -> fastq ?
CLOSED: [2023-05-24 Wed 23:19]
******** DONE Test bamtofastq sur le BAM originel: idem
CLOSED: [2023-05-23 Tue 01:16]
Dans ~/recherche/code/bisonex/Bamscissors.jl
samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -hf 0x2 -o 63003856_chr22.bam
On trie bien par nom de read
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
Conversion en fastq. Attention à bien prendre le *fichier trié* !
samtools bam2fq -1 63003856_chr22_init1.fq.gz -2 63003856_chr22_init2.fq.gz -n 63003856_chr22_sorted.bam
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 2592110 reads
On envoie sur le mésocentre
scp 63003856_chr22_init{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors
On démarre un job avec 24 coeurs pour aller vite
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_init1.fq.gz 63003856_chr22_init2.fq.gz | samtools sort -@24 -o output.bam -
Dans /home/alex/recherche/bisonex/code/BamScissors.jl/aligned
❯ scp meso:/Work/Users/apraga/bisonex/tests/bamscissors/output.bam*
******** DONE Avec samtools fastq
CLOSED: [2023-05-23 Tue 01:16]
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -o 63003856_chr22_sorted.bam
scp 63003856_chr22_{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
mesocentre
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz | samtools sort -@24 -o output.bam
scp meso:/Work/Users/apraga/bisonex/tests/bamscissors/output.bam\* aligned/
******** DONE En rajoutant .fna dans le doserrie (+samtools fastq): idem
CLOSED: [2023-05-23 Tue 01:16]
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef* .
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef* .
bwa mem -t 24 genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz | samtools sort -@24 -o output.bam
******** DONE Test picard : idem
CLOSED: [2023-05-23 Tue 01:16]
Dans bisonex/code/BamScissors.jl
❯ picard SamToFastq -I 63003856_chr22_sorted.bam -F testpicard1.fastq -F2 testpicard2.fastq
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef testpicard1.fastq.gz testpicard2.fastq.gz | samtools sort -@24 -o testpicard.bam -
******** DONE Il ne manque pas des reads dans les fastq :
CLOSED: [2023-05-23 Tue 01:26]
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools flagstat aligned/output.bam
1306273 + 0 in total (QC-passed reads + QC-failed reads)
1296055 + 0 primary
0 + 0 secondary
10218 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
1304871 + 0 mapped (99.89% : N/A)
1294653 + 0 primary mapped (99.89% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools flagstat 63003856_chr22.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
2592110 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2592110 + 0 properly paired (100.00% : N/A)
2592110 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
❯ echo $(zcat 63003856_chr22_init2.fq.gz | wc -l)/4 | bc
1296055
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0 took 2s
❯ echo $(zcat 63003856_chr22_init1.fq.gz | wc -l)/4 | bc
1296055
❯ samtools flagstat 63003856_chr22_sorted.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
2592110 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2592110 + 0 properly paired (100.00% : N/A)
2592110 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
******** DONE Avec sort -O BAM idem
CLOSED: [2023-05-23 Tue 01:21]
******** DONE Singleton ou non mapped ? nonVirtPosition
CLOSED: [2023-05-23 Tue 01:22]
❯ samtools bam2fq -1 63003856_chr22_init1.fq.gz -2 63003856_chr22_init2.fq.gz -0 both -s single.fq.gz -n 63003856_chr22_sorted.bam
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ wc -l 63003856_chr22_init* both single.fq.gz
403540 63003856_chr22_init1.fq.gz
404788 63003856_chr22_init2.fq.gz
0 both
0 single.fq.gz
808328 total
******** Problème d'aligner car les reads sont bien dans le .fastq ?
Ex:
zgrep "A00853:477:HMLWYDSX3:2:2624:2826:18630" 63003856_chr22_{1,2}.fq.gz
63003856_chr22_1.fq.gz:@A00853:477:HMLWYDSX3:2:2624:2826:18630
63003856_chr22_2.fq.gz:@A00853:477:HMLWYDSX3:2:2624:2826:18630
******* DONE Refaire la manip sur bam chr22 non modifié + mail Alexis
CLOSED: [2023-05-23 Tue 23:27]
[1]_samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
Les reads sont bien tous mappé
samtools flagstat 63003856_chr22.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
[2] samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
flagstat idem
[3] samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 2592110 reads
Nombre de reads ok (= la moitié) et les séquences ne sont pas identiques pour le premier read (= on n'a pas 2x la même chose)
echo $(zcat 63003856_chr22_1.fq.gz | wc -l)/4 | bc
1296055
echo $(zcat 63003856_chr22_2.fq.gz | wc -l)/4 | bc
1296055
[4]
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz -o wtf.bam
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1752492 sequences (240000014 bp)...
[M::mem_pestat]
# candidate unique pairs for (FF, FR, RF, RR): (12, 702821, 0, 18)
[M::mem_pestat] analyzing insert size distribution for orientation FF...
[M::mem_pestat] (25, 50, 75) percentile: (54, 197, 268)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 696)
[M::mem_pestat] mean and std.dev: (163.92, 129.71)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 910)
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (129, 175, 233)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 441)
[M::mem_pestat] mean and std.dev: (185.56, 75.37)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 545)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation RR...
[M::mem_pestat] (25, 50, 75) percentile: (56, 192, 239)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 605)
[M::mem_pestat] mean and std.dev: (158.00, 110.30)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 788)
[M::mem_pestat] skip orientation FF
[M::mem_pestat] skip orientation RR
[M::process] read 839618 sequences (114952336 bp)...
[M::mem_process_seqs] Processed 1752492 reads in 375.714 CPU sec, 17.645 real sec
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 336379, 0, 5)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (128, 174, 232)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 440)
[M::mem_pestat] mean and std.dev: (184.73, 74.63)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 544)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 839618 reads in 183.039 CPU sec, 7.961 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t 24 -o wtf.bam /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz
[main] Real time: 38.278 sec; CPU: 565.821 sec
Bon nombre de reads pourtant
samtools flagstat wtf.bam
2611059 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
18949 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2611058 + 0 mapped (100.00% : N/A)
2592109 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2590970 + 0 properly paired (99.96% : N/A)
2592108 + 0 with itself and mate mapped
1 + 0 singletons (0.00% : N/A)
458 + 0 with mate mapped to a different chr
63 + 0 with mate mapped to a different chr (mapQ>=5)
$ samtools sort -@24 -o wtf_sorted.bam wtf.sam
[bam_sort_core] merging from 0 files and 24 in-memory blocks...
samtools flagstat wtf.bam
2611059 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
18949 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2611058 + 0 mapped (100.00% : N/A)
2592109 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2590970 + 0 properly paired (99.96% : N/A)
2592108 + 0 with itself and mate mapped
1 + 0 singletons (0.00% : N/A)
458 + 0 with mate mapped to a different chr
63 + 0 with mate mapped to a different chr (mapQ>=5)
Effectivement, ce n'est pas un problème d'IGV
❯ samtools mpileup 63003856_chr22.bam -r NC_000022.11:42213078-42213078
[mpileup] 1 samples in 1 input files
NC_000022.11 42213078 N 19 TTtTTtTTTtTTtTtTttt kkk_kkkkFkFF_FkFkQk
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools mpileup aligned/wtf_sorted.bam -r NC_000022.11:42213078-42213078
[mpileup] 1 samples in 1 input files
NC_000022.11 42213078 N 5 TTtTT _FkFF
******* DONE Regarder où ont été aligné les reads (nouveau run)
CLOSED: [2023-05-24 Wed 23:18]
******** DONE Préparation
CLOSED: [2023-05-24 Wed 21:59]
On relance le pipeline pour avoir un BAM propre
On garde les reads non mappé à partsir de la sortie d'applybqsr
#+begin_src sh
NXF_OPTS=-D"user.name=apraga" nextflow run main.nf -c nextflow.config -profile standard,helios --input="/Work/Projects/bisonex/centogene/fastq/2200467051_63003856/63003856_S135_R{1,2}_001.fastq.gz" --outdir=out -bg
cd out/63003856_S135_R/preprocessing/applybqsr/
samtools view 63003856_S135_R.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
make run BG= READS="out/63003856_S135_R/preprocessing/applybqsr/63003856_chr22_{1,2}.fq.gz"
cd out/63003856_chr22/preprocessing/mapped/
samtools index 63003856_chr22.bam
samtools mpileup 63003856_chr22.bam -r NC_000022.11:42213078-42213078
#+end_src
On récupère les 2 bam dans
#+begin_src
cd /home/alex/recherche/bisonex/code/BamScissors.jl/
rsync -avz meso:/Work/Users/apraga/bisonex/out/63003856_chr22/preprocessing/mapped data
rsync -avz meso:/Work/Users/apraga/bisonex/out/63003856_S135_R/preprocessing/applybqsr/ data/init/
#+end_src
******** Vérification que le reads est ailleurs
On cherche un read manquant dans le second alignement
#+begin_src sh
samtools view data/init/63003856_chr22.bam | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
#+end_src
#+RESULTS:
: A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NC_000022.11 42212845 0 151M = 42212883 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT ACC+FBCDCBBBAEAEDEEBBCCCECACBAEBEBDCCBCBFDCCCCFACEBEBCEEDCCCCFDCAEDCACBCEBBCFEACCFBDCACDCBCEBDBBCFEEDCCCFAFEACECCCECAEEDCADCBEDC7BEBCCCFBAFDCECCFBEAACA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
: A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NC_000022.11 42212883 0 151M = 42212845 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT AADECCCBDCBFCE<?CDEEEEBDEACDEAC;:BFBCBCDCCBEAEACAEFCCEAFBCBCCDEECBDBCECBEECCEACDEEBBFGDEFGCCFFFFCFCCEFBFDCFCDAAEBEE:CECBABBEBEE;DBFCCCDBCDBCCBBC?@BEEDA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
#+begin_src sh
samtools view data/mapped/63003856_chr22.bam | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
#+end_src
#+RESULTS:
: A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NW_014040930.1 115017 0 151M = 115055 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT ACC+FBCDCBBBAEAEDEEBBCCCECACBAEBEBDCCBCBFDCCCCFACEBEBCEEDCCCCFDCAEDCACBCEBBCFEACCFBDCACDCBCEBDBBCFEEDCCCFAFEACECCCECAEEDCADCBEDC7BEBCCCFBAFDCECCFBEAACA NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
: A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NW_014040930.1 115055 0 151M = 115017 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT AADECCCBDCBFCE<?CDEEEEBDEACDEAC;:BFBCBCDCCBEAEACAEFCCEAFBCBC
CDEECBDBCECBEECCEACDEEBBFGDEFGCCFFFFCFCCEFBFDCFCDAAEBEE:CECBAB
BEBEE;DBFCCCDBCDBCCBBC?@BEEDA NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
Effectivement, on aligne sur une zonne supprimée !
******* DONE Corriger la qualité: non
CLOSED: [2023-05-24 Wed 22:19]
******** DONE Comparaison avec le fastq de référénce : qualité !!
CLOSED: [2023-05-24 Wed 22:17]
#+begin_src sh
cd /Work/Users/apraga/bisonex/work/6e/8548fc90263830bf677f36585f11dc
zgrep -A 3 "A00853:477:HMLWYDSX3:1:1413:4390:28573" 63003856_chr22_1.fq.gz
#+end_src
@A00853:477:HMLWYDSX3:1:1413:4390:28573
AGGGTTACCACCACCACCCTGACAGGAGATATTCTAGGAGTACTCAAGAGCATCAGGGGATGGCTGGTAGCCTAGAAGGAACCACAAGGCCCAATGTCTTGGTTAGTCAAACCAATGAATTAGCTAGCAGGGGCCTTCTGAACAAAAGCAT
+
ADEEB@?CBBCCBDCBDCCCFBD;EEBEBBABCEC:EEBEAADCFCDFBFECCFCFFFFCCGFEDGFBBEEDCAECCEEBCECBDBCEEDCCBCBFAECCFEACAEAEBCCDCBCBFB:;CAEDCAEDBEEEEDC?<ECFBCDBCCCEDAA
#+begin_src
zgrep -A 3 "A00853:477:HMLWYDSX3:1:1413:4390:28573" /Work/Projects/bisonex/centogene/fastq/2200467051_63003856/63003856_S135_R1_001.fastq.gz
#+end_src
#+RESULTS:
: @A00853:477:HMLWYDSX3:1:1413:4390:28573 1:N:0:ATTCCACACA+TAGGCGATTG
AGGGTTACCACCACCACCCTGACAGGAGATATTCTAGGAGTACTCAAGAGCATCAGGGGATGGCTGGTAGCCTAGAAGGAACCACAAGGCCCAATGTCTTGGTTAGTCAAACCAATGAATTAGCTAGCAGGGGCCTTCTGAACAAAAGCAT
: +
: FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFF
******** DONE Regarder la qualité après bwa mem vs applybqsr: différente
CLOSED: [2023-05-24 Wed 22:18]
Sur le mésocentre, dans /Work/Users/apraga/bisonex/out/63003856_S135_R/preprocessing
$ samtools view mapped/63003856_S135_R.bam NC_000022.11 | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NC_000022.11 42212845 0 151M = 42212883 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT FFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NC_000022.11 42212883 0 151M = 42212845 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT FFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
samtools view applybqsr/63003856_S135_R.bam NC_000022.11 | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NC_000022.11 42212845 0 151M = 42212883 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT ACC+FBCDCBBBAEAEDEEBBCCCECACBAEBEBDCCBCBFDCCCCFACEBEBCEEDCCCCFDCAEDCACBCEBBCFEACCFBDCACDCBCEBDBBCFEEDCCCFAFEACECCCECAEEDCADCBEDC7BEBCCCFBAFDCECCFBEAACA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NC_000022.11 42212883 0 151M = 42212845 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT AADECCCBDCBFCE<?CDEEEEBDEACDEAC;:BFBCBCDCCBEAEACAEFCCEAFBCBCCDEECBDBCECBEECCEACDEEBBFGDEFGCCFFFFCFCCEFBFDCFCDAAEBEE:CECBABBEBEE;DBFCCCDBCDBCCBBC?@BEEDA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
******** DONE Réaligner à partir de la sortie de bwa mem
CLOSED: [2023-05-24 Wed 22:32]
#+begin_src sh
cd out/63003856_S135_R/preprocessing/mapped/
samtools view 63003856_S135_R.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
#+end_src
ON vérifie la qualité
#+begin_src
zgrep -A 3 "A00853:477:HMLWYDSX3:1:1413:4390:28573" 63003856_chr22_1.fq.gz
#+end_src
#+RESULTS:
: @A00853:477:HMLWYDSX3:1:1413:4390:28573
: AGGGTTACCACCACCACCCTGACAGGAGATATTCTAGGAGTACTCAAGAGCATCAGGGGATGGCTGGTAGCCTAGAAGGAACCACAAGGCCCAATGTCTTGGTTAGTCAAACCAATGAATTAGCTAGCAGGGGCCTTCTGAACAAAAGCAT
: +
: FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFF
#+begin_src
NXF_OPTS=-D"user.name=apraga" nextflow run main.nf -c nextflow.config -profile standard,helios -resume --input="out/63003856_S135_R/preprocessing/mapped/63003856_chr22_{1,2}.fq.gz" --outdir=out/63003856_chr22-from-mapped
#+end_src
Puis ::
#+begin_src
cd /Work/Users/apraga/bisonex/out/63003856_chr22-from-mapped/63003856_chr22/preprocessing/mapped
samtools view 63003856_chr22.bam | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
#+end_src
#+RESULTS:
: A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NW_014040930.1 115017 0 151M = 115055 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT FFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
: A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NW_014040930.1 115055 0 151M = 115017 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT FFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
******* DONE Aligner sur génome de référence limité au chromosome 22
CLOSED: [2023-05-24 Wed 23:18]
******** KILL Test données non modifiées
CLOSED: [2023-05-24 Wed 23:18]
/Work/Users/apraga/bisonex/tests/bamscissors
#+begin_src
cd /Work/Groups/bisonex/data/genome/GRCh38.p13/
mkdir chr22/
samtools faidx genomeRef.fna NC_000022.11 > chr22/chr22.fna
cd chr22
samtools faidx chr22.fna
bwa index chr22.fna
#+end_src
#+begin_src
cd /Work/Users/apraga/bisonex/tests/bamscissors
ln -s ../ ../out/63003856_S135_R/preprocessing/applybqsr/63003856_chr22_{1,2}.fq.gz .
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/chr22/chr22.fna 63003856_chr22_1.fq.gz 63003856_chr22_1.fq.gz -o smallref.sam
#+end_src
******** DONE Test données modifiées: ok
CLOSED: [2023-05-24 Wed 23:18]
Données dans data/init
#+begin_src sh
time julia insertVariant.jl
rsync -avz data/init/*.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
#+end_src
#+begin_src
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/chr22/chr22.fna 63003856_chr22_1.fq.gz 63003856_chr22_1.fq.gz | samtools sort -@24 - -o smallref.bam
#+end_src
#+begin_src
rsync -avz meso:/Work/Users/apraga/bisonex/tests/bamscissors/smallref.bam mapped/
#+end_src
****** DONE Test haplotypecaller 1 variant
CLOSED: [2023-05-29 Mon 15:38]
****** DONE Test haplotypecaller tous les variants
****** DONE Comprendre pourquoi la répartiton ne suit pas la loi normale
CLOSED: [2023-06-01 Thu 21:44]
Certains hétérozygote soint à 0.01 ou 1...
******* DONE augmenter le nombre d'échantillions: idem
CLOSED: [2023-05-31 Wed 22:24]
******* DONE Vérifier le nombre de reads marqué vs édité
CLOSED: [2023-06-01 Thu 21:44]
******* DONE vérifier que 100 appel à rand(d, 1)[1] est semblable à un appel de rand(d, 100)
CLOSED: [2023-05-31 Wed 22:24]
julia> df = vcat(DataFrame(:y => z, :type => "z"), DataFrame(:y => y, :type => "y"));
julia> y = [rand(d, 1)[1] for x in 1:1000];
julia> z = rand(d,1000);
julia> df = vcat(DataFrame(:y => z, :type => "z"), DataFrame(:y => y, :type => "y"));
draw(data(df)*histogram(bins=100)*mapping(:y, color=:ty
pe,dodge=:type))
****** DONE Améliorer les performances
CLOSED: [2023-06-02 Fri 23:39]
#+begin_src julia
@time include("xamscissors.jl")
#+end_src
430s pour chromosome 22. Majorité dans l'édition de reads:
******* DONE Inserér tous les variants d'un reads d'un coup
CLOSED: [2023-06-01 Thu 23:09]
Ne change rien
******* DONE Test avec -t4: idem
CLOSED: [2023-06-01 Thu 23:17]
******* DONE Test mésocentre : idem
CLOSED: [2023-06-01 Thu 23:40]
348s
******* Changer la structure de données des
Dataframe -> dict = les performances horribles ont disparuse
****** KILL Refaire le test avec la nouvelle version
CLOSED: [2023-06-12 Mon 23:27]
******* DONE Génération des données
CLOSED: [2023-06-02 Fri 23:40]
Mésocentre
#+begin_src sh
cd /Work/Users/apraga/bisonex/out/63003856_S135_R/preprocessing/mapped
samtools view 63003856_S135_R.bam NC_000022.11 -o 63003856_S135_R_chr22.bam
samtools index 63003856_S135_R_chr22.bam
cp 63003856_S135_R_chr22.bam* /Work/Users/apraga/bisonex/tests/xamscissors/
cd /Work/Users/apraga/bisonex/tests/xamscissors
#+end_src
On génère les données
#+begin_src julia
using XAMScissors
insertSNV("./63003856_S135_R_chr22.bam", "snvs_chr22.csv", "out")
#+end_src
Puis
#+begin_src sh
julia xamscissors.jl
#+end_src
******* DONE Run
CLOSED: [2023-06-03 Sat 18:26]
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/runInserted.nf -profile standard,helios --input="tests/xamscissors/out/inserted_{1,2}.fq.gz"
******* DONE Après haplotypecaller : ok
CLOSED: [2023-06-03 Sat 18:27] SCHEDULED: <2023-06-03 Sat>
******* KILL Après filtre vep
CLOSED: [2023-06-12 Mon 23:27] SCHEDULED: <2023-06-03 Sat>
***** KILL PHase 3 : tous les SNV
CLOSED: [2023-06-12 Mon 23:28] SCHEDULED: <2023-06-03 Sat>
****** DONE Générer les données
CLOSED: [2023-06-03 Sat 20:16] SCHEDULED: <2023-06-03 Sat>
#+begin_src julia
using XAMScissors
insertSNV("../../out/63003856_S135_R/preprocessing/mapped/63003856_S135_R.bam", "snvs.csv", "out")
#+end_src
temps d'exécution 73min
#+begin_src sh
nohup bash -c 'time julia xamscissors.jl' &
xamscissors-63003856/*.fq.gz /Work/Groups/bisonex/data/xamscissors/
#+end_src
****** DONE Regénérer avec @time pour avoir les performaces
CLOSED: [2023-06-03 Sat 21:45]
markReads 6.265202 seconds (1.36 M allocations: 137.090 MiB, 1.00% gc time, 9.79% compilation time)
editReads 1327.701623 seconds (1.03 G allocations: 81.996 GiB, 0.59% gc time, 0.03% compilation time)
samtools index 117.743727 seconds (53 allocations: 1.883 KiB)
samtools sort 2820.074930 seconds (66 allocations: 2.789 KiB)
bam2fastq 134.148952 seconds (794 allocations: 40.539 KiB, 0.01% compilation time)
real 73m33.273s
user 77m38.194s
sys 1m26.684s
[bam_sort_core] merging from 60 files and 1 in-memory blocks...
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 126905130 reads
real 73m6.934s
user 77m25.397s
sys 1m21.339s
****** KILL Après haplotypecaller 556/590, majorité = échec alignement
CLOSED: [2023-06-12 Mon 23:27] SCHEDULED: <2023-06-04 Sun>
Haplotypecaller 556 found over 590
Amongst 34 missed variant, 2 have a mapping quality > 0
2×7 DataFrame
Row │ chrom pos ref alt zygosity meanQual stdQual
│ String15 Int64 String1 String1 String7 Float64 Float64
─────┼─────────────────────────────────────────────────────────────────────────
1 │ NC_000017.11 39672244 G A het 60.0 0.0
2 │ NC_000001.11 155235252 A G het 0.258065 2.48868
NC_000017.11 39672244 G A het => ok, problème de représentation car 2 variant côte à cote
NC_000001.11 155235252 A G het => peu de reads alternatifs (9/93 donc ok)
Position: chromoe 1 et 6 surtout
34×7 DataFrame
Row │ chrom pos ref alt zygosity
│ String15 Int64 String1 String1 String7
─────┼──────────────────────────────────────────────────────
1 │ NC_000001.11 153817496 A T het
2 │ NC_000001.11 155235252 A G het
3 │ NC_000001.11 155236268 G A het
4 │ NC_000001.11 155290591 C T het
5 │ NC_000001.11 155291918 G A het
6 │ NC_000001.11 155294358 G T het
7 │ NC_000002.12 149010343 C T het
8 │ NC_000006.12 32039426 T A het
9 │ NC_000006.12 32040110 G T het
10 │ NC_000006.12 32040723 G A het
11 │ NC_000006.12 32041006 C T het
12 │ NC_000006.12 32041147 G A het
13 │ NC_000006.12 33443054 G T het
14 │ NC_000006.12 33451815 C T het
15 │ NC_000006.12 170283230 C A het
16 │ NC_000006.12 170283754 G A het
17 │ NC_000006.12 170285637 T C het
18 │ NC_000006.12 170289678 A C het
19 │ NC_000010.11 87961118 G A het
20 │ NC_000012.12 2449086 C G het
21 │ NC_000015.10 74343027 C T het
22 │ NC_000016.10 16163078 G A het
23 │ NC_000016.10 21262032 C G het
24 │ NC_000016.10 21962506 C T homo
25 │ NC_000017.11 7513122 C T het
26 │ NC_000017.11 7513752 C T het
27 │ NC_000017.11 39672244 G A het
28 │ NC_000017.11 46018710 C T het
29 │ NC_000019.10 54144058 G A het
30 │ NC_000021.9 43063074 A G het
31 │ NC_000021.9 43426167 C T het
32 │ NC_000022.11 18918421 A G het
33 │ NC_000022.11 42087168 T A homo
34 │ NC_000022.11 42213078 T G het
****** DONE Voir où est l'alignement alternatif: sur NW_ (zone supprimée)
CLOSED: [2023-06-04 Sun 22:15] SCHEDULED: <2023-06-04 Sun>
ex chr15 74343027
A00853:477:HMLWYDSX3:2:2444:22354:28870
#+begin_src
cd /Work/Groups/bisonex/data/xamscissors
zgrep -A4 "A00853:477:HMLWYDSX3:2:2444:22354:28870" *.fq.gz
#+end_src
63003856_xamscissors_1.fq.gz:@A00853:477:HMLWYDSX3:2:2444:22354:28870
63003856_xamscissors_1.fq.gz:CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTTGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC
63003856_xamscissors_1.fq.gz:+
63003856_xamscissors_1.fq.gz:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFF:FFFFFFFFFF::FFFFFFFFFFF:FFFFFFFFFFFFFF:FFFFFFF,FFFFFF,FFFFFFFFFFFF:FF::FF
63003856_xamscissors_2.fq.gz:@A00853:477:HMLWYDSX3:2:2444:22354:28870
63003856_xamscissors_2.fq.gz:GACAGAAAGGAAGTGTTCACCACGATTACCGTGGCATCCTCTACAGACTCCTGGGAGACAGCAAGATGTCCTTCGAGGACATCAAGGCCAACGTCACAGAGATGCCGGCAGGAGGGGTGGACACGGTG
63003856_xamscissors_2.fq.gz:+
63003856_xamscissors_2.fq.gz:FF:FFF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFF:F:FF:FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFF,:FFF,FFFFFF:FFFFFFFFFFFFF
******* DONE Avec BLAT: sur _fix
CLOSED: [2023-06-04 Sun 21:07]
1er =
ACTIONS QUERY SCORE START END QSIZE IDENTITY CHROM STRAND START END SPAN
--------------------------------------------------------------------------------------------
------------------
browser details YourSeq 124 1 128 128 98.5% chr15_ML143370v1_fix + 172243 172370 128 What is chrom_fix?
browser details YourSeq 124 1 128 128 98.5% chr15 + 74342974 74343101 128
browser details YourSeq 23 1 25 128 96.0% chr19 - 33396097 33396121 25
Second
--------------------------------------------------------------------------------------------------------------
browser details YourSeq 126 1 128 128 99.3% chr15_ML143370v1_fix - 172243 172370 128 What is chrom_fix?
browser details YourSeq 126 1 128 128 99.3% chr15 - 74342974 74343101 128
browser details YourSeq 23 104 128 128 96.0% chr19 + 33396097 33396121 25
******* DONE Bwa mem à la main GRCh38.p13 : on est dans une zone NW
CLOSED: [2023-06-04 Sun 21:51]
On met les 2 reads dans des fichiers séparés puis
#+begin_src sh
cd /Work/Users/apraga/bisonex/tests/xamscissors/align
bwa mem /Work/Groups/bisonex/data/genome/GRCh38.p13/bwa/genomeRef test1.fq test2.fq
#+end_src
A00853:477:HMLWYDSX3:2:2444:22354:28870 97 NW_021160016.1 172243 0 128M = 172243 128 CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTTGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFF:FFFFFFFFFF::FFFFFFFFFFF:FFFFFFFFFFFFFF:FFFFFFF,FFFFFF,FFFFFFFFFFFF:FF::FF NM:i:2 MD:Z:22A30C7MC:Z:128M AS:i:118 XS:i:118 XA:Z:NC_000015.10,+74342974,128M,2;
A00853:477:HMLWYDSX3:2:2444:22354:28870 145 NW_021160016.1 172243 0 128M = 172243 -128 CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTCGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC FFFFFFFFFFFFF:FFFFFF,FFF:,FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFF:FF:F:FFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FF:FFF:FF NM:i:1 MD:Z:22A105 MC:Z:128M AS:i:123 XS:i:123 XA:Z:NC_000015.10,-74342974,128M,1;
******* DONE GRCh38.p14: idem
CLOSED: [2023-06-04 Sun 21:51]
A00853:477:HMLWYDSX3:2:2444:22354:28870 97 NW_021160016.1 172243 0 128M = 172243 128 CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTTGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFF:FFFFFFFFFF::FFFFFFFFFFF:FFFFFFFFFFFFFF:FFFFFFF,FFFFFF,FFFFFFFFFFFF:FF::FF NM:i:2 MD:Z:22A30C7MC:Z:128M AS:i:118 XS:i:118 XA:Z:NC_000015.10,+74342974,128M,2;
A00853:477:HMLWYDSX3:2:2444:22354:28870 145 NW_021160016.1 172243 0 128M = 172243 -128 CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTCGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC FFFFFFFFFFFFF:FFFFFF,FFF:,FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFF:FF:F:FFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FF:FFF:FF NM:i:1 MD:Z:22A105 MC:Z:128M AS:i:123 XS:i:123 XA:Z:NC_000015.10,-74342974,128M,1;
******* DONE GRCh38 : ok
CLOSED: [2023-06-04 Sun 22:15]
bwa mem /Work/Projects/bisonex/data/genome/GRCh38/GCA_000001405.15_GRCh38_full_analysis_set.fna test1.fq test2.fq
****** DONE Vérifier que les reads ont la même qualité sur les fichiers d'origine: oui
CLOSED: [2023-06-04 Sun 21:07]
****** DONE Supprimer les NW_ ?
CLOSED: [2023-06-10 Sat 10:40] SCHEDULED: <2023-06-04 Sun>
@A00853:477:HMLWYDSX3:3:2114:14742:8860
CAGGCCAGCCGCTCAGCCCGCTCCTTTCACCCTCTGCAGGAGAGCCTCGTGGCAGGCCAGTGGAGGGACATGATGGACTACATGCTCCAAGGGGTGGCGCAGCCGAGCATGGAAGAGGGCTCTGGACAGCTCCTGGAAGGGCACTTGCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@A00853:477:HMLWYDSX3:3:2114:14742:8860
CTTTTGCTTGTCCCCAGGACGCACCTCAGGGTGGTGAAGCAAAAAAACCACGGCCCAGGAGAGGGTGGGTGCTGTGGTCTCAGTGCCACCGATCAGGAGGTCCACTGCAGCCATGTGCAAGTGCCCTTCCAGGAGCTGTCCAGAGCCCTCT
+
FFFFFFFFFFFFFFFFFFFFFFF:FFF:FFFFFFFFFFFFF,FFFFFFFFFFFF:F:FFFF:FFFFF,,FFF:FFFFFFFFFF,FFFFFFF,FFFFFFFFFFF,FFFFFFFFF:FFFF,F:FFFFF:FFFFFFFFF:FFFF,FFFFFFFFF
****** DONE Supprimer NW_ et NT_
***** TODO Phase 2 : chr22, vaf variable :T2T:
SCHEDULED: <2023-06-16 Fri>
CLOSED: [2023-06-12 Mon 23:27]
***** TODO Phase 3 : tous SNV, vaf variable :T2T:
SCHEDULED: <2023-06-16 Fri>
CLOSED: [2023-06-12 Mon 23:27]
**** TODO Test Indel
*** Divers
**** DONE Vérifier nombre de reads fastq - bam
CLOSED: [2022-10-09 Sun 22:31]
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 HG003 :hg003:hg38:
CLOSED: [2023-04-16 Sun 00:20]
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run main.nf -profile standard,helios --input /Work/Groups/bisonex/data/giab/GRCh38/HG003_{1,2}.fq.gz -bg
#+end_src
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/compareVCF.nf -profile standard,helios -resume --outdir=compareHG003 --test.id=HG003 --test.query=out/HG003_1/variantCalling/haplotypecaller/HG003_1.vcf.gz --test.compare=vcfeval,happy --test.capture=data/AgilentSureSelectv05_hg38.bed
#+end_src
vcfeval
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
5.000 36745 36473 486 3988 0.9869 0.9021 0.9426
None 36748 36476 495 3985 0.9866 0.9022 0.9425
$ zcat NA12878.snp_roc.tsv.gz | tail -n 1 | awk '{print $7 $6}'
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 2731 2290 441 3092 208 577 62 53 0.838521 0.917296 0.186611 0.876141 NaN NaN 1.505145 1.888993
INDEL PASS 2731 2290 441 3092 208 577 62 53 0.838521 0.917296 0.186611 0.876141 NaN NaN 1.505145 1.888993
SNP ALL 37997 34481 3516 36861 306 2074 33 13 0.907466 0.991204 0.056265 0.947488 2.611269 2.565915 1.555780 1.621727
SNP PASS 37997 34481 3516 36861 306 2074 33 13 0.907466 0.991204 0.056265 0.947488 2.611269 2.5659
**** DONE HG004 :hg38:
CLOSED: [2023-04-16 Sun 00:20]
#+begin_src sh
NXF_OPTS=-D"user.name=${USER}" nextflow run main.nf -profile standard,helios --input /Work/Groups/bisonex/data/giab/GRCh38/HG004_{1,2}.fq.gz -bg
#+end_src
vcfeval
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
6.000 36938 36678 421 4040 0.9887 0.9014 0.9430
None 36942 36682 432 4036 0.9884 0.9015 0.9429
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 2787 2388 399 3183 195 580 53 38 0.856835 0.925086 0.182218 0.889654 NaN NaN 1.507834 1.848649
INDEL PASS 2787 2388 399 3183 195 580 53 38 0.856835 0.925086 0.182218 0.889654 NaN NaN 1.507834 1.848649
SNP ALL 38185 34560 3625 36921 254 2107 46 7 0.905067 0.992704 0.057068 0.946862 2.589175 2.553546 1.632595 1.653534
SNP PASS 38185 34560 3625 36921 254 2107 46 7 0.905067 0.992704 0.057068 0.946862 2.589175 2.553546 1.632595 1.653534
**** DONE Résumer résultats pour Paul + article :resultats:hg38:
CLOSED: [2023-04-06 Thu 21:41] SCHEDULED: <2023-04-02 Sun>
**** DONE Plot : ashkenazim trio :hg38:
CLOSED: [2023-04-18 Tue 21:27] SCHEDULED: <2023-04-16 Sun>
/Entered on/ [2023-04-16 Sun 17:29]
**** TODO NA12878 :na12878:T2T:
SCHEDULED: <2023-06-15 Thu>
*** KILL Platinum genome
CLOSED: [2023-06-14 Wed 22:37]
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 Insilico :centogene:
*** TODO tous les variants centogène
**** DONE Extraire liste des SNVs
CLOSED: [2023-04-22 Sat 17:32] SCHEDULED: <2023-04-17 Mon>
***** DONE Corriger manquant à la main
CLOSED: [2023-04-22 Sat 17:31]
La sortie est sauvegardé dans git-annex : variants_success.csv
***** DONE Automatique
CLOSED: [2023-04-22 Sat 17:31]
**** DONE Convert SNVs : transcript -> génomique
CLOSED: [2023-06-03 Sat 17:16]
***** DONE Variant_recoder
CLOSED: [2023-04-26 Wed 21:21] SCHEDULED: <2023-04-22 Sat>
****** KILL Haskell: 160 manquant : recoded-success.csv
CLOSED: [2023-04-25 Tue 18:32]
La liste des variants a été générée en Haskel l et nettoyée à la main.
On générer une liste de variant pour variant_rec oder et on soumet tout d'un coup.
[[file:~/recherche/bisonex/parsevariants/app/Main.hs][parsevariant]]
#+begin_src haskell
recodeVariant = do
prepareVariantRecod er "variant_success.csv" "renamed.csv"
runVariantRecoder "renamed.csv" "recoded.json"
#+end_src
#+RESULTS:
: <interactive>:4:3-19: error:
: Variable not in scope: runVariantRecoder :: String -> String -> t
: gh
Problème : 160 n'ont pas pu être lu sur 820, probablement à cause du numéro mineur de transcrit
La sortie est sauvegardé dans git-annex : variants-recoded-raw.json.
****** KILL Julia
CLOSED: [2023-04-25 Tue 18:32]
On regénère la liste de variant et on passe à Julia pour préparer l'appel en parallèle à variant recoder
[[file:~/recherche/bisonex/parsevariants/variantRecoder.jl][variantRecoder.jl]]
#+begin_src julia
setupVariantRecoder(unique(init), n)
#+end_src
Puis
#+begin_src sh
parallel -a parallel-recoder.sh --jobs 10
#+end_src
On récupère les résultats
#+begin_src julia
(fails, success) = mergeVariantRecoder(n)
CSV.write(fSuccess, success)
CSV.write(fFailures, fails)
#+end_src
Certains variants ne sont pas trouvé, donc on prépare un nouveau job en enlevant les versionrs mineures des transcrits
#+begin_src julia
# Cleanup json and txt
if isfile(fSuccess) && isfile(fFailures)
foreach(rm, variantRecoderInput())
foreach(rm, variantRecoderOutput())
end
redoFails(fFailures)
#+end_src
Puis
#+begin_src sh
parallel -a parallel-recoder.sh --jobs 3
#+end_src
Il manque encore 70 transcrits
***** DONE Julia avec mobidetails: recode-failures-mobidetails.csv
CLOSED: [2023-04-25 Tue 18:58]
Nouvelle stratégie : on essaie une fois variant recoder.
Pour tous les échecs, on utilise mobidetails (~170).
Si l'ID n'est pas trouvé, on incrémente le numéro de version 2 fois
***** DONE Reste une dizaine à corriger à la main
CLOSED: [2023-04-26 Wed 21:21]
- [X] certains transcrits ont juste été supprimé
- [X] Erreur de parsing, manque souvent un -
#+begin_src julia
lastTryMobidetails("recoded-failures-mobidetails.csv")
#+end_src
***** DONE Fusionner données
CLOSED: [2023-04-26 Wed 22:35]
#+begin_src julia
function mergeAllGenomic()
dNew = mergeAll("recoded-success.csv",
"recoded-failures-mobidetails.csv",
"recoded-failures-mobidetails-redo.csv")
dInit = @chain DataFrame(CSV.File("variant_success.csv")) begin
@transform :transcript = :transcript .* ":" .* :coding .* :codingPos .* :codingChange
@select :file :transcript :classification :zygosity
@rename :classificationCentogene = :classification
end
dTmp = outerjoin(dInit, dNew, on = :transcript)
CSV.write("variant_genomic.csv", dTmp)
end
fSuccess = "recoded-success.csv"
fFailures = "recoded-failures.csv"
# variantRecoder(fSuccess, fFailures)
# mobidetailsOnFailures(fFailures)
# lastTryMobidetails("recoded-failures-mobidetails.csv")
mergeAllGenomic()
#+end_src
***** DONE Formatter donner pour simuscop
CLOSED: [2023-04-28 Fri 11:55] SCHEDULED: <2023-04-26 Wed>
**** TODO Extraire liste des CNVs
SCHEDULED: <2023-04-17 Mon>
**** TODO Simuscop :simuscop:
***** DONE Entrainer le modèle sur 63003856/
CLOSED: [2023-04-29 Sat 19:56]
Relancer le modèle pour être sûr
***** DONE Générer fastq avec simuscop (del et ins seulement) 20x
CLOSED: [2023-04-28 Fri 23:35] SCHEDULED: <2023-04-22 Sat>
****** DONE Génerer un profile avec bed de centogène
CLOSED: [2023-04-28 Fri 11:54] SCHEDULED: <2023-04-22 Sat>
NA12878 mais à refaire avec un vrai séquencage
Voir [[*Centogène][Bed Centogène]] pour choix
****** DONE Générer les données en 20x
CLOSED: [2023-04-28 Fri 11:54] SCHEDULED: <2023-04-22 Sat>
capture de centogene
****** DONE Regénérer en supprimant les doublons
CLOSED: [2023-04-28 Fri 17:28]
***** DONE Quelle couverture ?
CLOSED: [2023-04-29 Sat 18:26]
ex sur chr11:16,014,966 où on a 11 reads dans la simulation contre 200 !
****** 200 est la plus proche
#+attr_html: :width 500px
[[./simuscop-200-chr1-1.png]]
#+attr_html: :width 500px
[[./simuscop-200-chr1-2.png]]
****** DONE 20x
CLOSED: [2023-04-29 Sat 15:38]
****** DONE 50x
CLOSED: [2023-04-29 Sat 15:38]
****** DONE 100x
CLOSED: [2023-04-29 Sat 15:39]
****** DONE 200x
CLOSED: [2023-04-29 Sat 15:39]
***** DONE Reads mal centrés sur des petits exons seuls
CLOSED: [2023-04-29 Sat 19:56] SCHEDULED: <2023-04-29 Sat>
Capture ok : [[https://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr1%3A153817168%2D153817824&hgsid=296556270_F4fkENLPXHXidi2oALXls2jxNH9l][UCSC]] (track noire)
Mais mauvaise répartitiopn
#+attr_html: :width 800px
[[./simuscop-error.png]]
À tester
- Problème de profile ?
- mauvais patient ?
- mauvaise génération ? -> comparer avec ceux donnés sur github
- nom des chromosomes ?
****** DONE [#A] Tester sur exon 6 GATAD2B pour NC_000001.11:g.153817496A>T
CLOSED: [2023-04-29 Sat 19:56] SCHEDULED: <2023-04-29 Sat>
******* DONE Configuration + Profile 63003856.profile: idem, mal centré
CLOSED: [2023-04-29 Sat 19:18]
Téléchargement des données
#+begin_src sh :dir ~/code/bisonex/test-simuscop
scp meso:/Work/Projects/bisonex/data/genome/GRCh38.p14/genomeRef.fna .
scp meso:Work/Projects/bisonex/data/simuscop/*.profile .
scp -r meso:/Work/Projects/bisonex/data/genome/GRCh38.p13/bwa .
#+end_src
On récupère l'exon (NB: org-mode ne lance pas le code...)
#+begin_src julia
using CSV,DataFramesMeta
d = CSV.read("VCGS_Exome_Covered_Targets_hg38_40.1MB_renamed.bed", header=false, delim="\t", DataFrame)
@subset d :Column1 .== "NC_000001.11" :Column2 .<= 153817496 :Column3 .>= 153817496
#+end_src
NC_000001.11 153817371 153817542
Génération du bed
#+begin_src sh :dir ~/code/bisonex/test-simuscop
echo -e "NC_000001.11\t153817371\t153817542" > gatad2b-exon6.bed
#+end_src
#+RESULTS:
Génération d'un variant
#+begin_src sh :dir ~/code/bisonex/test-simuscop
echo -e "s\tsingle\tNC_000001.11\t153817496\tA\tT\thet"> variant.txt
#+end_src
#+RESULTS:
Génération du fichier de config
#+begin_src sh :dir ~/code/bisonex/test-simuscop
cat > config_wes.txt << EOL
ref = genomeRef.fna
profile = ./63003856.profile
variation = ./variant.txt
target = ./gatad2b-exon6.bed
layout = PE
threads = 1
name = single
output = test-gatad2b
coverage = 20
EOL
#+end_src
#+RESULTS:
On démarre la simulation
#+begin_src sh :dir ~/code/bisonex/test-simuscop
simuReads config_wes.txt
#+end_src
#+RESULTS:
Alignement
#+begin_src sh :dir ~/code/bisonex/test-simuscop
bwa mem -R '@RG\tID:sample\tSM:sample\tPL:ILLUMINA\tPM:Miseq\tCN:lol\tLB:definition_to_add' bwa/genomeRef test-gatad2b/single_1.fq test-gatad2b/single_2.fq | samtools sort -o single.bam
#+end_src
#+RESULTS:
******* DONE Profile github HiSeq2000
CLOSED: [2023-04-29 Sat 19:56]
#+begin_src sh :dir ~/code/bisonex/test-simuscop :result file
wget https://raw.githubusercontent.com/qasimyu/simuscop/master/testData/Illumina_HiSeq2000.profile
#+end_src
#+RESULTS:
#+begin_src sh :dir ~/code/bisonex/test-simuscop
cat > config_wes.txt << EOL
ref = genomeRef.fna
profile = ./Illumina_HiSeq2000.profile
variation = ./variant.txt
target = ./gatad2b-exon6.bed
layout = PE
threads = 1
name = single
output = test-gatad2b-hiseq2000
coverage = 20
EOL
simuReads config_wes.txt
bwa mem -R '@RG\tID:sample\tSM:sample\tPL:ILLUMINA\tPM:Miseq\tCN:lol\tLB:definition_to_add' bwa/genomeRef test-gatad2b-hiseq2000/single_1.fq test-gatad2b-hiseq2000/single_2.fq | samtools sort -o single-hiseq2000.bam
samtools index single-hiseq2000.bam
#+end_src
#+RESULTS:
******* KILL Tester exemple sur github
CLOSED: [2023-04-29 Sat 19:56]
#+begin_src sh
git clone https://github.com/qasimyu/simuscop/
cd simuscop
simuReads configFiles/config_test_wes.txt
#+end_src
******* KILL Centrer la fenêtre sur les zones de capture
CLOSED: [2023-04-30 Sun 13:28] SCHEDULED: <2023-04-29 Sat>
1000bp par défaut, ce qui est plus grand que les zones de captures...
Changer fragzip ne fonctionne pas
Si on rajoute un offset sur l'exon: 200bp, est encore plus allongé
NC_000001.11 153817371 153817542 ->
NC_000001.11 153817171 153817742
Si on désactive les target ?
Regarder les target sur le chromosome 1
#+begin_src sh :dir ~/code/bisonex/test-simuscop :results silent
scp meso:/Work/Projects/bisonex/data/simuscop/VCGS_Exome_Covered_Targets_hg38_40.1MB_renamed.bed .
#+end_src
#+begin_src sh :dir ~/code/bisonex/test-simuscop :results silent
head -n 100 VCGS_Exome_Covered_Targets_hg38_40.1MB_renamed.bed > 100exons.bed
echo -e "s\tsingle\tNC_000001.11\t153817496\tA\tT\thet"> variant.txt
cat > config_wes.txt << EOL
ref = genomeRef.fna
profile = ./63003856.profile
variation = ./variant.txt
layout = PE
threads = 4
target = 100exons.bed
name = single
output = test-gatad2b
coverage = 200
EOL
./simuscop/bin/simuReads config_wes.txt
bwa mem bwa/genomeRef test-gatad2b/single_1.fq test-gatad2b/single_2.fq | samtools sort -o single.bam
samtools index single.bam
#+end_src
***** KILL Vérifier tous les variants sont retrouvés en 200x: hg38
CLOSED: [2023-06-12 Mon 23:25]
****** DONE Après alignement
CLOSED: [2023-04-29 Sat 18:27] SCHEDULED: <2023-04-28 Fri>
******* DONE SNV: avec doublons
CLOSED: [2023-04-28 Fri 18:12]
On utilise [[file:~/recherche/bisonex/simuscop/checkBam.jl][checkBam.jl]]
#+begin_src julia
d = prepareVariant("../parsevariants/variant_genomic.csv")
root = "/home/alex/code/bisonex/simuscop-centogene/cento"
bam = root * "/preprocessing/applybqsr/cento.bam"
bai = root * "/preprocessing/recalibrated/cento.bam.bai"
snv = getSNV(d, bam, bai)
#+end_src
Nombreux faux homozygouteS
Vérification avec checkFalseHemizygous(snv) : nombreux doublons dans le fichier pour simuscop...
******* DONE SNV sans doublons
CLOSED: [2023-04-29 Sat 18:27]
******** DONE 18 faux homozygote mais avec peu de reads
CLOSED: [2023-04-29 Sat 18:27]
julia> @subset snv :refCount .== 0 :altCount .> 0 :zygosity .== "heterozygous"
18×10 DataFrame
Row │ chrom pos variant variantType zygosity ref alt refCount altCount readsCount
│ SubStrin…? Int64 SubStrin…? String? String15 SubStrin… SubStrin… Int64 Int64 Int64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ NC_000022.11 42213078 g.42213078T>G snv heterozygous T G 0 1 1
2 │ NC_000012.12 101680427 g.101680427C>A snv heterozygous C A 0 3 3
3 │ NC_000014.9 105385684 g.105385684G>C snv heterozygous G C 0 4 4
4 │ NC_000011.10 125978299 g.125978299C>T snv heterozygous C T 0 3 3
5 │ NC_000023.11 77998618 g.77998618C>T snv heterozygous C T 0 2 2
6 │ NC_000015.10 66703292 g.66703292C>T snv heterozygous C T 0 3 3
7 │ NC_000010.11 87961118 g.87961118G>A snv heterozygous G A 0 3 3
8 │ NC_000012.12 112477719 g.112477719A>G snv heterozygous A G 0 2 2
9 │ NC_000020.11 6778406 g.6778406C>T snv heterozygous C T 0 3 3
10 │ NC_000023.11 68192943 g.68192943G>A snv heterozygous G A 0 2 2
11 │ NC_000004.12 987858 g.987858C>T snv heterozygous C T 0 3 4
12 │ NC_000015.10 66435145 g.66435145G>A snv heterozygous G A 0 1 2
13 │ NC_000002.12 47809595 g.47809595C>T snv heterozygous C T 0 2 2
14 │ NC_000003.12 136477305 g.136477305C>G snv heterozygous C G 0 4 4
15 │ NC_000005.10 157285458 g.157285458C>T snv heterozygous C T 0 3 3
16 │ NC_000012.12 23604413 g.23604413T>G snv heterozygous T G 0 5 5
17 │ NC_000019.10 52219703 g.52219703C>T snv heterozygous C T 0 1 1
18 │ NC_000016.10 88856757 g.88856757C>T snv heterozygous C T 0 8 8
******** DONE 8 non retrouvé => probablement hors de la zjone de capture
CLOSED: [2023-04-28 Fri 19:49]
julia> @subset snv :refCount .== 0 :altCount .== 0
8×10 DataFrame
Row │ chrom pos variant variantType zygosity ref alt refCount altCount readsCount
│ SubStrin…? Int64 SubStrin…? String? String15 SubStrin… SubStrin… Int64 Int64 Int64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ NC_000015.10 74343027 g.74343027C>T snv heterozygous C T 0 0 0
2 │ NC_000011.10 20638345 g.20638345A>G snv heterozygous A G 0 0 0
3 │ NC_000004.12 139370252 g.139370252C>T snv heterozygous C T 0 0 2
4 │ NC_000017.11 61966475 g.61966475G>T snv heterozygous G T 0 0 0
5 │ NC_000019.10 54144058 g.54144058G>A snv heterozygous G A 0 0 0
6 │ NC_000023.11 77635947 g.77635947A>G snv hemizygous A G 0 0 0
7 │ NC_000005.10 1258495 g.1258495G>A snv heterozygous G A 0 0 0
8 │ NC_000012.12 2449086 g.2449086C>G snv heterozygous C G 0 0 0
****** KILL Après haplotypecaller
CLOSED: [2023-06-12 Mon 23:24]
******* KILL 20x
CLOSED: [2023-04-29 Sat 15:39]
Manque 183 sur 766
[[file:~/recherche/bisonex/simuscop/checkVCF.jl][checkVCF.jl]]
#+begin_src julia
@subset leftjoin(d2, dHaplo2, on=:genomic) ismissing.(:Column1)
#+end_src
Problème de profondeur ?
Ex: chr13 nombre de 101081606
NC_000011.10 16014966 g.16014966G>A
1 read sur 11 pour allèle alternative
Sur le patient de référence, 202 reads!
Celui-ci n'est pas le fichier de capture (ni dans le bam !)
ex: NC_000015.10 74343027 g.74343027C>T
Pour les autres, on devrait les retrouver...
Vérifier le nombre de reads sur 63003856
Vérifier la paramétrisation du modèle également
******* DONE [#B] 200x
CLOSED: [2023-05-18 Thu 11:04] SCHEDULED: <2023-04-30 Sun>
120 manquants (99 sans doublon)!
On vérifie dans IGV (vcf + bam après alignement) :
******** snv NC_000015.10 74343027
- rien d'appelé
- pas une région répétée
- base quality (voir [[*Phred score][Phred score]] ) à 37 donc ok
- variant retrouvé à 26/42
- Bam après aplybqsr: base qualità 35 donc ok
chr15 également à 89318565, variant retrouvé à 25/33 avec basequal de 37
Sans oublier de charger les instructions avx
#+begin_src sh
module load gcc@11.3.0/gcc-12.1.0
#+end_src
On coupe le .bam par chromosome pour débugger (sur le mesocentre)
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/simuscop-centogene-200x/cento/testing :results silent
ln -s ../preprocessing/applybqsr/cento.bam .
ln -s ../preprocessing/recalibrated/cento.bam.bai .
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
On doit lancer à la main (org-mode ne connait pas le chemin de samtools)
samtools view -b cento.bam NC_000015.10 > cento_chr15.bam
samtools index cento_chr15.bam
Puis on se restreint au chronmosome 15
samtools faidx genomeRef.fna NC_000015.10 > genomeRef_chr15.fa
samtools faidx genomeRef_chr15.fa
gatk CreateSequenceDictionary -R genomeRef_chr15.fa -O genomeRef_chr15.dict
On restreint au chromosome 15 avec l'option -L (dure = 1min)
gatk --java-options "-Xmx3072M" HaplotypeCaller --input cento_chr15.bam \
--output test.vcf.gz --reference genomeRef.fna --dbsnp dbSNP.gz --tmp-dir . --max-mnp-distance 2 -L NC_000015.10
******** DONE Tutorial haplotycaller
CLOSED: [2023-05-01 Mon 19:58]
Procédure : https://gatk.broadinstitute.org/hc/en-us/articles/360043491652-When-HaplotypeCaller-and-Mutect2-do-not-call-an-expected-variant
********* DONE Supprimer --max-mnp-distance = 2: idem
CLOSED: [2023-04-30 Sun 15:42]
********* DONE --debug &> run.log : Non appelé...
CLOSED: [2023-04-30 Sun 15:52]
********* DONE --linked-de-bruijn-graph: idem
CLOSED: [2023-04-30 Sun 15:55]
********* DONE --recover-all-dangling-branches
CLOSED: [2023-04-30 Sun 16:01]
********* DONE --min-pruning 0 : plus mais pas celui là
CLOSED: [2023-04-30 Sun 15:59]
********* DONE --bam-output
CLOSED: [2023-04-30 Sun 16:50]
********** DONE : rien !
CLOSED: [2023-04-30 Sun 16:08]
********** DONE + --recover-all-dangling-branches : rien !
CLOSED: [2023-04-30 Sun 16:08]
********* DONE Données filtrées ? apparement non
CLOSED: [2023-04-30 Sun 16:41]
183122 read(s) filtered by: MappingQualityReadFilter
3674 read(s) filtered by: NotDuplicateReadFilter
********** DONE --disable-read-filter MappingQualityReadFilter: idem
CLOSED: [2023-04-30 Sun 16:34]
On a bien - 0 read(s) filtered by: MappingQualityAvailableReadFilter
********** DONE --disable-read-filter NotDuplicateReadFilter: idem
CLOSED: [2023-04-30 Sun 16:40]
********* DONE Essayer freebayes : idem
CLOSED: [2023-04-30 Sun 16:22]
freebayes -f genomeRef.fna -r NC_000015.10 cento_chr15.bam > freebayes-test-chr15.vcf
********* DONE Avec toutes les options : idem
--linked-de-bruijn-graph --recover-all-dangling-branches --min-pruning 0 --bam-output debug.bam
CLOSED: [2023-04-30 Sun 16:50]
********* DONE Vérifier qu'on regarde le même bam : oui
CLOSED: [2023-04-30 Sun 16:50]
********* DONE Désactiver dbSNP : idem
CLOSED: [2023-04-30 Sun 16:52]
********* DONE Changer kmer size : idem
CLOSED: [2023-04-30 Sun 16:56]
par exemple[[https://gatk.broadinstitute.org/hc/en-us/community/posts/360075653152-REAL-Variant-not-called-by-HaplotypeCaller][forum gatk]] --kmer-size 18 --kmer-size 22
********* DONE --adaptive-pruning true
CLOSED: [2023-05-01 Mon 19:57]
******** DONE Mapping quality : est à 0 !!!!
CLOSED: [2023-05-01 Mon 19:58]
******* KILL Comparer VCF avec vcfeval :haplotypecaller:
CLOSED: [2023-06-12 Mon 23:24]
On prépare les données en julia
#+begin_src ~/recherche/bisonex/simuscop
julia --project=. toVCF.jl
#+end_src
Puis on export sur le mésocentre
#+begin_src
scp variants_for_vcfeval.tsv.gz* meso:centogene_variants/
#+end_src
#+begin_src
z bis
cd simuscop-200x
rtg vcfeval -b ~/centogene_variants/variants_for_vcfeval.tsv.gz -c cento/variantCalling/haplotypecaller/cento.vcf.gz -o compare-haplotypecaller -t /Work/Groups/bisonex/data/giab/GRCh38/genomeRef.sdf
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
82.000 540 540 60 45 0.9000 0.9231 0.9114
None 546 546 329 39 0.6240 0.9333 0.7479
******* KILL Comparer avec hap.py :haplotypecaller:
CLOSED: [2023-06-12 Mon 23:24]
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/checkInserted.nf -profile standard,helios --outdir=compare-simuscop-200x --query=out/simuscop-centogene-200x/cento/callVariant/haplotypecaller/cento.vcf.gz --truth=centogene_variants/variants_for_vcfeval.tsv.gz --id=simuscop-200x-check
******* DONE Méthode naïve 549/585
CLOSED: [2023-05-04 Thu 21:57]
Haplotypecaller: Nb reference SNV 692 vs found 585
Variant calling, filter technical: reference SNV 692 vs found 521
****** KILL Avant annotation
CLOSED: [2023-06-12 Mon 23:25] SCHEDULED: <2023-04-28 Fri>
#+begin_src
cd cento/variantCalling
bgzip filter-technical.vcf
tabix -p vcf filter-technical.vcf.gz -f
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
12.000 519 519 55 66 0.9042 0.8872 0.8956
None 519 519 55 66 0.9042 0.8872 0.8956
******* DONE Méthode naïve 521/585
CLOSED: [2023-05-04 Thu 21:57]
Haplotypecaller: Nb reference SNV 692 vs found 585
Variant calling, filter technical: reference SNV 692 vs found 521
******* KILL Comparer avec hap.py
CLOSED: [2023-06-12 Mon 23:24]
****** KILL Après filtre annotation
CLOSED: [2023-06-12 Mon 23:25]
******* DONE Méthode naïve : 493/585
CLOSED: [2023-05-04 Thu 22:09]
******* KILL Comparer avec hap.py
CLOSED: [2023-06-12 Mon 23:25]
******* KILL VCf eval
CLOSED: [2023-06-12 Mon 23:25]
cd cento/annotation/
bgzip postvep-filter.vcf
tabix postvep-filter.vcf.gz
cd ../..
rtg vcfeval -b ~/centogene_variants/variants_for_vcfeval.tsv.gz -c cento/annotation/postvep-filter.vcf.gz -o compare-vepfilter -t /Work/Groups/bisonex/data/giab/GRCh38/genomeRef.sdf
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
12.000 491 491 50 94 0.9076 0.8393 0.8721
None 491 491 50 94 0.9076 0.8393 0.8721
***** TODO Vérifier tous les variants sont retrouvés en 200x: hg38 :T2T:
****** DONE Après alignement
CLOSED: [2023-04-29 Sat 18:27] SCHEDULED: <2023-04-28 Fri>
******* DONE SNV: avec doublons
CLOSED: [2023-04-28 Fri 18:12]
On utilise [[file:~/recherche/bisonex/simuscop/checkBam.jl][checkBam.jl]]
#+begin_src julia
d = prepareVariant("../parsevariants/variant_genomic.csv")
root = "/home/alex/code/bisonex/simuscop-centogene/cento"
bam = root * "/preprocessing/applybqsr/cento.bam"
bai = root * "/preprocessing/recalibrated/cento.bam.bai"
snv = getSNV(d, bam, bai)
#+end_src
Nombreux faux homozygouteS
Vérification avec checkFalseHemizygous(snv) : nombreux doublons dans le fichier pour simuscop...
******* DONE SNV sans doublons
CLOSED: [2023-04-29 Sat 18:27]
******** DONE 18 faux homozygote mais avec peu de reads
CLOSED: [2023-04-29 Sat 18:27]
julia> @subset snv :refCount .== 0 :altCount .> 0 :zygosity .== "heterozygous"
18×10 DataFrame
Row │ chrom pos variant variantType zygosity ref alt refCount altCount readsCount
│ SubStrin…? Int64 SubStrin…? String? String15 SubStrin… SubStrin… Int64 Int64 Int64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ NC_000022.11 42213078 g.42213078T>G snv heterozygous T G 0 1 1
2 │ NC_000012.12 101680427 g.101680427C>A snv heterozygous C A 0 3 3
3 │ NC_000014.9 105385684 g.105385684G>C snv heterozygous G C 0 4 4
4 │ NC_000011.10 125978299 g.125978299C>T snv heterozygous C T 0 3 3
5 │ NC_000023.11 77998618 g.77998618C>T snv heterozygous C T 0 2 2
6 │ NC_000015.10 66703292 g.66703292C>T snv heterozygous C T 0 3 3
7 │ NC_000010.11 87961118 g.87961118G>A snv heterozygous G A 0 3 3
8 │ NC_000012.12 112477719 g.112477719A>G snv heterozygous A G 0 2 2
9 │ NC_000020.11 6778406 g.6778406C>T snv heterozygous C T 0 3 3
10 │ NC_000023.11 68192943 g.68192943G>A snv heterozygous G A 0 2 2
11 │ NC_000004.12 987858 g.987858C>T snv heterozygous C T 0 3 4
12 │ NC_000015.10 66435145 g.66435145G>A snv heterozygous G A 0 1 2
13 │ NC_000002.12 47809595 g.47809595C>T snv heterozygous C T 0 2 2
14 │ NC_000003.12 136477305 g.136477305C>G snv heterozygous C G 0 4 4
15 │ NC_000005.10 157285458 g.157285458C>T snv heterozygous C T 0 3 3
16 │ NC_000012.12 23604413 g.23604413T>G snv heterozygous T G 0 5 5
17 │ NC_000019.10 52219703 g.52219703C>T snv heterozygous C T 0 1 1
18 │ NC_000016.10 88856757 g.88856757C>T snv heterozygous C T 0 8 8
******** DONE 8 non retrouvé => probablement hors de la zjone de capture
CLOSED: [2023-04-28 Fri 19:49]
julia> @subset snv :refCount .== 0 :altCount .== 0
8×10 DataFrame
Row │ chrom pos variant variantType zygosity ref alt refCount altCount readsCount
│ SubStrin…? Int64 SubStrin…? String? String15 SubStrin… SubStrin… Int64 Int64 Int64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ NC_000015.10 74343027 g.74343027C>T snv heterozygous C T 0 0 0
2 │ NC_000011.10 20638345 g.20638345A>G snv heterozygous A G 0 0 0
3 │ NC_000004.12 139370252 g.139370252C>T snv heterozygous C T 0 0 2
4 │ NC_000017.11 61966475 g.61966475G>T snv heterozygous G T 0 0 0
5 │ NC_000019.10 54144058 g.54144058G>A snv heterozygous G A 0 0 0
6 │ NC_000023.11 77635947 g.77635947A>G snv hemizygous A G 0 0 0
7 │ NC_000005.10 1258495 g.1258495G>A snv heterozygous G A 0 0 0
8 │ NC_000012.12 2449086 g.2449086C>G snv heterozygous C G 0 0 0
****** KILL Après haplotypecaller
CLOSED: [2023-06-12 Mon 23:24]
******* KILL 20x
CLOSED: [2023-04-29 Sat 15:39]
Manque 183 sur 766
[[file:~/recherche/bisonex/simuscop/checkVCF.jl][checkVCF.jl]]
#+begin_src julia
@subset leftjoin(d2, dHaplo2, on=:genomic) ismissing.(:Column1)
#+end_src
Problème de profondeur ?
Ex: chr13 nombre de 101081606
NC_000011.10 16014966 g.16014966G>A
1 read sur 11 pour allèle alternative
Sur le patient de référence, 202 reads!
Celui-ci n'est pas le fichier de capture (ni dans le bam !)
ex: NC_000015.10 74343027 g.74343027C>T
Pour les autres, on devrait les retrouver...
Vérifier le nombre de reads sur 63003856
Vérifier la paramétrisation du modèle également
******* DONE [#B] 200x
CLOSED: [2023-05-18 Thu 11:04] SCHEDULED: <2023-04-30 Sun>
120 manquants (99 sans doublon)!
On vérifie dans IGV (vcf + bam après alignement) :
******** snv NC_000015.10 74343027
- rien d'appelé
- pas une région répétée
- base quality (voir [[*Phred score][Phred score]] ) à 37 donc ok
- variant retrouvé à 26/42
- Bam après aplybqsr: base qualità 35 donc ok
chr15 également à 89318565, variant retrouvé à 25/33 avec basequal de 37
Sans oublier de charger les instructions avx
#+begin_src sh
module load gcc@11.3.0/gcc-12.1.0
#+end_src
On coupe le .bam par chromosome pour débugger (sur le mesocentre)
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/simuscop-centogene-200x/cento/testing :results silent
ln -s ../preprocessing/applybqsr/cento.bam .
ln -s ../preprocessing/recalibrated/cento.bam.bai .
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
On doit lancer à la main (org-mode ne connait pas le chemin de samtools)
samtools view -b cento.bam NC_000015.10 > cento_chr15.bam
samtools index cento_chr15.bam
Puis on se restreint au chronmosome 15
samtools faidx genomeRef.fna NC_000015.10 > genomeRef_chr15.fa
samtools faidx genomeRef_chr15.fa
gatk CreateSequenceDictionary -R genomeRef_chr15.fa -O genomeRef_chr15.dict
On restreint au chromosome 15 avec l'option -L (dure = 1min)
gatk --java-options "-Xmx3072M" HaplotypeCaller --input cento_chr15.bam \
--output test.vcf.gz --reference genomeRef.fna --dbsnp dbSNP.gz --tmp-dir . --max-mnp-distance 2 -L NC_000015.10
******** DONE Tutorial haplotycaller
CLOSED: [2023-05-01 Mon 19:58]
Procédure : https://gatk.broadinstitute.org/hc/en-us/articles/360043491652-When-HaplotypeCaller-and-Mutect2-do-not-call-an-expected-variant
********* DONE Supprimer --max-mnp-distance = 2: idem
CLOSED: [2023-04-30 Sun 15:42]
********* DONE --debug &> run.log : Non appelé...
CLOSED: [2023-04-30 Sun 15:52]
********* DONE --linked-de-bruijn-graph: idem
CLOSED: [2023-04-30 Sun 15:55]
********* DONE --recover-all-dangling-branches
CLOSED: [2023-04-30 Sun 16:01]
********* DONE --min-pruning 0 : plus mais pas celui là
CLOSED: [2023-04-30 Sun 15:59]
********* DONE --bam-output
CLOSED: [2023-04-30 Sun 16:50]
********** DONE : rien !
CLOSED: [2023-04-30 Sun 16:08]
********** DONE + --recover-all-dangling-branches : rien !
CLOSED: [2023-04-30 Sun 16:08]
********* DONE Données filtrées ? apparement non
CLOSED: [2023-04-30 Sun 16:41]
183122 read(s) filtered by: MappingQualityReadFilter
3674 read(s) filtered by: NotDuplicateReadFilter
********** DONE --disable-read-filter MappingQualityReadFilter: idem
CLOSED: [2023-04-30 Sun 16:34]
On a bien - 0 read(s) filtered by: MappingQualityAvailableReadFilter
********** DONE --disable-read-filter NotDuplicateReadFilter: idem
CLOSED: [2023-04-30 Sun 16:40]
********* DONE Essayer freebayes : idem
CLOSED: [2023-04-30 Sun 16:22]
freebayes -f genomeRef.fna -r NC_000015.10 cento_chr15.bam > freebayes-test-chr15.vcf
********* DONE Avec toutes les options : idem
--linked-de-bruijn-graph --recover-all-dangling-branches --min-pruning 0 --bam-output debug.bam
CLOSED: [2023-04-30 Sun 16:50]
********* DONE Vérifier qu'on regarde le même bam : oui
CLOSED: [2023-04-30 Sun 16:50]
********* DONE Désactiver dbSNP : idem
CLOSED: [2023-04-30 Sun 16:52]
********* DONE Changer kmer size : idem
CLOSED: [2023-04-30 Sun 16:56]
par exemple[[https://gatk.broadinstitute.org/hc/en-us/community/posts/360075653152-REAL-Variant-not-called-by-HaplotypeCaller][forum gatk]] --kmer-size 18 --kmer-size 22
********* DONE --adaptive-pruning true
CLOSED: [2023-05-01 Mon 19:57]
******** DONE Mapping quality : est à 0 !!!!
CLOSED: [2023-05-01 Mon 19:58]
******* KILL Comparer VCF avec vcfeval :haplotypecaller:
CLOSED: [2023-06-12 Mon 23:24]
On prépare les données en julia
#+begin_src ~/recherche/bisonex/simuscop
julia --project=. toVCF.jl
#+end_src
Puis on export sur le mésocentre
#+begin_src
scp variants_for_vcfeval.tsv.gz* meso:centogene_variants/
#+end_src
#+begin_src
z bis
cd simuscop-200x
rtg vcfeval -b ~/centogene_variants/variants_for_vcfeval.tsv.gz -c cento/variantCalling/haplotypecaller/cento.vcf.gz -o compare-haplotypecaller -t /Work/Groups/bisonex/data/giab/GRCh38/genomeRef.sdf
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
82.000 540 540 60 45 0.9000 0.9231 0.9114
None 546 546 329 39 0.6240 0.9333 0.7479
******* KILL Comparer avec hap.py :haplotypecaller:
CLOSED: [2023-06-12 Mon 23:24]
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/checkInserted.nf -profile standard,helios --outdir=compare-simuscop-200x --query=out/simuscop-centogene-200x/cento/callVariant/haplotypecaller/cento.vcf.gz --truth=centogene_variants/variants_for_vcfeval.tsv.gz --id=simuscop-200x-check
******* DONE Méthode naïve 549/585
CLOSED: [2023-05-04 Thu 21:57]
Haplotypecaller: Nb reference SNV 692 vs found 585
Variant calling, filter technical: reference SNV 692 vs found 521
****** KILL Avant annotation
CLOSED: [2023-06-12 Mon 23:25] SCHEDULED: <2023-04-28 Fri>
#+begin_src
cd cento/variantCalling
bgzip filter-technical.vcf
tabix -p vcf filter-technical.vcf.gz -f
#+end_src
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
12.000 519 519 55 66 0.9042 0.8872 0.8956
None 519 519 55 66 0.9042 0.8872 0.8956
******* DONE Méthode naïve 521/585
CLOSED: [2023-05-04 Thu 21:57]
Haplotypecaller: Nb reference SNV 692 vs found 585
Variant calling, filter technical: reference SNV 692 vs found 521
******* KILL Comparer avec hap.py
CLOSED: [2023-06-12 Mon 23:24]
****** KILL Après filtre annotation
CLOSED: [2023-06-12 Mon 23:25]
******* DONE Méthode naïve : 493/585
CLOSED: [2023-05-04 Thu 22:09]
******* KILL Comparer avec hap.py
CLOSED: [2023-06-12 Mon 23:25]
******* KILL VCf eval
CLOSED: [2023-06-12 Mon 23:25]
cd cento/annotation/
bgzip postvep-filter.vcf
tabix postvep-filter.vcf.gz
cd ../..
rtg vcfeval -b ~/centogene_variants/variants_for_vcfeval.tsv.gz -c cento/annotation/postvep-filter.vcf.gz -o compare-vepfilter -t /Work/Groups/bisonex/data/giab/GRCh38/genomeRef.sdf
Threshold True-pos-baseline True-pos-call False-pos False-neg Precision Sensitivity F-measure
----------------------------------------------------------------------------------------------------
12.000 491 491 50 94 0.9076 0.8393 0.8721
None 491 491 50 94 0.9076 0.8393 0.8721
**** KILL NEAT : trop lent :neat:
CLOSED: [2023-04-29 Sat 22:06]
***** KILL Génération fastq sur exno 5 GATAD2B
CLOSED: [2023-04-29 Sat 22:06]
Trop lent : pour 1 exon : 1500 secondes !
#+begin_src sh
samtools faidx genomeRef.fna NC_000001.11 | save -f genomeRef_chr1.fna
python gen_reads.py -r ../test-simuscop/genomeRef_chr1.fna -o lol -tr ../test-simuscop/gatad2b-exon6.bed -R 147 --pe 150 10
#+end_src
**** KILL ReSeq : exome avec exons comme fasta mais ne gère pas des exons trop petits :reseq:
CLOSED: [2023-04-30 Sun 19:44] SCHEDULED: <2023-04-29 Sat>
#+begin_quote
Can I simulate exome sequencing? Yes. You need to use a reference that only contains the exons as individual scaffolds. Using --refBiasFile you can specify the coverage of individual exons. To simulate intron contamination you can add the whole reference to the reference containing the exons and strongly reduce the coverage for these scaffolds using --refBiasFile.
#+end_quote
Par contre, rapide
***** DONE Fasta pour exons seuls
CLOSED: [2023-04-30 Sun 19:25]
Depuis le GFF
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_genomic.gff.gz
#+end_src
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
gunzip -c GCF_000001405.39_GRCh38.p13_genomic.gff.gz | grep -w "exon" > exons.gff
#+end_src
On génère les exons
#+begin_src sh :dir ~/code/bisonex/test-reseq
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed exons.gff -fo exons.fna
#+end_src
A tester avec un profile déjà fait :
https://github.com/schmeing/ReSeq-profiles/tree/master/profiles
On cherche l'exons qui nous intéresse
NC_000001.11 g.153817496 A>T
N'y est pas ??
****** DONE On test sur les 2 premiers : exec
CLOSED: [2023-04-30 Sun 18:39]
#+begin_src
head exons.fa -n 2 > 2exons.fna
#+end_src
#+begin_src sh
../ReSeq/bin/reseq illuminaPE -j 32 -R exons.fa -s Ec-Hi2000-TruSeq.reseq --ipfIterations 0 -1 reseq-sim_1.fq reseq_sim_2.fq
#+end_src
#+begin_quote
error: All reference sequences are too short for simulating. They should have at least 1991 bases
#+end_quote
#+begin_src sh
grep '^>NC_000001.10' exons.fa | sed 's/:/,/;s/-/,/;s/^>//' > exons.csv
#+end_src
****** DONE Sur 200 premiers exons du chr1
CLOSED: [2023-04-30 Sun 19:17]
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
head -n200 exons.fna > exons-200.fna
bwa index exons-200.fna
#+end_src
Simulation avec 30x
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
../ReSeq/bin/reseq illuminaPE -R exons-200.fna -s Ec-Hi2000-TruSeq.reseq --ipfIterations 0 -1 reseq1.fq -2 reseq2.fq -c 30
#+end_src
Attention, pour l'alignement, il faut le nfa complet ! Sinon erreur du type
Erreurs:::sam_hdr_create] Duplicated sequence "NC_000001.10:762970-763155" in file "-"
Et pas de bam avec
samtools sort: failed to change sort order header to 'coordinate'
#+begin_src
bwa mem ../test-simuscop/bwa/genomeRef.fna reseq1.fq reseq2.fq | samtools sort -o reseq.bam
#+end_src
Manque des exons et l'allure ne correspond pas...
****** DONE Utiliser le fichier de capture : exons trop petits
CLOSED: [2023-04-30 Sun 19:25]
Comme pour ART
Trop court avec
echo -e "NC_000001.11\t153817371\t153817542" > gatad2b-exon6.bed
Donc on ajoute 1000 de chaque côté
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
echo -e "NC_000001.11\t153816371\t153818542" > gatad2b-exon6.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6.bed -fo gatad2b-exon6.fna
bwa index gatad2b-exon6.bed
../ReSeq/bin/reseq illuminaPE -R gatad2b-exon6.fna -s Ec-Hi2000-TruSeq.reseq --ipfIterations 0 -1 reseq1.fq -2 reseq2.fq -c 30
bwa mem ../test-simuscop/bwa/genomeRef.fna reseq1.fq reseq2.fq | samtools sort -o reseq.bam
samtools index reseq.bam
#+end_src
***** KILL Sur le chromosome 15 puis trier à la main sur les zones de capture ?
CLOSED: [2023-04-30 Sun 19:44]
#+begin_src sh :dir ~/code/bisonex/test-reseq :results silent
samtools faidx ../test-simuscop/genomeRef.fna NC_000015.10 > chr15.fna
../ReSeq/bin/reseq illuminaPE -R chr15.fna -s Ec-Hi2000-TruSeq.reseq --ipfIterations 0 -1 reseq1.fq -2 reseq2.fq -c 30
#+end_src
**** DONE ART : fonctionne très mal en targeted
CLOSED: [2023-04-30 Sun 11:49]
***** DONE Génération de reads
CLOSED: [2023-04-30 Sun 11:49]
****** DONE Avec seulement les exons en séquence
CLOSED: [2023-04-30 Sun 10:24]
head -n6 exons.fa | save three-exons.fna
../art_bin_MountRainier/art_illumina -ss HS25 -i three-exons.fna -o ./paired_end_com -l 150 -f 10 -p -m 500 -s 10 -sam
Le sam n'est pas visible sur igv mais si on aligne avec bwa mem, on a quelques reads
****** DONE Extraire une zone de capture dans le fasta
CLOSED: [2023-04-30 Sun 11:49]
NC_000001.11 g.153817496 A>T
******* DONE Essai 1: ne dépasse pas la zone
CLOSED: [2023-04-30 Sun 10:49]
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153817371\t153817542" > gatad2b-exon6.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6.bed -fo gatad2b-exon6.fa
#+end_src
-ss HS25 : nom du profile illumina
-l 150 : reads de 150
-f 10 : coverage de 10
-p : paired end
-m 500 : longueur moyenne des fragment d'ADN
-s 10 : déviation standard
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6.fa -o ./gatad2b -l 150 -f 100 -p -m 500 -s 10
#+end_src
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.b
am
samtools index gatad2b.bam
#+end_src
#+attr_html: :width 800px
[[./art-capture-1.png]]
******* Avec offset
50bp idem
NC_000001.11 153817371 153817542 -> NC_000001.11 153817321 153817592
On essaie 1000
NC_000001.11 153817371 153817542 -> NC_000001.11 153816371 153818542 ->
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153816371\t153818542" > gatad2b-exon6-offset.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6-offset.bed -fo gatad2b-exon6-offset.fa
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6-offset.fa -o ./gatad2b -l 150 -f 100 -p -m 500 -s 10
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.bam
samtools index gatad2b.bam
#+end_src
mieux mais trop large
#+attr_html: :width 800px
[[./art-exon6-offset1000.png]]
Sur les vraies données, on a une large de 500bp environ
#+attr_html: :width 800px
[[./illumina-ref-exon6.png]]
Ici l'exon fait 250bp donc on rajouter 125bp de chaque côté
NC_000001.11 153817371 153817542
NC_000001.11 153817246 153817667
Résulats incohérents :on a 2 colonnes séparées !
Essai avec +200bp de chaque côt
NC_000001.11 153817371 153817542
NC_000001.11 153817171 153817742
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153817171\t153817742"> gatad2b-exon6-offset.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6-offset.bed -fo gatad2b-exon6-offset.fa
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6-offset.fa -o ./gatad2b -l 150 -f 100 -p -m 500 -s 50
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.bam
samtools index gatad2b.bam
#+end_src
En changeant la longueur des reads
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153817171\t153817742"> gatad2b-exon6-offset.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6-offset.bed -fo gatad2b-exon6-offset.fa
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6-offset.fa -o ./gatad2b -l 126 -f 100 -p -m 500 -s 50
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.bam
samtools index gatad2b.bam
#+end_src
Idem, fait 2 "tas" séparés de chaque côté de l'exon
+/- 100
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
echo -e "NC_000001.11\t153817271\t153817642" > gatad2b-exon6-offset.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6-offset.bed -fo gatad2b-exon6-offset.fa
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6-offset.fa -o ./gatad2b -l 125 -f 100 -p -m 500 -s 30
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.bam
samtools index gatad2b.bam
#+end_src
**** KILL NGSNG : plante
CLOSED: [2023-05-01 Mon 09:27]
https://github.com/RAHenriksen/NGSNGS
***** KILL Essai avec fasta
CLOSED: [2023-05-01 Mon 09:27]
#+begin_src sh :dir ~/code/bisonex/test-ngsngs :results silent
echo -e "NC_000001.11\t153817371\t153817542" > gatad2b-exon6.bed
bedtools getfasta -fi ../test-simuscop/genomeRef.fna -bed gatad2b-exon6.bed -fo gatad2b-exon6.fa
#+end_src
-ss HS25 : nom du profile illumina
-l 150 : reads de 150
-f 10 : coverage de 10
-p : paired end
-m 500 : longueur moyenne des fragment d'ADN
-s 10 : déviation standard
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
../art_bin_MountRainier/art_illumina -ss HS25 -i gatad2b-exon6.fa -o ./gatad2b -l 150 -f 100 -p -m 500 -s 10
#+end_src
#+begin_src sh :dir ~/code/bisonex/test-art :results silent
bwa mem ../test-simuscop/bwa/genomeRef gatad2b1.fq gatad2b2.fq | samtools sort -o gatad2b.b
am
samtools index gatad2b.bam
#+end_src
**** KILL Script maison :generate:
CLOSED: [2023-05-13 Sat 18:29] SCHEDULED: <2023-05-01 Mon>
***** KILL SNV
CLOSED: [2023-05-13 Sat 18:29] SCHEDULED: <2023-05-01 Mon>
****** DONE Script python: ok seulement pour reads corrigé, trop long sinon
CLOSED: [2023-05-01 Mon 13:32]
Même sur un seul chromosome (15)...
****** DONE Couper les données avec bedtools intersect ?
CLOSED: [2023-05-01 Mon 20:33]
-wa pour avoir les reads qui corresponds
-v pour ceux qui n'intersecte pass
ok sur un exemple
******* DONE Test simple : ok
CLOSED: [2023-05-01 Mon 13:43]
#+attr_html: :width 500px
[[./test-intersect-chr15.png]]
#+attr_html: :width 500px
[[./test-nointersect-chr15.png]]
******* DONE Chromosome 15 : vérifier BAM
CLOSED: [2023-05-01 Mon 20:33] SCHEDULED: <2023-05-01 Mon>
#+begin_src
samtools view -b `63003856_S135.bam` NC_000015.10 > `63003856_S135_chr15.bam`
#+end_src
On génère un python avec les dépendances
#+begin_src
nix-build
#+end_src
Puis on lance un script julia qui va couper le bam en 2, lancer le script python sur l'intersection et fusionner le résultat
#+begin_src
julia -Jbisonex.so --project=. insertVariants.jl `63003856_S135_chr15.bam` test_new.bam
#+end_src
NB: pour accélerer l'exécution, générer une sysimage :
#+begin_src julia
(1.7) pkg> activate .
julia> create_sysimage(; sysimage_path="bisonex.so")
#+end_src
3 variants: Ok sur le nombre de reads et les variants
Ok pour homozygote
#+attr_html: :width 500px
[[./check-snv-chr15.png]]
******* DONE Chromosome 15 Vérifier VAF avec checkBam.jl: ok
julia> @subset snv :chrom .== "NC_000015.10"
26×10 DataFrame
Row │ chrom pos variant variantType zygosity ref alt refCount altCount readsCount
│ SubStrin…? Int64 SubStrin…? String? String15 SubStrin… SubStrin… Int64 Int64 Int64
─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ NC_000015.10 74343027 g.74343027C>T snv heterozygous C T 61 58 120
2 │ NC_000015.10 75400778 g.75400778C>G snv heterozygous C G 108 79 187
3 │ NC_000015.10 89327201 g.89327201C>T snv heterozygous C T 243 241 486
4 │ NC_000015.10 48767448 g.48767448A>C snv heterozygous A C 72 70 142
5 │ NC_000015.10 75411685 g.75411685T>C snv heterozygous T C 79 81 160
6 │ NC_000015.10 66703292 g.66703292C>T snv heterozygous C T 70 60 130
7 │ NC_000015.10 89325639 g.89325639G>A snv heterozygous G A 257 267 524
8 │ NC_000015.10 89330184 g.89330184G>A snv heterozygous G A 258 287 548
9 │ NC_000015.10 89330184 g.89330184G>A snv heterozygous G A 258 287 548
10 │ NC_000015.10 89325639 g.89325639G>A snv heterozygous G A 257 267 524
11 │ NC_000015.10 42401752 g.42401752G>A snv homozygous G A 61 212 273
12 │ NC_000015.10 89327201 g.89327201C>T snv heterozygous C T 243 241 486
13 │ NC_000015.10 38339896 g.38339896G>A snv heterozygous G A 56 86 144
14 │ NC_000015.10 26869324 g.26869324A>T snv heterozygous A T 62 49 113
15 │ NC_000015.10 66435145 g.66435145G>A snv heterozygous G A 98 95 193
16 │ NC_000015.10 60514655 g.60514655G>A snv heterozygous G A 94 99 194
17 │ NC_000015.10 42410947 g.42410947A>G snv heterozygous A G 153 123 276
18 │ NC_000015.10 75430368 g.75430368C>T snv heterozygous C T 80 62 142
19 │ NC_000015.10 25375494 g.25375494T>C snv heterozygous T C 103 104 207
20 │ NC_000015.10 60497497 g.60497497C>A snv heterozygous C A 61 65 126
21 │ NC_000015.10 74891539 g.74891539C>T snv heterozygous C T 118 124 242
22 │ NC_000015.10 48488433 g.48488433A>G snv heterozygous A G 367 122 492
23 │ NC_000015.10 89318565 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
**** HOLD Bamsurgeon :bamsurgeon:
***** HOLD 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
****** HOLD Corriger run nextflow pour éviter les erreurs
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]
***** HOLD Variants cento
****** STRT SNV
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>
******* HOLD Corrigier position pour avoir une bonne VAF
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 [[id:966a298c-948a-4694-a6f5-c326b1046a05][XAMscissors.jl]] :xamscissors:
***** TODO Test SNV
****** DONE Phase 1 : chr22, VAF=1
CLOSED: [2023-05-29 Mon 15:36]
******* 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
******* KILL Tester SNV chromosome 22
CLOSED: [2023-05-29 Mon 15:36] SCHEDULED: <2023-05-20 Sat>
****** KILL PHase 2 : chr22, VAF variable
CLOSED: [2023-06-12 Mon 23:27] SCHEDULED: <2023-06-03 Sat>
******* DONE de nombreux reads sont perdus -> ok sur un SNV en alignant sur chromosome 22
CLOSED: [2023-05-29 Mon 15:38]
Problème dansr le BAM car sans les insertion
******** DONE Filtrer les reads sans pair : idem
CLOSED: [2023-05-23 Tue 00:01]
FOUND:Found 16662 unpaired mates
at htsjdk.samtools.SAMUtils.processValidationError(SAMUtils.java:470)
at picard.sam.SamToFastq.doWork(SamToFastq.java:224)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:308)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)
ex: chr22:42213078
On essaie
samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -hf 0x2 -o 63003856_chr22.bam
Ne rale pas...
SCHEDULED: <2023-05-22 Mon>
******** DONE Problème de la conversion BAM -> fastq ?
CLOSED: [2023-05-24 Wed 23:19]
********* DONE Test bamtofastq sur le BAM originel: idem
CLOSED: [2023-05-23 Tue 01:16]
Dans ~/recherche/code/bisonex/Bamscissors.jl
samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -hf 0x2 -o 63003856_chr22.bam
On trie bien par nom de read
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
Conversion en fastq. Attention à bien prendre le *fichier trié* !
samtools bam2fq -1 63003856_chr22_init1.fq.gz -2 63003856_chr22_init2.fq.gz -n 63003856_chr22_sorted.bam
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 2592110 reads
On envoie sur le mésocentre
scp 63003856_chr22_init{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors
On démarre un job avec 24 coeurs pour aller vite
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_init1.fq.gz 63003856_chr22_init2.fq.gz | samtools sort -@24 -o output.bam -
Dans /home/alex/recherche/bisonex/code/BamScissors.jl/aligned
❯ scp meso:/Work/Users/apraga/bisonex/tests/bamscissors/output.bam*
********* DONE Avec samtools fastq
CLOSED: [2023-05-23 Tue 01:16]
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -o 63003856_chr22_sorted.bam
scp 63003856_chr22_{1,2}.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
mesocentre
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz | samtools sort -@24 -o output.bam
scp meso:/Work/Users/apraga/bisonex/tests/bamscissors/output.bam\* aligned/
********* DONE En rajoutant .fna dans le doserrie (+samtools fastq): idem
CLOSED: [2023-05-23 Tue 01:16]
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef* .
ln -s /Work/Projects/bisonex/data/genome/GRCh38.p13/genomeRef* .
bwa mem -t 24 genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz | samtools sort -@24 -o output.bam
********* DONE Test picard : idem
CLOSED: [2023-05-23 Tue 01:16]
Dans bisonex/code/BamScissors.jl
❯ picard SamToFastq -I 63003856_chr22_sorted.bam -F testpicard1.fastq -F2 testpicard2.fastq
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef testpicard1.fastq.gz testpicard2.fastq.gz | samtools sort -@24 -o testpicard.bam -
********* DONE Il ne manque pas des reads dans les fastq :
CLOSED: [2023-05-23 Tue 01:26]
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools flagstat aligned/output.bam
1306273 + 0 in total (QC-passed reads + QC-failed reads)
1296055 + 0 primary
0 + 0 secondary
10218 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
1304871 + 0 mapped (99.89% : N/A)
1294653 + 0 primary mapped (99.89% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools flagstat 63003856_chr22.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
2592110 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2592110 + 0 properly paired (100.00% : N/A)
2592110 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
❯ echo $(zcat 63003856_chr22_init2.fq.gz | wc -l)/4 | bc
1296055
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0 took 2s
❯ echo $(zcat 63003856_chr22_init1.fq.gz | wc -l)/4 | bc
1296055
❯ samtools flagstat 63003856_chr22_sorted.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
2592110 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2592110 + 0 properly paired (100.00% : N/A)
2592110 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
********* DONE Avec sort -O BAM idem
CLOSED: [2023-05-23 Tue 01:21]
********* DONE Singleton ou non mapped ? nonVirtPosition
CLOSED: [2023-05-23 Tue 01:22]
❯ samtools bam2fq -1 63003856_chr22_init1.fq.gz -2 63003856_chr22_init2.fq.gz -0 both -s single.fq.gz -n 63003856_chr22_sorted.bam
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ wc -l 63003856_chr22_init* both single.fq.gz
403540 63003856_chr22_init1.fq.gz
404788 63003856_chr22_init2.fq.gz
0 both
0 single.fq.gz
808328 total
********* Problème d'aligner car les reads sont bien dans le .fastq ?
Ex:
zgrep "A00853:477:HMLWYDSX3:2:2624:2826:18630" 63003856_chr22_{1,2}.fq.gz
63003856_chr22_1.fq.gz:@A00853:477:HMLWYDSX3:2:2624:2826:18630
63003856_chr22_2.fq.gz:@A00853:477:HMLWYDSX3:2:2624:2826:18630
******** DONE Refaire la manip sur bam chr22 non modifié + mail Alexis
CLOSED: [2023-05-23 Tue 23:27]
[1]_samtools view ~/code/bisonex/out/63003856/preprocessing/mapped/63003856_S135.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
Les reads sont bien tous mappé
samtools flagstat 63003856_chr22.bam
2609214 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
17104 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2609214 + 0 mapped (100.00% : N/A)
[2] samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
flagstat idem
[3] samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 2592110 reads
Nombre de reads ok (= la moitié) et les séquences ne sont pas identiques pour le premier read (= on n'a pas 2x la même chose)
echo $(zcat 63003856_chr22_1.fq.gz | wc -l)/4 | bc
1296055
echo $(zcat 63003856_chr22_2.fq.gz | wc -l)/4 | bc
1296055
[4]
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz -o wtf.bam
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1752492 sequences (240000014 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (12, 702821, 0, 18)
[M::mem_pestat] analyzing insert size distribution for orientation FF...
[M::mem_pestat] (25, 50, 75) percentile: (54, 197, 268)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 696)
[M::mem_pestat] mean and std.dev: (163.92, 129.71)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 910)
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (129, 175, 233)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 441)
[M::mem_pestat] mean and std.dev: (185.56, 75.37)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 545)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation RR...
[M::mem_pestat] (25, 50, 75) percentile: (56, 192, 239)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 605)
[M::mem_pestat] mean and std.dev: (158.00, 110.30)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 788)
[M::mem_pestat] skip orientation FF
[M::mem_pestat] skip orientation RR
[M::process] read 839618 sequences (114952336 bp)...
[M::mem_process_seqs] Processed 1752492 reads in 375.714 CPU sec, 17.645 real sec
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 336379, 0, 5)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (128, 174, 232)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 440)
[M::mem_pestat] mean and std.dev: (184.73, 74.63)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 544)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 839618 reads in 183.039 CPU sec, 7.961 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t 24 -o wtf.bam /Work/Projects/bisonex/data/genome/GRCh38.p13/bwa/genomeRef 63003856_chr22_1.fq.gz 63003856_chr22_2.fq.gz
[main] Real time: 38.278 sec; CPU: 565.821 sec
Bon nombre de reads pourtant
samtools flagstat wtf.bam
2611059 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
18949 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2611058 + 0 mapped (100.00% : N/A)
2592109 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2590970 + 0 properly paired (99.96% : N/A)
2592108 + 0 with itself and mate mapped
1 + 0 singletons (0.00% : N/A)
458 + 0 with mate mapped to a different chr
63 + 0 with mate mapped to a different chr (mapQ>=5)
$ samtools sort -@24 -o wtf_sorted.bam wtf.sam
[bam_sort_core] merging from 0 files and 24 in-memory blocks...
samtools flagstat wtf.bam
2611059 + 0 in total (QC-passed reads + QC-failed reads)
2592110 + 0 primary
0 + 0 secondary
18949 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
2611058 + 0 mapped (100.00% : N/A)
2592109 + 0 primary mapped (100.00% : N/A)
2592110 + 0 paired in sequencing
1296055 + 0 read1
1296055 + 0 read2
2590970 + 0 properly paired (99.96% : N/A)
2592108 + 0 with itself and mate mapped
1 + 0 singletons (0.00% : N/A)
458 + 0 with mate mapped to a different chr
63 + 0 with mate mapped to a different chr (mapQ>=5)
Effectivement, ce n'est pas un problème d'IGV
❯ samtools mpileup 63003856_chr22.bam -r NC_000022.11:42213078-42213078
[mpileup] 1 samples in 1 input files
NC_000022.11 42213078 N 19 TTtTTtTTTtTTtTtTttt kkk_kkkkFkFF_FkFkQk
bisonex/code/BamScissors.jl on bamscissors [!?] via ஃ v1.9.0
❯ samtools mpileup aligned/wtf_sorted.bam -r NC_000022.11:42213078-42213078
[mpileup] 1 samples in 1 input files
NC_000022.11 42213078 N 5 TTtTT _FkFF
******** DONE Regarder où ont été aligné les reads (nouveau run)
CLOSED: [2023-05-24 Wed 23:18]
********* DONE Préparation
CLOSED: [2023-05-24 Wed 21:59]
On relance le pipeline pour avoir un BAM propre
On garde les reads non mappé à partsir de la sortie d'applybqsr
#+begin_src sh
NXF_OPTS=-D"user.name=apraga" nextflow run main.nf -c nextflow.config -profile standard,helios --input="/Work/Projects/bisonex/centogene/fastq/2200467051_63003856/63003856_S135_R{1,2}_001.fastq.gz" --outdir=out -bg
cd out/63003856_S135_R/preprocessing/applybqsr/
samtools view 63003856_S135_R.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
make run BG= READS="out/63003856_S135_R/preprocessing/applybqsr/63003856_chr22_{1,2}.fq.gz"
cd out/63003856_chr22/preprocessing/mapped/
samtools index 63003856_chr22.bam
samtools mpileup 63003856_chr22.bam -r NC_000022.11:42213078-42213078
#+end_src
On récupère les 2 bam dans
#+begin_src
cd /home/alex/recherche/bisonex/code/BamScissors.jl/
rsync -avz meso:/Work/Users/apraga/bisonex/out/63003856_chr22/preprocessing/mapped data
rsync -avz meso:/Work/Users/apraga/bisonex/out/63003856_S135_R/preprocessing/applybqsr/ data/init/
#+end_src
********* Vérification que le reads est ailleurs
On cherche un read manquant dans le second alignement
#+begin_src sh
samtools view data/init/63003856_chr22.bam | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
#+end_src
#+RESULTS:
: A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NC_000022.11 42212845 0 151M = 42212883 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT ACC+FBCDCBBBAEAEDEEBBCCCECACBAEBEBDCCBCBFDCCCCFACEBEBCEEDCCCCFDCAEDCACBCEBBCFEACCFBDCACDCBCEBDBBCFEEDCCCFAFEACECCCECAEEDCADCBEDC7BEBCCCFBAFDCECCFBEAACA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
: A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NC_000022.11 42212883 0 151M = 42212845 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT AADECCCBDCBFCE<?CDEEEEBDEACDEAC;:BFBCBCDCCBEAEACAEFCCEAFBCBCCDEECBDBCECBEECCEACDEEBBFGDEFGCCFFFFCFCCEFBFDCFCDAAEBEE:CECBABBEBEE;DBFCCCDBCDBCCBBC?@BEEDA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
#+begin_src sh
samtools view data/mapped/63003856_chr22.bam | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
#+end_src
#+RESULTS:
: A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NW_014040930.1 115017 0 151M = 115055 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT ACC+FBCDCBBBAEAEDEEBBCCCECACBAEBEBDCCBCBFDCCCCFACEBEBCEEDCCCCFDCAEDCACBCEBBCFEACCFBDCACDCBCEBDBBCFEEDCCCFAFEACECCCECAEEDCADCBEDC7BEBCCCFBAFDCECCFBEAACA NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
: A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NW_014040930.1 115055 0 151M = 115017 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT AADECCCBDCBFCE<?CDEEEEBDEACDEAC;:BFBCBCDCCBEAEACAEFCCEAFBCBCCDEECBDBCECBEECCEACDEEBBFGDEFGCCFFFFCFCCEFBFDCFCDAAEBEE:CECBABBEBEE;DBFCCCDBCDBCCBBC?@BEEDA NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
Effectivement, on aligne sur une zonne supprimée !
******** DONE Corriger la qualité: non
CLOSED: [2023-05-24 Wed 22:19]
********* DONE Comparaison avec le fastq de référénce : qualité !!
CLOSED: [2023-05-24 Wed 22:17]
#+begin_src sh
cd /Work/Users/apraga/bisonex/work/6e/8548fc90263830bf677f36585f11dc
zgrep -A 3 "A00853:477:HMLWYDSX3:1:1413:4390:28573" 63003856_chr22_1.fq.gz
#+end_src
@A00853:477:HMLWYDSX3:1:1413:4390:28573
AGGGTTACCACCACCACCCTGACAGGAGATATTCTAGGAGTACTCAAGAGCATCAGGGGATGGCTGGTAGCCTAGAAGGAACCACAAGGCCCAATGTCTTGGTTAGTCAAACCAATGAATTAGCTAGCAGGGGCCTTCTGAACAAAAGCAT
+
ADEEB@?CBBCCBDCBDCCCFBD;EEBEBBABCEC:EEBEAADCFCDFBFECCFCFFFFCCGFEDGFBBEEDCAECCEEBCECBDBCEEDCCBCBFAECCFEACAEAEBCCDCBCBFB:;CAEDCAEDBEEEEDC?<ECFBCDBCCCEDAA
#+begin_src
zgrep -A 3 "A00853:477:HMLWYDSX3:1:1413:4390:28573" /Work/Projects/bisonex/centogene/fastq/2200467051_63003856/63003856_S135_R1_001.fastq.gz
#+end_src
#+RESULTS:
: @A00853:477:HMLWYDSX3:1:1413:4390:28573 1:N:0:ATTCCACACA+TAGGCGATTG
AGGGTTACCACCACCACCCTGACAGGAGATATTCTAGGAGTACTCAAGAGCATCAGGGGATGGCTGGTAGCCTAGAAGGAACCACAAGGCCCAATGTCTTGGTTAGTCAAACCAATGAATTAGCTAGCAGGGGCCTTCTGAACAAAAGCAT
: +
: FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFF
********* DONE Regarder la qualité après bwa mem vs applybqsr: différente
CLOSED: [2023-05-24 Wed 22:18]
Sur le mésocentre, dans /Work/Users/apraga/bisonex/out/63003856_S135_R/preprocessing
$ samtools view mapped/63003856_S135_R.bam NC_000022.11 | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NC_000022.11 42212845 0 151M = 42212883 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT FFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NC_000022.11 42212883 0 151M = 42212845 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT FFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
samtools view applybqsr/63003856_S135_R.bam NC_000022.11 | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NC_000022.11 42212845 0 151M = 42212883 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT ACC+FBCDCBBBAEAEDEEBBCCCECACBAEBEBDCCBCBFDCCCCFACEBEBCEEDCCCCFDCAEDCACBCEBBCFEACCFBDCACDCBCEBDBBCFEEDCCCFAFEACECCCECAEEDCADCBEDC7BEBCCCFBAFDCECCFBEAACA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NC_000022.11 42212883 0 151M = 42212845 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT AADECCCBDCBFCE<?CDEEEEBDEACDEAC;:BFBCBCDCCBEAEACAEFCCEAFBCBCCDEECBDBCECBEECCEACDEEBBFGDEFGCCFFFFCFCCEFBFDCFCDAAEBEE:CECBABBEBEE;DBFCCCDBCDBCCBBC?@BEEDA MC:Z:151M MD:Z:151 PG:Z:MarkDuplicates RG:Z:sample NM:i:0 AS:i:151 XS:i:151
********* DONE Réaligner à partir de la sortie de bwa mem
CLOSED: [2023-05-24 Wed 22:32]
#+begin_src sh
cd out/63003856_S135_R/preprocessing/mapped/
samtools view 63003856_S135_R.bam NC_000022.11 -f 0x2 -o 63003856_chr22.bam
samtools sort -n 63003856_chr22.bam -o 63003856_chr22_sorted.bam
samtools fastq -1 63003856_chr22_1.fq.gz -2 63003856_chr22_2.fq.gz -0 /dev/null -s /dev/null -n 63003856_chr22_sorted.bam
#+end_src
ON vérifie la qualité
#+begin_src
zgrep -A 3 "A00853:477:HMLWYDSX3:1:1413:4390:28573" 63003856_chr22_1.fq.gz
#+end_src
#+RESULTS:
: @A00853:477:HMLWYDSX3:1:1413:4390:28573
: AGGGTTACCACCACCACCCTGACAGGAGATATTCTAGGAGTACTCAAGAGCATCAGGGGATGGCTGGTAGCCTAGAAGGAACCACAAGGCCCAATGTCTTGGTTAGTCAAACCAATGAATTAGCTAGCAGGGGCCTTCTGAACAAAAGCAT
: +
: FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFF
#+begin_src
NXF_OPTS=-D"user.name=apraga" nextflow run main.nf -c nextflow.config -profile standard,helios -resume --input="out/63003856_S135_R/preprocessing/mapped/63003856_chr22_{1,2}.fq.gz" --outdir=out/63003856_chr22-from-mapped
#+end_src
Puis ::
#+begin_src
cd /Work/Users/apraga/bisonex/out/63003856_chr22-from-mapped/63003856_chr22/preprocessing/mapped
samtools view 63003856_chr22.bam | rg "A00853:477:HMLWYDSX3:1:1413:4390:28573"
#+end_src
#+RESULTS:
: A00853:477:HMLWYDSX3:1:1413:4390:28573 163 NW_014040930.1 115017 0 151M = 115055 189 CCCAGGGGCCCCAGTGGGGATTTTCTAATAGAGACCCAATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACT FFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
: A00853:477:HMLWYDSX3:1:1413:4390:28573 83 NW_014040930.1 115055 0 151M = 115017 -189 ATGCTTTTGTTCAGAAGGCCCCTGCTAGCTAATTCATTGGTTTGACTAACCAAGACATTGGGCCTTGTGGTTCCTTCTAGGCTACCAGCCATCCCCTGATGCTCTTGAGTACTCCTAGAATATCTCCTGTCAGGGTGGTGGTGGTAACCCT FFFFFFFFFFFFFF::FFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:151 MC:Z:151M AS:i:151 XS:i:151 RG:Z:sample
******** DONE Aligner sur génome de référence limité au chromosome 22
CLOSED: [2023-05-24 Wed 23:18]
********* KILL Test données non modifiées
CLOSED: [2023-05-24 Wed 23:18]
/Work/Users/apraga/bisonex/tests/bamscissors
#+begin_src
cd /Work/Groups/bisonex/data/genome/GRCh38.p13/
mkdir chr22/
samtools faidx genomeRef.fna NC_000022.11 > chr22/chr22.fna
cd chr22
samtools faidx chr22.fna
bwa index chr22.fna
#+end_src
#+begin_src
cd /Work/Users/apraga/bisonex/tests/bamscissors
ln -s ../ ../out/63003856_S135_R/preprocessing/applybqsr/63003856_chr22_{1,2}.fq.gz .
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/chr22/chr22.fna 63003856_chr22_1.fq.gz 63003856_chr22_1.fq.gz -o smallref.sam
#+end_src
********* DONE Test données modifiées: ok
CLOSED: [2023-05-24 Wed 23:18]
Données dans data/init
#+begin_src sh
time julia insertVariant.jl
rsync -avz data/init/*.fq.gz meso:/Work/Users/apraga/bisonex/tests/bamscissors/
#+end_src
#+begin_src
srun -c 24 -p smp -t 1:00:00 --pty bash
bwa mem -t 24 /Work/Projects/bisonex/data/genome/GRCh38.p13/chr22/chr22.fna 63003856_chr22_1.fq.gz 63003856_chr22_1.fq.gz | samtools sort -@24 - -o smallref.bam
#+end_src
#+begin_src
rsync -avz meso:/Work/Users/apraga/bisonex/tests/bamscissors/smallref.bam mapped/
#+end_src
******* DONE Test haplotypecaller 1 variant
CLOSED: [2023-05-29 Mon 15:38]
******* DONE Test haplotypecaller tous les variants
******* DONE Comprendre pourquoi la répartiton ne suit pas la loi normale
CLOSED: [2023-06-01 Thu 21:44]
Certains hétérozygote soint à 0.01 ou 1...
******** DONE augmenter le nombre d'échantillions: idem
CLOSED: [2023-05-31 Wed 22:24]
******** DONE Vérifier le nombre de reads marqué vs édité
CLOSED: [2023-06-01 Thu 21:44]
******** DONE vérifier que 100 appel à rand(d, 1)[1] est semblable à un appel de rand(d, 100)
CLOSED: [2023-05-31 Wed 22:24]
julia> df = vcat(DataFrame(:y => z, :type => "z"), DataFrame(:y => y, :type => "y"));
julia> y = [rand(d, 1)[1] for x in 1:1000];
julia> z = rand(d,1000);
julia> df = vcat(DataFrame(:y => z, :type => "z"), DataFrame(:y => y, :type => "y"));
draw(data(df)*histogram(bins=100)*mapping(:y, color=:type,dodge=:type))
******* DONE Améliorer les performances
CLOSED: [2023-06-02 Fri 23:39]
#+begin_src julia
@time include("xamscissors.jl")
#+end_src
430s pour chromosome 22. Majorité dans l'édition de reads:
******** DONE Inserér tous les variants d'un reads d'un coup
CLOSED: [2023-06-01 Thu 23:09]
Ne change rien
******** DONE Test avec -t4: idem
CLOSED: [2023-06-01 Thu 23:17]
******** DONE Test mésocentre : idem
CLOSED: [2023-06-01 Thu 23:40]
348s
******** Changer la structure de données des
Dataframe -> dict = les performances horribles ont disparuse
******* KILL Refaire le test avec la nouvelle version
CLOSED: [2023-06-12 Mon 23:27]
******** DONE Génération des données
CLOSED: [2023-06-02 Fri 23:40]
Mésocentre
#+begin_src sh
cd /Work/Users/apraga/bisonex/out/63003856_S135_R/preprocessing/mapped
samtools view 63003856_S135_R.bam NC_000022.11 -o 63003856_S135_R_chr22.bam
samtools index 63003856_S135_R_chr22.bam
cp 63003856_S135_R_chr22.bam* /Work/Users/apraga/bisonex/tests/xamscissors/
cd /Work/Users/apraga/bisonex/tests/xamscissors
#+end_src
On génère les données
#+begin_src julia
using XAMScissors
insertSNV("./63003856_S135_R_chr22.bam", "snvs_chr22.csv", "out")
#+end_src
Puis
#+begin_src sh
julia xamscissors.jl
#+end_src
******** DONE Run
CLOSED: [2023-06-03 Sat 18:26]
NXF_OPTS=-D"user.name=${USER}" nextflow run workflows/runInserted.nf -profile standard,helios --input="tests/xamscissors/out/inserted_{1,2}.fq.gz"
******** DONE Après haplotypecaller : ok
CLOSED: [2023-06-03 Sat 18:27] SCHEDULED: <2023-06-03 Sat>
******** KILL Après filtre vep
CLOSED: [2023-06-12 Mon 23:27] SCHEDULED: <2023-06-03 Sat>
****** KILL PHase 3 : tous les SNV
CLOSED: [2023-06-12 Mon 23:28] SCHEDULED: <2023-06-03 Sat>
******* DONE Générer les données
CLOSED: [2023-06-03 Sat 20:16] SCHEDULED: <2023-06-03 Sat>
#+begin_src julia
using XAMScissors
insertSNV("../../out/63003856_S135_R/preprocessing/mapped/63003856_S135_R.bam", "snvs.csv", "out")
#+end_src
temps d'exécution 73min
#+begin_src sh
nohup bash -c 'time julia xamscissors.jl' &
xamscissors-63003856/*.fq.gz /Work/Groups/bisonex/data/xamscissors/
#+end_src
******* DONE Regénérer avec @time pour avoir les performaces
CLOSED: [2023-06-03 Sat 21:45]
markReads 6.265202 seconds (1.36 M allocations: 137.090 MiB, 1.00% gc time, 9.79% compilation time)
editReads 1327.701623 seconds (1.03 G allocations: 81.996 GiB, 0.59% gc time, 0.03% compilation time)
samtools index 117.743727 seconds (53 allocations: 1.883 KiB)
samtools sort 2820.074930 seconds (66 allocations: 2.789 KiB)
bam2fastq 134.148952 seconds (794 allocations: 40.539 KiB, 0.01% compilation time)
real 73m33.273s
user 77m38.194s
sys 1m26.684s
[bam_sort_core] merging from 60 files and 1 in-memory blocks...
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 126905130 reads
real 73m6.934s
user 77m25.397s
sys 1m21.339s
******* KILL Après haplotypecaller 556/590, majorité = échec alignement
CLOSED: [2023-06-12 Mon 23:27] SCHEDULED: <2023-06-04 Sun>
Haplotypecaller 556 found over 590
Amongst 34 missed variant, 2 have a mapping quality > 0
2×7 DataFrame
Row │ chrom pos ref alt zygosity meanQual stdQual
│ String15 Int64 String1 String1 String7 Float64 Float64
─────┼─────────────────────────────────────────────────────────────────────────
1 │ NC_000017.11 39672244 G A het 60.0 0.0
2 │ NC_000001.11 155235252 A G het 0.258065 2.48868
NC_000017.11 39672244 G A het => ok, problème de représentation car 2 variant côte à cote
NC_000001.11 155235252 A G het => peu de reads alternatifs (9/93 donc ok)
Position: chromoe 1 et 6 surtout
34×7 DataFrame
Row │ chrom pos ref alt zygosity
│ String15 Int64 String1 String1 String7
─────┼──────────────────────────────────────────────────────
1 │ NC_000001.11 153817496 A T het
2 │ NC_000001.11 155235252 A G het
3 │ NC_000001.11 155236268 G A het
4 │ NC_000001.11 155290591 C T het
5 │ NC_000001.11 155291918 G A het
6 │ NC_000001.11 155294358 G T het
7 │ NC_000002.12 149010343 C T het
8 │ NC_000006.12 32039426 T A het
9 │ NC_000006.12 32040110 G T het
10 │ NC_000006.12 32040723 G A het
11 │ NC_000006.12 32041006 C T het
12 │ NC_000006.12 32041147 G A het
13 │ NC_000006.12 33443054 G T het
14 │ NC_000006.12 33451815 C T het
15 │ NC_000006.12 170283230 C A het
16 │ NC_000006.12 170283754 G A het
17 │ NC_000006.12 170285637 T C het
18 │ NC_000006.12 170289678 A C het
19 │ NC_000010.11 87961118 G A het
20 │ NC_000012.12 2449086 C G het
21 │ NC_000015.10 74343027 C T het
22 │ NC_000016.10 16163078 G A het
23 │ NC_000016.10 21262032 C G het
24 │ NC_000016.10 21962506 C T homo
25 │ NC_000017.11 7513122 C T het
26 │ NC_000017.11 7513752 C T het
27 │ NC_000017.11 39672244 G A het
28 │ NC_000017.11 46018710 C T het
29 │ NC_000019.10 54144058 G A het
30 │ NC_000021.9 43063074 A G het
31 │ NC_000021.9 43426167 C T het
32 │ NC_000022.11 18918421 A G het
33 │ NC_000022.11 42087168 T A homo
34 │ NC_000022.11 42213078 T G het
******* DONE Voir où est l'alignement alternatif: sur NW_ (zone supprimée)
CLOSED: [2023-06-04 Sun 22:15] SCHEDULED: <2023-06-04 Sun>
ex chr15 74343027
A00853:477:HMLWYDSX3:2:2444:22354:28870
#+begin_src
cd /Work/Groups/bisonex/data/xamscissors
zgrep -A4 "A00853:477:HMLWYDSX3:2:2444:22354:28870" *.fq.gz
#+end_src
63003856_xamscissors_1.fq.gz:@A00853:477:HMLWYDSX3:2:2444:22354:28870
63003856_xamscissors_1.fq.gz:CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTTGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC
63003856_xamscissors_1.fq.gz:+
63003856_xamscissors_1.fq.gz:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFF:FFFFFFFFFF::FFFFFFFFFFF:FFFFFFFFFFFFFF:FFFFFFF,FFFFFF,FFFFFFFFFFFF:FF::FF
63003856_xamscissors_2.fq.gz:@A00853:477:HMLWYDSX3:2:2444:22354:28870
63003856_xamscissors_2.fq.gz:GACAGAAAGGAAGTGTTCACCACGATTACCGTGGCATCCTCTACAGACTCCTGGGAGACAGCAAGATGTCCTTCGAGGACATCAAGGCCAACGTCACAGAGATGCCGGCAGGAGGGGTGGACACGGTG
63003856_xamscissors_2.fq.gz:+
63003856_xamscissors_2.fq.gz:FF:FFF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFF:F:FF:FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFF,:FFF,FFFFFF:FFFFFFFFFFFFF
******** DONE Avec BLAT: sur _fix
CLOSED: [2023-06-04 Sun 21:07]
1er =
ACTIONS QUERY SCORE START END QSIZE IDENTITY CHROM STRAND START END SPAN
--------------------------------------------------------------------------------------------------------------
browser details YourSeq 124 1 128 128 98.5% chr15_ML143370v1_fix + 172243 172370 128 What is chrom_fix?
browser details YourSeq 124 1 128 128 98.5% chr15 + 74342974 74343101 128
browser details YourSeq 23 1 25 128 96.0% chr19 - 33396097 33396121 25
Second
--------------------------------------------------------------------------------------------------------------
browser details YourSeq 126 1 128 128 99.3% chr15_ML143370v1_fix - 172243 172370 128 What is chrom_fix?
browser details YourSeq 126 1 128 128 99.3% chr15 - 74342974 74343101 128
browser details YourSeq 23 104 128 128 96.0% chr19 + 33396097 33396121 25
******** DONE Bwa mem à la main GRCh38.p13 : on est dans une zone NW
CLOSED: [2023-06-04 Sun 21:51]
On met les 2 reads dans des fichiers séparés puis
#+begin_src sh
cd /Work/Users/apraga/bisonex/tests/xamscissors/align
bwa mem /Work/Groups/bisonex/data/genome/GRCh38.p13/bwa/genomeRef test1.fq test2.fq
#+end_src
A00853:477:HMLWYDSX3:2:2444:22354:28870 97 NW_021160016.1 172243 0 128M = 172243 128 CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTTGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFF:FFFFFFFFFF::FFFFFFFFFFF:FFFFFFFFFFFFFF:FFFFFFF,FFFFFF,FFFFFFFFFFFF:FF::FF NM:i:2 MD:Z:22A30C7MC:Z:128M AS:i:118 XS:i:118 XA:Z:NC_000015.10,+74342974,128M,2;
A00853:477:HMLWYDSX3:2:2444:22354:28870 145 NW_021160016.1 172243 0 128M = 172243 -128 CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTCGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC FFFFFFFFFFFFF:FFFFFF,FFF:,FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFF:FF:F:FFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FF:FFF:FF NM:i:1 MD:Z:22A105 MC:Z:128M AS:i:123 XS:i:123 XA:Z:NC_000015.10,-74342974,128M,1;
******** DONE GRCh38.p14: idem
CLOSED: [2023-06-04 Sun 21:51]
A00853:477:HMLWYDSX3:2:2444:22354:28870 97 NW_021160016.1 172243 0 128M = 172243 128 CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTTGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFF:FFFFFFFFFF::FFFFFFFFFFF:FFFFFFFFFFFFFF:FFFFFFF,FFFFFF,FFFFFFFFFFFF:FF::FF NM:i:2 MD:Z:22A30C7MC:Z:128M AS:i:118 XS:i:118 XA:Z:NC_000015.10,+74342974,128M,2;
A00853:477:HMLWYDSX3:2:2444:22354:28870 145 NW_021160016.1 172243 0 128M = 172243 -128 CACCGTGTCCACCCCTCCTGCCGGCATCTCTGTGACGTTGGCCTTGATGTCCTCGAAGGACATCTTGCTGTCTCCCAGGAGTCTGTAGAGGATGCCACGGTAATCGTGGTGAACACTTCCTTTCTGTC FFFFFFFFFFFFF:FFFFFF,FFF:,FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFF:FF:F:FFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FF:FFF:FF NM:i:1 MD:Z:22A105 MC:Z:128M AS:i:123 XS:i:123 XA:Z:NC_000015.10,-74342974,128M,1;
******** DONE GRCh38 : ok
CLOSED: [2023-06-04 Sun 22:15]
bwa mem /Work/Projects/bisonex/data/genome/GRCh38/GCA_000001405.15_GRCh38_full_analysis_set.fna test1.fq test2.fq
******* DONE Vérifier que les reads ont la même qualité sur les fichiers d'origine: oui
CLOSED: [2023-06-04 Sun 21:07]
******* DONE Supprimer les NW_ ?
CLOSED: [2023-06-10 Sat 10:40] SCHEDULED: <2023-06-04 Sun>
@A00853:477:HMLWYDSX3:3:2114:14742:8860
CAGGCCAGCCGCTCAGCCCGCTCCTTTCACCCTCTGCAGGAGAGCCTCGTGGCAGGCCAGTGGAGGGACATGATGGACTACATGCTCCAAGGGGTGGCGCAGCCGAGCATGGAAGAGGGCTCTGGACAGCTCCTGGAAGGGCACTTGCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@A00853:477:HMLWYDSX3:3:2114:14742:8860
CTTTTGCTTGTCCCCAGGACGCACCTCAGGGTGGTGAAGCAAAAAAACCACGGCCCAGGAGAGGGTGGGTGCTGTGGTCTCAGTGCCACCGATCAGGAGGTCCACTGCAGCCATGTGCAAGTGCCCTTCCAGGAGCTGTCCAGAGCCCTCT
+
FFFFFFFFFFFFFFFFFFFFFFF:FFF:FFFFFFFFFFFFF,FFFFFFFFFFFF:F:FFFF:FFFFF,,FFF:FFFFFFFFFF,FFFFFFF,FFFFFFFFFFF,FFFFFFFFF:FFFF,F:FFFFF:FFFFFFFFF:FFFF,FFFFFFFFF
******* DONE Supprimer NW_ et NT_
****** TODO Phase 2 : chr22, vaf variable :T2T:
SCHEDULED: <2023-06-16 Fri>
CLOSED: [2023-06-12 Mon 23:27]
****** TODO Phase 3 : tous SNV, vaf variable :T2T:
SCHEDULED: <2023-06-16 Fri>
CLOSED: [2023-06-12 Mon 23:27]
***** TODO Test Indel
**** Divers
***** DONE Vérifier nombre de reads fastq - bam
CLOSED: [2022-10-09 Sun 22:31]
*** TODO Liste varants "clinically relevent" (Clinge - CT-R d)
SCHEDULED: <2023-06-25 Sun>
[cite:@wilcox2021]