misc/wgsim.c

   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 }