1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
|
# Claire Lemaitre
# 16/10/2015
# Creates a small dataset to test MindTheGap find
# With :
# - several chromosomes
# - all types of variants = solo SNP, multi SNP, homozygous insertion fuzzy or not, heterozygous insertion fuzzy or not, homozygous insertion with a SNP close before or after the insertion, homozygous deletion
# Here, positions are 0-based
# 1. Generates random sequences
~/workspace/divers_scripts/gener_alea 1000 5
mv alea.seq init.fasta
init.fasta = reference.fasta
cp init.fasta reference.fasta
# 2. Put 13 SNPs
cp init.fasta initSnp.fasta
liste des SNPs : (0-based)
Seq0 100 T -> C
Seq0 815 C -> A
Seq1 205 G -> C # 2 snps proches
Seq1 218 T -> A
Seq2 319 T -> C # 2 snp proches (k=31)
Seq2 343 C -> A
Seq2 378 G -> C
Seq3 255 A -> T
Seq3 510 C -> A
Seq3 765 G -> A
Seq4 256 C -> T
Seq4 511 A -> G
Seq4 840 C -> T
# 3. Put 3 deletions (homozygous)
cp initSnp.fasta initSnpAndDel.fasta
dans initSnpAndDel.fasta :
Seq0 296 -> 410 clean seq = TAGCTTGAGAGTGCGTATCTCACCGATCCCCTGGCTATGCTCCGCGATTCACTAGTAGTTTCACGCCGACAGAGCGAAACCGTGATAGGTCATCATGCCGGTCTGCAGTCACGT
Seq1 739 -> 847 fuzzy repeat 2 CT debut del seq = CTGTTGGGAAGGAATTGCAATACTCTCCGAACCAGCTTAGGGCCCCCCGCCGCCGCAATTCGAGCGTTATGCCCGGAGCATTTGCACGATGCCATTAAACTATATCAA
Seq4 884 -> 928 fuzzy repeat 2, 1 lettre de chaque côté : A et T seq = AGGGACCTAGACGCAACAGTAACCGCCTCGGAGTAAGCCCTGGT
# 4. Put 5 homozygous insertions
~/workspace/divers_scripts/gener_alea 3000 1 : pour piocher des séquences à insérer
cp initSnpAndDel.fasta allele1.fasta
dans allele1.fasta :
Seq2 834 fuzzy rep 1 (T) seq=TGCACGCTGCAGGATTGGAACCACAATGTACGCCGATCCAAGCAGTAGTGGTTCATTGTATAAGTATCCTCCCTTGATTGGTCGAATATTAGGCATGCCCCGGGAGCATGTGGGCTCGAGCCACGGAGAGCAACTAATCGCGCATAAAACAAATACCTCATGGTTTTTGTGCGGAAAACCGTTGGGTGGACCATCAGCGGTTGTGAT
Seq2 535 clean seq= TAACGTTCGCTGAACATCGACTCCGGTGACGACATACGATTCAAGAAGAGAGTGACTCTGTAGGATAACATCCCGCAACGCCTAATCCATCCAGCCTGGCACCATGTATAAAGGGCGTCAGGTATGTTAACGAGACTATT
Seq4 820 clean before SNP seq = GGATGGCGGCCGGAGAGCGCTGCAATCGCATGGCTCGGGA
Seq4 600 fuzzy rep 3 (1+2 GG + A) seq=GGTGTATTCCTGGGTTGAGTGGCAGGTTTCTCTTAATTCTTCCCTAAGTAGCTCCGA
Seq3 780 clean after SNP seq = TTTGCAGCACTAGCCGTTCCTTGACATCTGCGGCCAACTTGTGCCTGAACCTGGAGTTTCGACAGCGTGGCGCTCTGGCCTAGTTCTTCGCTGGCACCTGGAAGAGCCGCC
# 5. Put 3 heterozygous insertions
cp allele1.fasta allele2.fasta
dans allele2.fasta
Seq1 341 clean seq = ATGGTTTATAGAACCCGGGCGTTCATGTCCGTCAGAACGATCTTGGCACGGTAGCCCCTGGTCCAGAGAGCCAAGGTGACTCAGCCCCACGATGGTGGTCTAGAGCGAAATAACCCTCGCCGAGA
Seq4 348 fuzzy 2 CA seq= CAGTCTTAACCTTAAGACCGTTCATTGATAAAACTTGCTCACGCTCTAGATGGCGTGAAGCGAAACCTAGGAAAAAGTTTTGCAGATAATTAGATTATGCGCGATACTCCGCCGTGTGTT
Seq0 122 clean after a SNP seq= ATCTAAGCTGTGACCTTGTGGCCGAGGCGCTTTTCACGCCTACATTAACTCCTGGGAAGCTCTCTGCTCTAGTTTCAGTGCACATCTCCAGGTGAGCAACCCTGGCAAGCAGCCCCTTCCTGTAGAAATTACTTAGC
../build/MindTheGap find -in full_test/allele1.fasta,full_test/allele2.fasta -ref full_test/reference.fasta -abundance-min 1 -out full-test1
## Pb MTG : deletion position apres la deletion (premiere base deletée)
insertion HOM avant/apres un SNP, position décalée de 1 ??? vérifier les kmers breakpoints (seq4 821
insertion HET apres SNP non détectée, détéctable ?
2 backups d'ou viennent ces trous ??? vérifier les positions
MSNP seq2 319 raté : non mais affiché avec position 367 = 343 + 24 au lieu de 343 - 24
Backup position 305 = début MSNP en position 319 (right-kmer : position 305)
Backup position 289 = position 337 (right kmer n'existe pas dans ref. existe dans allele1 finit en position 319
Si on fait tourner que sur initSnp -> aucun soucis, idem sur initSnpAndDel, allele1 seul, allele2 seul : pas de backup. Si -snp-only : pas de backup. Vient de hétéro ? non -> si homo-only : 2 backup
# Ajout de print pour debugger : séquences de 01 : dans FindBreakpoints.hpp, lignes ~400
sequence Seq2
1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111000000000000000011000000000000000000000000000000000000011110000000000000000000000000000000111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111000000000000000000000000000000111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111000000000000000010000000000001111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
2 FP successifs (ces kmers n'existent pas dans les allele1 ou allele2.fasta) en position 305 et 306
=> ces FPs existent seulement avec -in full_test/allele1.fasta,full_test/allele2.fasta --> explique pourquoi tout se passe bien avec -in initSnp ou initSnpAndDel.fasta
=> avec -snp-only : ne sort pas de backup même s'il y en a (pas demandé)
=> position 367 : à cause de ces FPs, le SNP 343 est trouvé en partant de la droite, puis trouve un 0 à gauche, donc essaie un autre SNP (sans se souvenir qu'il va dépasser la taille du trou), donc trouve le SNP en position 319, MAIS reste le bug de position 343+24 au lieu de 343-24 : est-ce que ça avait été testé ?
Le premier backup est quand même renvoyé (la position est : position du premier 1 après le trou), car ne sait pas qu'il a finalement résolu le pb en partant de la droite.
Le deuxième backup : gap_stretch_size = 18446744073709551598 : car avec SNP 319 est allé plus loin que le trou, donc taille négative !
--> résolution pb de position pour SNPs proches (bug dans FindSNP.hpp : FindMultiSNPrev<span>::update() -> begin_pos-= nb_kmer_val au lieu de +=)
=> 2 solutions :
- soit faire attention dans MultipleSNP de ne pas dépasser la taille du trou : tant pis on ne trouve pas le SNP 319 du tout --> DONE
- soit ne pas renvoyer de backup du tout : NON cf. plus bas
pour full_test : ajouter le SNP 367, voir si même avec les FPs, on retrouve 3 SNPs sur 4 des multiples (du coup test si partir de la droite marche bien)
attention : les FPs peuvent dépendre des fonctions de hash, et ce test n'ura pas toujours les mêmes résus...
Nouveau jeu avec les FP après le muti-SNp --> pour tester multi-SNP avant FP
reference-bis.fasta seq2 292 A -> T
--> si autorise multiSNP à dépasser la taille : 2 SNP FP en trop --> option n°1 validée (et implémentée)
UPDATE 27/05/2016 :
Not detected : - Seq0 122 HET clean after SNP # not detected for the moment : not possible due to update after SNPs only touching gap-stretch-size ?
- Seq2 319 MSNP T -> C # 2 close snp, note : due to bloom-FP can be missed
Now reported positions are 1-based and points to the nucleotide of the variant (SNP), or just before (at the left) : deletions, insertions.
# 6. Generate reads :
l=5000
1000 reads -> c=20x
~/Bin/samtools-0.1.18/misc/wgsim -e 0.01 -d 200 -s 20 -N 500 -1 100 -2 100 -r 0 -R 0 allele1.fasta allele1_r1.fq allele1_r2.fq
~/Bin/samtools-0.1.18/misc/wgsim -e 0.01 -d 200 -s 20 -N 500 -1 100 -2 100 -r 0 -R 0 allele2.fasta allele2_r1.fq allele2_r2.fq
cat allele1_r1.fq allele2_r1.fq > reads_r1.fastq
cat allele1_r2.fq allele2_r2.fq > reads_r2.fastq
../build/bin/MindTheGap find -in full_test/reads_r1.fastq,full_test/reads_r2.fastq -ref full_test/reference.fasta -abundance-min auto -out test-reads
Avec les reads le SNP 320 est trouvé (pas de pb de bloom-FP)
cp reads_r1.fq ../../data/reads_r1.fastq
cp reads_r2.fq ../../data/reads_r2.fastq
cp reference.fasta ../../data/reference.fasta
# 7. Create Gold files for automated tests
../../build/bin/MindTheGap find -in ../../data/reads_r1.fastq,../../data/reads_r2.fastq -ref ../../data/reference.fasta -out gold -nb-cores 1 > gold_find.output
../../build/bin/MindTheGap fill -graph gold.h5 -bkpt gold.breakpoints -out gold -nb-cores 1 > gold_fill.output
# 08/02/2022 : rajoute des petites insertions 1-2 bp
ou ajout de seq5 et seq6 depuis : Projets/mindTheGap/test-small-indels (990 premiers nt de chr1 et chr2, 6/10 HOM, 4/10 HET seulement dans allele 1)
# change les param de simulation des reads : augmente couverture (+diminue tx erreurs)
~/Bin/samtools-0.1.18/misc/wgsim -e 0.001 -d 200 -s 20 -N 1000 -1 100 -2 100 -r 0 -R 0 allele1.fasta allele1_r1.fq allele1_r2.fq
~/Bin/samtools-0.1.18/misc/wgsim -e 0.001 -d 200 -s 20 -N 1000 -1 100 -2 100 -r 0 -R 0 allele2.fasta allele2_r1.fq allele2_r2.fq
cat allele1_r1.fq allele2_r1.fq > reads_r1.fastq
cat allele1_r2.fq allele2_r2.fq > reads_r2.fastq
ATTENTION : changements dans les résultats : othervariants.vcf : perd 2 snps (Seq1 206, 219), gagne 1 del (Seq0 297) mais en perd une autre (Seq1 740) ; ne change pas les résultats des grandes insertions
RQ : rate 2 petites insertions de taille 2 (Seq6 : pos 500 et 900)
# List of all mutations :
Seq0 101
Seq0 123
Seq0 816
Seq1 206
Seq1 219
Seq1 342
Seq1 740
Seq2 320
Seq2 344
Seq2 379
Seq2 535
Seq2 834
Seq3 256
Seq3 511
Seq3 766
Seq3 781
Seq4 257
Seq4 349
Seq4 512
Seq4 600
Seq4 841
Seq4 884
Seq4 821
|
