* samtools-0.1.6-13 (r479)

 * merge: optionally use file names as RG tags

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 3b49020a81a9ccca89c45745f60605df3cc89661..7aed741ca9c5df15df51c0d72af743f23388ef0f 100644 (file)
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -25,14 +25,13 @@ char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
 /*
   P(<b1,b2>) = \theta \sum_{i=1}^{N-1} 1/i
   P(D|<b1,b2>) = \sum_{k=1}^{N-1} p_k 1/2 [(k/N)^n_2(1-k/N)^n_1 + (k/N)^n1(1-k/N)^n_2]
-  p_k = i/k / \sum_{i=1}^{N-1} 1/i
+  p_k = 1/k / \sum_{i=1}^{N-1} 1/i
  */
 static void cal_het(bam_maqcns_t *aa)
 {
        int k, n1, n2;
        double sum_harmo; // harmonic sum
        double poly_rate;
-       double p1 = 0.0, p3 = 0.0; // just for testing
 
        free(aa->lhet);
        aa->lhet = (double*)calloc(256 * 256, sizeof(double));
@@ -42,7 +41,7 @@ static void cal_het(bam_maqcns_t *aa)
        for (n1 = 0; n1 < 256; ++n1) {
                for (n2 = 0; n2 < 256; ++n2) {
                        long double sum = 0.0;
-                       double lC = lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1); // \binom{n1+n2}{n1}
+                       double lC = aa->is_soap? 0 : lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1); // \binom{n1+n2}{n1}
                        for (k = 1; k <= aa->n_hap - 1; ++k) {
                                double pk = 1.0 / k / sum_harmo;
                                double log1 = log((double)k/aa->n_hap);
@@ -50,8 +49,6 @@ static void cal_het(bam_maqcns_t *aa)
                                sum += pk * 0.5 * (expl(log1*n2) * expl(log2*n1) + expl(log1*n1) * expl(log2*n2));
                        }
                        aa->lhet[n1<<8|n2] = lC + logl(sum);
-                       if (n1 == 17 && n2 == 3) p3 = lC + logl(expl(logl(0.5) * 20));
-                       if (n1 == 19 && n2 == 1) p1 = lC + logl(expl(logl(0.5) * 20));
                }
        }
        poly_rate = aa->het_rate * sum_harmo;
@@ -65,16 +62,18 @@ static void cal_coef(bam_maqcns_t *aa)
        long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256];
        double *lC;
 
-       lC = (double*)calloc(256 * 256, sizeof(double));
        // aa->lhet will be allocated and initialized 
        free(aa->fk); free(aa->coef);
+       aa->coef = 0;
        aa->fk = (double*)calloc(256, sizeof(double));
-       aa->coef = (double*)calloc(256*256*64, sizeof(double));
        aa->fk[0] = fk2[0] = 1.0;
        for (n = 1; n != 256; ++n) {
                aa->fk[n] = pow(aa->theta, n) * (1.0 - aa->eta) + aa->eta;
                fk2[n] = aa->fk[n>>1]; // this is an approximation, assuming reads equally likely come from both strands
        }
+       if (aa->is_soap) return;
+       aa->coef = (double*)calloc(256*256*64, sizeof(double));
+       lC = (double*)calloc(256 * 256, sizeof(double));
        for (n = 1; n != 256; ++n)
                for (k = 1; k <= n; ++k)
                        lC[n<<8|k] = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
@@ -183,56 +182,73 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
                for (j = 0; j != 4; ++j) b->c[j] = (int)(254.0 * b->c[j] / c + 0.5);
                for (j = c = 0; j != 4; ++j) c += b->c[j];
        }
