apraga/org - Change F4SKQVRCFPVCVXG7JGG3UONESIH5TMSQ2Q7DI2U7VX4R4TBP5SAAC

Bisonex save

Created by Alexis Praga on June 13, 2023

F4SKQVRCFPVCVXG7JGG3UONESIH5TMSQ2Q7DI2U7VX4R4TBP5SAAC

Dependencies

In channels

main

Change contents

Replacement in projects/bisonex.org at line 2 [5.35]

B:BD[3.8203] → [4.29:8221]

B:BD[4.8221] → [2.73:8265]

 be explained by their divergence in
 sequencing strategy that producing different length of reads (all BGI platforms
 were 100 base pair read length while all Illumina platforms were 150 base pair
 read length). The read length effects, as a key factor between two platforms,
 would bring alignment bias and error which are higher for short reads and
 ultimately affect the variants calling especially the INDELs identification
#+end_quote
*** Débugger variant calling (haplotypecaller)
https://gatk.broadinstitute.org/hc/en-us/articles/360043491652-When-HaplotypeCaller-and-Mutect2-do-not-call-an-expected-variant
https://gatk.broadinstitute.org/hc/en-us/articles/360035891111-Expected-variant-at-a-specific-site-was-not-called
*** Hap.py
Format de sortie :
#+begin_src r
vcf_field_names(vcf, tag = "FORMAT")
#+end_src
#+RESULTS:
: FORMAT BD    1      String  Decision for call (TP/FP/FN/N)
: FORMAT BK    1      String  Sub-type for decision (match/mismatch type)
: FORMAT BVT   1      String  High-level variant type (SNP|INDEL).
: FORMAT BLT   1      String  High-level location type (het|homref|hetalt|homa
am = genotype mismatch
lm = allele/haplotype mismatch
. = non vu
**** On vérifie que am = genotype mismatch
référence  = T/T
high-confidence = T/C
notre = C/C
#+begin_src sh
bcftools filter -i 'POS=19196584'  /Work/Groups/bisonex/data/giab/GRCh38/HG001_GRCh38_1_22_v4.2.1_benchmark.vcf.gz | grep -v '#'
bcftools filter -i 'POS=19196584'  ../out/NA12878_NIST7035-dbsnp/variantCalling/haplotypecaller/NA12878_NIST.vcf.gz | grep -v '#'
#+end_src
#+RESULTS:
: NC_000022.11    19196584        .       T       C       50      PASS    platforms=5;platformnames=Illumina,PacBio,10X,Ion,Solid;datasets=5;datasetnames=HiSeqPE300x,CCS15kb_20kb,10XChromiumLR,IonExome,SolidSE75bp;callsets=7;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbDV,CCS15kb_20kbGATK4,HiSeqPE300xfreebayes,10XLRGATK,IonExomeTVC,SolidSE75GATKHC;datasetsmissingcall=CGnormal;callable=CS_HiSeqPE300xGATK_callable,CS_CCS15kb_20kbDV_callable,CS_10XLRGATK_callable,CS_CCS15kb_20kbGATK4_callable,CS_HiSeqPE300xfreebayes_callable GT:PS:DP:ADALL:AD:GQ    0/1:.:781:109,123:138,150:348
: NC_000022.11    19196584        rs1061325       T       C       59.32   PASS    AC=2;AF=1;AN=2;DB;DP=2;ExcessHet=0;FS=0;MLEAC=1;MLEAF=0.5;MQ=60;QD=29.66;SOR=2.303      GT:AD:DP:GQ:PL  1/1:0,2:2:6:71,6,0
**** On vérifie que lm = allele/haplotype mismatch
référence  = CAA/CAA
high-confidence = CA/CA
notre = C/CA
#+begin_src sh
 bcftools filter -i 'POS=31277416'  /Work/Groups/bisonex/data/giab/GRCh38/HG001_GRCh38_1_22_v4.2.1_benchmark.vcf.gz | grep -v '#'
 bcftools filter -i 'POS=31277416'  ../out/NA12878_NIST7035-dbsnp/variantCalling/haplotypecaller/NA12878_NIST.vcf.gz | grep -v '#'
#+end_src
#+RESULTS:
: NC_000022.11    31277416        .       CA      C       50      PASS    platforms=3;platformnames=Illumina,PacBio,10X;datasets=3;datasetnames=HiSeqPE300x,CCS15kb_20kb,10XChromiumLR;callsets=4;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbDV,10XLRGATK,HiSeqPE300xfreebayes;datasetsmissingcall=CCS15kb_20kb,CGnormal,IonExome,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable;difficultregion=GRCh38_AllHomopolymers_gt6bp_imperfectgt10bp_slop5,GRCh38_SimpleRepeat_imperfecthomopolgt10_slop5  GT:PS:DP:ADALL:AD:GQ    1/1:.:465:16,229:0,190:129
: NC_000022.11    31277416        rs57244615      CAA     C,CA    389.02  PASS    AC=1,1;AF=0.5,0.5;AN=2;BaseQRankSum=0.37;DB;DP=37;ExcessHet=0;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;MQRankSum=0;QD=13.41;ReadPosRankSum=-0.651;SOR=0.572    GT:AD:DP:GQ:PL  1/2:5,10,14:29:64:406,202,313,64,0,88
*** Génération de reads
Biblio récente
https://www.biorxiv.org/content/10.1101/2022.03.29.486262v1.full.pdf
Parmi ceux qui gèrent les variations
- *simuscop* reads non centré sur les zones de capture
- *NEAT: exome* mais trop lent en pratique
- *Reseq* exome
- gensim : pas d'exome
- pIRS : non plus
- varsim : non plus
  ...
  Temps de calcul selon l'article de reseq https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02265-7
  #+begin_quote
  Due to ReSeq’s effective parallelization, its elapsed times are low for this benchmark with 48 virtual CPUs (Additional file 1: Figure S34b,e). In contrast, the single-threaded processes implemented in perl or python have strikingly high elapsed times. This is well visible in Hs-HiX-TruSeq and applies to the training of pIRS (over a week), NEAT (several days), and BEAR (half a week) as well as the simulation of NEAT (close to 2 weeks) and BEAR (several weeks).
Biblio : https://www.nature.com/articles/s41437-022-00577-3
  #+end_quote
Divers
- Liste ancienne : https://www.biostars.org/p/128762/
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02265-7
* Idées
** Validation analytique
mail Yannis : données patients +/- simulées
*** Utiliser données GCAT et uploader le notre ?
https://www.nature.com/articles/ncomms7275
*** [#A] Variant calling : Genome in a bottle : NA12878 + autres
Résumé : https://www.nist.gov/programs-projects/genome-bottle
Manuscript : https://www.nature.com/articles/s41587-019-0054-x.epdf?author_access_token=E_1bL0MtBBwZr91xEsy6B9RgN0jAjWel9jnR3ZoTv0OLNnFBR7rUIZNDXq0DIKdg3w6KhBF8Rz2RWQFFc0St45kC6CZs3cDYc87HNHovbWSOubJHDa9CeJV-pN0BW_mQ0n7cM13KF2JRr_wAAn524w%3D%3D
Article comparant les variant calling : https://www.biorxiv.org/content/10.1101/2020.12.11.422022v1.full.pdf
**** KILL Tester le séquencage aussi
CLOSED: [2023-01-30 lun. 18:30]
Depuis un fastq correspondant à Illumina  https://github.com/genome-in-a-bottle/giab_data_indexes
   puis on compare le VCF avec les "high confidence"
On séquence directement NA12878 -> inutile pour le pipeline seul
**** TODO Tester seul la partie bioinformatique
   Tout résumé ici : https://www.nist.gov/programs-projects/genome-bottle
- methode https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/NA12878/analysis/Illumina_PlatinumGenomes_NA12877_NA12878_09162015/IlluminaPlatinumGenomes-user-guide.pdf
- vcf
     https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/NA12878_HG001/latest/GRCh38/
NB: à quoi correspond https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/NA12878/analysis/Illumina_PlatinumGenomes_NA12877_NA12878_09162015/hg38/2.0.1/NA12878/ ??
   Article comparant les variant calling : https://www.biorxiv.org/content/10.1101/2020.12.11.422022v1.full.pdf
   Article pour vcfeval : https://www.nature.com/articles/s41587-019-0054-x
   La version 4 ajoute 273 gènes "clinically relevant" https://www.biorxiv.org/content/10.1101/2021.06.07.444885v3.full.pdf
   Ajout des zones "difficiles"
   https://www.biorxiv.org/content/10.1101/2020.07.24.212712v5.full.pdf
*** [#B] Pipeline : générer patient avec tous les variants retrouvés à Centogene
Comparaison de génération ADN (2019)
https://academic.oup.com/bfg/article/19/1/49/5680294
**** SimuSCop (exome)
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-020-03665-5
https://github.com/qasimyu/simuscop
1. Crééer un modèle depuis bam + vcf : Setoprofile
2. Génerer données NGS
** Annotation :
*** Comparaison vep / snpeff et annovar
* Changement nouvelle version
- Dernière version du génome (la version "prête à l'emploi" est seulement GRCh38 sans les version patchées)
* Notes
** Nextflow
*** afficher les résultats d'un process/workflow
#+begin_src
lol.out.view()
#+end_src
Attention, ne fonctionne pas si plusieurs sortie:
#+begin_src
lol.out[0].view()
#+end_src
ou si /a/ est le nom de la sortie
#+begin_src
lol.out.a.view()
#+end_src
** Quelle version du génome ?
Il y a 2 notations pour les chrosome: Refseq (NC_0001) ou chr1, chr2...
dbSNP utilise Refseq
pour le fasta, 2 solutions
- refseq : "https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/${genome}_latest/refseq_identifiers/${fna}.gz"
  -> nécessite d'indexer le fichier (long !)
- chromosome https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/
  -> nécessite d'annoter les chromosomes pour corriger
 (avec le fichier gff)
  On utilise la version chromosome donc on annote dbSNP (à faire)
** Performances
Ordinateur de Carine (WSL2) : 4h dont 1h15 alignement (parallélisé) et 1h15 haplotypecaller (séquentiel)
** Chromosomes NC, NT, NW
Correspondance :
https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&chromInfoPage=
Signification
https://genome.ucsc.edu/FAQ/FAQdownloads.html#downloadAlt
- alt = séquences alternatives (utilisables)
- fix = patch (correction ou amélioration)
- random = séquence connue sur un chromosome mais non encore utilisée
** Pipelines prêt-à-l’emploi nextflow
Problème : nécessite singularity ou docker (ou conda)
Potentiellement utilisable avec nix...
** Validation : Quelles données de référence ?
Discussion avec Alexis
- Platinum genomes = génome seul
*** [[https://github.com/genome-in-a-bottle/giab_data_indexes][Genome in a bottle]]
  - NA12878 :
    - Illumina HiSeq Exome : fastq + capture
    - Illumina TruSeq Exome : bam, pas de capture
      ici ww
  - HG002,3,4
    - Illumina Whole Exome  : bam. le kit de capture est "Agilent SureSelect Human All Exon V5 kit" selon [[https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/README.txt][README]]. On il faut les régions [[https://kb.10xgenomics.com/hc/en-us/articles/115004150923-Where-can-I-find-the-Agilent-Target-BED-files-][selon ce site]]
      Un autre fichier est disponible (capture ???)
    https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/wex_Agilent_SureSelect_v05_b37.baits.slop50.merged.list
  "target region" +/- 50bp
    testé sur chr311780-312086 : ok
Autres technologies non adaptées au pipeline (vu avec Alexis)
*** [[https://www.illumina.com/platinumgenomes.html][Platinum genome
]] Que du génome « sequenced to 50x depth on a HiSeq 2000 system”
Genome possible
** Zone de capture
GIAB fourni le .bed pour l'exome . INfo : https://support.illumina.com/sequencing/sequencing_kits/nextera-rapid-capture-exome-kit/downloads.html
** Centogène
https://www.twistbioscience.com/node/23906
Bed non fourni pour exactement cette capture
On prend https://www.twistbioscience.com/resources/data-files/twist-alliance-vcgs-exome-401mb-bed-files
qui content la majeure partie
* Données :data:
** DONE Remplacer bam par fastq sur mesocentre
CLOSED: [2023-04-16 Sun 16:33]
Commande
*** DONE Supprimer les fastq non "paired"
CLOSED: [2023-04-16 Sun 16:33]
nushell
Liste des fastq avec "paired-end" manquant
#+begin_src nu
ls **/*.fastq.gz | get name | path basename | split column "_" | get column1 | uniq -u | save single.txt
#+end_src
#+RESULTS:
: 62907927
: 62907970
: 62899606
: 62911287
: 62913201
: 62914084
: 62915905
: 62921595
: 62923065
: 62925220
: 62926503
: 62926502
: 62926500
: 62926499
: 62926498
: 62931719
: 62943423
: 62943400
: 62948290
: 62949205
: 62949206
: 62949118
: 62951284
: 62960792
: 62960785
: 62960787
: 62960617
: 62962561
: 62962692
: 62967473
: 62972194
: 62979102
On vérifie
#+begin_src nu
open single.txt  | lines | each {|e| ls $"fastq/*_($in)/*" | get 0  }
open single.txt  | lines | each {|e| ls $"fastq/*_($in)/*" | get 0.name }  | path basename | split column "_" | get column1 | uniq -c
#+end_src
On met tous dans un dossier (pas de suppression )
#+begin_src
open single.txt  | lines | each {|e| ls $"fastq/*_($in)/*" | get 0  }  | each {|e| ^mv $e.name bad-fastq/}
#+end_src
On vérifie que les dossiier sont videsj
 open single.txt  | lines | each {|e| ls $"fastq/*_($in)" | get 0.name } | ^ls -l $in
 Puis on supprime
 open single.txt  | lines | each {|e| ls $"fastq/*_($in)" | get 0.name } | ^rm -r $in
*** DONE Supprimer bam qui ont des fastq
CLOSED: [2023-04-16 Sun 16:33]
On liste les identifiants des fastq et bam dans un tableau avec leur type :
#+begin_src
let fastq = (ls fastq/*/*.fastq.gz | get name | parse "{dir}/{full_id}/{id}_{R}_001.fastq.gz"  | select dir id | uniq )
let bam = (ls bam/*/*.bam | get name | parse "{dir}/{full_id}/{id}_{S}.bqrt.bam"  | select dir id)
#+end_src
On groupe les résultat par identifiant (résultats = liste de records qui doit être convertie en table)
et on trie ceux qui n'ont qu'un fastq ou un bam
#+begin_src
let single = ( $bam | append $fastq | group-by id | transpose id files | get files | where {|x| ($x | length) == 1})
#+end_src
On convertit en table et on récupère seulement les bam
#+begin_src
$single | reduce {|it, acc| $acc | append $it} | where dir == bam | get id | each {|e| ^ls $"bam/*_($e)/*.bam"}
#+end_src
#+RESULTS:
: bam/2100656174_62913201/62913201_S52.bqrt.bam
: bam/2100733271_62925220/62925220_S33.bqrt.bam
: bam/2100738763_62926502/62926502_S108.bqrt.bam
: bam/2100746726_62926498/62926498_S105.bqrt.bam
: bam/2100787936_62931955/62931955_S4.bqrt.bam
: bam/2200066374_62948290/62948290_S130.bqrt.bam
: bam/2200074722_62948298/62948298_S131.bqrt.bam
: bam/2200074990_62948306/62948306_S218.bqrt.bam
: bam/2200214581_62967331/62967331_S267.bqrt.bam
: bam/2200225399_62972187/62972187_S85.bqrt.bam
: bam/2200293962_62979117/62979117_S63.bqrt.bam
: bam/2200423985_62999352/62999352_S1.bqrt.bam
: bam/2200495073_63010427/63010427_S20.bqrt.bam
: bam/2200511274_63012586/63012586_S114.bqrt.bam
: bam/2200669188_63036688/63036688_S150.bqrt.bam
* Nouveau workflow :workflow:
** TODO Bases de données
*** KILL Nix pour télécharger les données brutes
**** Conclusion
Non viable sur cluster car en dehors de /nix/store
On peut utiliser des symlink mais trop compliqué
**** KILL Axel au lieu de curl pour gérer les timeout?
CLOSED: [2022-08-19 Fri 15:18]
*** DONE Tester patch de @pennae pour gros fichiers
SCHEDULED: <2022-08-19 Fri>
*** KILL Télécharger les données avec nextflow: hg38
CLOSED: [2023-06-12 Mon 23:29]
**** DONE Genome de référence
**** DONE dbSNP
**** DONE VEP 20G
CLOSED: [2023-06-12 Mon 23:29]
Ajout vérification checksum -> à vérifier
**** DONE transcriptome (spip)
CLOSED: [2023-06-12 Mon 23:29]
Rajouter checksum manuel
**** KILL Refseq
**** KILL OMIM
CLOSED: [2023-06-12 Mon 23:29]
codé, à vérifier
**** KILL ACMG incidental
CLOSED: [2023-06-12 Mon 23:29]
*** TODO Télécharge données :TIT:
SCHEDULED: <2023-06-15 Thu>
**** DONE fasta
CLOSED: [2023-06-12 Mon 23:30]
**** DONE fasta : compatibilité hg38
CLOSED: [2023-06-12 Mon 23:30]
**** DONE fasta index
CLOSED: [2023-06-13 Tue 00:07]
**** DONE fasta index: compatibilité hg38
CLOSED: [2023-06-13 Tue 00:07]
**** DONE fasta dictionnaire
CLOSED: [2023-06-13 Tue 00:07]
**** DONE dbSNP
CLOSED: [2023-06-12 Mon 23:30]
**** DONE dbSNP : compatibilité hg38
CLOSED: [2023-06-12 Mon 23:30]
*** HOLD Processing bases de données
**** DONE dbSNP common
**** DONE Seulement les ID dans dbSNP common !
CLOSED: [2022-11-19 Sat 21:42]
172G au lieu de 253M...
**** HOLD common dbSNP not clinvar patho
***** DONE Conclusion partielle
CLOSED: [2022-12-12 Mon 22:25]
- vcfeval : prometteur mais n'arrive pas à traiter toutes les régions
- isec : trop de problèmes avec
- classif clinvar directement dans dbSNP: le plus simple
  Et ça permet de rattraper quelques erreurs dans le script d'Alexis
***** KILL Utiliser directement le numéro dbSNP dans clinvar ? Non
CLOSED: [2022-11-20 Sun 19:51]
Ex: chr20
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f 'rs%INFO/RS \n' -i 'INFO/RS != "." & INFO/CLNSIG="Pathogenic"' clinvar_chr20.vcf.gz | sort > ID_clinvar_patho.txt
bcftools query -f '%ID\n' dbSNP_common_chr20.vcf.gz | sort > ID_of_common_snp.txt
comm -23 ID_of_common_snp.txt ID_clinvar_patho.txt > ID_of_common_snp_not_clinvar_patho.txt
wc -l ID_of_common_snp_not_clinvar_patho.txt
# sort ID
#+end_src
#+RESULTS:
: 518846 ID_of_common_snp_not_clinvar_patho.txt
Version d'alexis
#+begin_src sh :dir ~/code/bisonex/test_isec
snp=dbSNP_common_chr20.vcf.gz
clinvar=clinvar_chr20_notremapped.vcf.gz
python ../script/pythonScript/clinvar_sbSNP.py \
    --clinvar $clinvar \
    --chrm_name_table ../

[3.8203]

[2.8265]

 be explained by their divergence in
 sequencing strategy that producing different length of reads (all BGI platforms
 were 100 base pair read length while all Illumina platforms were 150 base pair
 read length). The read length effects, as a key factor between two platforms,
 would bring alignment bias and error which are higher for short reads and
 ultimately affect the variants calling especially the INDELs identification
#+end_quote
*** Débugger variant calling (haplotypecaller)
https://gatk.broadinstitute.org/hc/en-us/articles/360043491652-When-HaplotypeCaller-and-Mutect2-do-not-call-an-expected-variant
https://gatk.broadinstitute.org/hc/en-us/articles/360035891111-Expected-variant-at-a-specific-site-was-not-called
*** Hap.py
Format de sortie :
#+begin_src r
vcf_field_names(vcf, tag = "FORMAT")
#+end_src
#+RESULTS:
: FORMAT BD    1      String  Decision for call (TP/FP/FN/N)
: FORMAT BK    1      String  Sub-type for decision (match/mismatch type)
: FORMAT BVT   1      String  High-level variant type (SNP|INDEL).
: FORMAT BLT   1      String  High-level location type (het|homref|hetalt|homa
am = genotype mismatch
lm = allele/haplotype mismatch
. = non vu
**** On vérifie que am = genotype mismatch
référence  = T/T
high-confidence = T/C
notre = C/C
#+begin_src sh
bcftools filter -i 'POS=19196584'  /Work/Groups/bisonex/data/giab/GRCh38/HG001_GRCh38_1_22_v4.2.1_benchmark.vcf.gz | grep -v '#'
bcftools filter -i 'POS=19196584'  ../out/NA12878_NIST7035-dbsnp/variantCalling/haplotypecaller/NA12878_NIST.vcf.gz | grep -v '#'
#+end_src
#+RESULTS:
: NC_000022.11    19196584        .       T       C       50      PASS    platforms=5;platformnames=Illumina,PacBio,10X,Ion,Solid;datasets=5;datasetnames=HiSeqPE300x,CCS15kb_20kb,10XChromiumLR,IonExome,SolidSE75bp;callsets=7;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbDV,CCS15kb_20kbGATK4,HiSeqPE300xfreebayes,10XLRGATK,IonExomeTVC,SolidSE75GATKHC;datasetsmissingcall=CGnormal;callable=CS_HiSeqPE300xGATK_callable,CS_CCS15kb_20kbDV_callable,CS_10XLRGATK_callable,CS_CCS15kb_20kbGATK4_callable,CS_HiSeqPE300xfreebayes_callable GT:PS:DP:ADALL:AD:GQ    0/1:.:781:109,123:138,150:348
: NC_000022.11    19196584        rs1061325       T       C       59.32   PASS    AC=2;AF=1;AN=2;DB;DP=2;ExcessHet=0;FS=0;MLEAC=1;MLEAF=0.5;MQ=60;QD=29.66;SOR=2.303      GT:AD:DP:GQ:PL  1/1:0,2:2:6:71,6,0
**** On vérifie que lm = allele/haplotype mismatch
référence  = CAA/CAA
high-confidence = CA/CA
notre = C/CA
#+begin_src sh
 bcftools filter -i 'POS=31277416'  /Work/Groups/bisonex/data/giab/GRCh38/HG001_GRCh38_1_22_v4.2.1_benchmark.vcf.gz | grep -v '#'
 bcftools filter -i 'POS=31277416'  ../out/NA12878_NIST7035-dbsnp/variantCalling/haplotypecaller/NA12878_NIST.vcf.gz | grep -v '#'
#+end_src
#+RESULTS:
: NC_000022.11    31277416        .       CA      C       50      PASS    platforms=3;platformnames=Illumina,PacBio,10X;datasets=3;datasetnames=HiSeqPE300x,CCS15kb_20kb,10XChromiumLR;callsets=4;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbDV,10XLRGATK,HiSeqPE300xfreebayes;datasetsmissingcall=CCS15kb_20kb,CGnormal,IonExome,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable;difficultregion=GRCh38_AllHomopolymers_gt6bp_imperfectgt10bp_slop5,GRCh38_SimpleRepeat_imperfecthomopolgt10_slop5  GT:PS:DP:ADALL:AD:GQ    1/1:.:465:16,229:0,190:129
: NC_000022.11    31277416        rs57244615      CAA     C,CA    389.02  PASS    AC=1,1;AF=0.5,0.5;AN=2;BaseQRankSum=0.37;DB;DP=37;ExcessHet=0;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;MQRankSum=0;QD=13.41;ReadPosRankSum=-0.651;SOR=0.572    GT:AD:DP:GQ:PL  1/2:5,10,14:29:64:406,202,313,64,0,88
*** Génération de reads
Biblio récente
https://www.biorxiv.org/content/10.1101/2022.03.29.486262v1.full.pdf
Parmi ceux qui gèrent les variations
- *simuscop* reads non centré sur les zones de capture
- *NEAT: exome* mais trop lent en pratique
- *Reseq* exome
- gensim : pas d'exome
- pIRS : non plus
- varsim : non plus
  ...
  Temps de calcul selon l'article de reseq https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02265-7
  #+begin_quote
  Due to ReSeq’s effective parallelization, its elapsed times are low for this benchmark with 48 virtual CPUs (Additional file 1: Figure S34b,e). In contrast, the single-threaded processes implemented in perl or python have strikingly high elapsed times. This is well visible in Hs-HiX-TruSeq and applies to the training of pIRS (over a week), NEAT (several days), and BEAR (half a week) as well as the simulation of NEAT (close to 2 weeks) and BEAR (several weeks).
Biblio : https://www.nature.com/articles/s41437-022-00577-3
  #+end_quote
Divers
- Liste ancienne : https://www.biostars.org/p/128762/
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02265-7
* Idées
** Validation analytique
mail Yannis : données patients +/- simulées
*** Utiliser données GCAT et uploader le notre ?
https://www.nature.com/articles/ncomms7275
*** [#A] Variant calling : Genome in a bottle : NA12878 + autres
Résumé : https://www.nist.gov/programs-projects/genome-bottle
Manuscript : https://www.nature.com/articles/s41587-019-0054-x.epdf?author_access_token=E_1bL0MtBBwZr91xEsy6B9RgN0jAjWel9jnR3ZoTv0OLNnFBR7rUIZNDXq0DIKdg3w6KhBF8Rz2RWQFFc0St45kC6CZs3cDYc87HNHovbWSOubJHDa9CeJV-pN0BW_mQ0n7cM13KF2JRr_wAAn524w%3D%3D
Article comparant les variant calling : https://www.biorxiv.org/content/10.1101/2020.12.11.422022v1.full.pdf
**** KILL Tester le séquencage aussi
CLOSED: [2023-01-30 lun. 18:30]
Depuis un fastq correspondant à Illumina  https://github.com/genome-in-a-bottle/giab_data_indexes
   puis on compare le VCF avec les "high confidence"
On séquence directement NA12878 -> inutile pour le pipeline seul
**** TODO Tester seul la partie bioinformatique
   Tout résumé ici : https://www.nist.gov/programs-projects/genome-bottle
- methode https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/NA12878/analysis/Illumina_PlatinumGenomes_NA12877_NA12878_09162015/IlluminaPlatinumGenomes-user-guide.pdf
- vcf
     https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/NA12878_HG001/latest/GRCh38/
NB: à quoi correspond https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/NA12878/analysis/Illumina_PlatinumGenomes_NA12877_NA12878_09162015/hg38/2.0.1/NA12878/ ??
   Article comparant les variant calling : https://www.biorxiv.org/content/10.1101/2020.12.11.422022v1.full.pdf
   Article pour vcfeval : https://www.nature.com/articles/s41587-019-0054-x
   La version 4 ajoute 273 gènes "clinically relevant" https://www.biorxiv.org/content/10.1101/2021.06.07.444885v3.full.pdf
   Ajout des zones "difficiles"
   https://www.biorxiv.org/content/10.1101/2020.07.24.212712v5.full.pdf
*** [#B] Pipeline : générer patient avec tous les variants retrouvés à Centogene
Comparaison de génération ADN (2019)
https://academic.oup.com/bfg/article/19/1/49/5680294
**** SimuSCop (exome)
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-020-03665-5
https://github.com/qasimyu/simuscop
1. Crééer un modèle depuis bam + vcf : Setoprofile
2. Génerer données NGS
** Annotation :
*** Comparaison vep / snpeff et annovar
* Changement nouvelle version
- Dernière version du génome (la version "prête à l'emploi" est seulement GRCh38 sans les version patchées)
* Notes
** Nextflow
*** afficher les résultats d'un process/workflow
#+begin_src
lol.out.view()
#+end_src
Attention, ne fonctionne pas si plusieurs sortie:
#+begin_src
lol.out[0].view()
#+end_src
ou si /a/ est le nom de la sortie
#+begin_src
lol.out.a.view()
#+end_src
** Quelle version du génome ?
- T2T: notation chromose = chR1,2 : ok genome, clinvar, dbSNP
- GRCh38: notation chromose = NC_... : ok genome, clinvar, dbSNP
** Performances
Ordinateur de Carine (WSL2) : 4h dont 1h15 alignement (parallélisé) et 1h15 haplotypecaller (séquentiel)
** Chromosomes NC, NT, NW
Correspondance :
https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&chromInfoPage=
Signification
https://genome.ucsc.edu/FAQ/FAQdownloads.html#downloadAlt
- alt = séquences alternatives (utilisables)
- fix = patch (correction ou amélioration)
- random = séquence connue sur un chromosome mais non encore utilisée
** Pipelines prêt-à-l’emploi nextflow
Problème : nécessite singularity ou docker (ou conda)
Potentiellement utilisable avec nix...
** Validation : Quelles données de référence ?
Discussion avec Alexis
- Platinum genomes = génome seul
*** [[https://github.com/genome-in-a-bottle/giab_data_indexes][Genome in a bottle]]
  - NA12878 :
    - Illumina HiSeq Exome : fastq + capture
    - Illumina TruSeq Exome : bam, pas de capture
      ici ww
  - HG002,3,4
    - Illumina Whole Exome  : bam. le kit de capture est "Agilent SureSelect Human All Exon V5 kit" selon [[https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/README.txt][README]]. On il faut les régions [[https://kb.10xgenomics.com/hc/en-us/articles/115004150923-Where-can-I-find-the-Agilent-Target-BED-files-][selon ce site]]
      Un autre fichier est disponible (capture ???)
    https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/analysis/OsloUniversityHospital_Exome_GATK_jointVC_11242015/wex_Agilent_SureSelect_v05_b37.baits.slop50.merged.list
  "target region" +/- 50bp
    testé sur chr311780-312086 : ok
Autres technologies non adaptées au pipeline (vu avec Alexis)
*** [[https://www.illumina.com/platinumgenomes.html][Platinum genome
]] Que du génome « sequenced to 50x depth on a HiSeq 2000 system”
Genome possible
** Zone de capture
GIAB fourni le .bed pour l'exome . INfo : https://support.illumina.com/sequencing/sequencing_kits/nextera-rapid-capture-exome-kit/downloads.html
** Centogène
https://www.twistbioscience.com/node/23906
Bed non fourni pour exactement cette capture
On prend https://www.twistbioscience.com/resources/data-files/twist-alliance-vcgs-exome-401mb-bed-files
qui content la majeure partie
* Données :data:
** DONE Remplacer bam par fastq sur mesocentre
CLOSED: [2023-04-16 Sun 16:33]
Commande
*** DONE Supprimer les fastq non "paired"
CLOSED: [2023-04-16 Sun 16:33]
nushell
Liste des fastq avec "paired-end" manquant
#+begin_src nu
ls **/*.fastq.gz | get name | path basename | split column "_" | get column1 | uniq -u | save single.txt
#+end_src
#+RESULTS:
: 62907927
: 62907970
: 62899606
: 62911287
: 62913201
: 62914084
: 62915905
: 62921595
: 62923065
: 62925220
: 62926503
: 62926502
: 62926500
: 62926499
: 62926498
: 62931719
: 62943423
: 62943400
: 62948290
: 62949205
: 62949206
: 62949118
: 62951284
: 62960792
: 62960785
: 62960787
: 62960617
: 62962561
: 62962692
: 62967473
: 62972194
: 62979102
On vérifie
#+begin_src nu
open single.txt  | lines | each {|e| ls $"fastq/*_($in)/*" | get 0  }
open single.txt  | lines | each {|e| ls $"fastq/*_($in)/*" | get 0.name }  | path basename | split column "_" | get column1 | uniq -c
#+end_src
On met tous dans un dossier (pas de suppression )
#+begin_src
open single.txt  | lines | each {|e| ls $"fastq/*_($in)/*" | get 0  }  | each {|e| ^mv $e.name bad-fastq/}
#+end_src
On vérifie que les dossiier sont videsj
 open single.txt  | lines | each {|e| ls $"fastq/*_($in)" | get 0.name } | ^ls -l $in
 Puis on supprime
 open single.txt  | lines | each {|e| ls $"fastq/*_($in)" | get 0.name } | ^rm -r $in
*** DONE Supprimer bam qui ont des fastq
CLOSED: [2023-04-16 Sun 16:33]
On liste les identifiants des fastq et bam dans un tableau avec leur type :
#+begin_src
let fastq = (ls fastq/*/*.fastq.gz | get name | parse "{dir}/{full_id}/{id}_{R}_001.fastq.gz"  | select dir id | uniq )
let bam = (ls bam/*/*.bam | get name | parse "{dir}/{full_id}/{id}_{S}.bqrt.bam"  | select dir id)
#+end_src
On groupe les résultat par identifiant (résultats = liste de records qui doit être convertie en table)
et on trie ceux qui n'ont qu'un fastq ou un bam
#+begin_src
let single = ( $bam | append $fastq | group-by id | transpose id files | get files | where {|x| ($x | length) == 1})
#+end_src
On convertit en table et on récupère seulement les bam
#+begin_src
$single | reduce {|it, acc| $acc | append $it} | where dir == bam | get id | each {|e| ^ls $"bam/*_($e)/*.bam"}
#+end_src
#+RESULTS:
: bam/2100656174_62913201/62913201_S52.bqrt.bam
: bam/2100733271_62925220/62925220_S33.bqrt.bam
: bam/2100738763_62926502/62926502_S108.bqrt.bam
: bam/2100746726_62926498/62926498_S105.bqrt.bam
: bam/2100787936_62931955/62931955_S4.bqrt.bam
: bam/2200066374_62948290/62948290_S130.bqrt.bam
: bam/2200074722_62948298/62948298_S131.bqrt.bam
: bam/2200074990_62948306/62948306_S218.bqrt.bam
: bam/2200214581_62967331/62967331_S267.bqrt.bam
: bam/2200225399_62972187/62972187_S85.bqrt.bam
: bam/2200293962_62979117/62979117_S63.bqrt.bam
: bam/2200423985_62999352/62999352_S1.bqrt.bam
: bam/2200495073_63010427/63010427_S20.bqrt.bam
: bam/2200511274_63012586/63012586_S114.bqrt.bam
: bam/2200669188_63036688/63036688_S150.bqrt.bam
* Nouveau workflow :workflow:
** TODO Bases de données
*** KILL Nix pour télécharger les données brutes
**** Conclusion
Non viable sur cluster car en dehors de /nix/store
On peut utiliser des symlink mais trop compliqué
**** KILL Axel au lieu de curl pour gérer les timeout?
CLOSED: [2022-08-19 Fri 15:18]
*** DONE Tester patch de @pennae pour gros fichiers
SCHEDULED: <2022-08-19 Fri>
*** KILL Télécharger les données avec nextflow: hg38
CLOSED: [2023-06-12 Mon 23:29]
**** DONE Genome de référence
**** DONE dbSNP
**** DONE VEP 20G
CLOSED: [2023-06-12 Mon 23:29]
Ajout vérification checksum -> à vérifier
**** DONE transcriptome (spip)
CLOSED: [2023-06-12 Mon 23:29]
Rajouter checksum manuel
**** KILL Refseq
**** KILL OMIM
CLOSED: [2023-06-12 Mon 23:29]
codé, à vérifier
**** KILL ACMG incidental
CLOSED: [2023-06-12 Mon 23:29]
*** TODO Télécharge données :TIT:
SCHEDULED: <2023-06-15 Thu>
**** DONE fasta
CLOSED: [2023-06-12 Mon 23:30]
***** DONE compatibilité hg38
CLOSED: [2023-06-12 Mon 23:30]
**** DONE fasta index
CLOSED: [2023-06-13 Tue 00:07]
***** DONE compatibilité hg38
CLOSED: [2023-06-13 Tue 00:07]
**** DONE fasta dictionnaire
CLOSED: [2023-06-13 Tue 00:07]
**** DONE dbSNP
CLOSED: [2023-06-12 Mon 23:30]
***** DONE compatibilité hg38
CLOSED: [2023-06-12 Mon 23:30]
**** TODO commonSNP
***** TODO compatibilité hg38
cd /Work/Groups/bisonex/data/dbsnp/GRCh38.p14
❯ ga@mesointeractive GRCh38.p14]$ zgrep -c '^NC' dbSNP_common.vcf.gz
21340485
[apraga@mesointeractive GRCh38.p14]$ pwd
[apraga@mesointeractive GRCh38.p14]$ zgrep -c '^NC'
dbSNP_common.vcf.gz                     ID_of_common_snp_not_clinvar_patho.txt
dbSNP_common.vcf.gz.tbi                 ID_of_common_snp.txt
[apraga@mesointeractive dbsnp]$ cd chm13v2.0/
[apraga@mesointeractive chm13v2.0]$ ls
chm13v2.0_dbSNPv155.vcf.gz      dbSNP_common.vcf.gz.tbi                 versions.yml
chm13v2.0_dbSNPv155.vcf.gz.tbi  ID_of_common_snp_not_clinvar_patho.txt
dbSNP_common.vcf.gz             ID_of_common_snp.txt
[apraga@mesointeractive chm13v2.0]$ zgrep -c '^chr' dbSNP_common.vcf.gz
19433713
[apraga@mesointeractive chm13v2.0] $
❯ man tmux
**** TODO commonSNP non patho
***** TODO compatibilité hg38
*** HOLD Processing bases de données
**** DONE dbSNP common
**** DONE Seulement les ID dans dbSNP common !
CLOSED: [2022-11-19 Sat 21:42]
172G au lieu de 253M...
**** HOLD common dbSNP not clinvar patho
***** DONE Conclusion partielle
CLOSED: [2022-12-12 Mon 22:25]
- vcfeval : prometteur mais n'arrive pas à traiter toutes les régions
- isec : trop de problèmes avec
- classif clinvar directement dans dbSNP: le plus simple
  Et ça permet de rattraper quelques erreurs dans le script d'Alexis
***** KILL Utiliser directement le numéro dbSNP dans clinvar ? Non
CLOSED: [2022-11-20 Sun 19:51]
Ex: chr20
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f 'rs%INFO/RS \n' -i 'INFO/RS != "." & INFO/CLNSIG="Pathogenic"' clinvar_chr20.vcf.gz | sort > ID_clinvar_patho.txt
bcftools query -f '%ID\n' dbSNP_common_chr20.vcf.gz | sort > ID_of_common_snp.txt
comm -23 ID_of_common_snp.txt ID_clinvar_patho.txt > ID_of_common_snp_not_clinvar_patho.txt
wc -l ID_of_common_snp_not_clinvar_patho.txt
# sort ID
#+end_src
#+RESULTS:
: 518846 ID_of_common_snp_not_clinvar_patho.txt
Version d'alexis
#+begin_src sh :dir ~/code/bisonex/test_isec
snp=dbSNP_common_chr20.vcf.gz
clinvar=clinvar_chr20_notremapped.vcf.gz
python ../script/pythonScript/clinvar_sbSNP.py \
    --clinvar $clinvar \
    --chrm_name_table ../

Replacement in projects/bisonex.org at line 48 [5.35]

B:BD[6.89] → [6.89:786]

B:BD[6.786] → [2.17156:24651]

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@?CBB
CCBDCBDCCCFBD;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

[6.89]

[2.24651]

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