1 /* The MIT License
2
3 Copyright (c) 2008 Genome Research Ltd (GRL).
4
5 Permission is hereby granted, free of charge, to any person obtaining
6 a copy of this software and associated documentation files (the
7 "Software"), to deal in the Software without restriction, including
8 without limitation the rights to use, copy, modify, merge, publish,
9 distribute, sublicense, and/or sell copies of the Software, and to
10 permit persons to whom the Software is furnished to do so, subject to
11 the following conditions:
12
13 The above copyright notice and this permission notice shall be
14 included in all copies or substantial portions of the Software.
15
16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
20 BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
21 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
22 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
23 SOFTWARE.
24 */
25
26 /* Contact: Heng Li <lh3@sanger.ac.uk> */
27
28 /* This program is separated from maq's read simulator with Colin
29 * Hercus' modification to allow longer indels. Colin is the chief
30 * developer of novoalign. */
31
32 #include <stdlib.h>
33 #include <math.h>
34 #include <time.h>
35 #include <assert.h>
36 #include <stdio.h>
37 #include <unistd.h>
38 #include <stdint.h>
39 #include <ctype.h>
40 #include <string.h>
41
42 #define PACKAGE_VERSION "0.2.3"
43
44 const uint8_t nst_nt4_table[256] = {
45 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
46 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
47 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5 /*'-'*/, 4, 4,
48 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
49 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
50 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
51 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
52 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
53 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
54 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
55 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
56 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
57 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
58 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
59 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
60 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
61 };
62
63 const int nst_color_space_table[] = { 4, 0, 0, 1, 0, 2, 3, 4, 0, 3, 2, 4, 1, 4, 4, 4};
64
65 /* Simple normal random number generator, copied from genran.c */
66
ran_normal()67 double ran_normal()
68 {
69 static int iset = 0;
70 static double gset;
71 double fac, rsq, v1, v2;
72 if (iset == 0) {
73 do {
74 v1 = 2.0 * drand48() - 1.0;
75 v2 = 2.0 * drand48() - 1.0;
76 rsq = v1 * v1 + v2 * v2;
77 } while (rsq >= 1.0 || rsq == 0.0);
78 fac = sqrt(-2.0 * log(rsq) / rsq);
79 gset = v1 * fac;
80 iset = 1;
81 return v2 * fac;
82 } else {
83 iset = 0;
84 return gset;
85 }
86 }
87
88 /* FASTA parser, copied from seq.c */
89
90 typedef struct {
91 int l, m; /* length and maximum buffer size */
92 unsigned char *s; /* sequence */
93 } seq_t;
94
95 #define INIT_SEQ(seq) (seq).s = 0; (seq).l = (seq).m = 0
96
97 static int SEQ_BLOCK_SIZE = 512;
98
seq_set_block_size(int size)99 void seq_set_block_size(int size)
100 {
101 SEQ_BLOCK_SIZE = size;
102 }
103
seq_read_fasta(FILE * fp,seq_t * seq,char * locus,char * comment)104 int seq_read_fasta(FILE *fp, seq_t *seq, char *locus, char *comment)
105 {
106 int c, l, max;
107 char *p;
108
109 c = 0;
110 while (!feof(fp) && fgetc(fp) != '>');
111 if (feof(fp)) return -1;
112 p = locus;
113 while (!feof(fp) && (c = fgetc(fp)) != ' ' && c != '\t' && c != '\n')
114 if (c != '\r') *p++ = c;
115 *p = '\0';
116 if (comment) {
117 p = comment;
118 if (c != '\n') {
119 while (!feof(fp) && ((c = fgetc(fp)) == ' ' || c == '\t'));
120 if (c != '\n') {
121 *p++ = c;
122 while (!feof(fp) && (c = fgetc(fp)) != '\n')
123 if (c != '\r') *p++ = c;
124 }
125 }
126 *p = '\0';
127 } else if (c != '\n') while (!feof(fp) && fgetc(fp) != '\n');
128 l = 0; max = seq->m;
129 while (!feof(fp) && (c = fgetc(fp)) != '>') {
130 if (isalpha(c) || c == '-' || c == '.') {
131 if (l + 1 >= max) {
132 max += SEQ_BLOCK_SIZE;
133 seq->s = (unsigned char*)realloc(seq->s, sizeof(char) * max);
134 }
135 seq->s[l++] = (unsigned char)c;
136 }
137 }
138 if (c == '>') ungetc(c,fp);
139 seq->s[l] = 0;
140 seq->m = max; seq->l = l;
141 return l;
142 }
143
144 /* Error-checking open, copied from utils.c */
145
146 #define xopen(fn, mode) err_xopen_core(__func__, fn, mode)
147
err_xopen_core(const char * func,const char * fn,const char * mode)148 FILE *err_xopen_core(const char *func, const char *fn, const char *mode)
149 {
150 FILE *fp = 0;
151 if (strcmp(fn, "-") == 0)
152 return (strstr(mode, "r"))? stdin : stdout;
153 if ((fp = fopen(fn, mode)) == 0) {
154 fprintf(stderr, "[%s] fail to open file '%s'. Abort!\n", func, fn);
155 abort();
156 }
157 return fp;
158 }
159
160 /* wgsim */
161
162 enum muttype_t {NOCHANGE = 0, INSERT = 0x1000, SUBSTITUTE = 0xe000, DELETE = 0xf000};
163 typedef unsigned short mut_t;
164 static mut_t mutmsk = (mut_t)0xf000;
165
166 typedef struct {
167 int l, m; /* length and maximum buffer size */
168 mut_t *s; /* sequence */
169 } mutseq_t;
170
171 static double ERR_RATE = 0.02;
172 static double MUT_RATE = 0.001;
173 static double INDEL_FRAC = 0.1;
174 static double INDEL_EXTEND = 0.3;
175 static int IS_SOLID = 0;
176 static int SHOW_MM_INFO = 1;
177
maq_mut_diref(const seq_t * seq,int is_hap,mutseq_t * hap1,mutseq_t * hap2)178 void maq_mut_diref(const seq_t *seq, int is_hap, mutseq_t *hap1, mutseq_t *hap2)
179 {
180 int i, deleting = 0;
181 mutseq_t *ret[2];
182
183 ret[0] = hap1; ret[1] = hap2;
184 ret[0]->l = seq->l; ret[1]->l = seq->l;
185 ret[0]->m = seq->m; ret[1]->m = seq->m;
186 ret[0]->s = (mut_t *)calloc(seq->m, sizeof(mut_t));
187 ret[1]->s = (mut_t *)calloc(seq->m, sizeof(mut_t));
188 for (i = 0; i != seq->l; ++i) {
189 int c;
190 c = ret[0]->s[i] = ret[1]->s[i] = (mut_t)nst_nt4_table[(int)seq->s[i]];
191 if (deleting) {
192 if (drand48() < INDEL_EXTEND) {
193 if (deleting & 1) ret[0]->s[i] |= DELETE;
194 if (deleting & 2) ret[1]->s[i] |= DELETE;
195 continue;
196 } else deleting = 0;
197 }
198 if (c < 4 && drand48() < MUT_RATE) { // mutation
199 if (drand48() >= INDEL_FRAC) { // substitution
200 double r = drand48();
201 c = (c + (int)(r * 3.0 + 1)) & 3;
202 if (is_hap || drand48() < 0.333333) { // hom
203 ret[0]->s[i] = ret[1]->s[i] = SUBSTITUTE|c;
204 } else { // het
205 ret[drand48()<0.5?0:1]->s[i] = SUBSTITUTE|c;
206 }
207 } else { // indel
208 if (drand48() < 0.5) { // deletion
209 if (is_hap || drand48() < 0.333333) { // hom-del
210 ret[0]->s[i] = ret[1]->s[i] = DELETE;
211 deleting = 3;
212 } else { // het-del
213 deleting = drand48()<0.5?1:2;
214 ret[deleting-1]->s[i] = DELETE;
215 }
216 } else { // insertion
217 int num_ins = 0, ins = 0;
218 do {
219 num_ins++;
220 ins = (ins << 2) | (int)(drand48() * 4.0);
221 } while (num_ins < 4 && drand48() < INDEL_EXTEND);
222
223 if (is_hap || drand48() < 0.333333) { // hom-ins
224 ret[0]->s[i] = ret[1]->s[i] = (num_ins << 12) | (ins << 4) | c;
225 } else { // het-ins
226 ret[drand48()<0.5?0:1]->s[i] = (num_ins << 12) | (ins << 4) | c;
227 }
228 }
229 }
230 }
231 }
232 }
maq_print_mutref(const char * name,const seq_t * seq,mutseq_t * hap1,mutseq_t * hap2)233 void maq_print_mutref(const char *name, const seq_t *seq, mutseq_t *hap1, mutseq_t *hap2)
234 {
235 int i;
236 for (i = 0; i != seq->l; ++i) {
237 int c[3];
238 c[0] = nst_nt4_table[(int)seq->s[i]];
239 c[1] = hap1->s[i]; c[2] = hap2->s[i];
240 if (c[0] >= 4) continue;
241 if ((c[1] & mutmsk) != NOCHANGE || (c[2] & mutmsk) != NOCHANGE) {
242 printf("%s\t%d\t", name, i+1);
243 if (c[1] == c[2]) { // hom
244 if ((c[1]&mutmsk) == SUBSTITUTE) { // substitution
245 printf("%c\t%c\t-\n", "ACGTN"[c[0]], "ACGTN"[c[1]&0xf]);
246 } else if ((c[1]&mutmsk) == DELETE) { // del
247 printf("%c\t-\t-\n", "ACGTN"[c[0]]);
248 } else if (((c[1] & mutmsk) >> 12) <= 5) { // ins
249 printf("-\t");
250 int n = (c[1]&mutmsk) >> 12, ins = c[1] >> 4;
251 while(n > 0) {
252 putchar("ACGTN"[ins & 0x3]);
253 n--;
254 }
255 printf("\t-\n");
256 } else assert(0);
257 } else { // het
258 if ((c[1]&mutmsk) == SUBSTITUTE || (c[2]&mutmsk) == SUBSTITUTE) { // substitution
259 printf("%c\t%c\t+\n", "ACGTN"[c[0]], "XACMGRSVTWYHKDBN"[1<<(c[1]&0x3)|1<<(c[2]&0x3)]);
260 } else if ((c[1]&mutmsk) == DELETE) {
261 printf("%c\t-\t+\n", "ACGTN"[c[0]]);
262 } else if ((c[2]&mutmsk) == DELETE) {
263 printf("%c\t-\t+\n", "ACGTN"[c[0]]);
264 } else if (((c[1] & mutmsk) >> 12) <= 4) { // ins1
265 printf("-\t");
266 int n = (c[1]&mutmsk) >> 12, ins = c[1] >> 4;
267 while (n > 0) {
268 putchar("ACGTN"[ins & 0x3]);
269 n--;
270 }
271 printf("\t+\n");
272 } else if (((c[2] & mutmsk) >> 12) <= 5) { // ins2
273 printf("-\t");
274 int n = (c[2]&mutmsk) >> 12, ins = c[2] >> 4;
275 while (n > 0) {
276 putchar("ACGTN"[ins & 0x3]);
277 ins >>= 2;
278 n--;
279 }
280 printf("\t+\n");
281 } else assert(0);
282 }
283 }
284 }
285 }
286
wgsim_core(FILE * fpout1,FILE * fpout2,FILE * fp_fa,int is_hap,uint64_t N,int dist,int std_dev,int size_l,int size_r)287 void wgsim_core(FILE *fpout1, FILE *fpout2, FILE *fp_fa, int is_hap, uint64_t N, int dist, int std_dev, int size_l, int size_r)
288 {
289 seq_t seq;
290 mutseq_t rseq[2];
291 uint64_t tot_len, ii;
292 int i, l, n_ref;
293 char name[256], *qstr;
294 int size[2], Q;
295 uint8_t *tmp_seq[2];
296 mut_t *target;
297
298 INIT_SEQ(seq);
299 srand48(time(0));
300 seq_set_block_size(0x1000000);
301 l = size_l > size_r? size_l : size_r;
302 qstr = (char*)calloc(l+1, 1);
303 tmp_seq[0] = (uint8_t*)calloc(l+2, 1);
304 tmp_seq[1] = (uint8_t*)calloc(l+2, 1);
305 size[0] = size_l; size[1] = size_r;
306
307 Q = (ERR_RATE == 0.0)? 'I' : (int)(-10.0 * log(ERR_RATE) / log(10.0) + 0.499) + 33;
308
309 tot_len = n_ref = 0;
310 while ((l = seq_read_fasta(fp_fa, &seq, name, 0)) >= 0) {
311 tot_len += l;
312 ++n_ref;
313 }
314 fprintf(stderr, "[wgsim_core] %d sequences, total length: %llu\n", n_ref, (long long)tot_len);
315 rewind(fp_fa);
316
317 while ((l = seq_read_fasta(fp_fa, &seq, name, 0)) >= 0) {
318 uint64_t n_pairs = (uint64_t)((long double)l / tot_len * N + 0.5);
319 if (l < dist + 3 * std_dev) {
320 fprintf(stderr, "[wgsim_core] kkip sequence '%s' as it is shorter than %d!\n", name, dist + 3 * std_dev);
321 continue;
322 }
323
324 // generate mutations and print them out
325 maq_mut_diref(&seq, is_hap, rseq, rseq+1);
326 maq_print_mutref(name, &seq, rseq, rseq+1);
327
328 for (ii = 0; ii != n_pairs; ++ii) { // the core loop
329 double ran;
330 int d, pos, s[2], is_flip = 0;
331 int n_sub[2], n_indel[2], n_err[2], ext_coor[2], j, k;
332 FILE *fpo[2];
333
334 do { // avoid boundary failure
335 ran = ran_normal();
336 ran = ran * std_dev + dist;
337 d = (int)(ran + 0.5);
338 pos = (int)((l - d + 1) * drand48());
339 } while (pos < 0 || pos >= seq.l || pos + d - 1 >= seq.l);
340
341 // flip or not
342 if (drand48() < 0.5) {
343 fpo[0] = fpout1; fpo[1] = fpout2;
344 s[0] = size[0]; s[1] = size[1];
345 } else {
346 fpo[1] = fpout1; fpo[0] = fpout2;
347 s[1] = size[0]; s[0] = size[1];
348 is_flip = 1;
349 }
350
351 // generate the read sequences
352 target = rseq[drand48()<0.5?0:1].s; // haplotype from which the reads are generated
353 n_sub[0] = n_sub[1] = n_indel[0] = n_indel[1] = n_err[0] = n_err[1] = 0;
354
355 #define __gen_read(x, start, iter) do { \
356 for (i = (start), k = 0, ext_coor[x] = -10; i >= 0 && i < seq.l && k < s[x]; iter) { \
357 int c = target[i], mut_type = c & mutmsk; \
358 if (ext_coor[x] < 0) { \
359 if (mut_type != NOCHANGE && mut_type != SUBSTITUTE) continue; \
360 ext_coor[x] = i; \
361 } \
362 if (mut_type == DELETE) ++n_indel[x]; \
363 else if (mut_type == NOCHANGE || mut_type == SUBSTITUTE) { \
364 tmp_seq[x][k++] = c & 0xf; \
365 if (mut_type == SUBSTITUTE) ++n_sub[x]; \
366 } else { \
367 int n, ins; \
368 ++n_indel[x]; \
369 tmp_seq[x][k++] = c & 0xf; \
370 for (n = mut_type>>12, ins = c>>4; n > 0 && k < s[x]; --n, ins >>= 2) \
371 tmp_seq[x][k++] = ins & 0x3; \
372 } \
373 } \
374 if (k != s[x]) ext_coor[x] = -10; \
375 } while (0)
376
377 if (!IS_SOLID) {
378 __gen_read(0, pos, ++i);
379 __gen_read(1, pos + d - 1, --i);
380 for (k = 0; k < s[1]; ++k) tmp_seq[1][k] = tmp_seq[1][k] < 4? 3 - tmp_seq[1][k] : 4; // complement
381 } else {
382 int c1, c2, c;
383 ++s[0]; ++s[1]; // temporarily increase read length by 1
384 if (is_flip) { // RR pair
385 __gen_read(0, pos + s[0], --i);
386 __gen_read(1, pos + d - 1, --i);
387 } else { // FF pair
388 __gen_read(0, pos, ++i);
389 __gen_read(1, pos + d - 1 - s[1], ++i);
390 ++ext_coor[0]; ++ext_coor[1];
391 }
392 // change to color sequence: (0,1,2,3) -> (A,C,G,T)
393 for (j = 0; j < 2; ++j) {
394 c1 = tmp_seq[j][0];
395 for (i = 1; i < s[j]; ++i) {
396 c2 = tmp_seq[j][i];
397 c = (c1 >= 4 || c2 >= 4)? 4 : nst_color_space_table[(1<<c1)|(1<<c2)];
398 tmp_seq[j][i-1] = c;
399 c1 = c2;
400 }
401 }
402 --s[0]; --s[1]; // change back
403 }
404 if (ext_coor[0] < 0 || ext_coor[1] < 0) { // fail to generate the read(s)
405 --ii;
406 continue;
407 }
408
409 // generate sequencing errors
410 for (j = 0; j < 2; ++j) {
411 for (i = 0; i < s[j]; ++i) {
412 int c = tmp_seq[j][i];
413 if (c >= 4) c = 4; // actually c should be never larger than 4 if everything is correct
414 else if (drand48() < ERR_RATE) {
415 c = (c + (int)(drand48() * 3.0 + 1)) & 3;
416 ++n_err[j];
417 }
418 tmp_seq[j][i] = c;
419 }
420 }
421
422 // print
423 for (j = 0; j < 2; ++j) {
424 for (i = 0; i < s[j]; ++i) qstr[i] = Q;
425 qstr[i] = 0;
426 if (SHOW_MM_INFO) {
427 fprintf(fpo[j], "@%s_%u_%u_%d:%d:%d_%d:%d:%d_%llx/%d\n", name, ext_coor[0]+1, ext_coor[1]+1,
428 n_err[0], n_sub[0], n_indel[0], n_err[1], n_sub[1], n_indel[1],
429 (long long)ii, j==0? is_flip+1 : 2-is_flip);
430 } else {
431 fprintf(fpo[j], "@%s_%u_%u_%llx/%d %d:%d:%d_%d:%d:%d\n", name, ext_coor[0]+1, ext_coor[1]+1,
432 (long long)ii, j==0? is_flip+1 : 2-is_flip,
433 n_err[0], n_sub[0], n_indel[0], n_err[1], n_sub[1], n_indel[1]);
434 }
435 for (i = 0; i < s[j]; ++i)
436 fputc("ACGTN"[(int)tmp_seq[j][i]], fpo[j]);
437 fprintf(fpo[j], "\n+\n%s\n", qstr);
438 }
439 }
440 free(rseq[0].s); free(rseq[1].s);
441 }
442 free(seq.s); free(qstr);
443 free(tmp_seq[0]); free(tmp_seq[1]);
444 }
445
simu_usage()446 static int simu_usage()
447 {
448 fprintf(stderr, "\n");
449 fprintf(stderr, "Program: wgsim (short read simulator)\n");
450 fprintf(stderr, "Version: %s\n", PACKAGE_VERSION);
451 fprintf(stderr, "Contact: Heng Li <lh3@sanger.ac.uk>\n\n");
452 fprintf(stderr, "Usage: wgsim [options] <in.ref.fa> <out.read1.fq> <out.read2.fq>\n\n");
453 fprintf(stderr, "Options: -e FLOAT base error rate [%.3f]\n", ERR_RATE);
454 fprintf(stderr, " -d INT outer distance between the two ends [500]\n");
455 fprintf(stderr, " -s INT standard deviation [50]\n");
456 fprintf(stderr, " -N INT number of read pairs [1000000]\n");
457 fprintf(stderr, " -1 INT length of the first read [70]\n");
458 fprintf(stderr, " -2 INT length of the second read [70]\n");
459 fprintf(stderr, " -r FLOAT rate of mutations [%.4f]\n", MUT_RATE);
460 fprintf(stderr, " -R FLOAT fraction of indels [%.2f]\n", INDEL_FRAC);
461 fprintf(stderr, " -X FLOAT probability an indel is extended [%.2f]\n", INDEL_EXTEND);
462 fprintf(stderr, " -c generate reads in color space (SOLiD reads)\n");
463 fprintf(stderr, " -C show mismatch info in comment rather than read name\n");
464 fprintf(stderr, " -h haplotype mode\n");
465 fprintf(stderr, "\n");
466 fprintf(stderr, "Note: For SOLiD reads, the first read is F3 and the second is R3.\n\n");
467 return 1;
468 }
469
main(int argc,char * argv[])470 int main(int argc, char *argv[])
471 {
472 int64_t N;
473 int dist, std_dev, c, size_l, size_r, is_hap = 0;
474 FILE *fpout1, *fpout2, *fp_fa;
475
476 N = 1000000; dist = 500; std_dev = 50;
477 size_l = size_r = 70;
478 while ((c = getopt(argc, argv, "e:d:s:N:1:2:r:R:hX:cC")) >= 0) {
479 switch (c) {
480 case 'd': dist = atoi(optarg); break;
481 case 's': std_dev = atoi(optarg); break;
482 case 'N': N = atoi(optarg); break;
483 case '1': size_l = atoi(optarg); break;
484 case '2': size_r = atoi(optarg); break;
485 case 'e': ERR_RATE = atof(optarg); break;
486 case 'r': MUT_RATE = atof(optarg); break;
487 case 'R': INDEL_FRAC = atof(optarg); break;
488 case 'X': INDEL_EXTEND = atof(optarg); break;
489 case 'c': IS_SOLID = 1; break;
490 case 'C': SHOW_MM_INFO = 0; break;
491 case 'h': is_hap = 1; break;
492 }
493 }
494 if (argc - optind < 3) return simu_usage();
495 fp_fa = (strcmp(argv[optind+0], "-") == 0)? stdin : xopen(argv[optind+0], "r");
496 fpout1 = xopen(argv[optind+1], "w");
497 fpout2 = xopen(argv[optind+2], "w");
498 wgsim_core(fpout1, fpout2, fp_fa, is_hap, N, dist, std_dev, size_l, size_r);
499
500 fclose(fpout1); fclose(fpout2); fclose(fp_fa);
501 return 0;
502 }
503