|
0
|
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
|
|
|
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
|
|
|
99 void seq_set_block_size(int size)
|
|
|
100 {
|
|
|
101 SEQ_BLOCK_SIZE = size;
|
|
|
102 }
|
|
|
103
|
|
|
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
|
|
|
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
|
|
|
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 }
|
|
|
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
|
|
|
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
|
|
|
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
|
|
|
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 }
|
