3WBY7ELOD4XA65WEV3UCLE6VUAGGF4T27WUYXN75QAOSYRPFEEMAC
WY5OHMYD4ZRUVYNHUQ63E34H6KFPR5RSYUN4EJR4OTREACIHROQAC
PAKCUZSF23LOLYHJJHFRJFXRYDNLGVWSQ3RZ3ZM5TCLFA2PCUIXAC
FXA3ZBV64FML7W47IPHTAJFJHN3J3XHVHFVNYED47XFSBIGMBKRQC
QWCRR5CXIUZYZADCPUIE35VYCH4IEIEPPCOP2RV6TETADPS6H6SAC
B2YQPU5UAXLACFYTD3FXVUP7AUC24GRHE6RR63CVFONRW6F7HMVQC
MPOLY5PJKJ74PWLI6PUOZ3KOIPLFLDAQO3WSTKRGA5MANYH5ZLUAC
DOTFSLPG6V4NKETA7D3MDB2MXSAA7CYXLRMYKAYURHGACX6JNP6AC
T3IPJM6TYF25RE2EQGASVSGIYJTJSJNZXBCEK3BNJ2LUYIBMTSSAC
MHIFI3P3R5PVLHHZRDZH3FPFN5IEDJRTWIVCNEUHNUKVSGYZW6YQC
XVXPOKRQXHHP3EPTGU3WCJWV7QXYDAKRCW6XWA2WH36ONKN2JA5AC
#+title: Bisonex
* Biblio
Comparaison WDL, Cromwell, nextflow
https://www.nature.com/articles/s41598-021-99288-8
Nextflow = bon compromis ?
* Changement nouvelle version
- Dernière version du génome (la version "prête à l'emploi" est seulement GRCh38 sans les version patchées)
* Notes
** 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...
* Données
** DONE Vérifier qualité données sur mesocentre
*** DONE BAM
picard ValidateSamFile
On regarde juste le code d'erreur (0 = pas d'erreur)
*** DONE Fastq
fastqc
Il faut ensuite extraire les zip and chercher les erreur dedans
** DONE Lister données sur mesocentre
<2022-12-23 ven.>
176 patients ok
Données manquantes :
- 5 avec juste 1 fastq sur les 2 sans bam
- 19 patients sans fastq ni bam
* Nouveau 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>
*** STRT Télécharger
- [X] Genome de référence
- [X] dbSNP
- [X] OMIM
- [X] VEP 20G
- [X] transcriptome (spip)
- [ ] Refseq
*** DONE Télécharger les données avec nextflow
CLOSED: [2022-09-13 Tue 21:37]
*** 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 ../database/RefSeq/refseq_to_number_only_consensual.txt \
--dbSNP $snp --output prod.txt
wc -l prod.txt
zgrep '^NC' dbSNP_common_chr20.vcf.gz | wc -l
#+end_src
#+RESULTS:
| 518832 | prod.txt |
| 518846 | |
***** KILL classification clinvar codée dbSNP ?
CLOSED: [2022-12-04 Sun 14:38]
Sur le chromosome 20
*Attention* CLNSIG a plusieurs champs (séparé par une virgule)
On y accède avec INFO/CLNSIG[*]
Ensuite, chaque item peut avoir plusieurs haploïdie (séparé par un |). IL faut donc utiliser une regexp
NB: *ne pas mettre la condition* dans une variable !!
Pour avoir les clinvar patho, on veut 5 mais pas 255 (= autre) pour la classification !`
Il faut également les likely patho et conflicting
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%INFO/CLNSIG\n' dbSNP_common_chr20.vcf.gz -i \
'INFO/CLNSIG[*]~"^5|" | INFO/CLNSIG[*]=="5" | INFO/CLNSIG[*]~"|5" | INFO/CLNSIG[*]~"^4|" | INFO/CLNSIG[*]=="4" | INFO/CLNSIG[*]~"|4" | INFO/CLNSIG[*]~"^12|" | INFO/CLNSIG[*]=="12" | INFO/CLNSIG[*]~"|12"' | sort
#+end_src
#+RESULTS:
| . | . | 12 | | | | | | | | |
| . | 12 | 0 | 2 | | | | | | | |
| 2 | 3 | 2 | 2 | 2 | 5 | . | | | | |
| . | 2 | 3 | 2 | 2 | 4 | | | | | |
| . | . | 3 | 12 | 3 | | | | | | |
| . | 5 | 2 | . | | | | | | | |
| . | . | . | 5 | 2 | 2 | | | | | |
| . | 9 | 9 | 9 | 5 | 5 | 2 | 3 | 2 | 3 | 2 |
Si on les exclut :
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%ID\n' dbSNP_common_chr20.vcf.gz -e \
'INFO/CLNSIG[*]~"^5|" | INFO/CLNSIG[*]=="5" | INFO/CLNSIG[*]~"|5" | INFO/CLNSIG[*]~"4" | INFO/CLNSIG[*]~"12"' | sort | uniq > common-notpatho.txt
#+end_src
#+RESULTS:
#+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 ../database/RefSeq/refseq_to_number_only_consensual.txt \
--dbSNP $snp --output tmp.txt
sort tmp.txt | uniq > common-notpatho-alexis.txt
wc -l common-notpatho-alexis.txt
#+end_src
#+RESULTS:
: 518832 common-notpatho-alexis.txt
On en a 6 de plus que la version d'Alexis mais quelques différences
Ceux d'Alexis qui manquent:
#+begin_src sh :dir ~/code/bisonex/test_isec
comm -23 common-notpatho-alexis.txt common-notpatho.txt > alexis-only.txt
cat alexis-only.txt
#+end_src
#+RESULTS:
| rs1064039 |
| rs3833341 |
| rs73598374 |
On les teste dans clinvar et dbSNP
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'ID=@alexis-only.txt' dbSNP_common_chr20.vcf.gz
bcftools query -f '%POS\n' -i 'ID=@alexis-only.txt' dbSNP_common_chr20.vcf.gz > alexis-only-pos.txt
while read -r line; do
bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'POS='$line clinvar_chr20.vcf.gz
done < alexis-only-pos.txt
# bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'POS=23637790' clinvar_chr20.vcf.gz
#+end_src
#+RESULTS:
| 764018 | A | ACAGGTCAAT,ACAGGT | .,5 | 2,. | |
| 23637790 | C | G,T | .,.,12 | | |
| 44651586 | C | A,G,T | .,.,.,5 | 2 | 2 |
| 764018 | A | ACAGGTCAAT | Benign | | |
| 23637790 | C | T | Benign | | |
| 44651586 | C | T | Benign | | |
On a donc une discordance entre clinvar et dbSNP.
On dirait qu'ils ont mal fait l'intersection avec clinvar.
Par exemple https://www.ncbi.nlm.nih.gov/snp/rs3833341#clinical_significance
Tu as l'impression qu'il y a un 1 clinvar bénin et 1 patho.
En cherchant par NM, tu vois qu'il est bénin sur clinvar car il y a d'autres soumissions ! https://www.ncbi.nlm.nih.gov/clinvar/variation/262235/
Confirmation sur nos bases de données :
$ bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'POS=764018' dbSNP_common_chr20.vcf.gz
764018 A ACAGGTCAAT,ACAGGT .,5|2,.
$ bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'POS=764018' clinvar_chr20.vcf.gz
764018 A ACAGGTCAAT Benign
***** KILL Corriger script alexi
CLOSED: [2022-12-04 Sun 13:03]
Gère clinvar patho, probablement patho ou conflicting !
***** HOLD Rtg tools
****** Test
1. Générer SDf file
#+begin_src sh
rtg format genomeRef.fna -o genomeRef.sdf
#+end_src
2. Pour les bases de donnés, il faut l'option --sample ALT sinon on a
#+begin_src
$ rtg vcfeval -b dbSNP_common.vcf.gz -c clinvar.vcf.gz -o test -t genomeRef.sdf/^C
VCF header does not contain a FORMAT field named GQ
Error: Record did not contain enough samples: NC_000001.11 10001 rs1570391677 A,C . PASS RS=1570391677;dbSNPBuildID=154;SSR=0;PSEUDOGENEINFO=DDX11L1:100287102;VC=SNV;R5;GNO;FREQ=KOREAN:0.9891,0.0109,.|SGDP_PRJ:0,1,.|dbGaP_PopFreq:1,.,0;COMMON
#+end_src
Essai intersection clinvar (patho ou non) dbSNP
- faux négatif = dbSNP common qui ne sont pas dans clinvar
- faux positif = clinvar qui ne sont pas dbSNP common
- vrai positif = clinvar qui sont dans dbSNP common
- vrai positif baseline = dbSNP common qui sont dans clinvar
On calcule le nombre de lignes
#+begin_src ssh
zgrep '^[^#]' /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz | wc -l
for i in *.vcf.gz; do echo $i; zgrep '^[^#]' $i | wc -l; done
#+end_src
| clinvar | 1493470 |
| fn.vcf.gz | 22330220 |
| fp.vcf.gz | 1222529 |
| tp-baseline.vcf.gz | 131040 |
| tp.vcf.gz | 136638 |
À noter qu'on ne retrouve pas tout clinvar...
1222529 + 131040 = 1353569 < 1493470
certains régions ne sont pas traitées :
#+begin_quote
Evaluation too complex (50002 unresolved paths, 34891 iterations) at reference region NC_000001.11:790930-790970. Variants in this region will not be included in results
#+end_quote
#+begin_src sh
grep 'not be included' vcfeval.log | wc -l
56192
#+end_src
Le total est quand même inférieur
On veut les clinvar non patho dans dbSNP soit les faux négatif (dbSNP common not contenu dans clinvar patho)
#+begin_src sh
bcftools filter -i 'INFO/CLNSIG="Pathogenic"' /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz -o /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar-patho.vcf.gz
tabix /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar-patho.vcf.gz
#+end_src
On lance le script (dbSNP common et clinvar = 9h)
#+begin_src sh
#!/bin/bash
#SBATCH --nodes=1
#SBATCH -p smp
#SBATCH --time=12:00:00
#SBATCH --mem=12G
dir=/Work/Groups/bisonex/data
dbSNP=$dir/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz
clinvar=$dir/clinvar/GRCh38/clinvar-patho.vcf.gz
genome=$dir/genome/GRCh38.p13/genomeRef.sdf
srun rtg vcfeval -b $dbSNP -c $clinvar -o common-not-patho -t $genome --sample ALT
#+end_src
****** HOLD Voir pour régions complexes non traitées
***** DONE bcftools isec : non
CLOSED: [2022-11-27 Sun 00:38]
#+begin_src sh
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -p common
#+end_src
On vérifie bien que les 2 fichiers commons on le même nombre de lignes
#+begin_src sh
$ grep -e '^NC' 0002.vcf | wc -l
74302
alex@gentoo ~/code/bisonex/data/common $ grep -e '^NC' 0003.vcf | wc -l
74302
#+end_src
****** DONE Impact option -n
CLOSED: [2022-10-23 Sun 13:56]
Mais en spécifiant -n =2:
#+begin_src sh
$ bedtools intersect -a dbSNP_common.vcf.gz -b clinvar.vcf.gz
74978
#+end_src
Si on ne regarde que les variants, on retrouve bien 74302
#+begin_src sh
rg "^NC" none_sorted.vcf | wc -l
#+end_src
NB : test fait avec
#+begin_src
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c none -n =2 -w 1 | sort > none.vcf
sort common/0003.vcf > common/0003_sorted.vcf
comm -13 common/0003_sorted.vcf none_sorted.vcf
#+end_src
****** DONE Géstion des duplicates: -c none
CLOSED: [2022-10-23 Sun 13:56]
Si on ne garde que ceux avec REF et ALT identiques
#+begin_src sh
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c none -n =2 -w 1 | wc -l
74978
#+end_src
Si on garde tout
#+begin_src sh
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c all -n =2 -w 1 | wc -l
137777
#+end_src
Pour regarder la différence :
#+begin_src sh
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c none -n =2 -w 1 | sort > none_sorted.vcf
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c all -n =2 -w 1 | sort > all_sorted.vcf
comm -13 none_sorted.vcf all_sorted.vcf | head
#+end_src
Sur un exemple,on a bien des variants différents
****** DONE Suppression des clinvar patho
CLOSED: [2022-10-23 Sun 18:55]
Semble faire le travail vu que dbSNP_commo a 23194960 lignes (donc ~80 000 de moins)
#+begin_src sh
$ bcftools isec -e 'INFO/CLNSIG="Pathogenic" & INFO/CLNSIG="Pathogenic/Likely_pathogenic"' -c none -n~10 dbSNP_common.vcf.gz clinvar.vcf.gz | wc -l
Note: -w option not given, printing list of sites...
23119984
#+end_src
Par contre, l'o'ption -w ou -p fait des ficher "data"...
Après un nouvel essai, plus de problème
#+begin_src
$ bcftools isec -e 'INFO/CLNSIG="Pathogenic" & INFO/CLNSIG="Pathogenic/Likely_pathogenic"' -c none -n=1 dbSNP_common.vcf.gz clinvar.vcf.gz -w 1 -o lol.vcf.gz
$ zcat lol.vcf.gz | wc -l
23120660
#+end_src
À noter le choix de l'option -n qui change entre "=1" et "~10"...
En effet "=1" = au moins 1 fichier et "~10" fait exactement dans le premier et non dans le second
#+begin_src
$ bcftools isec -e 'INFO/CLNSIG="Pathogenic" & INFO/CLNSIG="Pathogenic/Likely_pathogenic"' -c none -n~10 dbSNP_common.vcf.gz clinvar.vcf.gz -w 1 -o lol.vcf.gz
$ zcat lol.vcf.gz | wc -l
23120660
#+end_src
****** DONE Valider avec Alexis : bcftool isec
CLOSED: [2022-11-07 Mon 21:42 ]
****** DONE Pourquoi nombre de lignes différentes avec la version d'Alexis -> isec ne gère pas plusieurs ALT
CLOSED: [2022-11-26 Sat 23:36]
Grosse différence !
#+begin_src
$ wc -l ID_of_common_snp_not_clinvar_patho.txt
23119915 ID_of_common_snp_not_clinvar_patho.txt
$ wc -l /Work/Users/apraga/bisonex/database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
85820 /Work/Users/apraga/bisonex/database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
#+end_src
À noter que tout dbSNP = 23194960
******* Clinvar classe 4 ? Moins mais toujours trop
#+begin_src
$ zgrep '^NC' tmp.vcf.gz | wc -l
21081654
#+end_src
******* Comparer les ID et regarder ceux en plus
#+begin_src sh
bcftools isec -e 'INFO/CLNSIG="Pathogenic"' -c none -n~10 /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz -w 1 -o tmp.vcf.gz
zgrep -o -e 'rs[[:digit:]]\' tmp.vcf.gz | sort | id_sorted.txt
sort ../database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt > reference_sorted.txt
comm -23 id_sorted.txt reference_sorted.txt > unique1.txt
#+end_src
Par exemple
#+begin_src sh
zgrep rs1000000561 ../database/dbSNP/dbSNP_common.vcf.gz
#+end_src
NC_000002.12 136732859 rs1000000561 ACG A,ACGCG . PASS RS=1000000561;dbSNPBuildID=151;SSR=0;VC=INDEL;GNO;FREQ=ALSPAC:0.2506,0.7494,.|TOMMO:0.9971,0.002865,.|TWINSUK:0.2473,0.7527,.|dbGaP_PopFreq:0.993,0.006943,8.902e-05;COMMON
Attention, clinvar est en numéro de chromosomoe et dbSNP en NC...
Normalement, géré lors du calcul d'intersection !
Ce SNP n'est pas dans clinvar (vérifié dans UCSC)
******* Tester sur chromosome 20
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools view --regions NC_000020.11 ../database/dbSNP/dbSNP_common.vcf.gz -o dbSNP_common_chr20.vcf.gz
bcftools view --regions 20 ../database/clinvar/clinvar.vcf.gz -o clinvar_chr20.vcf.gz
tabix dbSNP_common_chr20.vcf.gz
tabix clinvar_chr20.vcf.gz
#+end_src
#+RESULTS:
Attention à bien renommer clinvar !
#+begin_src sh :dir ~/code/bisonex/test_isec
mv clinvar_chr20.vcf.gz clinvar_chr20_notremapped.vcf.gz
bcftools annotate --rename-chrs chromosome_mapping.txt clinvar_chr20_notremapped.vcf.gz -o clinvar_chr20.vcf.gz
#+end_src
#+RESULTS:
*ATTENTION*: sans indexer les vcf, les fichiers seront *VIDES*
*ATTENTION*: par défaut les filtres s'appliquent sur les 2. Cela est un problème si on joue sur l'inclusion et non l'exclusion
Attention: vérifier la conventdion de nommage des chromosomes
******** Test pathogene: ne prend pas en compte les multi-allèles ????
On teste l'intersection dbsnp et clinvar patho ainsi que le complémentaire
#+begin_src sh :dir ~/code/bisonex/test_isec
clinvar=clinvar_chr20_patho.vcf.gz
snp=dbSNP_common_chr20.vcf.gz
bcftools index $clinvar
bcftools index $snp
bcftools filter -i 'INFO/CLNSIG="Pathogenic"' clinvar_chr20.vcf.gz -o $clinvar
bcftools isec $snp $clinvar -p tmp
for i in tmp/*.vcf ; do echo $i; grep '^[^#]' $i | wc -l; done
#+end_src
#+RESULTS:
| tmp/0000.vcf |
| 518846 |
| tmp/0001.vcf |
| 0 |
| tmp/0002.vcf |
| 0 |
| tmp/0003.vcf |
| 0 |
Aucun clinvar patho... Clairement faux !
Autre méthode : on inclut tous les SNP et clinvar patho et on regarde ceux uniquement dans dbsnp
#+begin_src sh :dir ~/code/bisonex/test_isec
snp=dbSNP_common_chr20.vcf.gz
clinvar=clinvar_chr20.vcf.gz
bcftools isec -n=2 -i - -i 'INFO/CLNSIG="Pathogenic"' $snp $clinvar -p tmp
# grep '^[^#]' tmp/0000.vcf | wc -l
#+end_src
#+RESULTS:
Soit tout dbsnp donc rien
Note : on ne peut pas exclure les clinvar patho directement
#+begin_src sh :dir ~/code/bisonex/test_isec
snp=dbSNP_common_chr20.vcf.gz
clinvar=clinvar_chr20.vcf.gz
bcftools isec -i - -e 'INFO/CLNSIG="Pathogenic"' $snp $clinvar -p tmp
for i in tmp/*.vcf ; do echo $i; grep '^[^#]' $i | wc -l; done
#+end_src
Car on ne peut plus faire la différence !
Si on utilise la 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 ../database/RefSeq/refseq_to_number_only_consensual.txt \
--dbSNP $snp --output tmp.txt
sort tmp.txt > common-notpatho-alexis.txt
wc -l common-notpatho-alexis.txt
#+end_src
#+RESULTS:
: 518832 common-notpatho-alexis.txt
Si on cherche les clinvar patho (donc non présent dans la sortie)
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%ID\n' dbSNP_common_chr20.vcf.gz | sort > all.txt
sort common-notpatho-alexis.txt > alexis.txt
comm -23 all.txt alexis.txt > patho.txt
#+end_src
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%POS\n' -i 'ID=@patho.txt' dbSNP_common_chr20.vcf.gz -o pos.txt
for pos in $(cat pos.txt); do
bcftools query -f '%CHROM %POS %ID %REF %ALT\n' -i 'POS='$pos dbSNP_common_chr20.vcf.gz
bcftools query -f '%CHROM %POS %ID %REF %ALT %INFO/CLNSIG\n' -i 'POS='$pos clinvar_chr20.vcf.gz
echo "------"
done
#+end_src
#+RESULTS:
| NC_000020.11 | 3234173 | rs3827075 | T | A,C,G | |
| NC_000020.11 | 3234173 | 262001 | T | G | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 3234173 | 1072511 | T | TGGCGAAGC | Pathogenic |
| NC_000020.11 | 3234173 | 208613 | TGGCGAAGC | G | Pathogenic |
| NC_000020.11 | 3234173 | 1312 | TGGCGAAGC | T | Pathogenic |
| ------ | | | | | |
| NC_000020.11 | 4699605 | rs1799990 | A | G | |
| NC_000020.11 | 4699605 | 13397 | A | G | Benign/Likely_benign |
| ------ | | | | | |
| NC_000020.11 | 10652589 | rs1131695 | G | A,C,T | |
| NC_000020.11 | 10652589 | 163705 | G | . | Benign |
| NC_000020.11 | 10652589 | 143063 | G | A | Benign |
| NC_000020.11 | 10652589 | 234555 | G | C | Pathogenic |
| ------ | | | | | |
| NC_000020.11 | 10658574 | rs1801138 | G | A,T | |
| NC_000020.11 | 10658574 | 42481 | G | A | Benign |
| NC_000020.11 | 10658574 | 992651 | G | T | Likely_pathogenic |
| NC_000020.11 | 10658574 | 213550 | GC | A | Pathogenic |
| ------ | | | | | |
| NC_000020.11 | 10672794 | rs79338570 | G | A,C | |
| NC_000020.11 | 10672794 | 255557 | G | A | Benign/Likely_benign |
| NC_000020.11 | 10672794 | 594067 | G | C | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 10672794 | 1324603 | G | GGA | Likely_pathogenic |
| ------ | | | | | |
| NC_000020.11 | 18525868 | rs146917730 | C | T | |
| NC_000020.11 | 18525868 | 811603 | C | T | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 25390747 | rs373200654 | G | C | |
| NC_000020.11 | 25390747 | 338000 | G | C | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 32800145 | rs2424926 | C | G,T | |
| NC_000020.11 | 32800145 | 338173 | C | G | Benign |
| NC_000020.11 | 32800145 | 338174 | C | T | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 33412656 | rs35938843 | C | G,T | |
| NC_000020.11 | 33412656 | 220958 | C | T | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 45891622 | rs181943893 | G | A,C,T | |
| NC_000020.11 | 45891622 | 459632 | G | C | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 45891622 | 797035 | G | T | Likely_benign |
| NC_000020.11 | 45891622 | 1572689 | GCTA | G | Likely_benign |
| ------ | | | | | |
| NC_000020.11 | 54171651 | rs35873579 | G | A,T | |
| NC_000020.11 | 54171651 | 285894 | G | A | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 54171651 | 1373583 | G | C | Uncertain_significance |
| NC_000020.11 | 54171651 | 895614 | G | T | Benign/Likely_benign |
| ------ | | | | | |
| NC_000020.11 | 62172726 | rs36106901 | G | A | |
| NC_000020.11 | 62172726 | 981031 | G | A | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 63349782 | rs1044396 | G | A,C | |
| NC_000020.11 | 63349782 | 93427 | G | A | Benign |
| NC_000020.11 | 63349782 | 857384 | G | C | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 63414925 | rs1801545 | G | A,C,T | |
| NC_000020.11 | 63414925 | 194284 | G | A | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 63414925 | 129337 | G | C | Benign |
| NC_000020.11 | 63414925 | 851545 | GG | CA | Uncertain_significance |
| ------ | | | | | |
On a donc plusieurs problèmes :
1. isec devrait fonctionner au moins sur
| NC_000020.11 | 25390747 | rs373200654 | G | C | |
| NC_000020.11 | 25390747 | 338000 | G | C | Conflicting_interpretations_of_pathogenicity |
On teste juste sur cette ligne
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools filter -i 'POS=25390747' clinvar_chr20.vcf.gz -o clinvar_test.vcf.gz
bcftools filter -i 'POS=25390747' dbSNP_common_chr20.vcf.gz -o dbSNP_test.vcf.gz
#+end_src
On retrouve bien la ligne dans l'intersection...
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools filter -i 'POS=25390747' clinvar_chr20.vcf.gz -o clinvar_test.vcf.gz
bcftools index dbSNP_test.vcf.gz dbSNP_test.vcf.gz
bcftools index dbSNP_test.vcf.gz clinvar_test.vcf.gz
bcftools isec dbSNP_test.vcf.gz clinvar_test.vcf.gz -p test
#+end_src
#+RESULTS:
2. isec ne semble pas fonctionner sur en cas d'ALT multiples
| NC_000020.11 | 32800145 | rs2424926 | C | G,T | |
| NC_000020.11 | 32800145 | 338173 | C | G | Benign |
| NC_000020.11 | 32800145 | 338174 | C | T | Conflicting_interpretations_of_pathogenicity |
| | | | | | |
3. s'il y a plusieurs variantions à une position, il faut bien vérifier que tous ne sont pas patho.
La version d'Alexis le fait bien
| NC_000020.11 | 3234173 | rs3827075 | T | A,C,G | |
| NC_000020.11 | 3234173 | 262001 | T | G | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 3234173 | 1072511 | T | TGGCGAAGC | Pathogenic |
| NC_000020.11 | 3234173 | 208613 | TGGCGAAGC | G | Pathogenic |
| NC_000020.11 | 3234173 | 1312 | TGGCGAAGC | T | Pathogenic |
****** DONE Voir si isec gère les multiallélique (chr20) : non, impossible de faire marcher
CLOSED: [2022-11-27 Sun 00:37]
******* DONE chr20 en prenant un patho clinvar aussi dans dbSNP
CLOSED: [2022-11-27 Sun 00:37]
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools filter dbSNP_common_chr20.vcf.gz -i 'POS=10652589' -o test_dbsnp.vcf.gz
bcftools filter clinvar_chr20.vcf.gz -i 'POS=10652589' -o test_clinvar.vcf.gz
bcftools index test_dbsnp.vcf.gz
bcftools index test_clinvar.vcf.gz
#+end_src
#+RESULTS:
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools isec test_dbsnp.vcf.gz test_clinvar.vcf.gz -p tmp
grep '^[^#]' tmp/0002.vcf
grep '^[^#]' tmp/0003.vcf
#+end_src
#+RESULTS:
Même en biallélique, ne fonctionne pas.
Testé en modifiant test_dbsnp !
Fonctionne avec un variant par ligne
****** DONE isec en coupant les sites multialléliques: non
CLOSED: [2022-11-27 Sun 00:37]
******* DONE Exemple simple ok
CLOSED: [2022-11-27 Sun 00:34]
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools filter -i 'POS=10652589' dbSNP_common_chr20.vcf.gz -o dbsnp_mwi.vcf.gz
bcftools filter -i 'POS=10652589' clinvar_chr20.vcf.gz -o clinvar_mwi.vcf.gz
bcftools index -f dbsnp_mwi.vcf.gz
bcftools index -f clinvar_mwi.vcf.gz
bcftools isec dbsnp_mwi.vcf.gz clinvar_mwi.vcf.gz -n=2
#+end_src
#+RESULTS:
Même en biallélique, ne fonctionne pas.
Chr 20
Avec les fichiers du teste précédent
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools norm -m -any dbsnp_mwi.vcf.gz -o dbsnp_mwi_norm.vcf.gz
bcftools index dbsnp_mwi_norm.vcf.gz
bcftools isec dbsnp_mwi_norm.vcf.gz clinvar_mwi.vcf.gz -n=2
#+end_src
#+RESULTS:
| NC_000020.11 | 10652589 | G | A | 11 |
| NC_000020.11 | 10652589 | G | C | 11 |
******* TODO Sur dbSNP chr20 non
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools norm -m -any dbSNP_common_chr20 -o dbSNP_common_chr20_norm.vcf.gz
#+end_src
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools isec -i 'INFO/CLNSIG="Pathogenic"' dbSNP_common_chr20_norm.vcf.gz clinvar_chr20.vcf.gz -p tmp
#+end_src
#+RESULTS:
***** DONE Essai bedtools intersect
#+begin_src sh
bedtools intersect -a dbSNP_common.vcf.gz -b clinvar.vcf.gz
#+end_src
$ wc -l intersect.vcf
220206 intersect.vcf
** TODO Dépendences avec Nix
*** DONE GATK
CLOSED: [2022-10-21 Fri 21:59]
*** WAIT BioDBHTS
Contribuer pull request
*** DONE BioExtAlign
CLOSED: [2022-10-22 Sat 00:38]
*** WAIT BioBigFile
Revoir si on peut utliser kent dernière version
Contribuer pull request
*** HOLD rtg-tools
Convertir clinvar NC
*** DONE Spip
CLOSED: [2022-12-04 Sun 12:49]
Pas de pull request
*** DONE R + packages
CLOSED: [2022-11-19 Sat 21:05]
** DONE Exécution
CLOSED: [2022-09-13 Tue 21:37]
*** KILL test Bionix
*** KILL Implémenter execution avec Nix ?
Voir https://academic.oup.com/gigascience/article/9/11/giaa121/5987272?login=false
pour un exemple.
Probablement plus simple d’utiliser Nix pour gestion de l’environnement et snakemake pour l’exécution
Pas d’accès internet depuis le cluster
*** DONE nextflow
CLOSED: [2022-09-13 Tue 21:37]
** DONE Preprocessing avec nextflow
CLOSED: [2022-10-09 Sun 22:30]
*** DONE Map to reference
CLOSED: [2022-10-09 Sun 22:30]
*** DONE Mark duplicate
CLOSED: [2022-10-09 Sun 22:30]
*** DONE Recalibrate base quality score
CLOSED: [2022-10-09 Sun 22:30]
** DONE Variant calling avec Nextflow
CLOSED: [2022-11-19 Sat 21:34]
*** DONE Haplotype caller
CLOSED: [2022-10-09 Sun 22:40]
*** DONE Filter variants
CLOSED: [2022-10-09 Sun 22:40]
*** DONE Filter common snp not clinvar path
CLOSED: [2022-11-07 Mon 23:00]
Voir [[*common dbSNP not clinvar patho][common dbSNP not clinvar patho]]
*** DONE Filter variant only in consensual sequence
CLOSED: [2022-11-08 Tue 22:23]
*** DONE Filter technical variants
CLOSED: [2022-11-19 Sat 21:34]
** TODO Annotation avec nextflow
*** TODO VEP
*** TODO Spip
*** TODO Filtrer après VEP
On doit pouvoir se passer d'un script R avec bcftools
** STRT Tester version d'alexis avec Nix
*** DONE Ajouter clinvar
CLOSED: [2022-11-13 Sun 19:37]
*** DONE Alignement
CLOSED: [2022-11-13 Sun 12:52]
*** DONE Haplotype caller
CLOSED: [2022-11-13 Sun 13:00]
*** TODO Filter
- [X] depth
- [ ] comon snp not path
Problème avec liste des ID
**** TODO variant annotation
Besoin de vep
*** TODO Variant calling
* TODO Tests
** TODO Test de non régression avec version ALexis avec nix
*** DONE ID common snp
CLOSED: [2022-11-19 Sat 21:36]
#+begin_src
$ wc -l ID_of_common_snp.txt
23194290 ID_of_common_snp.txt
$ wc -l /Work/Users/apraga/bisonex/database/dbSNP/ID_of_common_snp.txt
23194290 /Work/Users/apraga/bisonex/database/dbSNP/ID_of_common_snp.txt
#+end_src
*** DONE ID common snp not clinvar patho
CLOSED: [2022-12-11 Sun 20:11]
**** DONE Vérification du problème
CLOSED: [2022-12-11 Sun 16:30]
Sur le J:
21155134 /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt.ref
Version de "non-régression"
21155076 database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
Nouvelle version
23193391 /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt
Si on enlève les doublons
$ sort database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt | uniq > old.txt
$ wc -l old.txt
21107097 old.txt
$ sort /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt | uniq > new.txt
$ wc -l new.txt
21174578 new.txt
$ sort /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt.ref | uniq > ref.txt
$ wc -l ref.txt
21107155 ref.txt
Si on regarde la différence
comm -23 ref.txt old.txt
rs1052692
rs1057518973
rs1057518973
rs11074121
rs112848754
rs12573787
rs145033890
rs147889095
rs1553904159
rs1560294695
rs1560296615
rs1560310926
rs1560325547
rs1560342418
rs1560356225
rs1578287542
...
On cherche le premier
bcftools query -i 'ID="rs1052692"' database/dbSNP/dbSNP_common.vcf.gz -f '%CHROM %POS %REF %ALT\n'
NC_000019.10 1619351 C A,T
Il est bien patho...
$ bcftools query -i 'POS=1619351' database/clinvar/clinvar.vcf.gz -f '%CHROM %POS %REF %ALT %INFO/CLNSIG\n'
19 1619351 C T Conflicting_interpretations_of_pathogenicity
On vérifie pour tous les autres
$ comm -23 ref.txt old.txt > tocheck.txt
On génère les régions à vérifier (chromosome number:position)
$ bcftools query -i 'ID=@tocheck.txt' database/dbSNP/dbSNP_common.vcf.gz -f '%CHROM\t%POS\n' > tocheck.pos
On génère le mapping inverse (chromosome number -> NC)
$ awk ' { t = $1; $1 = $2; $2 = t; print; } ' database/RefSeq/refseq_to_number_only_consensual.txt > mapping.txt
On remap clinvar
$ bcftools annotate --rename-chrs mapping.txt database/clinvar/clinvar.vcf.gz -o clinvar_remapped.vcf.gz
$ tabix clinvar_remapped.vcf.gz
Enfin, on cherche dans clinvar la classification
$ bcftools query -R tocheck.pos clinvar_remapped.vcf.gz -f '%CHROM %POS %INFO/CLNSIG\n'
$ bcftools query -R tocheck.pos database/dbSNP/dbSNP_common.vcf.gz -f '%CHROM %POS %ID \n' | grep '^NC'
#+RESULTS:
**** DONE Comprendre pourquoi la nouvelle version donne un résultat différent
CLOSED: [2022-12-11 Sun 20:11]
***** DONE Même version dbsnp et clinvar ?
CLOSED: [2022-12-10 Sat 23:02]
Clinvar différent !
$ bcftools stats clinvar.gz
clinvar (Alexis)
SN 0 number of samples: 0
SN 0 number of records: 1492828
SN 0 number of no-ALTs: 965
SN 0 number of SNPs: 1338007
SN 0 number of MNPs: 5562
SN 0 number of indels: 144580
SN 0 number of others: 3714
SN 0 number of multiallelic sites: 0
SN 0 number of multiallelic SNP sites: 0
clinvar (new)
SN 0 number of samples: 0
SN 0 number of records: 1493470
SN 0 number of no-ALTs: 965
SN 0 number of SNPs: 1338561
SN 0 number of MNPs: 5565
SN 0 number of indels: 144663
SN 0 number of others: 3716
SN 0 number of multiallelic sites: 0
SN 0 number of multiallelic SNP sites: 0
***** DONE Mettre à jour clinvar et dbnSNP pour travailler sur les mêm bases
CLOSED: [2022-12-11 Sun 12:10]
Problème persiste
***** DONE Supprimer la conversion en int du chromosome
CLOSED: [2022-12-10 Sat 19:29]
***** KILL Même NC ?
CLOSED: [2022-12-10 Sat 19:29]
$ zgrep "contig=<ID=NC_\(.*\)" clinvar/GRCh38/clinvar.vcf.gz > contig.clinvar
$ diff contig.txt contig.clinvar
< ##contig=<ID=NC_012920.1>
***** DONE Tester sur chromosome 19: ok
CLOSED: [2022-12-11 Sun 13:53]
On prépare les données
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
bcftools filter -i 'CHROM="NC_000019.10"' /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz -o dbSNP_common_19.vcf.gz
bcftools filter -i 'CHROM="NC_000019.10"' /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz -o clinvar_19.vcf.gz
bcftools filter -i 'CHROM="NC_000019.10"' /Work/Groups/bisonex/data-alexis/dbSNP/dbSNP_common.vcf.gz -o dbSNP_common_19_old.vcf.gz
bcftools filter -i 'CHROM="19"' /Work/Groups/bisonex/data-alexis/clinvar/clinvar.vcf.gz -o clinvar_19_old.vcf.gz
#+end_src
On récupère les 2 versions du script
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
git checkout regression ../../script/pythonScript/clinvar_sbSNP.py
cp ../../script/pythonScript/clinvar_sbSNP.py clinvar_sbSNP_old.py
git checkout HEAD ../../script/pythonScript/clinvar_sbSNP.py
#+end_src
#+RESULTS:
On compare
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
python ../../script/pythonScript/clinvar_sbSNP.py clinvar_sbSNP.py --clinvar clinvar_19.vcf.gz --dbSNP dbSNP_common_19.vcf.gz --output tmp.txt
sort tmp.txt | uniq > new.txt
table=/Work/Groups/bisonex/data-alexis/RefSeq/refseq_to_number_only_consensual.txt
python clinvar_sbSNP_old.py --clinvar clinvar_19_old.vcf.gz --dbSNP dbSNP_common_19_old.vcf.gz --output tmp_old.txt --chrm_name_table $table
sort tmp_old.txt | uniq > old.txt
wc -l old.txt new.txt
#+end_src
#+RESULTS:
| 535155 | old.txt |
| 535194 | new.txt |
| 1070349 | total |
Si on prend le premier manquant dans new, il est conflicting patho donc il ne devrait pas y être...
$ bcftools query -i 'ID="rs10418277"' dbSNP
_common_19.vcf.gz -f '%CHROM %POS %REF %ALT\n'
NC_000019.10 54939682 C G,T
$ bcftools query -i 'ID="rs10418277"' dbSNP_common_19_old.vcf.gz -f '%CHROM %POS %REF %ALT\n'
NC_000019.10 54939682 C G,T
$ bcftools query -i 'POS=54939682' clinvar_19.vcf.gz -f '%POS %REF %ALT %INFO/CLNSIG\n'
54939682 C G Conflicting_interpretations_of_pathogenicity
54939682 C T Benign
$ bcftools query -i 'POS=54939682' clinvar_19_old.vcf.gz -f '%POS %REF %ALT %INFO/CLNSIG\n'
54939682 C G Conflicting_interpretations_of_pathogenicity
54939682 C T Benign
$ grep rs10418277 *.txt
new.txt:rs10418277
tmp.txt:rs10418277
Le problème venait de la POS qui n'était plus convertie en int (suppression de la ligne par erreur ??)
On vérifie
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
python ../../script/pythonScript/clinvar_sbSNP.py --clinvar clinvar_19.vcf.gz --dbSNP dbSNP_common_19.vcf.gz --output tmp.txt
sort tmp.txt | uniq > new.txt
table=/Work/Groups/bisonex/data-alexis/RefSeq/refseq_to_number_only_consensual.txt
python clinvar_sbSNP_old.py --clinvar clinvar_19_old.vcf.gz --dbSNP dbSNP_common_19_old.vcf.gz --output tmp_old.txt --chrm_name_table $table
sort tmp_old.txt | uniq > old.txt
wc -l old.txt new.txt
diff old.txt new.txt
#+end_src
#+RESULTS:
| 535155 | old.txt |
| 535155 | new.txt |
| 1070310 | total |
***** DONE Tester sur chromosome 19 et 20: ok
CLOSED: [2022-12-11 Sun 15:56]
On prépare les données
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
bcftools filter -i 'CHROM="NC_000019.10" | CHROM="NC_000020.11"' /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz -o dbSNP_common_19_20.vcf.gz
bcftools filter -i 'CHROM="NC_000019.10" | CHROM="NC_000020.11"' /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz -o clinvar_19_20.vcf.gz
bcftools filter -i 'CHROM="NC_000019.10" | CHROM="NC_000020.11"' /Work/Groups/bisonex/data-alexis/dbSNP/dbSNP_common.vcf.gz -o dbSNP_common_19_20_old.vcf.gz
bcftools filter -i 'CHROM="19" | CHROM="20"' /Work/Groups/bisonex/data-alexis/clinvar/clinvar.vcf.gz -o clinvar_19_20_old.vcf.gz
#+end_src
#+RESULTS:
On récupère les 2 versions du script
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
git checkout regression ../../script/pythonScript/clinvar_sbSNP.py
cp ../../script/pythonScript/clinvar_sbSNP.py clinvar_sbSNP_old.py
git checkout HEAD ../../script/pythonScript/clinvar_sbSNP.py
#+end_src
#+RESULTS:
On compare
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
python ../../script/pythonScript/clinvar_sbSNP.py clinvar_sbSNP.py --clinvar clinvar_19_20.vcf.gz --dbSNP dbSNP_common_19_20.vcf.gz --output tmp.txt
sort tmp.txt | uniq > new.txt
table=/Work/Groups/bisonex/data-alexis/RefSeq/refseq_to_number_only_consensual.txt
python clinvar_sbSNP_old.py --clinvar clinvar_19_20_old.vcf.gz --dbSNP dbSNP_common_19_20_old.vcf.gz --output tmp_old.txt --chrm_name_table $table
sort tmp_old.txt | uniq > old.txt
wc -l old.txt new.txt
#+end_src
***** DONE Regarder la répartition des différences
CLOSED: [2022-12-11 Sun 16:29]
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
sort /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt | uniq > notpatho.new
sort /Work/Groups/bisonex/data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt | uniq > notpatho.old
comm -23 notpatho.new notpatho.old > nopatho.diff
#+end_src
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
bcftools query -i 'ID=@nopatho.diff' /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz -f '%CHROM\n' | sort | uniq -c
#+end_src
On a principalement des coordonnées non consensuelles (non "NC_", voir notes)
#+RESULTS:
: 2 NC_000002.12
: 18 NC_000003.12
: 2 NC_000004.12
: 2 NC_000005.10
: 14 NC_000006.12
: 6 NC_000007.14
: 2 NC_000009.12
: 1 NC_000010.11
: 6 NC_000014.9
: 1 NC_000015.10
: 3 NC_000016.10
: 3 NC_000017.11
: 1 NC_000019.10
: 1 NC_000020.11
: 1 NC_000021.9
: 2 NC_000022.11
: 16018 NT_113793.3
: 17010 NT_113796.3
: 14 NT_113891.3
: 1 NT_167244.2
: 13 NT_167245.2
: 2 NT_167246.2
: 13 NT_167247.2
: 7 NT_167248.2
: 14 NT_167249.2
: 14857 NT_187361.1
: 92 NT_187367.1
: 1 NT_187369.1
: 13 NT_187381.1
: 54 NT_187383.1
: 6 NT_187499.1
: 46 NT_187502.1
: 13754 NT_187513.1
: 611 NT_187517.1
: 1 NT_187520.1
: 1 NT_187524.1
: 249 NT_187526.1
: 18 NT_187532.1
: 1 NT_187546.1
: 886 NT_187562.1
: 1 NT_187564.1
: 346 NT_187576.1
: 13 NT_187600.1
: 5 NT_187601.1
: 494 NT_187606.1
: 1 NT_187607.1
: 12 NT_187613.1
: 307 NT_187614.1
: 1 NT_187625.1
: 445 NT_187633.1
: 43 NT_187648.1
: 18 NT_187649.1
: 1 NT_187652.1
: 512 NT_187661.1
: 18 NT_187678.1
: 49 NT_187681.1
: 1 NT_187682.1
: 18 NT_187688.1
: 12 NT_187689.1
: 18 NT_187690.1
: 18 NT_187691.1
: 404 NT_187693.1
: 2 NW_003315952.3
: 1 NW_003315970.2
: 203 NW_003571054.1
: 322 NW_003571055.2
: 16 NW_003571056.2
: 16 NW_003571057.2
: 16 NW_003571058.2
: 16 NW_003571059.2
: 16 NW_003571060.1
: 213 NW_003571061.2
: 2 NW_009646201.1
: 322 NW_009646205.1
: 321 NW_009646206.1
: 371 NW_012132914.1
: 1 NW_012132915.1
: 13 NW_012132918.1
: 2 NW_013171801.1
: 1 NW_013171807.1
: 49 NW_015148966.1
: 14 NW_015495298.1
: 2 NW_015495299.1
: 1 NW_016107298.1
: 4 NW_017363813.1
: 2 NW_017852933.1
: 1 NW_018654722.1
: 38 NW_021160001.1
: 1 NW_021160003.1
: 1 NW_021160007.1
: 7 NW_021160017.1
***** DONE Regarder la différence avec la version sans les sites non consensuels: ok !
CLOSED: [2022-12-11 Sun 20:11]
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
sort /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt | uniq > notpatho.new
sort /Work/Groups/bisonex/data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt | uniq > notpatho.old
comm -13 notpatho.new notpatho.old > notpatho.diff
wc -l
#+end_src
#+RESULTS:
: 528 notpatho.diff
Il manque 528 variants
rs1057520103
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
bcftools query -i 'ID=@notpatho.diff' /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz -f '%CHROM\n' | sort | uniq -c
#+end_src
#+RESULTS:
: 528 NC_012920.1
Donc la nouvelle version fonctionne mieux !
ON vérifie bien qu'ils sont dans l'ancienne version et la nouvelle:
$ grep -w -f notpatho.diff /Work/Groups/bisonex/data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt | wc -l
528
$ grep -w -f notpatho.diff /Work/Groups/bisonex/data/d
bSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt
#+end_src
**** DONE Supprimer les sites non consensuels
CLOSED: [2022-12-11 Sun 19:51]
**** DONE Rajouter les mitochondries (vu avec Paul)
CLOSED: [2022-12-13 Tue 17:26]
Ok avec notre version générée. Sur le J: 21155134
$ wc -l dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt
21155065 dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt
$ wc -l ../data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
21155065 ../data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
La différence vient probablement d'une vieille version de clinvar
*** TODO alignement + variant:
**** filtre DP : différent
On en a en plus
$ grep '^NC' filter-depth.vcf | wc -l
86580
Alexis
$ zgrep '^NC' 63003856_S135_DP_over_30.vcf | wc -l
82033
Ne vient pas du filtre sur la profondeur
On a testé
bcftools filter -i 'FORMAT/AD[0:1]<=10' 63003856_S135_DP_over_30.vcf
cftools filter -i 'FORMAT/DP<=30' 63003856_S135_DP_over_30.vcf
Idem pour notre version. Rien ne sort.
Haplotypecaller a les même arguments ??
**** Après haplotypecaller : différent
$ zgrep '^NC' /Work/Groups/bisonex/ref-vcf/63003856_S135.vcf | wc -l
1506894
$ zgrep '^NC' /Work/Users/apraga/bisonex/out/63003856_S135/variantCalling/haplotypecaller/63003856_S135.vcf | wc -l
1631935
**** Aligneur
3 lignes de différences
$ samtools view -c /Work/Users/apraga/bisonex/out/63003856_S135/preprocessing/mapped/63003856_S135.bam
128077211
$ samtools view -c /Work/Groups/bisonex/ref_63003856_S135/63003856_S135.bam
128077207
**** Clean sam
corrigé bug : fichier nettoyé après markduplicates... On relance. En attendant :
$ samtools view -c work/3f/532cde42fabd8973384a817b7539de/sorted.bam
85710928
$ samtools view -c work/3f/532cde42fabd8973384a817b7539de/62903840_S29.bam
85710928
mais :
$ samtools view -c /Work/Groups/bisonex/ref_63003856_S135/63003856_S135_cleaned.bam 128077207
$ samtools view -c /Work/Groups/bisonex/ref_63003856_S135/63003856_S135.bam 128077207
*** TODO 63003856_S135
** Divers
*** DONE Vérifier nombre de reads fastq - bam
CLOSED: [2022-10-09 Sun 22:31]
** TODO Genome in a bottle ?
On n'a pas l'ADN.. séquencer à Centogène ?
* Améliorations
** TODO Quality score recalibration avec un ensemble de fichier
Voir GATK best practice
** TODO Utiliser T-to-T comme références
Semble compliqué avec les nouvelles bases de données
** TODO Macro excel
** TODO Utiliser le XML de clinvar
Extraction sous VCF possible avec
https://github.com/SeqOne/clinvcf
** Annotation
Liste complète
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9252745/
*** TODO Utilise une version allégée de GnomAD (une seule colonne)
*** TODO Digenisme (cf nomenclature omim)
C’est dans le nom de la maladie
* HOLD Implémenter d’autres pipeline
Voir https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-04407-x
** KILL GATK
CLOSED: [2022 -11-11 Fri 20:01]
https://broadinstitute.github.io/warp/docs/Pipelines/Exome_Germline_Single_Sample_Pipeline/README
A priori, respecte les bonnes pratiques
** KILL Essayer snmake avec bonne pratiques
https://github.com/snakemake-workflows/dna-seq-gatk-variant-calling/blob/main/.github/workflows/main.yml
Installer Mamba (micromamba ne fonctionne pas sous nix)
Ne fonctionne pas sous WSL2... MultiQC n’est pas assez à jour
Problèmes de versions...
** KILL Sarek
CLOSED: [2022-12-11 Sun 11:09]
*** Dépendences
**** Nix
#+begin_src sh
nix profile install nixpkgs#mosdepth nixpkgs#python3
nix-shell -p python310Packages.pyyaml --run "nextflow run nf-core/sarek -profile test --executor slurm --queue smp --outdir test -resume"
#+end_src
***** KILL derivation nix pour profile complet
CLOSED: [2022-12-11 Sun 11:09]
**** KILL Sans nix
CLOSED: [2022-09-24 Sat 10:20]
On utilise conda
#+begin_src sh
module unload nix
module load anaconda3@2021.05/gcc-12.1.0
module load nextflow@22.04.0/gcc-12.1.0
module load openjdk@11.0.14.1_1/gcc-12.1.0
nextflow run nf-core/sarek -profile conda,test --executor slurm --queue smp --outdir test -resume
#+end_src
Essai 1: erreurs de permissions, corrigé en relancant le programme
#+begin_quote
Failed to create Conda environment
command: conda create --mkdir --yes --quiet --prefix /Work/Users/apraga/test-sarek/work/conda/env-2d53b1db50de676670cf1a91ef0cf6db bioconda::tabix=1.11
status : 1
message:
NotWritableError: The current user does not have write permissions to a required path.
path: /Home/Users/apraga/.conda/pkgs/urls.txt
uid: 1696
gid: 513
If you feel that permissions on this path are set incorrectly, you can manually
change them by executing
$ sudo chown 1696:513 /Home/Users/apraga/.conda/pkgs/urls.txt
#+end_quote
Corrigé avec
#+begin_src sh
chown 1696:513 /Home/Users/apraga/.conda/pkgs/urls.txt
#+end_src
Mais problème de proxy
*** KILL Dérivation nix pour modules python
CLOSED: [2022-12-11 Sun 11:09]
*** KILL Lancer sarek en mode test
CLOSED: [2022-12-11 Sun 11:09]
#+begin_src sh
nix-shell -p python310Packages.pyyaml --run "nextflow run nf-core/sarek -profile test --executor slurm --queue smp --outdir test -resume"
#+end_src
*** KILL Lancer sarek sur données allégées
CLOSED: [2022-12-11 Sun 11:09]
#+RESULTS:
| e0761a7ba5d10de9e7e97fa331667963925531c0199575bcceafbb13c3147e3f | data-alexis-reference/genome/GRCh38_latest_genomic.fna |
| e0761a7ba5d10de9e7e97fa331667963925531c0199575bcceafbb13c3147e3f | data/genome/GRCh38.p13/genomeRef.fna |
Dict ok si on renome le ficdhier d'origine
#+begin_src sh :dir /ssh:meso:/Work/Groups/bisonex/
sed 's/UR:.*/UR:genomeRef.fna/' data-alexis-reference/genome/GRCh38_latest_genomic.dict > lol.dict
diff lol.dict data/genome/GRCh38.p13/genomeRef.dict
#+end_src
#+RESULTS:
****** DONE dbSNP et dbSNP common: ok
CLOSED: [2023-01-03 Tue 23:17]
sha256sum GCF_000001405.39.gz
452e1112b6339a9b19821c2a226a8a3ba946e92a47e03e6ae464ef8820ee130d GCF_000001405.39.gz
sha256sum data-alexis-reference/dbSNP/GCF_000001405.39.gz
452e1112b6339a9b19821c2a226a8a3ba946e92a47e03e6ae464ef8820ee130d data-alexis-reference/dbSNP/GCF_000001405.39.gz"
sha256sum dbSNP_common.vcf.gz
70dfd9be859c39916598d23b5744cc1fbda04add5840cd90a6d0cd005bd3075b dbSNP_common.vcf.gz
sha256sum data-alexis-reference/dbSNP/dbSNP_common.vcf.gz
70dfd9be859c39916598d23b5744cc1fbda04add5840cd90a6d0cd005bd3075b data-alexis-reference/dbSNP/dbSNP_common.vcf.gz
***** TODO Outils
| | Prod | Test |
| VCFtools | 0.1.17 | 0.1.16 |
| bcftools | 1.14 | 1.16 |
| samtools | 1.14 | 1.13 |
| gatk | 4.2.4.1 | 4.3.0.0 |
On a des versions plus vieilles sauf (le plus important) Gatk
**** KILL Gatk 4.3.0
CLOSED: [2023-01-04 Wed 19:16]
***** KILL Alignement
CLOSED: [2023-01-04 Wed 19:16]
****** DONE Brut
CLOSED: [2022-12-26 Mon 22:03]
Bam alexis
$ samtools view -c /Work/Groups/bisonex/ref_63003856_S135/63003856_S135.bam
128077207
Notre
9f/26cf3d] Cached process > preprocess:BWA_MEM (63003856_S135)
$ samtools view -c work/9f/26cf3deb07b425a3e851be2a7bd782/63003856_S135.bam
On vérifie la sortie
$ samtools view -c out/63003856_S135/preprocessing/mapped/63003856_S135.bam
128077211
Petite différence (< 1e-8) mais selon Alexis, bwa mem est non reproductible. d'autant qu'on utilise une version parallélisée
128077211
****** On vérifie les arguments: ok
#+begin_src
bwa mem \
-R '@RG\tID:sample\tSM:sample\tPL:ILLUMINA\tPM:Miseq\tCN:CHU_Minjoz\tLB:definition_to_add' \
-t 24 \
$INDEX \
63003856_S135_R1_001.fastq.gz 63003856_S135_R2_001.fastq.gz \
| samtools sort --threads 24 -o 63003856_S135.bam -
#+end_src
#+begin_src
bwa mem -R "@RG\tID:$sample\tSM:$sample\tPL:ILLUMINA\tPM:Miseq\tCN:CHU_Minjoz\tLB:definition_to_add" -v 2 -t `nproc` $genomeRef ${fastq[0]} ${fastq[1]} | $samtools sort -@ `nproc` -O BAM -o $tmpDir/$post_bwa
#+end_src
****** DONE Avec gatk 4.2 (version alexis) : idem
$ samtools view -c mapped/63003856_S135.bam
128077211
****** DONE Nettoyé
CLOSED: [2022-12-26 Mon 22:08]
[57/4b5b4c] Cached process > preprocess:GATK4_CLEANSAM (63003856_S135)
$ samtools view -c work/57/4b5b4c647b98bb7099c4d1ba24bd75/63003856_S135.bam
128077211
Et la sortie
$ samtools view -c out/63003856_S135/preprocessing/clean-sam/sorted.bam
128077211
contre
$ samtools view -c /Work/Groups/bisonex/ref_63003856_S135/63003856_S135_cleaned.bam
128077207
On regarde les stats en détails (de la version nettoyée)
$ samtools flagstat work/57/4b5b4c647b98bb7099c4d1ba24bd75/63003856_S135.bam
128077211 + 0 in total (QC-passed reads + QC-failed reads)
126905130 + 0 primary
0 + 0 secondary
1172081 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
127941054 + 0 mapped (99.89% : N/A)
126768973 + 0 primary mapped (99.89% : N/A)
126905130 + 0 paired in sequencing
63452565 + 0 read1
63452565 + 0 read2
125263664 + 0 properly paired (98.71% : N/A)
126676024 + 0 with itself and mate mapped
92949 + 0 singletons (0.07% : N/A)
979608 + 0 with mate mapped to a different chr
675398 + 0 with mate mapped to a different chr (mapQ>=5)
$ samtools flagstat /Work/Groups/bisonex/ref_63003856_S135/63003856_S135_cleaned.bam
128077207 + 0 in total (QC-passed reads + QC-failed reads)
126905130 + 0 primary
0 + 0 secondary
1172077 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
127941051 + 0 mapped (99.89% : N/A)
126768974 + 0 primary mapped (99.89% : N/A)
126905130 + 0 paired in sequencing
63452565 + 0 read1
63452565 + 0 read2
125263790 + 0 properly paired (98.71% : N/A)
126676026 + 0 with itself and mate mapped
92948 + 0 singletons (0.07% : N/A)
979618 + 0 with mate mapped to a different chr
675412 + 0 with mate mapped to a different chr (mapQ>=5)
****** DONE mark duplicate
CLOSED: [2022-12-26 Mon 22:27]
Alexis
$ samtools view -c /Work/Groups/bisonex/ref_63003856_S135/63003856_S135_marked_dup.bam
128077207
Nous (pas de sortie dans out/)
$ samtools view -c work/46/bd75b4547452af36ee2c6b45362922/63003856_S135
128077211
logique car on ne supprime pas de donné...
******* Arguments ok
#+begin_src
gatk --java-options "-Xmx3g" MarkDuplicates \
--INPUT sorted.bam \
--OUTPUT marked_dups.bam \
--METRICS_FILE marked_dups.bam.metrics \
--TMP_DIR . \
--REFERENCE_SEQUENCE genomeRef.fna \
#+end_src
#+begin_src
$gatk MarkDuplicates \
-I $tmpDir/$post_cleanSam \
-O $tmpDir/$post_markDuplicate \
-M $tmpDir/"$sample"_marked_dup.metrix \
--CREATE_INDEX true \
--VERBOSITY WARNING
#+end_src
****** KILL baserecalibrator
CLOSED: [2023-01-04 Wed 19:15]
#+begin_src
$gatk BaseRecalibrator \
-I $tmpDir/$post_markDuplicate \
-R $genomeRef \
--known-sites $dbsnpDir/dbSNP_common.vcf.gz \
-O $tmpDir/$post_BaseRecalibrator
#+end_src
$ cd /Work/Users/apraga/bisonex/out/63003856_S135/preprocessing/baserecalibrator
$ sed 's/sample/63003856_S135/' 63003856_S135.table > 63003856_S135.table2
Les fichiers n'ont pas le même nombre d'erreurs mais assez proches. Sur le premier table, 3 score
#:GATKTable:3:94:%d:%d:%d:;
#:GATKTable:Quantized:Quality quantization map
QualityScore Count QuantizedScore
11 298878631 11
25 542282996 25
34 12846268833 34
vs (référence0)
11 298877785 11
25 542282089 25
34 12846264839 34
******* options ok
#+begin_src
gatk --java-options "-Xmx3g" BaseRecalibrator \
--input marked_dups.bam \
--output 63003856_S135.table \
--reference genomeRef.fna \
\
--known-sites dbSNP_common.vcf.gz \
--tmp-dir . \
#+end_src
****** KILL applybqsr
CLOSED: [2023-01-04 Wed 19:15]
options ok
#+begin_src
gatk --java-options "-Xmx3g" ApplyBQSR \
--input marked_dups.bam \
--output 63003856_S135.bam \
--reference genomeRef.fna \
--bqsr-recal-file 63003856_S135.table \
\
--tmp-dir . \
#+end_src
#+begin_src
$gatk ApplyBQSR -R $genomeRef \
-I $tmpDir/$post_markDuplicate \
--bqsr-recal-file $tmpDir/$post_BaseRecalibrator \
-O $bamDir/$post_ApplyBQSR \
--verbosity WARNING
#+end_src
missing file
***** KILL Variant caling
CLOSED: [2023-01-04 Wed 19:16]
****** KILL haplotypecaller
CLOSED: [2023-01-04 Wed 19:15]
******** options ok
#+begin_src
gatk --java-options "-Xmx3g" HaplotypeCaller \
--input 63003856_S135.bam \
--output 63003856_S135.vcf.gz \
--reference genomeRef.fna \
--dbsnp dbSNP.gz \
\
\
--tmp-dir . \
--max-mnp-distance 2
#+end_src
#+begin_src
$gatk --java-options "-Xmx32g" HaplotypeCaller \
-R $genomeRef \
-I $bamDir/$post_ApplyBQSR \
-O $vcfDir/$post_haplotypecaller \
-D "$dbsnpDir"/GCF_000001405.39.gz \
--max-mnp-distance 2 \
--verbosity WARNING
#+end_src
******** DONE Nombres lignes gatk 4.3.0 : trop différent
CLOSED: [2023-01-04 Wed 19:11]
$ grep '^NC' out/63003856_S135/variantCalling/haplotypecaller/63003856_S135.vcf | wc -l
1631935
$ grep '^NC' /Work/Groups/bisonex/ref-vcf/63003856_S135 .vcf | wc -l
1506894
******** DONE Regarder les flags d'haplotypecaller : nombreuses différences...
CLOSED: [2023-01-04 Wed 19:02]
| --dbsnp /Work/Users/apraga/bisonex/work/08/fca52ac598f21a2812f866bd590792/dbSNP.gz | --dbsnp /mnt/j/bases_de_donnees/dbSNP/GCF_000001405.39.gz |
| --max-mnp-distance 2 | --max-mnp-distance 2 |
| --output 63003856_S135.vcf.gz | --output /mnt/j/working_directory_pipeline_analyse_exome/vcf/63003856_S135.vcf |
| --input 63003856_S135.bam | --input /mnt/j/working_directory_pipeline_analyse_exome/bam/63003856_S135_recalibrated_hg38.bam |
| --reference genomeRef.fna | --reference /mnt/j/bases_de_donnees/genome/GRCh38_latest_genomic.fna |
| --tmp-dir . | --verbosity WARNING |
| --use-posteriors-to-calculate-qual false | --use-posteriors-to-calculate-qual false |
| --dont-use-dragstr-priors false | --dont-use-dragstr-priors false |
| --use-new-qual-calculator true | --use-new-qual-calculator true |
| --annotate-with-num-discovered-alleles false | --annotate-with-num-discovered-alleles false |
| --heterozygosity 0.001 | --heterozygosity 0.001 |
| --indel-heterozygosity 1.25E-4 | --indel-heterozygosity 1.25E-4 |
| --heterozygosity-stdev 0.01 | --heterozygosity-stdev 0.01 |
| --standard-min-confidence-threshold-for-calling 30.0 | --standard-min-confidence-threshold-for-calling 30.0 |
| --max-alternate-alleles 6 | --max-alternate-alleles 6 |
| --max-genotype-count 1024 | --max-genotype-count 1024 |
| --sample-ploidy 2 | --sample-ploidy 2 |
| --num-reference-samples-if-no-call 0 | --num-reference-samples-if-no-call 0 |
| --genotype-assignment-method USE_PLS_TO_ASSIGN | --genotype-assignment-method USE_PLS_TO_ASSIGN |
| --contamination-fraction-to-filter 0.0 | --contamination-fraction-to-filter 0.0 |
| --output-mode EMIT_VARIANTS_ONLY | --output-mode EMIT_VARIANTS_ONLY |
| --all-site-pls false | --all-site-pls false |
| --flow-likelihood-parallel-threads 0 | --gvcf-gq-bands 1 |
| --flow-likelihood-optimized-comp false | --gvcf-gq-bands 2 |
| --flow-use-t0-tag false | --gvcf-gq-bands 3 |
| --flow-probability-threshold 0.003 | --gvcf-gq-bands 4 |
| --flow-remove-non-single-base-pair-indels false | --gvcf-gq-bands 5 |
| --flow-remove-one-zero-probs false | --gvcf-gq-bands 6 |
| --flow-quantization-bins 121 | --gvcf-gq-bands 7 |
| --flow-fill-empty-bins-value 0.001 | --gvcf-gq-bands 8 |
| --flow-symmetric-indel-probs false | --gvcf-gq-bands 9 |
| --flow-report-insertion-or-deletion false | --gvcf-gq-bands 10 |
| --flow-disallow-probs-larger-than-call false | --gvcf-gq-bands 11 |
| --flow-lump-probs false | --gvcf-gq-bands 12 |
| --flow-retain-max-n-probs-base-format false | --gvcf-gq-bands 13 |
| --flow-probability-scaling-factor 10 | --gvcf-gq-bands 14 |
| --flow-order-cycle-length 4 | --gvcf-gq-bands 15 |
| --flow-number-of-uncertain-flows-to-clip 0 | --gvcf-gq-bands 16 |
| --flow-nucleotide-of-first-uncertain-flow T | --gvcf-gq-bands 17 |
| --keep-boundary-flows false | --gvcf-gq-bands 18 |
| --gvcf-gq-bands 1 | --gvcf-gq-bands 19 |
| --gvcf-gq-bands 2 | --gvcf-gq-bands 20 |
| --gvcf-gq-bands 3 | --gvcf-gq-bands 21 |
| --gvcf-gq-bands 4 | --gvcf-gq-bands 22 |
| --gvcf-gq-bands 5 | --gvcf-gq-bands 23 |
| --gvcf-gq-bands 6 | --gvcf-gq-bands 24 |
| --gvcf-gq-bands 7 | --gvcf-gq-bands 25 |
| --gvcf-gq-bands 8 | --gvcf-gq-bands 26 |
| --gvcf-gq-bands 9 | --gvcf-gq-bands 27 |
| --gvcf-gq-bands 10 | --gvcf-gq-bands 28 |
| --gvcf-gq-bands 11 | --gvcf-gq-bands 29 |
| --gvcf-gq-bands 12 | --gvcf-gq-bands 30 |
| --gvcf-gq-bands 13 | --gvcf-gq-bands 31 |
| --gvcf-gq-bands 14 | --gvcf-gq-bands 32 |
| --gvcf-gq-bands 15 | --gvcf-gq-bands 33 |
| --gvcf-gq-bands 16 | --gvcf-gq-bands 34 |
| --gvcf-gq-bands 17 | --gvcf-gq-bands 35 |
| --gvcf-gq-bands 18 | --gvcf-gq-bands 36 |
| --gvcf-gq-bands 19 | --gvcf-gq-bands 37 |
| --gvcf-gq-bands 20 | --gvcf-gq-bands 38 |
| --gvcf-gq-bands 21 | --gvcf-gq-bands 39 |
| --gvcf-gq-bands 22 | --gvcf-gq-bands 40 |
| --gvcf-gq-bands 23 | --gvcf-gq-bands 41 |
| --gvcf-gq-bands 24 | --gvcf-gq-bands 42 |
| --gvcf-gq-bands 25 | --gvcf-gq-bands 43 |
| --gvcf-gq-bands 26 | --gvcf-gq-bands 44 |
| --gvcf-gq-bands 27 | --gvcf-gq-bands 45 |
| --gvcf-gq-bands 28 | --gvcf-gq-bands 46 |
| --gvcf-gq-bands 29 | --gvcf-gq-bands 47 |
| --gvcf-gq-bands 30 | --gvcf-gq-bands 48 |
| --gvcf-gq-bands 31 | --gvcf-gq-bands 49 |
| --gvcf-gq-bands 32 | --gvcf-gq-bands 50 |
| --gvcf-gq-bands 33 | --gvcf-gq-bands 51 |
| --gvcf-gq-bands 34 | --gvcf-gq-bands 52 |
| --gvcf-gq-bands 35 | --gvcf-gq-bands 53 |
| --gvcf-gq-bands 36 | --gvcf-gq-bands 54 |
| --gvcf-gq-bands 37 | --gvcf-gq-bands 55 |
| --gvcf-gq-bands 38 | --gvcf-gq-bands 56 |
| --gvcf-gq-bands 39 | --gvcf-gq-bands 57 |
| --gvcf-gq-bands 40 | --gvcf-gq-bands 58 |
| --gvcf-gq-bands 41 | --gvcf-gq-bands 59 |
| --gvcf-gq-bands 42 | --gvcf-gq-bands 60 |
| --gvcf-gq-bands 43 | --gvcf-gq-bands 70 |
| --gvcf-gq-bands 44 | --gvcf-gq-bands 80 |
| --gvcf-gq-bands 45 | --gvcf-gq-bands 90 |
| --gvcf-gq-bands 46 | --gvcf-gq-bands 99 |
| --gvcf-gq-bands 47 | --floor-blocks false |
| --gvcf-gq-bands 48 | --indel-size-to-eliminate-in-ref-model 10 |
| --gvcf-gq-bands 49 | --disable-optimizations false |
| --gvcf-gq-bands 50 | --dragen-mode false |
| --gvcf-gq-bands 51 | --apply-bqd false |
| --gvcf-gq-bands 52 | --apply-frd false |
| --gvcf-gq-bands 53 | --disable-spanning-event-genotyping false |
| --gvcf-gq-bands 54 | --transform-dragen-mapping-quality false |
| --gvcf-gq-bands 55 | --mapping-quality-threshold-for-genotyping 20 |
| --gvcf-gq-bands 56 | --max-effective-depth-adjustment-for-frd 0 |
| --gvcf-gq-bands 57 | --just-determine-active-regions false |
| --gvcf-gq-bands 58 | --dont-genotype false |
| --gvcf-gq-bands 59 | --do-not-run-physical-phasing false |
| --gvcf-gq-bands 60 | --do-not-correct-overlapping-quality false |
| --gvcf-gq-bands 70 | --use-filtered-reads-for-annotations false |
| --gvcf-gq-bands 80 | --adaptive-pruning false |
| --gvcf-gq-bands 90 | --do-not-recover-dangling-branches false |
| --gvcf-gq-bands 99 | --recover-dangling-heads false |
| --floor-blocks false | --kmer-size 10 |
| --indel-size-to-eliminate-in-ref-model 10 | --kmer-size 25 |
| --disable-optimizations false | --dont-increase-kmer-sizes-for-cycles false |
| --dragen-mode false | --allow-non-unique-kmers-in-ref false |
| --flow-mode NONE | --num-pruning-samples 1 |
| --apply-bqd false | --min-dangling-branch-length 4 |
| --apply-frd false | --recover-all-dangling-branches false |
| --disable-spanning-event-genotyping false | --max-num-haplotypes-in-population 128 |
| --transform-dragen-mapping-quality false | --min-pruning 2 |
| --mapping-quality-threshold-for-genotyping 20 | --adaptive-pruning-initial-error-rate 0.001 |
| --max-effective-depth-adjustment-for-frd 0 | --pruning-lod-threshold 2.302585092994046 |
| --just-determine-active-regions false | --pruning-seeding-lod-threshold 9.210340371976184 |
| --dont-genotype false | --max-unpruned-variants 100 |
| --do-not-run-physical-phasing false | --linked-de-bruijn-graph false |
| --do-not-correct-overlapping-quality false | --disable-artificial-haplotype-recovery false |
| --use-filtered-reads-for-annotations false | --enable-legacy-graph-cycle-detection false |
| --use-flow-aligner-for-stepwise-hc-filtering false | --debug-assembly false |
| --adaptive-pruning false | --debug-graph-transformations false |
| --do-not-recover-dangling-branches false | --capture-assembly-failure-bam false |
| --recover-dangling-heads false | --num-matching-bases-in-dangling-end-to-recover -1 |
| --kmer-size 10 | --error-correction-log-odds -Infinity |
| --kmer-size 25 | --error-correct-reads false |
| --dont-increase-kmer-sizes-for-cycles false | --kmer-length-for-read-error-correction 25 |
| --allow-non-unique-kmers-in-ref false | --min-observations-for-kmer-to-be-solid 20 |
| --num-pruning-samples 1 | --base-quality-score-threshold 18 |
| --min-dangling-branch-length 4 | --dragstr-het-hom-ratio 2 |
| --recover-all-dangling-branches false | --dont-use-dragstr-pair-hmm-scores false |
| --max-num-haplotypes-in-population 128 | --pair-hmm-gap-continuation-penalty 10 |
| --min-pruning 2 | --expected-mismatch-rate-for-read-disqualification 0.02 |
| --adaptive-pruning-initial-error-rate 0.001 | --pair-hmm-implementation FASTEST_AVAILABLE |
| --pruning-lod-threshold 2.302585092994046 | --pcr-indel-model CONSERVATIVE |
| --pruning-seeding-lod-threshold 9.210340371976184 | --phred-scaled-global-read-mismapping-rate 45 |
| --max-unpruned-variants 100 | --disable-symmetric-hmm-normalizing false |
| --linked-de-bruijn-graph false | --disable-cap-base-qualities-to-map-quality false |
| --disable-artificial-haplotype-recovery false | --enable-dynamic-read-disqualification-for-genotyping false |
| --enable-legacy-graph-cycle-detection false | --dynamic-read-disqualification-threshold 1.0 |
| --debug-assembly false | --native-pair-hmm-threads 4 |
| --debug-graph-transformations false | --native-pair-hmm-use-double-precision false |
| --capture-assembly-failure-bam false | --bam-writer-type CALLED_HAPLOTYPES |
| --num-matching-bases-in-dangling-end-to-recover -1 | --dont-use-soft-clipped-bases false |
| --error-correction-log-odds -Infinity | --min-base-quality-score 10 |
| --error-correct-reads false | --smith-waterman JAVA |
| --kmer-length-for-read-error-correction 25 | --emit-ref-confidence NONE |
| --min-observations-for-kmer-to-be-solid 20 | --force-call-filtered-alleles false |
| --likelihood-calculation-engine PairHMM | --soft-clip-low-quality-ends false |
| --base-quality-score-threshold 18 | --allele-informative-reads-overlap-margin 2 |
| --dragstr-het-hom-ratio 2 | --smith-waterman-dangling-end-match-value 25 |
| --dont-use-dragstr-pair-hmm-scores false | --smith-waterman-dangling-end-mismatch-penalty -50 |
| --pair-hmm-gap-continuation-penalty 10 | --smith-waterman-dangling-end-gap-open-penalty -110 |
| --expected-mismatch-rate-for-read-disqualification 0.02 | --smith-waterman-dangling-end-gap-extend-penalty -6 |
| --pair-hmm-implementation FASTEST_AVAILABLE | --smith-waterman-haplotype-to-reference-match-value 200 |
| --pcr-indel-model CONSERVATIVE | --smith-waterman-haplotype-to-reference-mismatch-penalty -150 |
| --phred-scaled-global-read-mismapping-rate 45 | --smith-waterman-haplotype-to-reference-gap-open-penalty -260 |
| --disable-symmetric-hmm-normalizing false | --smith-waterman-haplotype-to-reference-gap-extend-penalty -11 |
| --disable-cap-base-qualities-to-map-quality false | --smith-waterman-read-to-haplotype-match-value 10 |
| --enable-dynamic-read-disqualification-for-genotyping false | --smith-waterman-read-to-haplotype-mismatch-penalty -15 |
| --dynamic-read-disqualification-threshold 1.0 | --smith-waterman-read-to-haplotype-gap-open-penalty -30 |
| --native-pair-hmm-threads 4 | --smith-waterman-read-to-haplotype-gap-extend-penalty -5 |
| --native-pair-hmm-use-double-precision false | --min-assembly-region-size 50 |
| --flow-hmm-engine-min-indel-adjust 6 | --max-assembly-region-size 300 |
| --flow-hmm-engine-flat-insertion-penatly 45 | --active-probability-threshold 0.002 |
| --flow-hmm-engine-flat-deletion-penatly 45 | --max-prob-propagation-distance 50 |
| --pileup-detection false | --force-active false |
| --pileup-detection-enable-indel-pileup-calling false | --assembly-region-padding 100 |
| --num-artificial-haplotypes-to-add-per-allele 5 | --padding-around-indels 75 |
| --artifical-haplotype-filtering-kmer-size 10 | --padding-around-snps 20 |
| --pileup-detection-snp-alt-threshold 0.1 | --padding-around-strs 75 |
| --pileup-detection-indel-alt-threshold 0.5 | --max-extension-into-assembly-region-padding-legacy 25 |
| --pileup-detection-absolute-alt-depth 0.0 | --max-reads-per-alignment-start 50 |
| --pileup-detection-snp-adjacent-to-assembled-indel-range 5 | --enable-legacy-assembly-region-trimming false |
| --pileup-detection-bad-read-tolerance 0.0 | --interval-set-rule UNION |
| --pileup-detection-proper-pair-read-badness true | --interval-padding 0 |
| --pileup-detection-edit-distance-read-badness-threshold 0.08 | --interval-exclusion-padding 0 |
| --pileup-detection-chimeric-read-badness true | --interval-merging-rule ALL |
| --pileup-detection-template-mean-badness-threshold 0.0 | --read-validation-stringency SILENT |
| --pileup-detection-template-std-badness-threshold 0.0 | --seconds-between-progress-updates 10.0 |
| --bam-writer-type CALLED_HAPLOTYPES | --disable-sequence-dictionary-validation false |
| --dont-use-soft-clipped-bases false | --create-output-bam-index true |
| --override-fragment-softclip-check false | --create-output-bam-md5 false |
| --min-base-quality-score 10 | --create-output-variant-index true |
| --smith-waterman JAVA | --create-output-variant-md5 false |
| --emit-ref-confidence NONE | --max-variants-per-shard 0 |
| --force-call-filtered-alleles false | --lenient false |
| --reference-model-deletion-quality 30 | --add-output-sam-program-record true |
| --soft-clip-low-quality-ends false | --add-output-vcf-command-line true |
| --allele-informative-reads-overlap-margin 2 | --cloud-prefetch-buffer 40 |
| --smith-waterman-dangling-end-match-value 25 | --cloud-index-prefetch-buffer -1 |
| --smith-waterman-dangling-end-mismatch-penalty -50 | --disable-bam-index-caching false |
| --smith-waterman-dangling-end-gap-open-penalty -110 | --sites-only-vcf-output false |
| --smith-waterman-dangling-end-gap-extend-penalty -6 | --help false |
| --smith-waterman-haplotype-to-reference-match-value 200 | --version false |
| --smith-waterman-haplotype-to-reference-mismatch-penalty -150 | --showHidden false |
| --smith-waterman-haplotype-to-reference-gap-open-penalty -260 | --QUIET false |
| --smith-waterman-haplotype-to-reference-gap-extend-penalty -11 | --use-jdk-deflater false |
| --smith-waterman-read-to-haplotype-match-value 10 | --use-jdk-inflater false |
| --smith-waterman-read-to-haplotype-mismatch-penalty -15 | --gcs-max-retries 20 |
| --smith-waterman-read-to-haplotype-gap-open-penalty -30 | --gcs-project-for-requester-pays |
| --smith-waterman-read-to-haplotype-gap-extend-penalty -5 | --disable-tool-default-read-filters false |
| --flow-assembly-collapse-hmer-size 0 | --minimum-mapping-quality 20 |
| --flow-assembly-collapse-partial-mode false | --disable-tool-default-annotations false |
| --flow-filter-alleles false | --enable-all-annotations false |
| --flow-filter-alleles-qual-threshold 30.0 | --allow-old-rms-mapping-quality-annotation-data false |
| --flow-filter-alleles-sor-threshold 3.0 | Version="4.2.4.1",Date="December 3, 2022 at 1:20:38 AM CET"> |
| --flow-filter-lone-alleles false |
| --flow-filter-alleles-debug-graphs false |
| --min-assembly-region-size 50 |
| --max-assembly-region-size 300 |
| --active-probability-threshold 0.002 |
| --max-prob-propagation-distance 50 |
| --force-active false |
| --assembly-region-padding 100 |
| --padding-around-indels 75 |
| --padding-around-snps 20 |
| --padding-around-strs 75 |
| --max-extension-into-assembly-region-padding-legacy 25 |
| --max-reads-per-alignment-start 50 |
| --enable-legacy-assembly-region-trimming false |
| --interval-set-rule UNION |
| --interval-padding 0 |
| --interval-exclusion-padding 0 |
| --interval-merging-rule ALL |
| --read-validation-stringency SILENT |
| --seconds-between-progress-updates 10.0 |
| --disable-sequence-dictionary-validation false |
| --create-output-bam-index true |
| --create-output-bam-md5 false |
| --create-output-variant-index true |
| --create-output-variant-md5 false |
| --max-variants-per-shard 0 |
| --lenient false |
| --add-output-sam-program-record true |
| --add-output-vcf-command-line true |
| --cloud-prefetch-buffer 40 |
| --cloud-index-prefetch-buffer -1 |
| --disable-bam-index-caching false |
| --sites-only-vcf-output false |
| --help false |
| --version false |
| --showHidden false |
| --verbosity INFO |
| --QUIET false |
| --use-jdk-deflater false |
| --use-jdk-inflater false |
| --gcs-max-retries 20 |
| --gcs-project-for-requester-pays |
| --disable-tool-default-read-filters false |
| --minimum-mapping-quality 20 |
| --disable-tool-default-annotations false |
| --enable-all-annotations false |
| --allow-old-rms-mapping-quality-annotation-data false" |
| Version="4.3.0.0",Date="December 16, 2022 at 12:51:03 AM CET"> |
****** KILL [#B] filterDepth : 21 en trop
CLOSED: [2023-01-04 Wed 19:16]
Nouvelle version (correcton bug markdupicates)
#+begin_src sh
grep '^NC' out/63003856_S135/variantCalling/filter-depth.vcf |wc -l
82054
#+end_src
Alexis
#+begin_src sh
zgrep '^NC' /Work/Groups/bisonex/ref_63003856_S135/63003856_S135_DP_over_30.vcf | wc -l
82033
#+end_src
Ne vient pas du filtre sur la profondeur:
bcftools filter -i 'FORMAT/AD[0:1]<=10' 63003856_S135_DP_over_30.vcf
bcftools filter -i 'FORMAT/DP<=30' 63003856_S135_DP_over_30.vcf
Idem pour notre version. Rien ne sort.
On compare le nombre de lignes
#+begin_src sh
bgzip out/63003856_S135/variantCalling/filter-depth.vcf
tabix /Work/Groups/bisonex/ref_63003856_S135/63003856_S135_DP_over_30.vcf.gz
tabix out/63003856_S135/variantCalling/filter-depth.vcf.gz
bcftools isec out/63003856_S135/variantCalling/filter-depth.vcf.gz /Work/Groups/bisonex/ref_63003856_S135/63003856_S135_DP_over_30.vcf.gz -p compare-filter-depth
find compare-filter-depth/ -type f -exec wc -l {} \;
84763 compare-filter-depth/sites.txt
710 compare-filter-depth/0001.vcf
8 compare-filter-depth/README.txt
85339 compare-filter-depth/0002.vcf
85340 compare-filter-depth/0003.vcf
725 compare-filter-depth/0000.vcf
#+end_src
****** KILL exclude SNP + consensual : 34 en trop !!
CLOSED: [2023-01-04 Wed 19:16]
$ grep '^NC' out/63003856_S135/variantCalling/filter-polymorphisms.vcf | wc -l
8898
vs
$ grep '^NC' /Work/Groups/bisonex/ref_63003856_S135/63003856_S135_DP_over_30_not_SNP_consensual_sequence.vcf | wc -l
8864
****** KILL Filter technical variants
CLOSED: [2023-01-04 Wed 19:16]
**** Gatk 4.2.4 (même version qu'alexis)
***** TODO Variant calling
****** TODO haplotypecaller: mieux mais non identique !
******* DONE Nombres lignes gatk 4.2.2 : faible différence
CLOSED: [2023-01-04 Wed 19:18]
$ zgrep '^NC' 63003856_S135.vcf.gz | wc -l
1506931
$ grep '^NC' /Work/Groups/bisonex/ref-vcf/63003856_S135 .vcf | wc -l
1506894
******* DONE Flags la même version de gatk 4.2.2 : ok identique
CLOSED: [2023-01-04 Wed 19:09]
##GATKCommandLine=<ID=HaplotypeCaller,CommandLine="HaplotypeCaller
| ",Version="4.2.4.1",Date="January 4, 2023 at 1:46:41 AM CET"> | Version="4.2.4.1",Date="December 3, 2022 at 1:20:38 AM CET"> |
| --dbsnp /Work/Users/apraga/bisonex/work/5d/feb81028d262d7701bed0a759ff6f6/dbSNP.gz | --dbsnp /mnt/j/bases_de_donnees/dbSNP/GCF_000001405.39.gz |
| --max-mnp-distance 2 | --max-mnp-distance 2 |
| --output 63003856_S135.vcf.gz | --output /mnt/j/working_directory_pipeline_analyse_exome/vcf/63003856_S135.vcf |
| --input 63003856_S135.bam | --input /mnt/j/working_directory_pipeline_analyse_exome/bam/63003856_S135_recalibrated_hg38.bam |
| --reference genomeRef.fna | --reference /mnt/j/bases_de_donnees/genome/GRCh38_latest_genomic.fna |
| --tmp-dir . | --verbosity WARNING |
| --use-posteriors-to-calculate-qual false | --use-posteriors-to-calculate-qual false |
| --dont-use-dragstr-priors false | --dont-use-dragstr-priors false |
| --use-new-qual-calculator true | --use-new-qual-calculator true |
| --annotate-with-num-discovered-alleles false | --annotate-with-num-discovered-alleles false |
| --heterozygosity 0.001 | --heterozygosity 0.001 |
| --indel-heterozygosity 1.25E-4 | --indel-heterozygosity 1.25E-4 |
| --heterozygosity-stdev 0.01 | --heterozygosity-stdev 0.01 |
| --standard-min-confidence-threshold-for-calling 30.0 | --standard-min-confidence-threshold-for-calling 30.0 |
| --max-alternate-alleles 6 | --max-alternate-alleles 6 |
| --max-genotype-count 1024 | --max-genotype-count 1024 |
| --sample-ploidy 2 | --sample-ploidy 2 |
| --num-reference-samples-if-no-call 0 | --num-reference-samples-if-no-call 0 |
| --genotype-assignment-method USE_PLS_TO_ASSIGN | --genotype-assignment-method USE_PLS_TO_ASSIGN |
| --contamination-fraction-to-filter 0.0 | --contamination-fraction-to-filter 0.0 |
| --output-mode EMIT_VARIANTS_ONLY | --output-mode EMIT_VARIANTS_ONLY |
| --all-site-pls false | --all-site-pls false |
| --gvcf-gq-bands 1 | --gvcf-gq-bands 1 |
| --gvcf-gq-bands 2 | --gvcf-gq-bands 2 |
| --gvcf-gq-bands 3 | --gvcf-gq-bands 3 |
| --gvcf-gq-bands 4 | --gvcf-gq-bands 4 |
| --gvcf-gq-bands 5 | --gvcf-gq-bands 5 |
| --gvcf-gq-bands 6 | --gvcf-gq-bands 6 |
| --gvcf-gq-bands 7 | --gvcf-gq-bands 7 |
| --gvcf-gq-bands 8 | --gvcf-gq-bands 8 |
| --gvcf-gq-bands 9 | --gvcf-gq-bands 9 |
| --gvcf-gq-bands 10 | --gvcf-gq-bands 10 |
| --gvcf-gq-bands 11 | --gvcf-gq-bands 11 |
| --gvcf-gq-bands 12 | --gvcf-gq-bands 12 |
| --gvcf-gq-bands 13 | --gvcf-gq-bands 13 |
| --gvcf-gq-bands 14 | --gvcf-gq-bands 14 |
| --gvcf-gq-bands 15 | --gvcf-gq-bands 15 |
| --gvcf-gq-bands 16 | --gvcf-gq-bands 16 |
| --gvcf-gq-bands 17 | --gvcf-gq-bands 17 |
| --gvcf-gq-bands 18 | --gvcf-gq-bands 18 |
| --gvcf-gq-bands 19 | --gvcf-gq-bands 19 |
| --gvcf-gq-bands 20 | --gvcf-gq-bands 20 |
| --gvcf-gq-bands 21 | --gvcf-gq-bands 21 |
| --gvcf-gq-bands 22 | --gvcf-gq-bands 22 |
| --gvcf-gq-bands 23 | --gvcf-gq-bands 23 |
| --gvcf-gq-bands 24 | --gvcf-gq-bands 24 |
| --gvcf-gq-bands 25 | --gvcf-gq-bands 25 |
| --gvcf-gq-bands 26 | --gvcf-gq-bands 26 |
| --gvcf-gq-bands 27 | --gvcf-gq-bands 27 |
| --gvcf-gq-bands 28 | --gvcf-gq-bands 28 |
| --gvcf-gq-bands 29 | --gvcf-gq-bands 29 |
| --gvcf-gq-bands 30 | --gvcf-gq-bands 30 |
| --gvcf-gq-bands 31 | --gvcf-gq-bands 31 |
| --gvcf-gq-bands 32 | --gvcf-gq-bands 32 |
| --gvcf-gq-bands 33 | --gvcf-gq-bands 33 |
| --gvcf-gq-bands 34 | --gvcf-gq-bands 34 |
| --gvcf-gq-bands 35 | --gvcf-gq-bands 35 |
| --gvcf-gq-bands 36 | --gvcf-gq-bands 36 |
| --gvcf-gq-bands 37 | --gvcf-gq-bands 37 |
| --gvcf-gq-bands 38 | --gvcf-gq-bands 38 |
| --gvcf-gq-bands 39 | --gvcf-gq-bands 39 |
| --gvcf-gq-bands 40 | --gvcf-gq-bands 40 |
| --gvcf-gq-bands 41 | --gvcf-gq-bands 41 |
| --gvcf-gq-bands 42 | --gvcf-gq-bands 42 |
| --gvcf-gq-bands 43 | --gvcf-gq-bands 43 |
| --gvcf-gq-bands 44 | --gvcf-gq-bands 44 |
| --gvcf-gq-bands 45 | --gvcf-gq-bands 45 |
| --gvcf-gq-bands 46 | --gvcf-gq-bands 46 |
| --gvcf-gq-bands 47 | --gvcf-gq-bands 47 |
| --gvcf-gq-bands 48 | --gvcf-gq-bands 48 |
| --gvcf-gq-bands 49 | --gvcf-gq-bands 49 |
| --gvcf-gq-bands 50 | --gvcf-gq-bands 50 |
| --gvcf-gq-bands 51 | --gvcf-gq-bands 51 |
| --gvcf-gq-bands 52 | --gvcf-gq-bands 52 |
| --gvcf-gq-bands 53 | --gvcf-gq-bands 53 |
| --gvcf-gq-bands 54 | --gvcf-gq-bands 54 |
| --gvcf-gq-bands 55 | --gvcf-gq-bands 55 |
| --gvcf-gq-bands 56 | --gvcf-gq-bands 56 |
| --gvcf-gq-bands 57 | --gvcf-gq-bands 57 |
| --gvcf-gq-bands 58 | --gvcf-gq-bands 58 |
| --gvcf-gq-bands 59 | --gvcf-gq-bands 59 |
| --gvcf-gq-bands 60 | --gvcf-gq-bands 60 |
| --gvcf-gq-bands 70 | --gvcf-gq-bands 70 |
| --gvcf-gq-bands 80 | --gvcf-gq-bands 80 |
| --gvcf-gq-bands 90 | --gvcf-gq-bands 90 |
| --gvcf-gq-bands 99 | --gvcf-gq-bands 99 |
| --floor-blocks false | --floor-blocks false |
| --indel-size-to-eliminate-in-ref-model 10 | --indel-size-to-eliminate-in-ref-model 10 |
| --disable-optimizations false | --disable-optimizations false |
| --dragen-mode false | --dragen-mode false |
| --apply-bqd false | --apply-bqd false |
| --apply-frd false | --apply-frd false |
| --disable-spanning-event-genotyping false | --disable-spanning-event-genotyping false |
| --transform-dragen-mapping-quality false | --transform-dragen-mapping-quality false |
| --mapping-quality-threshold-for-genotyping 20 | --mapping-quality-threshold-for-genotyping 20 |
| --max-effective-depth-adjustment-for-frd 0 | --max-effective-depth-adjustment-for-frd 0 |
| --just-determine-active-regions false | --just-determine-active-regions false |
| --dont-genotype false | --dont-genotype false |
| --do-not-run-physical-phasing false | --do-not-run-physical-phasing false |
| --do-not-correct-overlapping-quality false | --do-not-correct-overlapping-quality false |
| --use-filtered-reads-for-annotations false | --use-filtered-reads-for-annotations false |
| --adaptive-pruning false | --adaptive-pruning false |
| --do-not-recover-dangling-branches false | --do-not-recover-dangling-branches false |
| --recover-dangling-heads false | --recover-dangling-heads false |
| --kmer-size 10 | --kmer-size 10 |
| --kmer-size 25 | --kmer-size 25 |
| --dont-increase-kmer-sizes-for-cycles false | --dont-increase-kmer-sizes-for-cycles false |
| --allow-non-unique-kmers-in-ref false | --allow-non-unique-kmers-in-ref false |
| --num-pruning-samples 1 | --num-pruning-samples 1 |
| --min-dangling-branch-length 4 | --min-dangling-branch-length 4 |
| --recover-all-dangling-branches false | --recover-all-dangling-branches false |
| --max-num-haplotypes-in-population 128 | --max-num-haplotypes-in-population 128 |
| --min-pruning 2 | --min-pruning 2 |
| --adaptive-pruning-initial-error-rate 0.001 | --adaptive-pruning-initial-error-rate 0.001 |
| --pruning-lod-threshold 2.302585092994046 | --pruning-lod-threshold 2.302585092994046 |
| --pruning-seeding-lod-threshold 9.210340371976184 | --pruning-seeding-lod-threshold 9.210340371976184 |
| --max-unpruned-variants 100 | --max-unpruned-variants 100 |
| --linked-de-bruijn-graph false | --linked-de-bruijn-graph false |
| --disable-artificial-haplotype-recovery false | --disable-artificial-haplotype-recovery false |
| --enable-legacy-graph-cycle-detection false | --enable-legacy-graph-cycle-detection false |
| --debug-assembly false | --debug-assembly false |
| --debug-graph-transformations false | --debug-graph-transformations false |
| --capture-assembly-failure-bam false | --capture-assembly-failure-bam false |
| --num-matching-bases-in-dangling-end-to-recover -1 | --num-matching-bases-in-dangling-end-to-recover -1 |
| --error-correction-log-odds -Infinity | --error-correction-log-odds -Infinity |
| --error-correct-reads false | --error-correct-reads false |
| --kmer-length-for-read-error-correction 25 | --kmer-length-for-read-error-correction 25 |
| --min-observations-for-kmer-to-be-solid 20 | --min-observations-for-kmer-to-be-solid 20 |
| --base-quality-score-threshold 18 | --base-quality-score-threshold 18 |
| --dragstr-het-hom-ratio 2 | --dragstr-het-hom-ratio 2 |
| --dont-use-dragstr-pair-hmm-scores false | --dont-use-dragstr-pair-hmm-scores false |
| --pair-hmm-gap-continuation-penalty 10 | --pair-hmm-gap-continuation-penalty 10 |
| --expected-mismatch-rate-for-read-disqualification 0.02 | --expected-mismatch-rate-for-read-disqualification 0.02 |
| --pair-hmm-implementation FASTEST_AVAILABLE | --pair-hmm-implementation FASTEST_AVAILABLE |
| --pcr-indel-model CONSERVATIVE | --pcr-indel-model CONSERVATIVE |
| --phred-scaled-global-read-mismapping-rate 45 | --phred-scaled-global-read-mismapping-rate 45 |
| --disable-symmetric-hmm-normalizing false | --disable-symmetric-hmm-normalizing false |
| --disable-cap-base-qualities-to-map-quality false | --disable-cap-base-qualities-to-map-quality false |
| --enable-dynamic-read-disqualification-for-genotyping false | --enable-dynamic-read-disqualification-for-genotyping false |
| --dynamic-read-disqualification-threshold 1.0 | --dynamic-read-disqualification-threshold 1.0 |
| --native-pair-hmm-threads 4 | --native-pair-hmm-threads 4 |
| --native-pair-hmm-use-double-precision false | --native-pair-hmm-use-double-precision false |
| --bam-writer-type CALLED_HAPLOTYPES | --bam-writer-type CALLED_HAPLOTYPES |
| --dont-use-soft-clipped-bases false | --dont-use-soft-clipped-bases false |
| --min-base-quality-score 10 | --min-base-quality-score 10 |
| --smith-waterman JAVA | --smith-waterman JAVA |
| --emit-ref-confidence NONE | --emit-ref-confidence NONE |
| --force-call-filtered-alleles false | --force-call-filtered-alleles false |
| --soft-clip-low-quality-ends false | --soft-clip-low-quality-ends false |
| --allele-informative-reads-overlap-margin 2 | --allele-informative-reads-overlap-margin 2 |
| --smith-waterman-dangling-end-match-value 25 | --smith-waterman-dangling-end-match-value 25 |
| --smith-waterman-dangling-end-mismatch-penalty -50 | --smith-waterman-dangling-end-mismatch-penalty -50 |
| --smith-waterman-dangling-end-gap-open-penalty -110 | --smith-waterman-dangling-end-gap-open-penalty -110 |
| --smith-waterman-dangling-end-gap-extend-penalty -6 | --smith-waterman-dangling-end-gap-extend-penalty -6 |
| --smith-waterman-haplotype-to-reference-match-value 200 | --smith-waterman-haplotype-to-reference-match-value 200 |
| --smith-waterman-haplotype-to-reference-mismatch-penalty -150 | --smith-waterman-haplotype-to-reference-mismatch-penalty -150 |
| --smith-waterman-haplotype-to-reference-gap-open-penalty -260 | --smith-waterman-haplotype-to-reference-gap-open-penalty -260 |
| --smith-waterman-haplotype-to-reference-gap-extend-penalty -11 | --smith-waterman-haplotype-to-reference-gap-extend-penalty -11 |
| --smith-waterman-read-to-haplotype-match-value 10 | --smith-waterman-read-to-haplotype-match-value 10 |
| --smith-waterman-read-to-haplotype-mismatch-penalty -15 | --smith-waterman-read-to-haplotype-mismatch-penalty -15 |
| --smith-waterman-read-to-haplotype-gap-open-penalty -30 | --smith-waterman-read-to-haplotype-gap-open-penalty -30 |
| --smith-waterman-read-to-haplotype-gap-extend-penalty -5 | --smith-waterman-read-to-haplotype-gap-extend-penalty -5 |
| --min-assembly-region-size 50 | --min-assembly-region-size 50 |
| --max-assembly-region-size 300 | --max-assembly-region-size 300 |
| --active-probability-threshold 0.002 | --active-probability-threshold 0.002 |
| --max-prob-propagation-distance 50 | --max-prob-propagation-distance 50 |
| --force-active false | --force-active false |
| --assembly-region-padding 100 | --assembly-region-padding 100 |
| --padding-around-indels 75 | --padding-around-indels 75 |
| --padding-around-snps 20 | --padding-around-snps 20 |
| --padding-around-strs 75 | --padding-around-strs 75 |
| --max-extension-into-assembly-region-padding-legacy 25 | --max-extension-into-assembly-region-padding-legacy 25 |
| --max-reads-per-alignment-start 50 | --max-reads-per-alignment-start 50 |
| --enable-legacy-assembly-region-trimming false | --enable-legacy-assembly-region-trimming false |
| --interval-set-rule UNION | --interval-set-rule UNION |
| --interval-padding 0 | --interval-padding 0 |
| --interval-exclusion-padding 0 | --interval-exclusion-padding 0 |
| --interval-merging-rule ALL | --interval-merging-rule ALL |
| --read-validation-stringency SILENT | --read-validation-stringency SILENT |
| --seconds-between-progress-updates 10.0 | --seconds-between-progress-updates 10.0 |
| --disable-sequence-dictionary-validation false | --disable-sequence-dictionary-validation false |
| --create-output-bam-index true | --create-output-bam-index true |
| --create-output-bam-md5 false | --create-output-bam-md5 false |
| --create-output-variant-index true | --create-output-variant-index true |
| --create-output-variant-md5 false | --create-output-variant-md5 false |
| --max-variants-per-shard 0 | --max-variants-per-shard 0 |
| --lenient false | --lenient false |
| --add-output-sam-program-record true | --add-output-sam-program-record true |
| --add-output-vcf-command-line true | --add-output-vcf-command-line true |
| --cloud-prefetch-buffer 40 | --cloud-prefetch-buffer 40 |
| --cloud-index-prefetch-buffer -1 | --cloud-index-prefetch-buffer -1 |
| --disable-bam-index-caching false | --disable-bam-index-caching false |
| --sites-only-vcf-output false | --sites-only-vcf-output false |
| --help false | --help false |
| --version false | --version false |
| --showHidden false | --showHidden false |
| --verbosity INFO | --QUIET false |
| --QUIET false | --use-jdk-deflater false |
| --use-jdk-deflater false | --use-jdk-inflater false |
| --use-jdk-inflater false | --gcs-max-retries 20 |
| --gcs-max-retries 20 | --gcs-project-for-requester-pays |
| --gcs-project-for-requester-pays | --disable-tool-default-read-filters false |
| --disable-tool-default-read-filters false | --minimum-mapping-quality 20 |
| --minimum-mapping-quality 20 | --disable-tool-default-annotations false |
| --disable-tool-default-annotations false | --enable-all-annotations false |
| --enable-all-annotations false | --allow-old-rms-mapping-quality-annotation-data false |
| --allow-old-rms-mapping-quality-annotation-data false | |
****** filter depth : Toujours la même différence...
$ grep '^NC' filter-depth.vcf | wc -l
82054
$ zgrep '^NC' /Work/Groups/bisonex/ref_63003856_S135/63003856_S135_DP_over_30.vcf.gz | wc -l
82033
Non lié à la profondeur : on teste avec
bcftools filter -i 'FORMAT/DP<=30' filter-depth.vcf
bcftools filter -i 'FORMAT/AD[0:1]<=10' filter-depth.vcf
****** Vérifier qu'en utilsant 2 filtres différents on a bien la même chose : oui
$ bcftools filter -e 'FORMAT/DP<=30' 63003856_S135.vcf.gz | bcftools filter -e 'FORMAT/AD[0:1]<=10' -o two-filters.vcf
$ grep '^NC' two-filters.vcf | wc -l
82054
***** Tester bwa en séquentiel
#+title: Bisonex
* Biblio
Comparaison WDL, Cromwell, nextflow
https://www.nature.com/articles/s41598-021-99288-8
Nextflow = bon compromis ?
* Idées
** Validation analytique
*** Génération de reads avec variants connus
Comparaison de génération ADN (2019)
https://academic.oup.com/bfg/article/19/1/49/5680294
**** SimuSCop
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
*** Utiliser données GCAT et uploader le notre ?
https://www.nature.com/articles/ncomms7275
*** Genome in a bottle : NA12878 + autres
2 versions :
1. Depuis un fastq correspondant à Illumina https://github.com/genome-in-a-bottle/giab_data_indexes
puis on compare le VCF avec les "high confidence" (Article : https://www.nature.com/articles/s41587-019-0054-x
)
- 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/data/NA12878/analysis/Illumina_PlatinumGenomes_NA12877_NA12878_09162015/hg38/2.0.1/NA12878/
2. On séquence directement NA12878
* 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...
* Données
** TODO Vérifier qualité données sur mesocentre
*** STRT BAM
picard ValidateSamFile
On regarde juste le code d'erreur (0 = pas d'erreur)
*** STRT Fastq
fastqc
Il faut ensuite extraire les zip and chercher les erreur dedans
** TODO Lister données sur mesocentre
* Nouveau 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>
*** STRT Télécharger
- [X] Genome de référence
- [X] dbSNP
- [X] OMIM
- [X] VEP 20G
- [X] transcriptome (spip)
- [ ] Refseq
*** DONE Télécharger les données avec nextflow
CLOSED: [2022-09-13 Tue 21:37]
*** 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 ../database/RefSeq/refseq_to_number_only_consensual.txt \
--dbSNP $snp --output prod.txt
wc -l prod.txt
zgrep '^NC' dbSNP_common_chr20.vcf.gz | wc -l
#+end_src
#+RESULTS:
| 518832 | prod.txt |
| 518846 | |
***** KILL classification clinvar codée dbSNP ?
CLOSED: [2022-12-04 Sun 14:38]
Sur le chromosome 20
*Attention* CLNSIG a plusieurs champs (séparé par une virgule)
On y accède avec INFO/CLNSIG[*]
Ensuite, chaque item peut avoir plusieurs haploïdie (séparé par un |). IL faut donc utiliser une regexp
NB: *ne pas mettre la condition* dans une variable !!
Pour avoir les clinvar patho, on veut 5 mais pas 255 (= autre) pour la classification !`
Il faut également les likely patho et conflicting
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%INFO/CLNSIG\n' dbSNP_common_chr20.vcf.gz -i \
'INFO/CLNSIG[*]~"^5|" | INFO/CLNSIG[*]=="5" | INFO/CLNSIG[*]~"|5" | INFO/CLNSIG[*]~"^4|" | INFO/CLNSIG[*]=="4" | INFO/CLNSIG[*]~"|4" | INFO/CLNSIG[*]~"^12|" | INFO/CLNSIG[*]=="12" | INFO/CLNSIG[*]~"|12"' | sort
#+end_src
#+RESULTS:
| . | . | 12 | | | | | | | | |
| . | 12 | 0 | 2 | | | | | | | |
| 2 | 3 | 2 | 2 | 2 | 5 | . | | | | |
| . | 2 | 3 | 2 | 2 | 4 | | | | | |
| . | . | 3 | 12 | 3 | | | | | | |
| . | 5 | 2 | . | | | | | | | |
| . | . | . | 5 | 2 | 2 | | | | | |
| . | 9 | 9 | 9 | 5 | 5 | 2 | 3 | 2 | 3 | 2 |
Si on les exclut :
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%ID\n' dbSNP_common_chr20.vcf.gz -e \
'INFO/CLNSIG[*]~"^5|" | INFO/CLNSIG[*]=="5" | INFO/CLNSIG[*]~"|5" | INFO/CLNSIG[*]~"4" | INFO/CLNSIG[*]~"12"' | sort | uniq > common-notpatho.txt
#+end_src
#+RESULTS:
#+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 ../database/RefSeq/refseq_to_number_only_consensual.txt \
--dbSNP $snp --output tmp.txt
sort tmp.txt | uniq > common-notpatho-alexis.txt
wc -l common-notpatho-alexis.txt
#+end_src
#+RESULTS:
: 518832 common-notpatho-alexis.txt
On en a 6 de plus que la version d'Alexis mais quelques différences
Ceux d'Alexis qui manquent:
#+begin_src sh :dir ~/code/bisonex/test_isec
comm -23 common-notpatho-alexis.txt common-notpatho.txt > alexis-only.txt
cat alexis-only.txt
#+end_src
#+RESULTS:
| rs1064039 |
| rs3833341 |
| rs73598374 |
On les teste dans clinvar et dbSNP
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'ID=@alexis-only.txt' dbSNP_common_chr20.vcf.gz
bcftools query -f '%POS\n' -i 'ID=@alexis-only.txt' dbSNP_common_chr20.vcf.gz > alexis-only-pos.txt
while read -r line; do
bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'POS='$line clinvar_chr20.vcf.gz
done < alexis-only-pos.txt
# bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'POS=23637790' clinvar_chr20.vcf.gz
#+end_src
#+RESULTS:
| 764018 | A | ACAGGTCAAT,ACAGGT | .,5 | 2,. | |
| 23637790 | C | G,T | .,.,12 | | |
| 44651586 | C | A,G,T | .,.,.,5 | 2 | 2 |
| 764018 | A | ACAGGTCAAT | Benign | | |
| 23637790 | C | T | Benign | | |
| 44651586 | C | T | Benign | | |
On a donc une discordance entre clinvar et dbSNP.
On dirait qu'ils ont mal fait l'intersection avec clinvar.
Par exemple https://www.ncbi.nlm.nih.gov/snp/rs3833341#clinical_significance
Tu as l'impression qu'il y a un 1 clinvar bénin et 1 patho.
En cherchant par NM, tu vois qu'il est bénin sur clinvar car il y a d'autres soumissions ! https://www.ncbi.nlm.nih.gov/clinvar/variation/262235/
Confirmation sur nos bases de données :
$ bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'POS=764018' dbSNP_common_chr20.vcf.gz
764018 A ACAGGTCAAT,ACAGGT .,5|2,.
$ bcftools query -f '%POS %REF %ALT %INFO/CLNSIG\n' -i 'POS=764018' clinvar_chr20.vcf.gz
764018 A ACAGGTCAAT Benign
***** KILL Corriger script alexi
CLOSED: [2022-12-04 Sun 13:03]
Gère clinvar patho, probablement patho ou conflicting !
***** HOLD Rtg tools
****** Test
1. Générer SDf file
#+begin_src sh
rtg format genomeRef.fna -o genomeRef.sdf
#+end_src
2. Pour les bases de donnés, il faut l'option --sample ALT sinon on a
#+begin_src
$ rtg vcfeval -b dbSNP_common.vcf.gz -c clinvar.vcf.gz -o test -t genomeRef.sdf/^C
VCF header does not contain a FORMAT field named GQ
Error: Record did not contain enough samples: NC_000001.11 10001 rs1570391677 A,C . PASS RS=1570391677;dbSNPBuildID=154;SSR=0;PSEUDOGENEINFO=DDX11L1:100287102;VC=SNV;R5;GNO;FREQ=KOREAN:0.9891,0.0109,.|SGDP_PRJ:0,1,.|dbGaP_PopFreq:1,.,0;COMMON
#+end_src
Essai intersection clinvar (patho ou non) dbSNP
- faux négatif = dbSNP common qui ne sont pas dans clinvar
- faux positif = clinvar qui ne sont pas dbSNP common
- vrai positif = clinvar qui sont dans dbSNP common
- vrai positif baseline = dbSNP common qui sont dans clinvar
On calcule le nombre de lignes
#+begin_src ssh
zgrep '^[^#]' /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz | wc -l
for i in *.vcf.gz; do echo $i; zgrep '^[^#]' $i | wc -l; done
#+end_src
| clinvar | 1493470 |
| fn.vcf.gz | 22330220 |
| fp.vcf.gz | 1222529 |
| tp-baseline.vcf.gz | 131040 |
| tp.vcf.gz | 136638 |
À noter qu'on ne retrouve pas tout clinvar...
1222529 + 131040 = 1353569 < 1493470
certains régions ne sont pas traitées :
#+begin_quote
Evaluation too complex (50002 unresolved paths, 34891 iterations) at reference region NC_000001.11:790930-790970. Variants in this region will not be included in results
#+end_quote
#+begin_src sh
grep 'not be included' vcfeval.log | wc -l
56192
#+end_src
Le total est quand même inférieur
On veut les clinvar non patho dans dbSNP soit les faux négatif (dbSNP common not contenu dans clinvar patho)
#+begin_src sh
bcftools filter -i 'INFO/CLNSIG="Pathogenic"' /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz -o /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar-patho.vcf.gz
tabix /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar-patho.vcf.gz
#+end_src
On lance le script (dbSNP common et clinvar = 9h)
#+begin_src sh
#!/bin/bash
#SBATCH --nodes=1
#SBATCH -p smp
#SBATCH --time=12:00:00
#SBATCH --mem=12G
dir=/Work/Groups/bisonex/data
dbSNP=$dir/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz
clinvar=$dir/clinvar/GRCh38/clinvar-patho.vcf.gz
genome=$dir/genome/GRCh38.p13/genomeRef.sdf
srun rtg vcfeval -b $dbSNP -c $clinvar -o common-not-patho -t $genome --sample ALT
#+end_src
****** HOLD Voir pour régions complexes non traitées
***** DONE bcftools isec : non
CLOSED: [2022-11-27 Sun 00:38]
#+begin_src sh
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -p common
#+end_src
On vérifie bien que les 2 fichiers commons on le même nombre de lignes
#+begin_src sh
$ grep -e '^NC' 0002.vcf | wc -l
74302
alex@gentoo ~/code/bisonex/data/common $ grep -e '^NC' 0003.vcf | wc -l
74302
#+end_src
****** DONE Impact option -n
CLOSED: [2022-10-23 Sun 13:56]
Mais en spécifiant -n =2:
#+begin_src sh
$ bedtools intersect -a dbSNP_common.vcf.gz -b clinvar.vcf.gz
74978
#+end_src
Si on ne regarde que les variants, on retrouve bien 74302
#+begin_src sh
rg "^NC" none_sorted.vcf | wc -l
#+end_src
NB : test fait avec
#+begin_src
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c none -n =2 -w 1 | sort > none.vcf
sort common/0003.vcf > common/0003_sorted.vcf
comm -13 common/0003_sorted.vcf none_sorted.vcf
#+end_src
****** DONE Géstion des duplicates: -c none
CLOSED: [2022-10-23 Sun 13:56]
Si on ne garde que ceux avec REF et ALT identiques
#+begin_src sh
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c none -n =2 -w 1 | wc -l
74978
#+end_src
Si on garde tout
#+begin_src sh
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c all -n =2 -w 1 | wc -l
137777
#+end_src
Pour regarder la différence :
#+begin_src sh
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c none -n =2 -w 1 | sort > none_sorted.vcf
bcftools isec dbSNP_common.vcf.gz clinvar.vcf.gz -c all -n =2 -w 1 | sort > all_sorted.vcf
comm -13 none_sorted.vcf all_sorted.vcf | head
#+end_src
Sur un exemple,on a bien des variants différents
****** DONE Suppression des clinvar patho
CLOSED: [2022-10-23 Sun 18:55]
Semble faire le travail vu que dbSNP_commo a 23194960 lignes (donc ~80 000 de moins)
#+begin_src sh
$ bcftools isec -e 'INFO/CLNSIG="Pathogenic" & INFO/CLNSIG="Pathogenic/Likely_pathogenic"' -c none -n~10 dbSNP_common.vcf.gz clinvar.vcf.gz | wc -l
Note: -w option not given, printing list of sites...
23119984
#+end_src
Par contre, l'o'ption -w ou -p fait des ficher "data"...
Après un nouvel essai, plus de problème
#+begin_src
$ bcftools isec -e 'INFO/CLNSIG="Pathogenic" & INFO/CLNSIG="Pathogenic/Likely_pathogenic"' -c none -n=1 dbSNP_common.vcf.gz clinvar.vcf.gz -w 1 -o lol.vcf.gz
$ zcat lol.vcf.gz | wc -l
23120660
#+end_src
À noter le choix de l'option -n qui change entre "=1" et "~10"...
En effet "=1" = au moins 1 fichier et "~10" fait exactement dans le premier et non dans le second
#+begin_src
$ bcftools isec -e 'INFO/CLNSIG="Pathogenic" & INFO/CLNSIG="Pathogenic/Likely_pathogenic"' -c none -n~10 dbSNP_common.vcf.gz clinvar.vcf.gz -w 1 -o lol.vcf.gz
$ zcat lol.vcf.gz | wc -l
23120660
#+end_src
****** DONE Valider avec Alexis : bcftool isec
CLOSED: [2022-11-07 Mon 21:42 ]
****** DONE Pourquoi nombre de lignes différentes avec la version d'Alexis -> isec ne gère pas plusieurs ALT
CLOSED: [2022-11-26 Sat 23:36]
Grosse différence !
#+begin_src
$ wc -l ID_of_common_snp_not_clinvar_patho.txt
23119915 ID_of_common_snp_not_clinvar_patho.txt
$ wc -l /Work/Users/apraga/bisonex/database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
85820 /Work/Users/apraga/bisonex/database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
#+end_src
À noter que tout dbSNP = 23194960
******* Clinvar classe 4 ? Moins mais toujours trop
#+begin_src
$ zgrep '^NC' tmp.vcf.gz | wc -l
21081654
#+end_src
******* Comparer les ID et regarder ceux en plus
#+begin_src sh
bcftools isec -e 'INFO/CLNSIG="Pathogenic"' -c none -n~10 /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz -w 1 -o tmp.vcf.gz
zgrep -o -e 'rs[[:digit:]]\' tmp.vcf.gz | sort | id_sorted.txt
sort ../database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt > reference_sorted.txt
comm -23 id_sorted.txt reference_sorted.txt > unique1.txt
#+end_src
Par exemple
#+begin_src sh
zgrep rs1000000561 ../database/dbSNP/dbSNP_common.vcf.gz
#+end_src
NC_000002.12 136732859 rs1000000561 ACG A,ACGCG . PASS RS=1000000561;dbSNPBuildID=151;SSR=0;VC=INDEL;GNO;FREQ=ALSPAC:0.2506,0.7494,.|TOMMO:0.9971,0.002865,.|TWINSUK:0.2473,0.7527,.|dbGaP_PopFreq:0.993,0.006943,8.902e-05;COMMON
Attention, clinvar est en numéro de chromosomoe et dbSNP en NC...
Normalement, géré lors du calcul d'intersection !
Ce SNP n'est pas dans clinvar (vérifié dans UCSC)
******* Tester sur chromosome 20
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools view --regions NC_000020.11 ../database/dbSNP/dbSNP_common.vcf.gz -o dbSNP_common_chr20.vcf.gz
bcftools view --regions 20 ../database/clinvar/clinvar.vcf.gz -o clinvar_chr20.vcf.gz
tabix dbSNP_common_chr20.vcf.gz
tabix clinvar_chr20.vcf.gz
#+end_src
#+RESULTS:
Attention à bien renommer clinvar !
#+begin_src sh :dir ~/code/bisonex/test_isec
mv clinvar_chr20.vcf.gz clinvar_chr20_notremapped.vcf.gz
bcftools annotate --rename-chrs chromosome_mapping.txt clinvar_chr20_notremapped.vcf.gz -o clinvar_chr20.vcf.gz
#+end_src
#+RESULTS:
*ATTENTION*: sans indexer les vcf, les fichiers seront *VIDES*
*ATTENTION*: par défaut les filtres s'appliquent sur les 2. Cela est un problème si on joue sur l'inclusion et non l'exclusion
Attention: vérifier la conventdion de nommage des chromosomes
******** Test pathogene: ne prend pas en compte les multi-allèles ????
On teste l'intersection dbsnp et clinvar patho ainsi que le complémentaire
#+begin_src sh :dir ~/code/bisonex/test_isec
clinvar=clinvar_chr20_patho.vcf.gz
snp=dbSNP_common_chr20.vcf.gz
bcftools index $clinvar
bcftools index $snp
bcftools filter -i 'INFO/CLNSIG="Pathogenic"' clinvar_chr20.vcf.gz -o $clinvar
bcftools isec $snp $clinvar -p tmp
for i in tmp/*.vcf ; do echo $i; grep '^[^#]' $i | wc -l; done
#+end_src
#+RESULTS:
| tmp/0000.vcf |
| 518846 |
| tmp/0001.vcf |
| 0 |
| tmp/0002.vcf |
| 0 |
| tmp/0003.vcf |
| 0 |
Aucun clinvar patho... Clairement faux !
Autre méthode : on inclut tous les SNP et clinvar patho et on regarde ceux uniquement dans dbsnp
#+begin_src sh :dir ~/code/bisonex/test_isec
snp=dbSNP_common_chr20.vcf.gz
clinvar=clinvar_chr20.vcf.gz
bcftools isec -n=2 -i - -i 'INFO/CLNSIG="Pathogenic"' $snp $clinvar -p tmp
# grep '^[^#]' tmp/0000.vcf | wc -l
#+end_src
#+RESULTS:
Soit tout dbsnp donc rien
Note : on ne peut pas exclure les clinvar patho directement
#+begin_src sh :dir ~/code/bisonex/test_isec
snp=dbSNP_common_chr20.vcf.gz
clinvar=clinvar_chr20.vcf.gz
bcftools isec -i - -e 'INFO/CLNSIG="Pathogenic"' $snp $clinvar -p tmp
for i in tmp/*.vcf ; do echo $i; grep '^[^#]' $i | wc -l; done
#+end_src
Car on ne peut plus faire la différence !
Si on utilise la 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 ../database/RefSeq/refseq_to_number_only_consensual.txt \
--dbSNP $snp --output tmp.txt
sort tmp.txt > common-notpatho-alexis.txt
wc -l common-notpatho-alexis.txt
#+end_src
#+RESULTS:
: 518832 common-notpatho-alexis.txt
Si on cherche les clinvar patho (donc non présent dans la sortie)
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%ID\n' dbSNP_common_chr20.vcf.gz | sort > all.txt
sort common-notpatho-alexis.txt > alexis.txt
comm -23 all.txt alexis.txt > patho.txt
#+end_src
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools query -f '%POS\n' -i 'ID=@patho.txt' dbSNP_common_chr20.vcf.gz -o pos.txt
for pos in $(cat pos.txt); do
bcftools query -f '%CHROM %POS %ID %REF %ALT\n' -i 'POS='$pos dbSNP_common_chr20.vcf.gz
bcftools query -f '%CHROM %POS %ID %REF %ALT %INFO/CLNSIG\n' -i 'POS='$pos clinvar_chr20.vcf.gz
echo "------"
done
#+end_src
#+RESULTS:
| NC_000020.11 | 3234173 | rs3827075 | T | A,C,G | |
| NC_000020.11 | 3234173 | 262001 | T | G | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 3234173 | 1072511 | T | TGGCGAAGC | Pathogenic |
| NC_000020.11 | 3234173 | 208613 | TGGCGAAGC | G | Pathogenic |
| NC_000020.11 | 3234173 | 1312 | TGGCGAAGC | T | Pathogenic |
| ------ | | | | | |
| NC_000020.11 | 4699605 | rs1799990 | A | G | |
| NC_000020.11 | 4699605 | 13397 | A | G | Benign/Likely_benign |
| ------ | | | | | |
| NC_000020.11 | 10652589 | rs1131695 | G | A,C,T | |
| NC_000020.11 | 10652589 | 163705 | G | . | Benign |
| NC_000020.11 | 10652589 | 143063 | G | A | Benign |
| NC_000020.11 | 10652589 | 234555 | G | C | Pathogenic |
| ------ | | | | | |
| NC_000020.11 | 10658574 | rs1801138 | G | A,T | |
| NC_000020.11 | 10658574 | 42481 | G | A | Benign |
| NC_000020.11 | 10658574 | 992651 | G | T | Likely_pathogenic |
| NC_000020.11 | 10658574 | 213550 | GC | A | Pathogenic |
| ------ | | | | | |
| NC_000020.11 | 10672794 | rs79338570 | G | A,C | |
| NC_000020.11 | 10672794 | 255557 | G | A | Benign/Likely_benign |
| NC_000020.11 | 10672794 | 594067 | G | C | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 10672794 | 1324603 | G | GGA | Likely_pathogenic |
| ------ | | | | | |
| NC_000020.11 | 18525868 | rs146917730 | C | T | |
| NC_000020.11 | 18525868 | 811603 | C | T | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 25390747 | rs373200654 | G | C | |
| NC_000020.11 | 25390747 | 338000 | G | C | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 32800145 | rs2424926 | C | G,T | |
| NC_000020.11 | 32800145 | 338173 | C | G | Benign |
| NC_000020.11 | 32800145 | 338174 | C | T | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 33412656 | rs35938843 | C | G,T | |
| NC_000020.11 | 33412656 | 220958 | C | T | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 45891622 | rs181943893 | G | A,C,T | |
| NC_000020.11 | 45891622 | 459632 | G | C | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 45891622 | 797035 | G | T | Likely_benign |
| NC_000020.11 | 45891622 | 1572689 | GCTA | G | Likely_benign |
| ------ | | | | | |
| NC_000020.11 | 54171651 | rs35873579 | G | A,T | |
| NC_000020.11 | 54171651 | 285894 | G | A | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 54171651 | 1373583 | G | C | Uncertain_significance |
| NC_000020.11 | 54171651 | 895614 | G | T | Benign/Likely_benign |
| ------ | | | | | |
| NC_000020.11 | 62172726 | rs36106901 | G | A | |
| NC_000020.11 | 62172726 | 981031 | G | A | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 63349782 | rs1044396 | G | A,C | |
| NC_000020.11 | 63349782 | 93427 | G | A | Benign |
| NC_000020.11 | 63349782 | 857384 | G | C | Conflicting_interpretations_of_pathogenicity |
| ------ | | | | | |
| NC_000020.11 | 63414925 | rs1801545 | G | A,C,T | |
| NC_000020.11 | 63414925 | 194284 | G | A | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 63414925 | 129337 | G | C | Benign |
| NC_000020.11 | 63414925 | 851545 | GG | CA | Uncertain_significance |
| ------ | | | | | |
On a donc plusieurs problèmes :
1. isec devrait fonctionner au moins sur
| NC_000020.11 | 25390747 | rs373200654 | G | C | |
| NC_000020.11 | 25390747 | 338000 | G | C | Conflicting_interpretations_of_pathogenicity |
On teste juste sur cette ligne
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools filter -i 'POS=25390747' clinvar_chr20.vcf.gz -o clinvar_test.vcf.gz
bcftools filter -i 'POS=25390747' dbSNP_common_chr20.vcf.gz -o dbSNP_test.vcf.gz
#+end_src
On retrouve bien la ligne dans l'intersection...
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools filter -i 'POS=25390747' clinvar_chr20.vcf.gz -o clinvar_test.vcf.gz
bcftools index dbSNP_test.vcf.gz dbSNP_test.vcf.gz
bcftools index dbSNP_test.vcf.gz clinvar_test.vcf.gz
bcftools isec dbSNP_test.vcf.gz clinvar_test.vcf.gz -p test
#+end_src
#+RESULTS:
2. isec ne semble pas fonctionner sur en cas d'ALT multiples
| NC_000020.11 | 32800145 | rs2424926 | C | G,T | |
| NC_000020.11 | 32800145 | 338173 | C | G | Benign |
| NC_000020.11 | 32800145 | 338174 | C | T | Conflicting_interpretations_of_pathogenicity |
| | | | | | |
3. s'il y a plusieurs variantions à une position, il faut bien vérifier que tous ne sont pas patho.
La version d'Alexis le fait bien
| NC_000020.11 | 3234173 | rs3827075 | T | A,C,G | |
| NC_000020.11 | 3234173 | 262001 | T | G | Conflicting_interpretations_of_pathogenicity |
| NC_000020.11 | 3234173 | 1072511 | T | TGGCGAAGC | Pathogenic |
| NC_000020.11 | 3234173 | 208613 | TGGCGAAGC | G | Pathogenic |
| NC_000020.11 | 3234173 | 1312 | TGGCGAAGC | T | Pathogenic |
****** DONE Voir si isec gère les multiallélique (chr20) : non, impossible de faire marcher
CLOSED: [2022-11-27 Sun 00:37]
******* DONE chr20 en prenant un patho clinvar aussi dans dbSNP
CLOSED: [2022-11-27 Sun 00:37]
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools filter dbSNP_common_chr20.vcf.gz -i 'POS=10652589' -o test_dbsnp.vcf.gz
bcftools filter clinvar_chr20.vcf.gz -i 'POS=10652589' -o test_clinvar.vcf.gz
bcftools index test_dbsnp.vcf.gz
bcftools index test_clinvar.vcf.gz
#+end_src
#+RESULTS:
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools isec test_dbsnp.vcf.gz test_clinvar.vcf.gz -p tmp
grep '^[^#]' tmp/0002.vcf
grep '^[^#]' tmp/0003.vcf
#+end_src
#+RESULTS:
Même en biallélique, ne fonctionne pas.
Testé en modifiant test_dbsnp !
Fonctionne avec un variant par ligne
****** DONE isec en coupant les sites multialléliques: non
CLOSED: [2022-11-27 Sun 00:37]
******* DONE Exemple simple ok
CLOSED: [2022-11-27 Sun 00:34]
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools filter -i 'POS=10652589' dbSNP_common_chr20.vcf.gz -o dbsnp_mwi.vcf.gz
bcftools filter -i 'POS=10652589' clinvar_chr20.vcf.gz -o clinvar_mwi.vcf.gz
bcftools index -f dbsnp_mwi.vcf.gz
bcftools index -f clinvar_mwi.vcf.gz
bcftools isec dbsnp_mwi.vcf.gz clinvar_mwi.vcf.gz -n=2
#+end_src
#+RESULTS:
Même en biallélique, ne fonctionne pas.
Chr 20
Avec les fichiers du teste précédent
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools norm -m -any dbsnp_mwi.vcf.gz -o dbsnp_mwi_norm.vcf.gz
bcftools index dbsnp_mwi_norm.vcf.gz
bcftools isec dbsnp_mwi_norm.vcf.gz clinvar_mwi.vcf.gz -n=2
#+end_src
#+RESULTS:
| NC_000020.11 | 10652589 | G | A | 11 |
| NC_000020.11 | 10652589 | G | C | 11 |
******* TODO Sur dbSNP chr20 non
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools norm -m -any dbSNP_common_chr20 -o dbSNP_common_chr20_norm.vcf.gz
#+end_src
#+begin_src sh :dir ~/code/bisonex/test_isec
bcftools isec -i 'INFO/CLNSIG="Pathogenic"' dbSNP_common_chr20_norm.vcf.gz clinvar_chr20.vcf.gz -p tmp
#+end_src
#+RESULTS:
***** DONE Essai bedtools intersect
#+begin_src sh
bedtools intersect -a dbSNP_common.vcf.gz -b clinvar.vcf.gz
#+end_src
$ wc -l intersect.vcf
220206 intersect.vcf
** TODO Dépendences avec Nix
*** DONE GATK
CLOSED: [2022-10-21 Fri 21:59]
*** WAIT BioDBHTS
Contribuer pull request
*** DONE BioExtAlign
CLOSED: [2022-10-22 Sat 00:38]
*** WAIT BioBigFile
Revoir si on peut utliser kent dernière version
Contribuer pull request
*** HOLD rtg-tools
Convertir clinvar NC
*** DONE simuscop
CLOSED: [2022-12-30 Fri 22:31]
*** DONE Spip
CLOSED: [2022-12-04 Sun 12:49]
Pas de pull request
*** DONE R + packages
CLOSED: [2022-11-19 Sat 21:05]
** DONE Exécution
CLOSED: [2022-09-13 Tue 21:37]
*** KILL test Bionix
*** KILL Implémenter execution avec Nix ?
Voir https://academic.oup.com/gigascience/article/9/11/giaa121/5987272?login=false
pour un exemple.
Probablement plus simple d’utiliser Nix pour gestion de l’environnement et snakemake pour l’exécution
Pas d’accès internet depuis le cluster
*** DONE nextflow
CLOSED: [2022-09-13 Tue 21:37]
** TODO Preprocessing avec nextflow
*** TODO Map to reference
**** TODO Sample ID dans header
/Work/Users/apraga/bisonex/out/63003856_S135/preprocessing/baserecalibrator
*** DONE Mark duplicate
CLOSED: [2022-10-09 Sun 22:30]
*** DONE Recalibrate base quality score
CLOSED: [2022-10-09 Sun 22:30]
** DONE Variant calling avec Nextflow
CLOSED: [2022-11-19 Sat 21:34]
*** DONE Haplotype caller
CLOSED: [2022-10-09 Sun 22:40]
*** DONE Filter variants
CLOSED: [2022-10-09 Sun 22:40]
*** DONE Filter common snp not clinvar path
CLOSED: [2022-11-07 Mon 23:00]
Voir [[*common dbSNP not clinvar patho][common dbSNP not clinvar patho]]
*** DONE Filter variant only in consensual sequence
CLOSED: [2022-11-08 Tue 22:23]
*** DONE Filter technical variants
CLOSED: [2022-11-19 Sat 21:34]
** TODO Annotation avec nextflow
*** TODO VEP
*** TODO Spip
*** TODO Filtrer après VEP
On doit pouvoir se passer d'un script R avec bcftools
** STRT Tester version d'alexis avec Nix
*** DONE Ajouter clinvar
CLOSED: [2022-11-13 Sun 19:37]
*** DONE Alignement
CLOSED: [2022-11-13 Sun 12:52]
*** DONE Haplotype caller
CLOSED: [2022-11-13 Sun 13:00]
*** TODO Filter
- [X] depth
- [ ] comon snp not path
Problème avec liste des ID
**** TODO variant annotation
Besoin de vep
*** TODO Variant calling
** TODO MAJ avec picard
Normalement, GATK inclut picard mais la dernière version utilise picard pour certains outils
https://gatk.broadinstitute.org/hc/en-us/articles/9570266920219--Tool-Documentation-Index
A compléter après validation
*** TODO markduplicates
La dernière version dans la documentation utilise picard !!
** TODO Bwa-mem2 au lieu de bwa mem
https://github.com/bwa-mem2/bwa-mem2
** TODO Parallélisation haplotypecaller
spark est en beta, ne pas utiliser
parallélisation du pauvre : se restreindre à un chromosome avec -L et paralléliser sur le nombre de chromosome
** KILL CRAM au lieu de SAM ?
CLOSED: [2022-12-30 Fri 20:38]
Version compressée de bam mais :
#+begin_quote
All GATK tools that take in mapped read data expect a BAM file as primary format. Some support the CRAM format, but we have observed performance issues when working directly from CRAM files, so in our own work we convert CRAM to BAM first, and we only use CRAM for archival purposes
#+end_quote
Source: https://gatk.broadinstitute.org/hc/en-us/articles/360035890791-SAM-or-BAM-or-CRAM-Mapped-sequence-data-formats
* TODO Tests
** TODO Test de non régression avec version ALexis avec nix
*** DONE ID common snp
CLOSED: [2022-11-19 Sat 21:36]
#+begin_src
$ wc -l ID_of_common_snp.txt
23194290 ID_of_common_snp.txt
$ wc -l /Work/Users/apraga/bisonex/database/dbSNP/ID_of_common_snp.txt
23194290 /Work/Users/apraga/bisonex/database/dbSNP/ID_of_common_snp.txt
#+end_src
*** DONE ID common snp not clinvar patho
CLOSED: [2022-12-11 Sun 20:11]
**** DONE Vérification du problème
CLOSED: [2022-12-11 Sun 16:30]
Sur le J:
21155134 /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt.ref
Version de "non-régression"
21155076 database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
Nouvelle version
23193391 /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt
Si on enlève les doublons
$ sort database/dbSNP/ID_of_common_snp_not_clinvar_patho.txt | uniq > old.txt
$ wc -l old.txt
21107097 old.txt
$ sort /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt | uniq > new.txt
$ wc -l new.txt
21174578 new.txt
$ sort /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt.ref | uniq > ref.txt
$ wc -l ref.txt
21107155 ref.txt
Si on regarde la différence
comm -23 ref.txt old.txt
rs1052692
rs1057518973
rs1057518973
rs11074121
rs112848754
rs12573787
rs145033890
rs147889095
rs1553904159
rs1560294695
rs1560296615
rs1560310926
rs1560325547
rs1560342418
rs1560356225
rs1578287542
...
On cherche le premier
bcftools query -i 'ID="rs1052692"' database/dbSNP/dbSNP_common.vcf.gz -f '%CHROM %POS %REF %ALT\n'
NC_000019.10 1619351 C A,T
Il est bien patho...
$ bcftools query -i 'POS=1619351' database/clinvar/clinvar.vcf.gz -f '%CHROM %POS %REF %ALT %INFO/CLNSIG\n'
19 1619351 C T Conflicting_interpretations_of_pathogenicity
On vérifie pour tous les autres
$ comm -23 ref.txt old.txt > tocheck.txt
On génère les régions à vérifier (chromosome number:position)
$ bcftools query -i 'ID=@tocheck.txt' database/dbSNP/dbSNP_common.vcf.gz -f '%CHROM\t%POS\n' > tocheck.pos
On génère le mapping inverse (chromosome number -> NC)
$ awk ' { t = $1; $1 = $2; $2 = t; print; } ' database/RefSeq/refseq_to_number_only_consensual.txt > mapping.txt
On remap clinvar
$ bcftools annotate --rename-chrs mapping.txt database/clinvar/clinvar.vcf.gz -o clinvar_remapped.vcf.gz
$ tabix clinvar_remapped.vcf.gz
Enfin, on cherche dans clinvar la classification
$ bcftools query -R tocheck.pos clinvar_remapped.vcf.gz -f '%CHROM %POS %INFO/CLNSIG\n'
$ bcftools query -R tocheck.pos database/dbSNP/dbSNP_common.vcf.gz -f '%CHROM %POS %ID \n' | grep '^NC'
#+RESULTS:
**** DONE Comprendre pourquoi la nouvelle version donne un résultat différent
CLOSED: [2022-12-11 Sun 20:11]
***** DONE Même version dbsnp et clinvar ?
CLOSED: [2022-12-10 Sat 23:02]
Clinvar différent !
$ bcftools stats clinvar.gz
clinvar (Alexis)
SN 0 number of samples: 0
SN 0 number of records: 1492828
SN 0 number of no-ALTs: 965
SN 0 number of SNPs: 1338007
SN 0 number of MNPs: 5562
SN 0 number of indels: 144580
SN 0 number of others: 3714
SN 0 number of multiallelic sites: 0
SN 0 number of multiallelic SNP sites: 0
clinvar (new)
SN 0 number of samples: 0
SN 0 number of records: 1493470
SN 0 number of no-ALTs: 965
SN 0 number of SNPs: 1338561
SN 0 number of MNPs: 5565
SN 0 number of indels: 144663
SN 0 number of others: 3716
SN 0 number of multiallelic sites: 0
SN 0 number of multiallelic SNP sites: 0
***** DONE Mettre à jour clinvar et dbnSNP pour travailler sur les mêm bases
CLOSED: [2022-12-11 Sun 12:10]
Problème persiste
***** DONE Supprimer la conversion en int du chromosome
CLOSED: [2022-12-10 Sat 19:29]
***** KILL Même NC ?
CLOSED: [2022-12-10 Sat 19:29]
$ zgrep "contig=<ID=NC_\(.*\)" clinvar/GRCh38/clinvar.vcf.gz > contig.clinvar
$ diff contig.txt contig.clinvar
< ##contig=<ID=NC_012920.1>
***** DONE Tester sur chromosome 19: ok
CLOSED: [2022-12-11 Sun 13:53]
On prépare les données
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
bcftools filter -i 'CHROM="NC_000019.10"' /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz -o dbSNP_common_19.vcf.gz
bcftools filter -i 'CHROM="NC_000019.10"' /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz -o clinvar_19.vcf.gz
bcftools filter -i 'CHROM="NC_000019.10"' /Work/Groups/bisonex/data-alexis/dbSNP/dbSNP_common.vcf.gz -o dbSNP_common_19_old.vcf.gz
bcftools filter -i 'CHROM="19"' /Work/Groups/bisonex/data-alexis/clinvar/clinvar.vcf.gz -o clinvar_19_old.vcf.gz
#+end_src
On récupère les 2 versions du script
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
git checkout regression ../../script/pythonScript/clinvar_sbSNP.py
cp ../../script/pythonScript/clinvar_sbSNP.py clinvar_sbSNP_old.py
git checkout HEAD ../../script/pythonScript/clinvar_sbSNP.py
#+end_src
#+RESULTS:
On compare
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
python ../../script/pythonScript/clinvar_sbSNP.py clinvar_sbSNP.py --clinvar clinvar_19.vcf.gz --dbSNP dbSNP_common_19.vcf.gz --output tmp.txt
sort tmp.txt | uniq > new.txt
table=/Work/Groups/bisonex/data-alexis/RefSeq/refseq_to_number_only_consensual.txt
python clinvar_sbSNP_old.py --clinvar clinvar_19_old.vcf.gz --dbSNP dbSNP_common_19_old.vcf.gz --output tmp_old.txt --chrm_name_table $table
sort tmp_old.txt | uniq > old.txt
wc -l old.txt new.txt
#+end_src
#+RESULTS:
| 535155 | old.txt |
| 535194 | new.txt |
| 1070349 | total |
Si on prend le premier manquant dans new, il est conflicting patho donc il ne devrait pas y être...
$ bcftools query -i 'ID="rs10418277"' dbSNP
_common_19.vcf.gz -f '%CHROM %POS %REF %ALT\n'
NC_000019.10 54939682 C G,T
$ bcftools query -i 'ID="rs10418277"' dbSNP_common_19_old.vcf.gz -f '%CHROM %POS %REF %ALT\n'
NC_000019.10 54939682 C G,T
$ bcftools query -i 'POS=54939682' clinvar_19.vcf.gz -f '%POS %REF %ALT %INFO/CLNSIG\n'
54939682 C G Conflicting_interpretations_of_pathogenicity
54939682 C T Benign
$ bcftools query -i 'POS=54939682' clinvar_19_old.vcf.gz -f '%POS %REF %ALT %INFO/CLNSIG\n'
54939682 C G Conflicting_interpretations_of_pathogenicity
54939682 C T Benign
$ grep rs10418277 *.txt
new.txt:rs10418277
tmp.txt:rs10418277
Le problème venait de la POS qui n'était plus convertie en int (suppression de la ligne par erreur ??)
On vérifie
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
python ../../script/pythonScript/clinvar_sbSNP.py --clinvar clinvar_19.vcf.gz --dbSNP dbSNP_common_19.vcf.gz --output tmp.txt
sort tmp.txt | uniq > new.txt
table=/Work/Groups/bisonex/data-alexis/RefSeq/refseq_to_number_only_consensual.txt
python clinvar_sbSNP_old.py --clinvar clinvar_19_old.vcf.gz --dbSNP dbSNP_common_19_old.vcf.gz --output tmp_old.txt --chrm_name_table $table
sort tmp_old.txt | uniq > old.txt
wc -l old.txt new.txt
diff old.txt new.txt
#+end_src
#+RESULTS:
| 535155 | old.txt |
| 535155 | new.txt |
| 1070310 | total |
***** DONE Tester sur chromosome 19 et 20: ok
CLOSED: [2022-12-11 Sun 15:56]
On prépare les données
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
bcftools filter -i 'CHROM="NC_000019.10" | CHROM="NC_000020.11"' /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz -o dbSNP_common_19_20.vcf.gz
bcftools filter -i 'CHROM="NC_000019.10" | CHROM="NC_000020.11"' /Work/Groups/bisonex/data/clinvar/GRCh38/clinvar.vcf.gz -o clinvar_19_20.vcf.gz
bcftools filter -i 'CHROM="NC_000019.10" | CHROM="NC_000020.11"' /Work/Groups/bisonex/data-alexis/dbSNP/dbSNP_common.vcf.gz -o dbSNP_common_19_20_old.vcf.gz
bcftools filter -i 'CHROM="19" | CHROM="20"' /Work/Groups/bisonex/data-alexis/clinvar/clinvar.vcf.gz -o clinvar_19_20_old.vcf.gz
#+end_src
#+RESULTS:
On récupère les 2 versions du script
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
git checkout regression ../../script/pythonScript/clinvar_sbSNP.py
cp ../../script/pythonScript/clinvar_sbSNP.py clinvar_sbSNP_old.py
git checkout HEAD ../../script/pythonScript/clinvar_sbSNP.py
#+end_src
#+RESULTS:
On compare
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
python ../../script/pythonScript/clinvar_sbSNP.py clinvar_sbSNP.py --clinvar clinvar_19_20.vcf.gz --dbSNP dbSNP_common_19_20.vcf.gz --output tmp.txt
sort tmp.txt | uniq > new.txt
table=/Work/Groups/bisonex/data-alexis/RefSeq/refseq_to_number_only_consensual.txt
python clinvar_sbSNP_old.py --clinvar clinvar_19_20_old.vcf.gz --dbSNP dbSNP_common_19_20_old.vcf.gz --output tmp_old.txt --chrm_name_table $table
sort tmp_old.txt | uniq > old.txt
wc -l old.txt new.txt
#+end_src
***** DONE Regarder la répartition des différences
CLOSED: [2022-12-11 Sun 16:29]
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
sort /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt | uniq > notpatho.new
sort /Work/Groups/bisonex/data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt | uniq > notpatho.old
comm -23 notpatho.new notpatho.old > nopatho.diff
#+end_src
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
bcftools query -i 'ID=@nopatho.diff' /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz -f '%CHROM\n' | sort | uniq -c
#+end_src
On a principalement des coordonnées non consensuelles (non "NC_", voir notes)
#+RESULTS:
: 2 NC_000002.12
: 18 NC_000003.12
: 2 NC_000004.12
: 2 NC_000005.10
: 14 NC_000006.12
: 6 NC_000007.14
: 2 NC_000009.12
: 1 NC_000010.11
: 6 NC_000014.9
: 1 NC_000015.10
: 3 NC_000016.10
: 3 NC_000017.11
: 1 NC_000019.10
: 1 NC_000020.11
: 1 NC_000021.9
: 2 NC_000022.11
: 16018 NT_113793.3
: 17010 NT_113796.3
: 14 NT_113891.3
: 1 NT_167244.2
: 13 NT_167245.2
: 2 NT_167246.2
: 13 NT_167247.2
: 7 NT_167248.2
: 14 NT_167249.2
: 14857 NT_187361.1
: 92 NT_187367.1
: 1 NT_187369.1
: 13 NT_187381.1
: 54 NT_187383.1
: 6 NT_187499.1
: 46 NT_187502.1
: 13754 NT_187513.1
: 611 NT_187517.1
: 1 NT_187520.1
: 1 NT_187524.1
: 249 NT_187526.1
: 18 NT_187532.1
: 1 NT_187546.1
: 886 NT_187562.1
: 1 NT_187564.1
: 346 NT_187576.1
: 13 NT_187600.1
: 5 NT_187601.1
: 494 NT_187606.1
: 1 NT_187607.1
: 12 NT_187613.1
: 307 NT_187614.1
: 1 NT_187625.1
: 445 NT_187633.1
: 43 NT_187648.1
: 18 NT_187649.1
: 1 NT_187652.1
: 512 NT_187661.1
: 18 NT_187678.1
: 49 NT_187681.1
: 1 NT_187682.1
: 18 NT_187688.1
: 12 NT_187689.1
: 18 NT_187690.1
: 18 NT_187691.1
: 404 NT_187693.1
: 2 NW_003315952.3
: 1 NW_003315970.2
: 203 NW_003571054.1
: 322 NW_003571055.2
: 16 NW_003571056.2
: 16 NW_003571057.2
: 16 NW_003571058.2
: 16 NW_003571059.2
: 16 NW_003571060.1
: 213 NW_003571061.2
: 2 NW_009646201.1
: 322 NW_009646205.1
: 321 NW_009646206.1
: 371 NW_012132914.1
: 1 NW_012132915.1
: 13 NW_012132918.1
: 2 NW_013171801.1
: 1 NW_013171807.1
: 49 NW_015148966.1
: 14 NW_015495298.1
: 2 NW_015495299.1
: 1 NW_016107298.1
: 4 NW_017363813.1
: 2 NW_017852933.1
: 1 NW_018654722.1
: 38 NW_021160001.1
: 1 NW_021160003.1
: 1 NW_021160007.1
: 7 NW_021160017.1
***** DONE Regarder la différence avec la version sans les sites non consensuels: ok !
CLOSED: [2022-12-11 Sun 20:11]
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
sort /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt | uniq > notpatho.new
sort /Work/Groups/bisonex/data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt | uniq > notpatho.old
comm -13 notpatho.new notpatho.old > notpatho.diff
wc -l
#+end_src
#+RESULTS:
: 528 notpatho.diff
Il manque 528 variants
rs1057520103
#+begin_src sh :dir /ssh:meso:/Work/Users/apraga/bisonex/tests/debug-commonsnp
PATH=$PATH:$HOME/.nix-profile/bin
bcftools query -i 'ID=@notpatho.diff' /Work/Groups/bisonex/data/dbSNP/GRCh38.p13/dbSNP_common.vcf.gz -f '%CHROM\n' | sort | uniq -c
#+end_src
#+RESULTS:
: 528 NC_012920.1
Donc la nouvelle version fonctionne mieux !
ON vérifie bien qu'ils sont dans l'ancienne version et la nouvelle:
$ grep -w -f notpatho.diff /Work/Groups/bisonex/data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt | wc -l
528
$ grep -w -f notpatho.diff /Work/Groups/bisonex/data/d
bSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt
#+end_src
**** DONE Supprimer les sites non consensuels
CLOSED: [2022-12-11 Sun 19:51]
**** DONE Rajouter les mitochondries (vu avec Paul)
CLOSED: [2022-12-13 Tue 17:26]
Ok avec notre version générée. Sur le J: 21155134
$ wc -l dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt
21155065 dbSNP/GRCh38.p13/ID_of_common_snp_not_clinvar_patho.txt
$ wc -l ../data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
21155065 ../data-alexis/dbSNP/ID_of_common_snp_not_clinvar_patho.txt
La différence vient probablement d'une vieille version de clinvar
*** TODO 63003856_S135
**** Notes
Actuellement
| Étapes | Production | Test |
|------------------+------------+-----------|
| bwa mem | 128077207 | 128077211 |
| haplotype caller | 1506894 | 1631935 |
| filterdepth | 82033 | 82054 |
#
**** STRT Version de non-régression: sans télécharger les bases de données
**** TODO Comparer les versions (log)
***** TODO bases de données
****** Genome de référénce:
fna ok
#+begin_src sh :dir /ssh:meso:/Work/Groups/bisonex/
sha256sum data-alexis-reference/genome/GRCh38_latest_genomic.fna
sha256sum data/genome/GRCh38.p13/genomeRef.fna
#+end_src
* Biblio
Comparaison WDL, Cromwell, nextflow
https://www.nature.com/articles/s41598-021-99288-8
Nextflow = bon compromis ?
* Idées
** Validation analytique
*** Génération de reads avec variants connus
Comparaison de génération ADN (2019)
https://academic.oup.com/bfg/article/19/1/49/5680294
**** SimuSCop
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
*** Utiliser données GCAT et uploader le notre ?
https://www.nature.com/articles/ncomms7275
*** Genome in a bottle : NA12878 + autres
2 versions :
1. Depuis un fastq correspondant à Illumina https://github.com/genome-in-a-bottle/giab_data_indexes
puis on compare le VCF avec les "high confidence" (Article : https://www.nature.com/articles/s41587-019-0054-x
)
- 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/data/NA12878/analysis/Illumina_PlatinumGenomes_NA12877_NA12878_09162015/hg38/2.0.1/NA12878/
2. On séquence directement NA12878
* 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()
** Divers
*** DONE Vérifier nombre de reads fastq - bam
CLOSED: [2022-10-09 Sun 22:31]
* Améliorations
** TODO Quality score recalibration avec un ensemble de fichier
Voir GATK best practice
** KILL Utiliser T-to-T comme références
CLOSED: [2023-01-01 Sun 21:35]
Semble compliqué avec les nouvelles bases de données
** TODO Macro excel
** TODO Utiliser le XML de clinvar
Extraction sous VCF possible avec
https://github.com/SeqOne/clinvcf
** Annotation
Liste complète
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9252745/
*** TODO Utilise une version allégée de GnomAD (une seule colonne)
*** TODO Digenisme (cf nomenclature omim)
C’est dans le nom de la maladie
* HOLD Implémenter d’autres pipeline
Voir https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-04407-x
** KILL GATK
CLOSED: [2022 -11-11 Fri 20:01]
https://broadinstitute.github.io/warp/docs/Pipelines/Exome_Germline_Single_Sample_Pipeline/README
A priori, respecte les bonnes pratiques
** KILL Essayer snmake avec bonne pratiques
https://github.com/snakemake-workflows/dna-seq-gatk-variant-calling/blob/main/.github/workflows/main.yml
Installer Mamba (micromamba ne fonctionne pas sous nix)
Ne fonctionne pas sous WSL2... MultiQC n’est pas assez à jour
Problèmes de versions...
** KILL Sarek
CLOSED: [2022-12-11 Sun 11:09]
*** Dépendences
**** Nix
#+begin_src sh
nix profile install nixpkgs#mosdepth nixpkgs#python3
nix-shell -p python310Packages.pyyaml --run "nextflow run nf-core/sarek -profile test --executor slurm --queue smp --outdir test -resume"
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()
***** KILL derivation nix pour profile complet
CLOSED: [2022-12-11 Sun 11:09]
**** KILL Sans nix
CLOSED: [2022-09-24 Sat 10:20]
On utilise conda
#+begin_src sh
module unload nix
module load anaconda3@2021.05/gcc-12.1.0
module load nextflow@22.04.0/gcc-12.1.0
module load openjdk@11.0.14.1_1/gcc-12.1.0
nextflow run nf-core/sarek -profile conda,test --executor slurm --queue smp --outdir test -resume
*** DONE simuscop
CLOSED: [2022-12-30 Fri 22:31]
** TODO MAJ avec picard
Normalement, GATK inclut picard mais la dernière version utilise picard pour certains outils
https://gatk.broadinstitute.org/hc/en-us/articles/9570266920219--Tool-Documentation-Index
A compléter après validation
*** TODO markduplicates
La dernière version dans la documentation utilise picard !!
** TODO Bwa-mem2 au lieu de bwa mem
https://github.com/bwa-mem2/bwa-mem2
** TODO Parallélisation haplotypecaller
spark est en beta, ne pas utiliser
parallélisation du pauvre : se restreindre à un chromosome avec -L et paralléliser sur le nombre de chromosome
** KILL CRAM au lieu de SAM ?
CLOSED: [2022-12-30 Fri 20:38]
Version compressée de bam mais :
Essai 1: erreurs de permissions, corrigé en relancant le programme
All GATK tools that take in mapped read data expect a BAM file as primary format. Some support the CRAM format, but we have observed performance issues when working directly from CRAM files, so in our own work we convert CRAM to BAM first, and we only use CRAM for archival purposes
Failed to create Conda environment
command: conda create --mkdir --yes --quiet --prefix /Work/Users/apraga/test-sarek/work/conda/env-2d53b1db50de676670cf1a91ef0cf6db bioconda::tabix=1.11
status : 1
message:
NotWritableError: The current user does not have write permissions to a required path.
path: /Home/Users/apraga/.conda/pkgs/urls.txt
uid: 1696
gid: 513
If you feel that permissions on this path are set incorrectly, you can manually
change them by executing
$ sudo chown 1696:513 /Home/Users/apraga/.conda/pkgs/urls.txt
Source: https://gatk.broadinstitute.org/hc/en-us/articles/360035890791-SAM-or-BAM-or-CRAM-Mapped-sequence-data-formats
** KILL Utiliser T-to-T comme références
CLOSED: [2023-01-01 Sun 21:35]
** TODO Preprocessing avec nextflow
*** TODO Map to reference
**** TODO Sample ID dans header
/Work/Users/apraga/bisonex/out/63003856_S135/preprocessing/baserecalibrator
*** TODO 63003856_S135
**** Notes
Actuellement
| Étapes | Production | Test |
|------------------+------------+-----------|
| bwa mem | 128077207 | 128077211 |
| haplotype caller | 1506894 | 1631935 |
| filterdepth | 82033 | 82054 |
#
**** STRT Version de non-régression: sans télécharger les bases de données
**** TODO Comparer les versions (log)
***** TODO bases de données
****** Genome de référénce:
fna ok
#+begin_src sh :dir /ssh:meso:/Work/Groups/bisonex/
sha256sum data-alexis-reference/genome/GRCh38_latest_genomic.fna
sha256sum data/genome/GRCh38.p13/genomeRef.fna
Corrigé avec
#+begin_src sh
chown 1696:513 /Home/Users/apraga/.conda/pkgs/urls.txt
#+end_src
Mais problème de proxy
*** KILL Dérivation nix pour modules python
CLOSED: [2022-12-11 Sun 11:09]
*** KILL Lancer sarek en mode test
CLOSED: [2022-12-11 Sun 11:09]
#+begin_src sh
nix-shell -p python310Packages.pyyaml --run "nextflow run nf-core/sarek -profile test --executor slurm --queue smp --outdir test -resume"
*** KILL Lancer sarek sur données allégées
CLOSED: [2022-12-11 Sun 11:09]