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 uint8_t nst_nt4_table[256] = {
  43         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  44         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  45         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 5 /*'-'*/, 4, 4,
  46         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  47         4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
  48         4, 4, 4, 4,  3, 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, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  52         4, 4, 4, 4,  4, 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 };
  60
  61 /* Simple normal random number generator, copied from genran.c */
  62
  63 double ran_normal()
  64 {
  65         static int iset = 0;
  66         static double gset;
  67         double fac, rsq, v1, v2;
  68         if (iset == 0) {
  69                 do {
  70                         v1 = 2.0 * drand48() - 1.0;
  71                         v2 = 2.0 * drand48() - 1.0;
  72                         rsq = v1 * v1 + v2 * v2;
  73                 } while (rsq >= 1.0 || rsq == 0.0);
  74                 fac = sqrt(-2.0 * log(rsq) / rsq);
  75                 gset = v1 * fac;
  76                 iset = 1;
  77                 return v2 * fac;
  78         } else {
  79                 iset = 0;
  80                 return gset;
  81         }
  82 }
  83
  84 /* FASTA parser, copied from seq.c */
  85
  86 typedef struct {
  87         int l, m; /* length and maximum buffer size */
  88         unsigned char *s; /* sequence */
  89 } seq_t;
  90
  91 #define INIT_SEQ(seq) (seq).s = 0; (seq).l = (seq).m = 0
  92
  93 static int SEQ_BLOCK_SIZE = 512;
  94
  95 void seq_set_block_size(int size)
  96 {
  97         SEQ_BLOCK_SIZE = size;
  98 }
  99
 100 int seq_read_fasta(FILE *fp, seq_t *seq, char *locus, char *comment)
 101 {
 102         int c, l, max;
 103         char *p;
 104
 105         c = 0;
 106         while (!feof(fp) && fgetc(fp) != '>');
 107         if (feof(fp)) return -1;
 108         p = locus;
 109         while (!feof(fp) && (c = fgetc(fp)) != ' ' && c != '\t' && c != '\n')
 110                 if (c != '\r') *p++ = c;
 111         *p = '\0';
 112         if (comment) {
 113                 p = comment;
 114                 if (c != '\n') {
 115                         while (!feof(fp) && ((c = fgetc(fp)) == ' ' || c == '\t'));
 116                         if (c != '\n') {
 117                                 *p++ = c;
 118                                 while (!feof(fp) && (c = fgetc(fp)) != '\n')
 119                                         if (c != '\r') *p++ = c;
 120                         }
 121                 }
 122                 *p = '\0';
 123         } else if (c != '\n') while (!feof(fp) && fgetc(fp) != '\n');
 124         l = 0; max = seq->m;
 125         while (!feof(fp) && (c = fgetc(fp)) != '>') {
 126                 if (isalpha(c) || c == '-' || c == '.') {
 127                         if (l + 1 >= max) {
 128                                 max += SEQ_BLOCK_SIZE;
 129                                 seq->s = (unsigned char*)realloc(seq->s, sizeof(char) * max);
 130                         }
 131                         seq->s[l++] = (unsigned char)c;
 132                 }
 133         }
 134         if (c == '>') ungetc(c,fp);
 135         seq->s[l] = 0;
 136         seq->m = max; seq->l = l;
 137         return l;
 138 }
 139
 140 /* Error-checking open, copied from utils.c */
 141
 142 #define xopen(fn, mode) err_xopen_core(__func__, fn, mode)
 143
 144 FILE *err_xopen_core(const char *func, const char *fn, const char *mode)
 145 {
 146         FILE *fp = 0;
 147         if (strcmp(fn, "-") == 0)
 148                 return (strstr(mode, "r"))? stdin : stdout;
 149         if ((fp = fopen(fn, mode)) == 0) {
 150                 fprintf(stderr, "[%s] fail to open file '%s'. Abort!\n", func, fn);
 151                 abort();
 152         }
 153         return fp;
 154 }
 155
 156 /* wgsim */
 157
 158 enum muttype_t {NOCHANGE = 0, INSERT = 0x1000, SUBSTITUTE = 0xe000, DELETE = 0xf000};
 159 typedef unsigned short mut_t;
 160 static mut_t mutmsk = (unsigned short)0xf000;
 161
 162 typedef struct
 163 {
 164         int l, m; /* length and maximum buffer size */
 165         mut_t *s; /* sequence */
 166 } mutseq_t;
 167
 168 static double ERR_RATE = 0.02;
 169 static double MUT_RATE = 0.001;
 170 static double INDEL_FRAC = 0.1;
 171 static double INDEL_EXTEND = 0.3;
 172
 173 void maq_mut_diref(const seq_t *seq, int is_hap, mutseq_t *hap1, mutseq_t *hap2)
 174 {
 175         int i, deleting = 0;
 176         mutseq_t *ret[2];
 177
 178         ret[0] = hap1; ret[1] = hap2;
 179         ret[0]->l = seq->l; ret[1]->l = seq->l;
 180         ret[0]->m = seq->m; ret[1]->m = seq->m;
 181         ret[0]->s = (mut_t *)calloc(seq->m, sizeof(mut_t));
 182         ret[1]->s = (mut_t *)calloc(seq->m, sizeof(mut_t));
 183         for (i = 0; i != seq->l; ++i) {
 184                 int c;
 185                 c = ret[0]->s[i] = ret[1]->s[i] = (mut_t)nst_nt4_table[(int)seq->s[i]];
 186         if (deleting) {
 187             if (drand48() < INDEL_EXTEND) {
 188                 if (deleting & 1) ret[0]->s[i] |= DELETE;
 189                 if (deleting & 2) ret[1]->s[i] |= DELETE;
 190                 continue;
 191             } else deleting = 0;
 192         }
 193                 if (c < 4 && drand48() < MUT_RATE) { // mutation
 194                         if (drand48() >= INDEL_FRAC) { // substitution
 195                                 double r = drand48();
 196                                 c = (c + (int)(r * 3.0 + 1)) & 3;
 197                                 if (is_hap || drand48() < 0.333333) { // hom
 198                                         ret[0]->s[i] = ret[1]->s[i] = SUBSTITUTE|c;
 199                                 } else { // het
 200                                         ret[drand48()<0.5?0:1]->s[i] = SUBSTITUTE|c;
 201                                 }
 202                         } else { // indel
 203                                 if (drand48() < 0.5) { // deletion
 204                                         if (is_hap || drand48() < 0.333333) { // hom-del
 205                                                 ret[0]->s[i] = ret[1]->s[i] = DELETE;
 206                         deleting = 3;
 207                                         } else { // het-del
 208                         deleting = drand48()<0.5?1:2;
 209                                                 ret[deleting-1]->s[i] = DELETE;
 210                                         }
 211                                 } else { // insertion
 212                     int num_ins = 0, ins = 0;
 213                     do {
 214                         num_ins++;
 215                         ins = (ins << 2) | (int)(drand48() * 4.0);
 216                     } while(num_ins < 4 && drand48() < INDEL_EXTEND);
 217
 218                                         if (is_hap || drand48() < 0.333333) { // hom-ins
 219                                                 ret[0]->s[i] = ret[1]->s[i] = (num_ins << 12) | (ins << 4) | c;
 220                                         } else { // het-ins
 221                                                 ret[drand48()<0.5?0:1]->s[i] = (num_ins << 12) | (ins << 4) | c;
 222                                         }
 223                                 }
 224                         }
 225                 }
 226         }
 227 }
 228 void maq_print_mutref(const char *name, const seq_t *seq, mutseq_t *hap1, mutseq_t *hap2)
 229 {
 230         int i;
 231         for (i = 0; i != seq->l; ++i) {
 232                 int c[3];
 233                 c[0] = nst_nt4_table[(int)seq->s[i]];
 234                 c[1] = hap1->s[i]; c[2] = hap2->s[i];
 235                 if (c[0] >= 4) continue;
 236                 if ((c[1] & mutmsk) != NOCHANGE || (c[1] & mutmsk) != NOCHANGE) {
 237                         printf("%s\t%d\t", name, i+1);
 238                         if (c[1] == c[2]) { // hom
 239                                 if ((c[1]&mutmsk) == SUBSTITUTE) { // substitution
 240                                         printf("%c\t%c\t-\n", "ACGTN"[c[0]], "ACGTN"[c[1]&0xf]);
 241                                 } else if ((c[1]&mutmsk) == DELETE) { // del
 242                                         printf("%c\t-\t-\n", "ACGTN"[c[0]]);
 243                                 } else if (((c[1] & mutmsk) >> 12) <= 5) { // ins
 244                                         printf("-\t");
 245                     int n = (c[1]&mutmsk) >> 12, ins = c[1] >> 4;
 246                     while(n > 0) {
 247                         putchar("ACGTN"[ins & 0x3]);
 248                         n--;
 249                     }
 250                     printf("\t-\n");
 251                                 }  else assert(0);
 252                         } else { // het
 253                                 if ((c[1]&mutmsk) == SUBSTITUTE || (c[2]&mutmsk) == SUBSTITUTE) { // substitution
 254                                         printf("%c\t%c\t+\n", "ACGTN"[c[0]], "XACMGRSVTWYHKDBN"[1<<(c[1]&0x3)|1<<(c[2]&0x3)]);
 255                                 } else if ((c[1]&mutmsk) == DELETE) {
 256                                         printf("%c\t-\t+\n", "ACGTN"[c[0]]);
 257                                 } else if ((c[2]&mutmsk) == DELETE) {
 258                                         printf("%c\t-\t+\n", "ACGTN"[c[0]]);
 259                                 } else if (((c[1] & mutmsk) >> 12) <= 4) { // ins1
 260                                         printf("-\t");
 261                     int n = (c[1]&mutmsk) >> 12, ins = c[1] >> 4;
 262                     while (n > 0) {
 263                         putchar("ACGTN"[ins & 0x3]);
 264                         n--;
 265                     }
 266                     printf("\t+\n");
 267                                 } else if (((c[2] & mutmsk) >> 12) <= 5) { // ins2
 268                                         printf("-\t");
 269                     int n = (c[2]&mutmsk) >> 12, ins = c[2] >> 4;
 270                     while (n > 0) {
 271                         putchar("ACGTN"[ins & 0x3]);
 272                         ins >>= 2;
 273                         n--;
 274                     }
 275                     printf("\t+\n");
 276                                 } else assert(0);
 277                         }
 278                 }
 279         }
 280 }
 281
 282 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)
 283 {
 284         seq_t seq;
 285     mutseq_t rseq[2];
 286         uint64_t tot_len, ii;
 287         int i, l, n_ref, k;
 288         char name[256], *qstr, *str;
 289         int size[2], Q;
 290         uint8_t *tmp_seq[2];
 291     mut_t *target;
 292
 293         INIT_SEQ(seq);
 294         srand48(time(0));
 295         seq_set_block_size(0x1000000);
 296         l = size_l > size_r? size_l : size_r;
 297         qstr = (char*)calloc(l+1, 1);
 298         str = (char*)calloc(l+1, 1);
 299         tmp_seq[0] = (uint8_t*)calloc(l+1, 1);
 300         tmp_seq[1] = (uint8_t*)calloc(l+1, 1);
 301         size[0] = size_l; size[1] = size_r;
 302
 303         Q = (int)(-10.0 * log(ERR_RATE) / log(10.0) + 0.499) + 33;
 304
 305         tot_len = n_ref = 0;
 306         while ((l = seq_read_fasta(fp_fa, &seq, name, 0)) >= 0) {
 307                 tot_len += l;
 308                 ++n_ref;
 309         }
 310         fprintf(stderr, "-- %d sequences, total length: %llu\n", n_ref, (long long)tot_len);
 311         rewind(fp_fa);
 312
 313         while ((l = seq_read_fasta(fp_fa, &seq, name, 0)) >= 0) {
 314                 uint64_t n_pairs = (uint64_t)((long double)l / tot_len * N + 0.5);
 315                 if (l < dist + 3 * std_dev) {
 316                         fprintf(stderr, "[wgsim_core] Skip sequence '%s' as it is shorter than %d!\n", name, dist + 3 * std_dev);
 317                         continue;
 318                 }
 319                 maq_mut_diref(&seq, is_hap, rseq, rseq+1);
 320                 maq_print_mutref(name, &seq, rseq, rseq+1);
 321                 for (ii = 0; ii != n_pairs; ++ii) {
 322                         double ran;
 323                         int d, pos, s[2], begin, end, is_flip = 0;
 324                         FILE *fpo[2];
 325                         do {
 326                                 ran = ran_normal();
 327                                 ran = ran * std_dev + dist;
 328                                 d = (int)(ran + 0.5);
 329                                 pos = (int)((l - d + 1) * drand48());
 330                         } while (pos < 0 || pos >= seq.l || pos + d - 1 >= seq.l);
 331                         // flip or not
 332                         if (drand48() < 0.5) {
 333                                 fpo[0] = fpout1; fpo[1] = fpout2;
 334                                 s[0] = size[0]; s[1] = size[1];
 335                         } else {
 336                                 fpo[1] = fpout1; fpo[0] = fpout2;
 337                                 s[1] = size[0]; s[0] = size[1];
 338                                 is_flip = 1;
 339                         }
 340                         // generate the read sequences
 341                         target = rseq[drand48()<0.5?0:1].s; // haploid from which the reads are generated
 342                         for (i = pos, k = 0, begin = 0; i < seq.l && k < s[0]; ++i) {
 343                                 int c = target[i];
 344                 int mut_type = c & mutmsk;
 345                                 if (mut_type == DELETE) continue; // deletion
 346                                 if (begin == 0) {
 347                                         begin = i;
 348                                         if (mut_type != NOCHANGE && mut_type != SUBSTITUTE) mut_type = NOCHANGE; // skip ins at the first base
 349                                 }
 350                 if(mut_type == NOCHANGE || mut_type == SUBSTITUTE) {
 351                     tmp_seq[0][k++] = c&0xf;
 352                     continue;
 353                 }
 354                 int n = mut_type >> 12, ins = c >> 4;
 355                 while (n > 0) {
 356                     tmp_seq[0][k++] = ins & 0x3;
 357                     ins >>= 2;
 358                     n--;
 359                     if(k == s[0]) break;
 360                 }
 361                                 tmp_seq[0][k++] = c&0xf;
 362                         }
 363                         for (i = pos + d - 1, k = 0, end = 0; i >= 0  && k < s[1]; --i) {
 364                                 int c = target[i];
 365                                 if ((c&mutmsk) == DELETE) continue; // deletion
 366                                 if (end == 0) end = i;
 367                 tmp_seq[1][k++] = c&0xf;
 368                 if((c&mutmsk) == NOCHANGE || (c&mutmsk) == SUBSTITUTE) continue;
 369                 int n = (c&mutmsk) >> 12, ins = c >> 4;
 370                 while (n > 0) {
 371                     if (k == s[1]) break;
 372                     tmp_seq[1][k++] = ins & 0x3;
 373                     ins >>= 2;
 374                     n--;
 375                 }
 376                         }
 377                         // start to print out
 378                         if (!is_flip) sprintf(str, "%s_%u_%u_%llx", name, begin+1, end+1, (long long)ii);
 379                         else sprintf(str, "%s_%u_%u_%llx", name, begin+1, end+1, (long long)ii);
 380                         // print forward read
 381                         fprintf(fpo[0], "@%s/%d\n", str, is_flip+1);
 382                         for (i = 0; i < s[0]; ++i) {
 383                                 int c = tmp_seq[0][i];
 384                                 if (c > 4) c = 4;
 385                                 if (c < 4) {
 386                                         if (drand48() < ERR_RATE)
 387                                                 c = (c + (int)(drand48() * 3.0 + 1)) & 3;
 388                                 }
 389                                 fputc("ACGTN"[c], fpo[0]);
 390                         }
 391                         for (i = 0; i < s[0]; ++i) qstr[i] = Q;
 392                         qstr[s[0]] = 0;
 393                         fprintf(fpo[0], "\n+\n%s\n", qstr);
 394                         // print reverse read
 395                         fprintf(fpo[1], "@%s/%d\n", str, 2-is_flip);
 396                         for (i = 0; i < s[1]; ++i) {
 397                                 int c = tmp_seq[1][i];
 398                                 if (c > 4) c = 4;
 399                                 if (c < 4) {
 400                                         c = 3 - c; // complement
 401                                         if (drand48() < ERR_RATE)
 402                                                 c = (c + (int)(drand48() * 3.0 + 1)) & 3;
 403                                 }
 404                                 fputc("ACGTN"[c], fpo[1]);
 405                         }
 406                         for (i = 0; i < s[1]; ++i) qstr[i] = Q;
 407                         qstr[s[1]] = 0;
 408                         fprintf(fpo[1], "\n+\n%s\n", qstr);
 409                 }
 410                 free(rseq[0].s); free(rseq[1].s);
 411         }
 412         free(seq.s);
 413         free(qstr); free(str);
 414         free(tmp_seq[0]); free(tmp_seq[1]);
 415 }
 416
 417 static int simu_usage()
 418 {
 419         fprintf(stderr, "\n");
 420         fprintf(stderr, "Usage:   wgsim [options] <in.ref.fa> <out.read1.fq> <out.read2.fq>\n\n");
 421         fprintf(stderr, "Options: -e FLOAT      base error rate [%.3f]\n", ERR_RATE);
 422         fprintf(stderr, "         -d INT        outer distance between the two ends [500]\n");
 423         fprintf(stderr, "         -s INT        standard deviation [50]\n");
 424         fprintf(stderr, "         -N INT        number of read pairs [1000000]\n");
 425         fprintf(stderr, "         -1 INT        length of the first read [70]\n");
 426         fprintf(stderr, "         -2 INT        length of the second read [70]\n");
 427         fprintf(stderr, "         -r FLOAT      rate of mutations [%.4f]\n", MUT_RATE);
 428         fprintf(stderr, "         -R FLOAT      fraction of indels [%.2f]\n", INDEL_FRAC);
 429         fprintf(stderr, "         -X FLOAT      probability an indel is extended [%.2f]\n", INDEL_EXTEND);
 430         fprintf(stderr, "         -h            haploid mode\n");
 431         fprintf(stderr, "\n");
 432         return 1;
 433 }
 434
 435 int main(int argc, char *argv[])
 436 {
 437         int64_t N;
 438         int dist, std_dev, c, size_l, size_r, is_hap = 0;
 439         FILE *fpout1, *fpout2, *fp_fa;
 440
 441         N = 1000000; dist = 500; std_dev = 50;
 442         size_l = size_r = 70;
 443         while ((c = getopt(argc, argv, "e:d:s:N:1:2:r:R:hX:")) >= 0) {
 444                 switch (c) {
 445                 case 'd': dist = atoi(optarg); break;
 446                 case 's': std_dev = atoi(optarg); break;
 447                 case 'N': N = atoi(optarg); break;
 448                 case '1': size_l = atoi(optarg); break;
 449                 case '2': size_r = atoi(optarg); break;
 450                 case 'e': ERR_RATE = atof(optarg); break;
 451                 case 'r': MUT_RATE = atof(optarg); break;
 452                 case 'R': INDEL_FRAC = atof(optarg); break;
 453                 case 'X': INDEL_EXTEND = atof(optarg); break;
 454                 case 'h': is_hap = 1; break;
 455                 }
 456         }
 457         if (argc - optind < 3) return simu_usage();
 458         fp_fa = (strcmp(argv[optind+0], "-") == 0)? stdin : xopen(argv[optind+0], "r");
 459         fpout1 = xopen(argv[optind+1], "w");
 460         fpout2 = xopen(argv[optind+2], "w");
 461         wgsim_core(fpout1, fpout2, fp_fa, is_hap, N, dist, std_dev, size_l, size_r);
 462
 463         fclose(fpout1); fclose(fpout2); fclose(fp_fa);
 464         return 0;
 465 }