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