- target = rseq[drand48()<0.5?0:1].s; // haploid from which the reads are generated
- for (i = pos, k = 0, begin = 0; i < seq.l && k < s[0]; ++i) {
- int c = target[i];
- int mut_type = c & mutmsk;
- if (mut_type == DELETE) continue; // deletion
- if (begin == 0) {
- begin = i;
- if (mut_type != NOCHANGE && mut_type != SUBSTITUTE) mut_type = NOCHANGE; // skip ins at the first base
+ target = rseq[drand48()<0.5?0:1].s; // haplotype from which the reads are generated
+ n_sub[0] = n_sub[1] = n_indel[0] = n_indel[1] = n_err[0] = n_err[1] = 0;
+
+#define __gen_read(x, start, iter) do { \
+ for (i = (start), k = 0, ext_coor[x] = -1; i >= 0 && i < seq.l && k < s[x]; iter) { \
+ int c = target[i], mut_type = c & mutmsk; \
+ if (ext_coor[x] < 0) { \
+ ext_coor[x] = i; \
+ if (mut_type != NOCHANGE && mut_type != SUBSTITUTE) { \
+ ext_coor[x] = -1; \
+ break; \
+ } \
+ } \
+ if (mut_type == DELETE) ++n_indel[x]; \
+ else if (mut_type == NOCHANGE || mut_type == SUBSTITUTE) { \
+ tmp_seq[x][k++] = c & 0xf; \
+ if (mut_type == SUBSTITUTE) ++n_sub[x]; \
+ } else { \
+ int n, ins; \
+ ++n_indel[x]; \
+ for (n = mut_type>>12, ins = c>>4; n > 0 && k < s[x]; --n, ins >>= 2) \
+ tmp_seq[x][k++] = ins & 0x3; \
+ } \
+ } \
+ if (k != s[x]) ext_coor[x] = -1; \
+ } while (0)
+
+ if (!IS_SOLID) {
+ __gen_read(0, pos, ++i);
+ __gen_read(1, pos + d - 1, --i);
+ for (k = 0; k < s[1]; ++k) tmp_seq[1][k] = tmp_seq[1][k] < 4? 3 - tmp_seq[1][k] : 4; // complement
+ } else {
+ int c1, c2, c;
+ ++s[0]; ++s[1]; // temporarily increase read length by 1
+ if (is_flip) { // RR pair
+ __gen_read(0, pos + s[0], --i);
+ __gen_read(1, pos + d - 1, --i);
+ } else { // FF pair
+ __gen_read(0, pos, ++i);
+ __gen_read(1, pos + d - 1 - s[1], ++i);
+ }
+ // change to color sequence: (0,1,2,3) -> (A,C,G,T)
+ for (j = 0; j < 2; ++j) {
+ c1 = tmp_seq[j][0];
+ for (i = 1; i < s[j]; ++i) {
+ c2 = tmp_seq[j][i];
+ c = (c1 >= 4 || c2 >= 4)? 4 : nst_color_space_table[(1<<c1)|(1<<c2)];
+ tmp_seq[j][i-1] = c;
+ c1 = c2;
+ }