From: Heng Li Date: Sun, 18 Oct 2009 00:57:26 +0000 (+0000) Subject: * samtools-0.1.6-13 (r479) X-Git-Url: https://git.donarmstrong.com/?p=samtools.git;a=commitdiff_plain;h=4c2687f34917ab1767b8c5ae60b8aba564d2c95a * 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 3b49020..7aed741 100644 --- 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() = \theta \sum_{i=1}^{N-1} 1/i P(D|) = \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 2a82aee..fa5489d 100644 --- 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 5bf1ed0..fff42f7 100644 --- 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] |\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 19a93cb..9884f3d 100644 --- 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] [...]\n\n"); + fprintf(stderr, "Usage: samtools merge [-nr] [-h inh.sam] [...]\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 [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 d4d2793..b67971f 100644 --- 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[]);