-       // generate likelihood
-       for (j = 0; j != 4; ++j) {
-               // homozygous
-               float tmp1, tmp3;
-               int tmp2, bar_e;
-               for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) {
-                       if (j == k) continue;
-                       tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k];
-               }
-               if (tmp2) {
-                       bar_e = (int)(tmp1 / tmp3 + 0.5);
-                       if (bar_e < 4) bar_e = 4; // should not happen
-                       if (bar_e > 63) bar_e = 63;
-                       p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
-               } else p[j<<2|j] = 0.0; // all the bases are j
-               // heterozygous
-               for (k = j + 1; k < 4; ++k) {
-                       for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) {
-                               if (i == j || i == k) continue;
-                               tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i];
+       if (!bm->is_soap) {
+               // generate likelihood
+               for (j = 0; j != 4; ++j) {
+                       // homozygous
+                       float tmp1, tmp3;
+                       int tmp2, bar_e;
+                       for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) {
+                               if (j == k) continue;
+                               tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k];
                        }
                        if (tmp2) {
                                bar_e = (int)(tmp1 / tmp3 + 0.5);
-                               if (bar_e < 4) bar_e = 4;
+                               if (bar_e < 4) bar_e = 4; // should not happen
                                if (bar_e > 63) bar_e = 63;
-                               p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
-                       } else p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]]; // all the bases are either j or k
+                               p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
+                       } else p[j<<2|j] = 0.0; // all the bases are j
+                       // heterozygous
+                       for (k = j + 1; k < 4; ++k) {
+                               for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) {
+                                       if (i == j || i == k) continue;
+                                       tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i];
+                               }
+                               if (tmp2) {
+                                       bar_e = (int)(tmp1 / tmp3 + 0.5);
+                                       if (bar_e < 4) bar_e = 4;
+                                       if (bar_e > 63) bar_e = 63;
+                                       p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
+                               } else p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]]; // all the bases are either j or k
+                       }
+                       //
+                       for (k = 0; k != 4; ++k)
+                               if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
                }
-               //
-               for (k = 0; k != 4; ++k)
-                       if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
-       }
 
