* samtools-0.1.7-6

[samtools.git] / misc / wgsim.c

diff --git a/misc/wgsim.c b/misc/wgsim.c
index 485dce1ac1da06cf4054c4d5a391871f929adb2a..7b5f095d910b06f3242419066ac8879b7cdb7454 100644 (file)
--- a/misc/wgsim.c
+++ b/misc/wgsim.c
@@ -39,7 +39,7 @@
 #include <ctype.h>
 #include <string.h>
 
-#define PACKAGE_VERSION "0.2.2"
+#define PACKAGE_VERSION "0.2.3"
 
 const uint8_t nst_nt4_table[256] = {
        4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
@@ -238,7 +238,7 @@ void maq_print_mutref(const char *name, const seq_t *seq, mutseq_t *hap1, mutseq
                c[0] = nst_nt4_table[(int)seq->s[i]];
                c[1] = hap1->s[i]; c[2] = hap2->s[i];
                if (c[0] >= 4) continue;
-               if ((c[1] & mutmsk) != NOCHANGE || (c[1] & mutmsk) != NOCHANGE) {
+               if ((c[1] & mutmsk) != NOCHANGE || (c[2] & mutmsk) != NOCHANGE) {
                        printf("%s\t%d\t", name, i+1);
                        if (c[1] == c[2]) { // hom
                                if ((c[1]&mutmsk) == SUBSTITUTE) { // substitution
@@ -304,7 +304,7 @@ void wgsim_core(FILE *fpout1, FILE *fpout2, FILE *fp_fa, int is_hap, uint64_t N,
        tmp_seq[1] = (uint8_t*)calloc(l+2, 1);
        size[0] = size_l; size[1] = size_r;
 
-       Q = (int)(-10.0 * log(ERR_RATE) / log(10.0) + 0.499) + 33;
+       Q = (ERR_RATE == 0.0)? 'I' : (int)(-10.0 * log(ERR_RATE) / log(10.0) + 0.499) + 33;
 
        tot_len = n_ref = 0;
        while ((l = seq_read_fasta(fp_fa, &seq, name, 0)) >= 0) {
@@ -366,6 +366,7 @@ void wgsim_core(FILE *fpout1, FILE *fpout2, FILE *fp_fa, int is_hap, uint64_t N,
                                        } else {                                                                                        \
                                                int n, ins;                                                                             \
                                                ++n_indel[x];                                                                   \
+                                               tmp_seq[x][k++] = c & 0xf;                                              \
                                                for (n = mut_type>>12, ins = c>>4; n > 0 && k < s[x]; --n, ins >>= 2) \
                                                        tmp_seq[x][k++] = ins & 0x3;                            \
                                        }                                                                                                       \