-       { // fix p[k<<2|k]
-               float max1, max2, min1, min2;
-               int max_k, min_k;
-               max_k = min_k = -1;
-               max1 = max2 = -1.0; min1 = min2 = 1e30;
-               for (k = 0; k < 4; ++k) {
-                       if (b->esum[k] > max1) {
-                               max2 = max1; max1 = b->esum[k]; max_k = k;
-                       } else if (b->esum[k] > max2) max2 = b->esum[k];
+               { // fix p[k<<2|k]
+                       float max1, max2, min1, min2;
+                       int max_k, min_k;
+                       max_k = min_k = -1;
+                       max1 = max2 = -1.0; min1 = min2 = 1e30;
+                       for (k = 0; k < 4; ++k) {
+                               if (b->esum[k] > max1) {
+                                       max2 = max1; max1 = b->esum[k]; max_k = k;
+                               } else if (b->esum[k] > max2) max2 = b->esum[k];
+                       }
+                       for (k = 0; k < 4; ++k) {
+                               if (p[k<<2|k] < min1) {
+                                       min2 = min1; min1 = p[k<<2|k]; min_k = k;
+                               } else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
+                       }
+                       if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
+                               p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
                }
-               for (k = 0; k < 4; ++k) {
-                       if (p[k<<2|k] < min1) {
-                               min2 = min1; min1 = p[k<<2|k]; min_k = k;
-                       } else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
+       } else { // apply the SOAP model
+               // generate likelihood
+               for (j = 0; j != 4; ++j) {
+                       float tmp;
+                       // homozygous
+                       for (k = 0, tmp = 0.0; k != 4; ++k)
+                               if (j != k) tmp += b->esum[k];
+                       p[j<<2|j] = tmp;
+                       // heterozygous
+                       for (k = j + 1; k < 4; ++k) {
+                               for (i = 0, tmp = 0.0; i != 4; ++i)
+                                       if (i != j && i != k) tmp += b->esum[i];
+                               p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp;
+                       }
                }
-               if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
-                       p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
        }
 
        // convert necessary information to glf1_t

diff --git a/bam_maqcns.h b/bam_maqcns.h
index 2a82aeee8057f5992570d0896a1621b9116057cd..fa5489d13caad7c5c20489091302bfb5b2a5799d 100644 (file)
--- a/bam_maqcns.h
+++ b/bam_maqcns.h
@@ -7,7 +7,7 @@ struct __bmc_aux_t;
 
 typedef struct {
        float het_rate, theta;
-       int n_hap, cap_mapQ;
+       int n_hap, cap_mapQ, is_soap;
 
        float eta, q_r;
        double *fk, *coef;

diff --git a/bam_plcmd.c b/bam_plcmd.c
index 5bf1ed095096b7cd91fc99d2256467b89e394d33..fff42f7c582b0f40947616723e3c3fc15ca630c9 100644 (file)
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -299,8 +299,9 @@ int bam_pileup(int argc, char *argv[])
        d->tid = -1; d->mask = BAM_DEF_MASK;
        d->c = bam_maqcns_init();
        d->ido = bam_maqindel_opt_init();
-       while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:S2")) >= 0) {
+       while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:S2a")) >= 0) {
                switch (c) {
+               case 'a': d->c->is_soap = 1; break;
                case 's': d->format |= BAM_PLF_SIMPLE; break;
                case 't': fn_list = strdup(optarg); break;
                case 'l': fn_pos = strdup(optarg); break;
@@ -327,6 +328,7 @@ int bam_pileup(int argc, char *argv[])
                fprintf(stderr, "Usage:  samtools pileup [options] <in.bam>|<in.sam>\n\n");
                fprintf(stderr, "Option: -s        simple (yet incomplete) pileup format\n");
                fprintf(stderr, "        -S        the input is in SAM\n");
+               fprintf(stderr, "        -a        use the SOAPsnp model for SNP calling\n");
                fprintf(stderr, "        -2        output the 2nd best call and quality\n");
                fprintf(stderr, "        -i        only show lines/consensus with indels\n");
                fprintf(stderr, "        -m INT    filtering reads with bits in INT [%d]\n", d->mask);

diff --git a/bam_sort.c b/bam_sort.c
index 19a93cb348c7123fc3392d348dbf674b36790306..9884f3d6273f680709ee0e44532a8371413fddf4 100644 (file)
--- a/bam_sort.c
+++ b/bam_sort.c
@@ -71,14 +71,15 @@ static void swap_header_text(bam_header_t *h1, bam_header_t *h2)
   @discussion Padding information may NOT correctly maintained. This
   function is NOT thread safe.
  */
-void bam_merge_core(int by_qname, const char *out, const char *headers, int n, char * const *fn)
+void bam_merge_core(int by_qname, const char *out, const char *headers, int n, char * const *fn, int add_RG)
 {
        bamFile fpout, *fp;
        heap1_t *heap;
        bam_header_t *hout = 0;
        bam_header_t *hheaders = NULL;
-       int i, j;
+       int i, j, *RG_len = 0;
        uint64_t idx = 0;
+       char **RG = 0;
 
        if (headers) {
                tamFile fpheaders = sam_open(headers);
@@ -93,6 +94,22 @@ void bam_merge_core(int by_qname, const char *out, const char *headers, int n, c
        g_is_by_qname = by_qname;
        fp = (bamFile*)calloc(n, sizeof(bamFile));
        heap = (heap1_t*)calloc(n, sizeof(heap1_t));
+       // prepare RG tag
+       if (add_RG) {
+               RG = (char**)calloc(n, sizeof(void*));
+               RG_len = (int*)calloc(n, sizeof(int));
+               for (i = 0; i != n; ++i) {
+                       int l = strlen(fn[i]);
+                       const char *s = fn[i];
+                       if (l > 4 && strcmp(s + l - 4, ".bam") == 0) l -= 4;
+                       for (j = l - 1; j >= 0; --j) if (s[j] == '/') break;
+                       ++j; l -= j;
+                       RG[i] = calloc(l + 1, 1);
+                       RG_len[i] = l;
+                       strncpy(RG[i], s + j, l);
+               }
+       }
+       // read the first
        for (i = 0; i != n; ++i) {
                heap1_t *h;
                bam_header_t *hin;
@@ -154,6 +171,8 @@ void bam_merge_core(int by_qname, const char *out, const char *headers, int n, c
        ks_heapmake(heap, n, heap);
        while (heap->pos != HEAP_EMPTY) {
                bam1_t *b = heap->b;
+               if (add_RG && bam_aux_get(b, "RG") == 0)
+                       bam_aux_append(b, "RG", 'Z', RG_len[heap->i] + 1, (uint8_t*)RG[heap->i]);
                bam_write1_core(fpout, &b->core, b->data_len, b->data);
                if ((j = bam_read1(fp[heap->i], b)) >= 0) {
                        heap->pos = ((uint64_t)b->core.tid<<32) | (uint32_t)b->core.pos<<1 | bam1_strand(b);
@@ -166,32 +185,38 @@ void bam_merge_core(int by_qname, const char *out, const char *headers, int n, c
                ks_heapadjust(heap, 0, n, heap);
        }
 
+       if (add_RG) {
+               for (i = 0; i != n; ++i) free(RG[i]);
+               free(RG); free(RG_len);
+       }
        for (i = 0; i != n; ++i) bam_close(fp[i]);
        bam_close(fpout);
        free(fp); free(heap);
 }
 int bam_merge(int argc, char *argv[])
 {
-       int c, is_by_qname = 0;
+       int c, is_by_qname = 0, add_RG = 0;
        char *fn_headers = NULL;
 
-       while ((c = getopt(argc, argv, "h:n")) >= 0) {
+       while ((c = getopt(argc, argv, "h:nr")) >= 0) {
                switch (c) {
+               case 'r': add_RG = 1; break;
                case 'h': fn_headers = strdup(optarg); break;
                case 'n': is_by_qname = 1; break;
                }
        }
        if (optind + 2 >= argc) {
                fprintf(stderr, "\n");
-               fprintf(stderr, "Usage:   samtools merge [-n] [-h inh.sam] <out.bam> <in1.bam> <in2.bam> [...]\n\n");
+               fprintf(stderr, "Usage:   samtools merge [-nr] [-h inh.sam] <out.bam> <in1.bam> <in2.bam> [...]\n\n");
                fprintf(stderr, "Options: -n       sort by read names\n");
+               fprintf(stderr, "         -r       attach RG tag (inferred from file names)\n");
                fprintf(stderr, "         -h FILE  copy the header in FILE to <out.bam> [in1.bam]\n\n");
                fprintf(stderr, "Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users\n");
                fprintf(stderr, "      must provide the correct header with -h, or uses Picard which properly maintains\n");
                fprintf(stderr, "      the header dictionary in merging.\n\n");
                return 1;
        }
-       bam_merge_core(is_by_qname, argv[optind], fn_headers, argc - optind - 1, argv + optind + 1);
+       bam_merge_core(is_by_qname, argv[optind], fn_headers, argc - optind - 1, argv + optind + 1, add_RG);
        free(fn_headers);
        return 0;
 }
@@ -288,7 +313,7 @@ void bam_sort_core_ext(int is_by_qname, const char *fn, const char *prefix, size
                        fns[i] = (char*)calloc(strlen(prefix) + 20, 1);
                        sprintf(fns[i], "%s.%.4d.bam", prefix, i);
                }
-               bam_merge_core(is_by_qname, fnout, 0, n, fns);
+               bam_merge_core(is_by_qname, fnout, 0, n, fns, 0);
                free(fnout);
                for (i = 0; i < n; ++i) {
                        unlink(fns[i]);

diff --git a/bamtk.c b/bamtk.c
index d4d27935bc934b5a1f1ed44f5723e1cf7c00d38b..b67971fd3ffec0b31b3ec8fc099de82ce9d82b62 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -9,7 +9,7 @@
 #endif
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.6-12 (r478)"
+#define PACKAGE_VERSION "0.1.6-13 (r479)"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);