DFLAGS= -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -D_USE_KNETFILE -D_CURSES_LIB=1
KNETFILE_O= knetfile.o
LOBJS= bgzf.o kstring.o bam_aux.o bam.o bam_import.o sam.o bam_index.o \
- bam_pileup.o bam_lpileup.o bam_md.o glf.o razf.o faidx.o bedidx.o \
+ bam_pileup.o bam_lpileup.o bam_md.o razf.o faidx.o bedidx.o \
$(KNETFILE_O) bam_sort.o sam_header.o bam_reheader.o kprobaln.o bam_cat.o
-AOBJS= bam_tview.o bam_maqcns.o bam_plcmd.o sam_view.o \
+AOBJS= bam_tview.o bam_plcmd.o sam_view.o \
bam_rmdup.o bam_rmdupse.o bam_mate.o bam_stat.o bam_color.o \
bamtk.o kaln.o bam2bcf.o bam2bcf_indel.o errmod.o sample.o \
cut_target.o phase.o bam2depth.o
lib:libbam.a
libbam.a:$(LOBJS)
- $(AR) -cru $@ $(LOBJS)
+ $(AR) -csru $@ $(LOBJS)
samtools:lib-recur $(AOBJS)
$(CC) $(CFLAGS) -o $@ $(AOBJS) -Lbcftools $(LIBPATH) libbam.a -lbcf $(LIBCURSES) -lm -lz
sam.o:sam.h bam.h
bam_import.o:bam.h kseq.h khash.h razf.h
bam_pileup.o:bam.h razf.h ksort.h
-bam_plcmd.o:bam.h faidx.h bam_maqcns.h glf.h bcftools/bcf.h bam2bcf.h
+bam_plcmd.o:bam.h faidx.h bcftools/bcf.h bam2bcf.h
bam_index.o:bam.h khash.h ksort.h razf.h bam_endian.h
bam_lpileup.o:bam.h ksort.h
-bam_tview.o:bam.h faidx.h bam_maqcns.h
-bam_maqcns.o:bam.h ksort.h bam_maqcns.h kaln.h
+bam_tview.o:bam.h faidx.h
bam_sort.o:bam.h ksort.h razf.h
bam_md.o:bam.h faidx.h
-glf.o:glf.h
sam_header.o:sam_header.h khash.h
bcf.o:bcftools/bcf.h
bam2bcf.o:bam2bcf.h errmod.h bcftools/bcf.h
+Beta Release 0.1.17 (6 July, 2011)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+With the maturity of `mpileup' and the lack of update in the `pileup' command,
+the `pileup' command is now formally dropped. Most of the pileup functionality,
+such as outputting mapping quality and read positions, have been added
+`mpileup'.
+
+Since this release, `bcftools view' is able to perform contrast SNP calling
+(option -T) for discovering de novo and/or somatic mutations between a pair of
+samples or in a family trio. Potential mutations are scored by a log likelihood
+ratio, which is very simple in math, but should be comparable to more
+sophisticated methods. Note that getting the score is only the very first step.
+A lot more need to be done to reduce systematical errors due to mapping and
+reference errors and structural variations.
+
+Other notable changes in samtools:
+
+ * Improved sorting order checking during indexing.
+
+ * Improved region parsing. Colons in reference sequence names are parsed
+ properly.
+
+ * Fixed an issue where mpileup does not apply BAQ for the first few reads when
+ a region is specified.
+
+ * Fixed an issue where `faidx' does not work with FASTA files with long lines.
+
+ * Bugfix: wrong SP genotype information in the BCF output.
+
+Other notable changes in bcftools:
+
+ * Output the ML esitmate of the allele count.
+
+ * Added the HWE plus F<0 filter to varFilter. For multiple samples, it
+ effectively filters false heterozygous calls around centromeres.
+
+ * For association mapping, perform both 1-degree and 2-degree test. The
+ 2-degree test is conservative but more robust to HWE violation.
+
+(0.1.17: 6 July 2011, r973:277)
+
+
+
Beta Release 0.1.16 (21 April, 2011)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@copyright Genome Research Ltd.
*/
-#define BAM_VERSION "0.1.16 (r963:234)"
+#define BAM_VERSION "0.1.17 (r973:277)"
#include <stdint.h>
#include <stdlib.h>
memcpy(b->gi[0].data, bc->PL, b->gi[0].len * bc->n);
if (bcr) {
uint16_t *dp = (uint16_t*)b->gi[1].data;
- uint8_t *sp = is_SP? b->gi[2].data : 0;
+ int32_t *sp = is_SP? b->gi[2].data : 0;
for (i = 0; i < bc->n; ++i) {
bcf_callret1_t *p = bcr + i;
dp[i] = p->depth < 0xffff? p->depth : 0xffff;
#include <string.h>
#include "bam.h"
#include "bam2bcf.h"
-#include "ksort.h"
#include "kaln.h"
#include "kprobaln.h"
#include "khash.h"
KHASH_SET_INIT_STR(rg)
+#include "ksort.h"
+KSORT_INIT_GENERIC(uint32_t)
+
#define MINUS_CONST 0x10000000
#define INDEL_WINDOW_SIZE 50
int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref,
const void *rghash)
{
- extern void ks_introsort_uint32_t(int, uint32_t*);
int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score1, *score2, max_ref2;
int N, K, l_run, ref_type, n_alt;
char *inscns = 0, *ref2, *query, **ref_sample;
if (n_types == 1 || (double)n_alt / n_tot < bca->min_frac || n_alt < bca->min_support) { // then skip
free(aux); return -1;
}
+ if (n_types >= 64) {
+ free(aux);
+ if (bam_verbose >= 2)
+ fprintf(stderr, "[%s] excessive INDEL alleles at position %d. Skip the position.\n", __func__, pos + 1);
+ return -1;
+ }
types = (int*)calloc(n_types, sizeof(int));
t = 0;
types[t++] = aux[0] - MINUS_CONST;
for (t = 0; t < n_types; ++t)
if (types[t] == 0) break;
ref_type = t; // the index of the reference type (0)
- assert(n_types < 64);
}
{ // calculate left and right boundary
left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;
return k == kh_end(h)? -1 : kh_value(h, k);
}
-int bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *begin, int *end)
+int bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *beg, int *end)
{
- char *s, *p;
- int i, l, k;
+ char *s;
+ int i, l, k, name_end;
khiter_t iter;
khash_t(s) *h;
bam_init_header_hash(header);
h = (khash_t(s)*)header->hash;
- l = strlen(str);
- p = s = (char*)malloc(l+1);
- /* squeeze out "," */
- for (i = k = 0; i != l; ++i)
- if (str[i] != ',' && !isspace(str[i])) s[k++] = str[i];
- s[k] = 0;
- for (i = 0; i != k; ++i) if (s[i] == ':') break;
- s[i] = 0;
- iter = kh_get(s, h, s); /* get the ref_id */
- if (iter == kh_end(h)) { // name not found
- *ref_id = -1; free(s);
- return -1;
- }
- *ref_id = kh_value(h, iter);
- if (i == k) { /* dump the whole sequence */
- *begin = 0; *end = 1<<29; free(s);
- return 0;
- }
- for (p = s + i + 1; i != k; ++i) if (s[i] == '-') break;
- *begin = atoi(p);
- if (i < k) {
- p = s + i + 1;
- *end = atoi(p);
- } else *end = 1<<29;
- if (*begin > 0) --*begin;
+ *ref_id = *beg = *end = -1;
+ name_end = l = strlen(str);
+ s = (char*)malloc(l+1);
+ // remove space
+ for (i = k = 0; i < l; ++i)
+ if (!isspace(str[i])) s[k++] = str[i];
+ s[k] = 0; l = k;
+ // determine the sequence name
+ for (i = l - 1; i >= 0; --i) if (s[i] == ':') break; // look for colon from the end
+ if (i >= 0) name_end = i;
+ if (name_end < l) { // check if this is really the end
+ int n_hyphen = 0;
+ for (i = name_end + 1; i < l; ++i) {
+ if (s[i] == '-') ++n_hyphen;
+ else if (!isdigit(s[i]) && s[i] != ',') break;
+ }
+ if (i < l || n_hyphen > 1) name_end = l; // malformated region string; then take str as the name
+ s[name_end] = 0;
+ iter = kh_get(s, h, s);
+ if (iter == kh_end(h)) { // cannot find the sequence name
+ iter = kh_get(s, h, str); // try str as the name
+ if (iter == kh_end(h)) {
+ if (bam_verbose >= 2) fprintf(stderr, "[%s] fail to determine the sequence name.\n", __func__);
+ free(s); return -1;
+ } else s[name_end] = ':', name_end = l;
+ }
+ } else iter = kh_get(s, h, str);
+ *ref_id = kh_val(h, iter);
+ // parse the interval
+ if (name_end < l) {
+ for (i = k = name_end + 1; i < l; ++i)
+ if (s[i] != ',') s[k++] = s[i];
+ s[k] = 0;
+ *beg = atoi(s + name_end + 1);
+ for (i = name_end + 1; i != k; ++i) if (s[i] == '-') break;
+ *end = i < k? atoi(s + i + 1) : 1<<29;
+ if (*beg > 0) --*beg;
+ } else *beg = 0, *end = 1<<29;
free(s);
- if (*begin > *end) {
- fprintf(stderr, "[bam_parse_region] invalid region.\n");
- return -1;
- }
- return 0;
+ return *beg <= *end? 0 : -1;
}
int32_t bam_aux2i(const uint8_t *s)
else if (str->s[5] == 'S') while (p < str->s + str->l) ((uint16_t*)s)[k++] = (uint16_t)strtol(p, &p, 0), ++p;
else if (str->s[5] == 'i') while (p < str->s + str->l) ((int32_t*)s)[k++] = (int32_t)strtol(p, &p, 0), ++p;
else if (str->s[5] == 'I') while (p < str->s + str->l) ((uint32_t*)s)[k++] = (uint32_t)strtol(p, &p, 0), ++p;
- else if (str->s[5] == 'f') while (p < str->s + str->l) ((float*)s)[k++] = (float)strtof(p, &p), ++p;
+ else if (str->s[5] == 'f') while (p < str->s + str->l) ((float*)s)[k++] = (float)strtod(p, &p), ++p;
else parse_error(fp->n_lines, "unrecognized array type");
s += Bsize * n; doff += size;
} else parse_error(fp->n_lines, "unrecognized type");
save_bin = save_tid = last_tid = last_bin = 0xffffffffu;
save_off = last_off = bam_tell(fp); last_coor = 0xffffffffu;
- n_mapped = n_unmapped = n_no_coor = off_end = 0;
+ n_mapped = n_unmapped = n_no_coor = off_end = 0;
off_beg = off_end = bam_tell(fp);
while ((ret = bam_read1(fp, b)) >= 0) {
if (c->tid < 0) ++n_no_coor;
- if (last_tid != c->tid) { // change of chromosomes
+ if (last_tid < c->tid || (last_tid >= 0 && c->tid < 0)) { // change of chromosomes
last_tid = c->tid;
last_bin = 0xffffffffu;
- } else if (last_coor > c->pos) {
+ } else if ((uint32_t)last_tid > (uint32_t)c->tid) {
+ fprintf(stderr, "[bam_index_core] the alignment is not sorted (%s): %d-th chr > %d-th chr\n",
+ bam1_qname(b), last_tid+1, c->tid+1);
+ return NULL;
+ } else if ((int32_t)c->tid >= 0 && last_coor > c->pos) {
fprintf(stderr, "[bam_index_core] the alignment is not sorted (%s): %u > %u in %d-th chr\n",
bam1_qname(b), last_coor, c->pos, c->tid+1);
- exit(1);
+ return NULL;
}
if (c->tid >= 0) insert_offset2(&idx->index2[b->core.tid], b, last_off);
if (c->bin != last_bin) { // then possibly write the binning index
if (bam_tell(fp) <= last_off) {
fprintf(stderr, "[bam_index_core] bug in BGZF/RAZF: %llx < %llx\n",
(unsigned long long)bam_tell(fp), (unsigned long long)last_off);
- exit(1);
+ return NULL;
}
if (c->flag & BAM_FUNMAP) ++n_unmapped;
else ++n_mapped;
++n_no_coor;
if (c->tid >= 0 && n_no_coor) {
fprintf(stderr, "[bam_index_core] the alignment is not sorted: reads without coordinates prior to reads with coordinates.\n");
- exit(1);
+ return NULL;
}
}
}
}
idx = bam_index_core(fp);
bam_close(fp);
+ if(idx == 0) {
+ fprintf(stderr, "[bam_index_build2] fail to index the BAM file.\n");
+ return -1;
+ }
if (_fnidx == 0) {
fnidx = (char*)calloc(strlen(fn) + 5, 1);
strcpy(fnidx, fn); strcat(fnidx, ".bai");
+++ /dev/null
-#include <math.h>
-#include <assert.h>
-#include "bam.h"
-#include "bam_maqcns.h"
-#include "ksort.h"
-#include "errmod.h"
-#include "kaln.h"
-KSORT_INIT_GENERIC(uint32_t)
-
-#define INDEL_WINDOW_SIZE 50
-#define INDEL_EXT_DEP 0.9
-
-typedef struct __bmc_aux_t {
- int max;
- uint32_t *info;
- uint16_t *info16;
- errmod_t *em;
-} bmc_aux_t;
-
-typedef struct {
- float esum[4], fsum[4];
- uint32_t c[4];
-} glf_call_aux_t;
-
-/*
- 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 = 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;
-
- free(aa->lhet);
- aa->lhet = (double*)calloc(256 * 256, sizeof(double));
- sum_harmo = 0.0;
- for (k = 1; k <= aa->n_hap - 1; ++k)
- sum_harmo += 1.0 / k;
- for (n1 = 0; n1 < 256; ++n1) {
- for (n2 = 0; n2 < 256; ++n2) {
- long double sum = 0.0;
- double lC = aa->errmod == BAM_ERRMOD_SOAP? 0 : lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1);
- for (k = 1; k <= aa->n_hap - 1; ++k) {
- double pk = 1.0 / k / sum_harmo;
- double log1 = log((double)k/aa->n_hap);
- double log2 = log(1.0 - (double)k/aa->n_hap);
- sum += pk * 0.5 * (expl(log1*n2) * expl(log2*n1) + expl(log1*n1) * expl(log2*n2));
- }
- aa->lhet[n1<<8|n2] = lC + logl(sum);
- }
- }
- poly_rate = aa->het_rate * sum_harmo;
- aa->q_r = -4.343 * log(2.0 * poly_rate / (1.0 - poly_rate));
-}
-
-/** initialize the helper structure */
-static void cal_coef(bam_maqcns_t *aa)
-{
- int k, n, q;
- long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256];
- double *lC;
-
- if (aa->errmod == BAM_ERRMOD_MAQ2) return; // no need to do the following
- // aa->lhet will be allocated and initialized
- free(aa->fk); free(aa->coef);
- aa->coef = 0;
- aa->fk = (double*)calloc(256, 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->errmod == BAM_ERRMOD_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);
- for (q = 1; q != 64; ++q) {
- double e = pow(10.0, -q/10.0);
- double le = log(e);
- double le1 = log(1.0-e);
- for (n = 1; n != 256; ++n) {
- double *coef = aa->coef + (q<<16|n<<8);
- sum_a[n+1] = 0.0;
- for (k = n; k >= 0; --k) { // a_k = \sum_{i=k}^n C^n_k \epsilon^k (1-\epsilon)^{n-k}
- sum_a[k] = sum_a[k+1] + expl(lC[n<<8|k] + k*le + (n-k)*le1);
- b[k] = sum_a[k+1] / sum_a[k];
- if (b[k] > 0.99) b[k] = 0.99;
- }
- for (k = 0; k != n; ++k) // log(\bar\beta_{nk}(\bar\epsilon)^{f_k})
- q_c[k] = -4.343 * fk2[k] * logl(b[k] / e);
- for (k = 1; k != n; ++k) q_c[k] += q_c[k-1]; // \prod_{i=0}^k c_i
- for (k = 0; k <= n; ++k) { // powl() in 64-bit mode seems broken on my Mac OS X 10.4.9
- tmp[k] = -4.343 * logl(1.0 - expl(fk2[k] * logl(b[k])));
- coef[k] = (k? q_c[k-1] : 0) + tmp[k]; // this is the final c_{nk}
- }
- }
- }
- free(lC);
-}
-
-bam_maqcns_t *bam_maqcns_init()
-{
- bam_maqcns_t *bm;
- bm = (bam_maqcns_t*)calloc(1, sizeof(bam_maqcns_t));
- bm->aux = (bmc_aux_t*)calloc(1, sizeof(bmc_aux_t));
- bm->het_rate = 0.001;
- bm->theta = 0.83f;
- bm->n_hap = 2;
- bm->eta = 0.03;
- bm->cap_mapQ = 60;
- bm->min_baseQ = 13;
- return bm;
-}
-
-void bam_maqcns_prepare(bam_maqcns_t *bm)
-{
- if (bm->errmod == BAM_ERRMOD_MAQ2) bm->aux->em = errmod_init(1. - bm->theta);
- cal_coef(bm); cal_het(bm);
-}
-
-void bam_maqcns_destroy(bam_maqcns_t *bm)
-{
- if (bm == 0) return;
- free(bm->lhet); free(bm->fk); free(bm->coef); free(bm->aux->info); free(bm->aux->info16);
- if (bm->aux->em) errmod_destroy(bm->aux->em);
- free(bm->aux); free(bm);
-}
-
-glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam_maqcns_t *bm)
-{
- glf_call_aux_t *b = 0;
- int i, j, k, w[8], c, n;
- glf1_t *g = (glf1_t*)calloc(1, sizeof(glf1_t));
- float p[16], min_p = 1e30;
- uint64_t rms;
-
- g->ref_base = ref_base;
- if (_n == 0) return g;
-
- // construct aux array
- if (bm->aux->max < _n) {
- bm->aux->max = _n;
- kroundup32(bm->aux->max);
- bm->aux->info = (uint32_t*)realloc(bm->aux->info, 4 * bm->aux->max);
- bm->aux->info16 = (uint16_t*)realloc(bm->aux->info16, 2 * bm->aux->max);
- }
- for (i = n = 0, rms = 0; i < _n; ++i) {
- const bam_pileup1_t *p = pl + i;
- uint32_t q, x = 0, qq;
- uint16_t y = 0;
- if (p->is_del || p->is_refskip || (p->b->core.flag&BAM_FUNMAP)) continue;
- q = (uint32_t)bam1_qual(p->b)[p->qpos];
- if (q < bm->min_baseQ) continue;
- x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
- y |= bam1_strand(p->b)<<4;
- if (p->b->core.qual < q) q = p->b->core.qual;
- c = p->b->core.qual < bm->cap_mapQ? p->b->core.qual : bm->cap_mapQ;
- rms += c * c;
- x |= q << 24;
- y |= q << 5;
- qq = bam1_seqi(bam1_seq(p->b), p->qpos);
- q = bam_nt16_nt4_table[qq? qq : ref_base];
- if (!p->is_del && !p->is_refskip && q < 4) x |= 1 << 21 | q << 16, y |= q;
- bm->aux->info16[n] = y;
- bm->aux->info[n++] = x;
- }
- rms = (uint8_t)(sqrt((double)rms / n) + .499);
- if (bm->errmod == BAM_ERRMOD_MAQ2) {
- errmod_cal(bm->aux->em, n, 4, bm->aux->info16, p);
- goto goto_glf;
- }
- ks_introsort(uint32_t, n, bm->aux->info);
- // generate esum and fsum
- b = (glf_call_aux_t*)calloc(1, sizeof(glf_call_aux_t));
- for (k = 0; k != 8; ++k) w[k] = 0;
- for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
- uint32_t info = bm->aux->info[j];
- if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
- k = info>>16&7;
- if (info>>24 > 0) {
- b->esum[k&3] += bm->fk[w[k]] * (info>>24);
- b->fsum[k&3] += bm->fk[w[k]];
- if (w[k] < 0xff) ++w[k];
- ++b->c[k&3];
- }
- }
- // rescale ->c[]
- for (j = c = 0; j != 4; ++j) c += b->c[j];
- if (c > 255) {
- 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];
- }
- if (bm->errmod == BAM_ERRMOD_MAQ) {
- // 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 (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;
- }
-
- { // 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;
- }
- } else if (bm->errmod == BAM_ERRMOD_SOAP) { // 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;
- }
- }
- }
-
-goto_glf:
- // convert necessary information to glf1_t
- g->ref_base = ref_base; g->max_mapQ = rms;
- g->depth = n > 16777215? 16777215 : n;
- for (j = 0; j != 4; ++j)
- for (k = j; k < 4; ++k)
- if (p[j<<2|k] < min_p) min_p = p[j<<2|k];
- g->min_lk = min_p > 255.0? 255 : (int)(min_p + 0.5);
- for (j = c = 0; j != 4; ++j)
- for (k = j; k < 4; ++k)
- g->lk[c++] = p[j<<2|k]-min_p > 255.0? 255 : (int)(p[j<<2|k]-min_p + 0.5);
-
- free(b);
- return g;
-}
-
-uint32_t glf2cns(const glf1_t *g, int q_r)
-{
- int i, j, k, p[10], ref4;
- uint32_t x = 0;
- ref4 = bam_nt16_nt4_table[g->ref_base];
- for (i = k = 0; i < 4; ++i)
- for (j = i; j < 4; ++j) {
- int prior = (i == ref4 && j == ref4? 0 : i == ref4 || j == ref4? q_r : q_r + 3);
- p[k] = (g->lk[k] + prior)<<4 | i<<2 | j;
- ++k;
- }
- for (i = 1; i < 10; ++i) // insertion sort
- for (j = i; j > 0 && p[j] < p[j-1]; --j)
- k = p[j], p[j] = p[j-1], p[j-1] = k;
- x = (1u<<(p[0]&3) | 1u<<(p[0]>>2&3)) << 28; // the best genotype
- x |= (uint32_t)g->max_mapQ << 16; // rms mapQ
- x |= ((p[1]>>4) - (p[0]>>4) < 256? (p[1]>>4) - (p[0]>>4) : 255) << 8; // consensus Q
- for (k = 0; k < 10; ++k)
- if ((p[k]&0xf) == (ref4<<2|ref4)) break;
- if (k == 10) k = 9;
- x |= (p[k]>>4) - (p[0]>>4) < 256? (p[k]>>4) - (p[0]>>4) : 255; // snp Q
- return x;
-}
-
-uint32_t bam_maqcns_call(int n, const bam_pileup1_t *pl, bam_maqcns_t *bm)
-{
- glf1_t *g;
- uint32_t x;
- if (n) {
- g = bam_maqcns_glfgen(n, pl, 0xf, bm);
- x = g->depth == 0? (0xfU<<28 | 0xfU<<24) : glf2cns(g, (int)(bm->q_r + 0.5));
- free(g);
- } else x = 0xfU<<28 | 0xfU<<24;
- return x;
-}
-
-/************** *****************/
-
-bam_maqindel_opt_t *bam_maqindel_opt_init()
-{
- bam_maqindel_opt_t *mi = (bam_maqindel_opt_t*)calloc(1, sizeof(bam_maqindel_opt_t));
- mi->q_indel = 40;
- mi->r_indel = 0.00015;
- mi->r_snp = 0.001;
- //
- mi->mm_penalty = 3;
- mi->indel_err = 4;
- mi->ambi_thres = 10;
- return mi;
-}
-
-void bam_maqindel_ret_destroy(bam_maqindel_ret_t *mir)
-{
- if (mir == 0) return;
- free(mir->s[0]); free(mir->s[1]); free(mir);
-}
-
-int bam_tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
-{
- int k, x = c->pos, y = 0, last_y = 0;
- *_tpos = c->pos;
- for (k = 0; k < c->n_cigar; ++k) {
- int op = cigar[k] & BAM_CIGAR_MASK;
- int l = cigar[k] >> BAM_CIGAR_SHIFT;
- if (op == BAM_CMATCH) {
- if (c->pos > tpos) return y;
- if (x + l > tpos) {
- *_tpos = tpos;
- return y + (tpos - x);
- }
- x += l; y += l;
- last_y = y;
- } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l;
- else if (op == BAM_CDEL || op == BAM_CREF_SKIP) {
- if (x + l > tpos) {
- *_tpos = is_left? x : x + l;
- return y;
- }
- x += l;
- }
- }
- *_tpos = x;
- return last_y;
-}
-
-#define MINUS_CONST 0x10000000
-
-bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
- int _n_types, int *_types)
-{
- int i, j, n_types, *types, left, right, max_rd_len = 0;
- bam_maqindel_ret_t *ret = 0;
- // if there is no proposed indel, check if there is an indel from the alignment
- if (_n_types == 0) {
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pl + i;
- if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
- }
- if (i == n) return 0; // no indel
- }
- { // calculate how many types of indels are available (set n_types and types)
- int m;
- uint32_t *aux;
- aux = (uint32_t*)calloc(n + _n_types + 1, 4);
- m = 0;
- aux[m++] = MINUS_CONST; // zero indel is always a type
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pl + i;
- if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0)
- aux[m++] = MINUS_CONST + p->indel;
- j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
- if (j > max_rd_len) max_rd_len = j;
- }
- if (_n_types) // then also add this to aux[]
- for (i = 0; i < _n_types; ++i)
- if (_types[i]) aux[m++] = MINUS_CONST + _types[i];
- ks_introsort(uint32_t, m, aux);
- // squeeze out identical types
- for (i = 1, n_types = 1; i < m; ++i)
- if (aux[i] != aux[i-1]) ++n_types;
- types = (int*)calloc(n_types, sizeof(int));
- j = 0;
- types[j++] = aux[0] - MINUS_CONST;
- for (i = 1; i < m; ++i) {
- if (aux[i] != aux[i-1])
- types[j++] = aux[i] - MINUS_CONST;
- }
- free(aux);
- }
- { // calculate left and right boundary
- left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;
- right = pos + INDEL_WINDOW_SIZE;
- if (types[0] < 0) right -= types[0];
- // in case the alignments stand out the reference
- for (i = pos; i < right; ++i)
- if (ref[i] == 0) break;
- right = i;
- }
- { // the core part
- char *ref2, *rs, *inscns = 0;
- int qr_snp, k, l, *score, *pscore, max_ins = types[n_types-1];
- qr_snp = (int)(-4.343 * log(mi->r_snp) + .499);
- if (max_ins > 0) { // get the consensus of inserted sequences
- int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int));
- // count occurrences
- for (i = 0; i < n_types; ++i) {
- if (types[i] <= 0) continue; // not insertion
- for (j = 0; j < n; ++j) {
- const bam_pileup1_t *p = pl + j;
- if (!(p->b->core.flag&BAM_FUNMAP) && p->indel == types[i]) {
- for (k = 1; k <= p->indel; ++k) {
- int c = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), p->qpos + k)];
- if (c < 4) ++inscns_aux[i*max_ins*4 + (k-1)*4 + c];
- }
- }
- }
- }
- // construct the consensus of inserted sequence
- inscns = (char*)calloc(n_types * max_ins, sizeof(char));
- for (i = 0; i < n_types; ++i) {
- for (j = 0; j < types[i]; ++j) {
- int max = 0, max_k = -1, *ia = inscns_aux + i*max_ins*4 + j*4;
- for (k = 0; k < 4; ++k) {
- if (ia[k] > max) {
- max = ia[k];
- max_k = k;
- }
- }
- inscns[i*max_ins + j] = max? 1<<max_k : 15;
- }
- }
- free(inscns_aux);
- }
- // calculate score
- ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1);
- rs = (char*)calloc(right - left + max_rd_len + types[n_types-1] + 2, 1);
- score = (int*)calloc(n_types * n, sizeof(int));
- pscore = (int*)calloc(n_types * n, sizeof(int));
- for (i = 0; i < n_types; ++i) {
- ka_param_t ap = ka_param_blast;
- ap.band_width = 2 * types[n_types - 1] + 2;
- ap.gap_end_ext = 0;
- // write ref2
- for (k = 0, j = left; j <= pos; ++j)
- ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
- if (types[i] <= 0) j += -types[i];
- else for (l = 0; l < types[i]; ++l)
- ref2[k++] = bam_nt16_nt4_table[(int)inscns[i*max_ins + l]];
- if (types[0] < 0) { // mask deleted sequences
- int jj, tmp = types[i] >= 0? -types[0] : -types[0] + types[i];
- for (jj = 0; jj < tmp && j < right && ref[j]; ++jj, ++j)
- ref2[k++] = 4;
- }
- for (; j < right && ref[j]; ++j)
- ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
- if (j < right) right = j;
- // calculate score for each read
- for (j = 0; j < n; ++j) {
- const bam_pileup1_t *p = pl + j;
- int qbeg, qend, tbeg, tend;
- if (p->b->core.flag & BAM_FUNMAP) continue;
- qbeg = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), left, 0, &tbeg);
- qend = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), right, 1, &tend);
- assert(tbeg >= left);
- for (l = qbeg; l < qend; ++l)
- rs[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), l)];
- {
- int x, y, n_acigar, ps;
- uint32_t *acigar;
- ps = 0;
- if (tend - tbeg + types[i] <= 0) {
- score[i*n+j] = -(1<<20);
- pscore[i*n+j] = 1<<20;
- continue;
- }
- acigar = ka_global_core((uint8_t*)ref2 + tbeg - left, tend - tbeg + types[i], (uint8_t*)rs, qend - qbeg, &ap, &score[i*n+j], &n_acigar);
- x = tbeg - left; y = 0;
- for (l = 0; l < n_acigar; ++l) {
- int op = acigar[l]&0xf;
- int len = acigar[l]>>4;
- if (op == BAM_CMATCH) {
- int k;
- for (k = 0; k < len; ++k)
- if (ref2[x+k] != rs[y+k] && ref2[x+k] < 4)
- ps += bam1_qual(p->b)[y+k] < qr_snp? bam1_qual(p->b)[y+k] : qr_snp;
- x += len; y += len;
- } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
- if (op == BAM_CINS && l > 0 && l < n_acigar - 1) ps += mi->q_indel * len;
- y += len;
- } else if (op == BAM_CDEL) {
- if (l > 0 && l < n_acigar - 1) ps += mi->q_indel * len;
- x += len;
- }
- }
- pscore[i*n+j] = ps;
- /*if (1) { // for debugging only
- fprintf(stderr, "id=%d, pos=%d, type=%d, j=%d, score=%d, psore=%d, %d, %d, %d, %d, %d, ",
- j, pos+1, types[i], j, score[i*n+j], pscore[i*n+j], tbeg, tend, qbeg, qend, mi->q_indel);
- for (l = 0; l < n_acigar; ++l) fprintf(stderr, "%d%c", acigar[l]>>4, "MIDS"[acigar[l]&0xf]);
- fprintf(stderr, "\n");
- for (l = 0; l < tend - tbeg + types[i]; ++l) fputc("ACGTN"[ref2[l+tbeg-left]], stderr);
- fputc('\n', stderr);
- for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[rs[l]], stderr);
- fputc('\n', stderr);
- }*/
- free(acigar);
- }
- }
- }
- { // get final result
- int *sum, max1, max2, max1_i, max2_i;
- // pick up the best two score
- sum = (int*)calloc(n_types, sizeof(int));
- for (i = 0; i < n_types; ++i)
- for (j = 0; j < n; ++j)
- sum[i] += -pscore[i*n+j];
- max1 = max2 = -0x7fffffff; max1_i = max2_i = -1;
- for (i = 0; i < n_types; ++i) {
- if (sum[i] > max1) {
- max2 = max1; max2_i = max1_i; max1 = sum[i]; max1_i = i;
- } else if (sum[i] > max2) {
- max2 = sum[i]; max2_i = i;
- }
- }
- free(sum);
- // write ret
- ret = (bam_maqindel_ret_t*)calloc(1, sizeof(bam_maqindel_ret_t));
- ret->indel1 = types[max1_i]; ret->indel2 = types[max2_i];
- ret->s[0] = (char*)calloc(abs(ret->indel1) + 2, 1);
- ret->s[1] = (char*)calloc(abs(ret->indel2) + 2, 1);
- // write indel sequence
- if (ret->indel1 > 0) {
- ret->s[0][0] = '+';
- for (k = 0; k < ret->indel1; ++k)
- ret->s[0][k+1] = bam_nt16_rev_table[(int)inscns[max1_i*max_ins + k]];
- } else if (ret->indel1 < 0) {
- ret->s[0][0] = '-';
- for (k = 0; k < -ret->indel1 && ref[pos + k + 1]; ++k)
- ret->s[0][k+1] = ref[pos + k + 1];
- } else ret->s[0][0] = '*';
- if (ret->indel2 > 0) {
- ret->s[1][0] = '+';
- for (k = 0; k < ret->indel2; ++k)
- ret->s[1][k+1] = bam_nt16_rev_table[(int)inscns[max2_i*max_ins + k]];
- } else if (ret->indel2 < 0) {
- ret->s[1][0] = '-';
- for (k = 0; k < -ret->indel2 && ref[pos + k + 1]; ++k)
- ret->s[1][k+1] = ref[pos + k + 1];
- } else ret->s[1][0] = '*';
- // write count
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pl + i;
- if (p->indel == ret->indel1) ++ret->cnt1;
- else if (p->indel == ret->indel2) ++ret->cnt2;
- else ++ret->cnt_anti;
- }
- { // write gl[]
- int tmp, seq_err = 0;
- double x = 1.0;
- tmp = max1_i - max2_i;
- if (tmp < 0) tmp = -tmp;
- for (j = 0; j < tmp + 1; ++j) x *= INDEL_EXT_DEP;
- seq_err = mi->q_indel * (1.0 - x) / (1.0 - INDEL_EXT_DEP);
- ret->gl[0] = ret->gl[1] = 0;
- for (j = 0; j < n; ++j) {
- int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j];
- //fprintf(stderr, "id=%d, %d, %d, %d, %d, %d\n", j, pl[j].b->core.pos+1, types[max1_i], types[max2_i], s1, s2);
- if (s1 > s2) ret->gl[0] += s1 - s2 < seq_err? s1 - s2 : seq_err;
- else ret->gl[1] += s2 - s1 < seq_err? s2 - s1 : seq_err;
- }
- }
- // write cnt_ref and cnt_ambi
- if (max1_i != 0 && max2_i != 0) {
- for (j = 0; j < n; ++j) {
- int diff1 = score[j] - score[max1_i * n + j];
- int diff2 = score[j] - score[max2_i * n + j];
- if (diff1 > 0 && diff2 > 0) ++ret->cnt_ref;
- else if (diff1 == 0 || diff2 == 0) ++ret->cnt_ambi;
- }
- }
- }
- free(score); free(pscore); free(ref2); free(rs); free(inscns);
- }
- { // call genotype
- int q[3], qr_indel = (int)(-4.343 * log(mi->r_indel) + 0.5);
- int min1, min2, min1_i;
- q[0] = ret->gl[0] + (ret->s[0][0] != '*'? 0 : 0) * qr_indel;
- q[1] = ret->gl[1] + (ret->s[1][0] != '*'? 0 : 0) * qr_indel;
- q[2] = n * 3 + (ret->s[0][0] == '*' || ret->s[1][0] == '*'? 1 : 1) * qr_indel;
- min1 = min2 = 0x7fffffff; min1_i = -1;
- for (i = 0; i < 3; ++i) {
- if (q[i] < min1) {
- min2 = min1; min1 = q[i]; min1_i = i;
- } else if (q[i] < min2) min2 = q[i];
- }
- ret->gt = min1_i;
- ret->q_cns = min2 - min1;
- // set q_ref
- if (ret->gt < 2) ret->q_ref = (ret->s[ret->gt][0] == '*')? 0 : q[1-ret->gt] - q[ret->gt] - qr_indel - 3;
- else ret->q_ref = (ret->s[0][0] == '*')? q[0] - q[2] : q[1] - q[2];
- if (ret->q_ref < 0) ret->q_ref = 0;
- }
- free(types);
- return ret;
-}
+++ /dev/null
-#ifndef BAM_MAQCNS_H
-#define BAM_MAQCNS_H
-
-#include "glf.h"
-
-#define BAM_ERRMOD_MAQ2 0
-#define BAM_ERRMOD_MAQ 1
-#define BAM_ERRMOD_SOAP 2
-
-struct __bmc_aux_t;
-
-typedef struct {
- float het_rate, theta;
- int n_hap, cap_mapQ, errmod, min_baseQ;
-
- float eta, q_r;
- double *fk, *coef;
- double *lhet;
- struct __bmc_aux_t *aux;
-} bam_maqcns_t;
-
-typedef struct {
- int q_indel; // indel sequencing error, phred scaled
- float r_indel; // indel prior
- float r_snp; // snp prior
- // hidden parameters, unchangeable from command line
- int mm_penalty, indel_err, ambi_thres;
-} bam_maqindel_opt_t;
-
-typedef struct {
- int indel1, indel2;
- int cnt1, cnt2, cnt_anti;
- int cnt_ref, cnt_ambi;
- char *s[2];
- //
- int gt, gl[2];
- int q_cns, q_ref;
-} bam_maqindel_ret_t;
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
- bam_maqcns_t *bam_maqcns_init();
- void bam_maqcns_prepare(bam_maqcns_t *bm);
- void bam_maqcns_destroy(bam_maqcns_t *bm);
- glf1_t *bam_maqcns_glfgen(int n, const bam_pileup1_t *pl, uint8_t ref_base, bam_maqcns_t *bm);
- uint32_t bam_maqcns_call(int n, const bam_pileup1_t *pl, bam_maqcns_t *bm);
- // return: cns<<28 | cns2<<24 | mapQ<<16 | cnsQ<<8 | cnsQ2
- uint32_t glf2cns(const glf1_t *g, int q_r);
-
- bam_maqindel_opt_t *bam_maqindel_opt_init();
- bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
- int _n_types, int *_types);
- void bam_maqindel_ret_destroy(bam_maqindel_ret_t*);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
#include <errno.h>
#include "sam.h"
#include "faidx.h"
-#include "bam_maqcns.h"
-#include "khash.h"
-#include "glf.h"
#include "kstring.h"
-typedef int *indel_list_t;
-KHASH_MAP_INIT_INT64(64, indel_list_t)
-
-#define BAM_PLF_SIMPLE 0x01
-#define BAM_PLF_CNS 0x02
-#define BAM_PLF_INDEL_ONLY 0x04
-#define BAM_PLF_GLF 0x08
-#define BAM_PLF_VAR_ONLY 0x10
-#define BAM_PLF_2ND 0x20
-#define BAM_PLF_RANBASE 0x40
-#define BAM_PLF_1STBASE 0x80
-#define BAM_PLF_ALLBASE 0x100
-#define BAM_PLF_READPOS 0x200
-#define BAM_PLF_NOBAQ 0x400
-
-typedef struct {
- bam_header_t *h;
- bam_maqcns_t *c;
- bam_maqindel_opt_t *ido;
- faidx_t *fai;
- khash_t(64) *hash;
- uint32_t format;
- int tid, len, last_pos;
- int mask;
- int capQ_thres, min_baseQ;
- int max_depth; // for indel calling, ignore reads with the depth too high. 0 for unlimited
- char *ref;
- glfFile fp_glf; // for glf output only
-} pu_data_t;
-
-char **__bam_get_lines(const char *fn, int *_n);
-void bam_init_header_hash(bam_header_t *header);
-int32_t bam_get_tid(const bam_header_t *header, const char *seq_name);
-
-static khash_t(64) *load_pos(const char *fn, bam_header_t *h)
-{
- char **list;
- int i, j, n, *fields, max_fields;
- khash_t(64) *hash;
- bam_init_header_hash(h);
- list = __bam_get_lines(fn, &n);
- hash = kh_init(64);
- max_fields = 0; fields = 0;
- for (i = 0; i < n; ++i) {
- char *str = list[i];
- int chr, n_fields, ret;
- khint_t k;
- uint64_t x;
- n_fields = ksplit_core(str, 0, &max_fields, &fields);
- if (n_fields < 2) continue;
- chr = bam_get_tid(h, str + fields[0]);
- if (chr < 0) {
- fprintf(stderr, "[load_pos] unknown reference sequence name: %s\n", str + fields[0]);
- continue;
- }
- x = (uint64_t)chr << 32 | (atoi(str + fields[1]) - 1);
- k = kh_put(64, hash, x, &ret);
- if (ret == 0) {
- fprintf(stderr, "[load_pos] position %s:%s has been loaded.\n", str+fields[0], str+fields[1]);
- continue;
- }
- kh_val(hash, k) = 0;
- if (n_fields > 2) {
- // count
- for (j = 2; j < n_fields; ++j) {
- char *s = str + fields[j];
- if ((*s != '+' && *s != '-') || !isdigit(s[1])) break;
- }
- if (j > 2) { // update kh_val()
- int *q, y, z;
- q = kh_val(hash, k) = (int*)calloc(j - 1, sizeof(int));
- q[0] = j - 2; z = j; y = 1;
- for (j = 2; j < z; ++j)
- q[y++] = atoi(str + fields[j]);
- }
- }
- free(str);
- }
- free(list); free(fields);
- return hash;
-}
-
static inline int printw(int c, FILE *fp)
{
char buf[16];
return 0;
}
-// an analogy to pileup_func() below
-static int glt3_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pu, void *data)
-{
- pu_data_t *d = (pu_data_t*)data;
- bam_maqindel_ret_t *r = 0;
- int rb, *proposed_indels = 0;
- glf1_t *g;
- glf3_t *g3;
-
- if (d->fai == 0) {
- fprintf(stderr, "[glt3_func] reference sequence is required for generating GLT. Abort!\n");
- exit(1);
- }
- if (d->hash) { // only output a list of sites
- khint_t k = kh_get(64, d->hash, (uint64_t)tid<<32|pos);
- if (k == kh_end(d->hash)) return 0;
- proposed_indels = kh_val(d->hash, k);
- }
- g3 = glf3_init1();
- if (d->fai && (int)tid != d->tid) {
- if (d->ref) { // then write the end mark
- g3->rtype = GLF3_RTYPE_END;
- glf3_write1(d->fp_glf, g3);
- }
- glf3_ref_write(d->fp_glf, d->h->target_name[tid], d->h->target_len[tid]); // write reference
- free(d->ref);
- d->ref = fai_fetch(d->fai, d->h->target_name[tid], &d->len);
- d->tid = tid;
- d->last_pos = 0;
- }
- rb = (d->ref && (int)pos < d->len)? d->ref[pos] : 'N';
- g = bam_maqcns_glfgen(n, pu, bam_nt16_table[rb], d->c);
- memcpy(g3, g, sizeof(glf1_t));
- g3->rtype = GLF3_RTYPE_SUB;
- g3->offset = pos - d->last_pos;
- d->last_pos = pos;
- glf3_write1(d->fp_glf, g3);
- if (pos < d->len) {
- int m = (!d->max_depth || d->max_depth>n) ? n : d->max_depth;
- if (proposed_indels)
- r = bam_maqindel(m, pos, d->ido, pu, d->ref, proposed_indels[0], proposed_indels+1);
- else r = bam_maqindel(m, pos, d->ido, pu, d->ref, 0, 0);
- }
- if (r) { // then write indel line
- int het = 3 * n, min;
- min = het;
- if (min > r->gl[0]) min = r->gl[0];
- if (min > r->gl[1]) min = r->gl[1];
- g3->ref_base = 0;
- g3->rtype = GLF3_RTYPE_INDEL;
- memset(g3->lk, 0, 10);
- g3->lk[0] = r->gl[0] - min < 255? r->gl[0] - min : 255;
- g3->lk[1] = r->gl[1] - min < 255? r->gl[1] - min : 255;
- g3->lk[2] = het - min < 255? het - min : 255;
- g3->offset = 0;
- g3->indel_len[0] = r->indel1;
- g3->indel_len[1] = r->indel2;
- g3->min_lk = min < 255? min : 255;
- g3->max_len = (abs(r->indel1) > abs(r->indel2)? abs(r->indel1) : abs(r->indel2)) + 1;
- g3->indel_seq[0] = strdup(r->s[0]+1);
- g3->indel_seq[1] = strdup(r->s[1]+1);
- glf3_write1(d->fp_glf, g3);
- bam_maqindel_ret_destroy(r);
- }
- free(g);
- glf3_destroy1(g3);
- return 0;
-}
-
static inline void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, const char *ref)
{
int j;
if (p->is_tail) putchar('$');
}
-static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pu, void *data)
-{
- pu_data_t *d = (pu_data_t*)data;
- bam_maqindel_ret_t *r = 0;
- int i, rb, rms_mapq = -1, *proposed_indels = 0;
- uint64_t rms_aux;
- uint32_t cns = 0;
-
- // if GLF is required, suppress -c completely
- if (d->format & BAM_PLF_GLF) return glt3_func(tid, pos, n, pu, data);
- // if d->hash is initialized, only output the sites in the hash table
- if (d->hash) {
- khint_t k = kh_get(64, d->hash, (uint64_t)tid<<32|pos);
- if (k == kh_end(d->hash)) return 0;
- proposed_indels = kh_val(d->hash, k);
- }
- // update d->ref if necessary
- if (d->fai && (int)tid != d->tid) {
- free(d->ref);
- d->ref = faidx_fetch_seq(d->fai, d->h->target_name[tid], 0, 0x7fffffff, &d->len);
- d->tid = tid;
- }
- rb = (d->ref && (int)pos < d->len)? d->ref[pos] : 'N';
- // when the indel-only mode is asked for, return if no reads mapped with indels
- if (d->format & BAM_PLF_INDEL_ONLY) {
- for (i = 0; i < n; ++i)
- if (pu[i].indel != 0) break;
- if (i == n) return 0;
- }
- // call the consensus and indel
- if (d->format & BAM_PLF_CNS) { // call consensus
- if (d->format & (BAM_PLF_RANBASE|BAM_PLF_1STBASE)) { // use a random base or the 1st base as the consensus call
- const bam_pileup1_t *p = (d->format & BAM_PLF_1STBASE)? pu : pu + (int)(drand48() * n);
- int q = bam1_qual(p->b)[p->qpos];
- int mapQ = p->b->core.qual < d->c->cap_mapQ? p->b->core.qual : d->c->cap_mapQ;
- uint32_t b = bam1_seqi(bam1_seq(p->b), p->qpos);
- cns = b<<28 | 0xf<<24 | mapQ<<16 | q<<8;
- } else if (d->format & BAM_PLF_ALLBASE) { // collapse all bases
- uint64_t rmsQ = 0;
- uint32_t b = 0;
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pu + i;
- int q = p->b->core.qual < d->c->cap_mapQ? p->b->core.qual : d->c->cap_mapQ;
- b |= bam1_seqi(bam1_seq(p->b), p->qpos);
- rmsQ += q * q;
- }
- rmsQ = (uint64_t)(sqrt((double)rmsQ / n) + .499);
- cns = b<<28 | 0xf<<24 | rmsQ<<16 | 60<<8;
- } else {
- glf1_t *g = bam_maqcns_glfgen(n, pu, bam_nt16_table[rb], d->c);
- cns = g->depth == 0? (0xfu<<28 | 0xf<<24) : glf2cns(g, (int)(d->c->q_r + .499));
- free(g);
- }
- }
- if ((d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)) && d->ref && pos < d->len) { // call indels
- int m = (!d->max_depth || d->max_depth>n) ? n : d->max_depth;
- if (proposed_indels) // the first element gives the size of the array
- r = bam_maqindel(m, pos, d->ido, pu, d->ref, proposed_indels[0], proposed_indels+1);
- else r = bam_maqindel(m, pos, d->ido, pu, d->ref, 0, 0);
- }
- // when only variant sites are asked for, test if the site is a variant
- if ((d->format & BAM_PLF_CNS) && (d->format & BAM_PLF_VAR_ONLY)) {
- if (!(bam_nt16_table[rb] != 15 && cns>>28 != 15 && cns>>28 != bam_nt16_table[rb])) { // not a SNP
- if (!(r && (r->gt == 2 || strcmp(r->s[r->gt], "*")))) { // not an indel
- if (r) bam_maqindel_ret_destroy(r);
- return 0;
- }
- }
- }
- // print the first 3 columns
- fputs(d->h->target_name[tid], stdout); putchar('\t');
- printw(pos+1, stdout); putchar('\t'); putchar(rb); putchar('\t');
- // print consensus information if required
- if (d->format & BAM_PLF_CNS) {
- putchar(bam_nt16_rev_table[cns>>28]); putchar('\t');
- printw(cns>>8&0xff, stdout); putchar('\t');
- printw(cns&0xff, stdout); putchar('\t');
- printw(cns>>16&0xff, stdout); putchar('\t');
- }
- // print pileup sequences
- printw(n, stdout); putchar('\t');
- for (i = 0; i < n; ++i)
- pileup_seq(pu + i, pos, d->len, d->ref);
- // finalize rms_mapq
- if (d->format & BAM_PLF_CNS) {
- for (i = rms_aux = 0; i < n; ++i) {
- const bam_pileup1_t *p = pu + i;
- int tmp = p->b->core.qual < d->c->cap_mapQ? p->b->core.qual : d->c->cap_mapQ;
- rms_aux += tmp * tmp;
- }
- rms_aux = (uint64_t)(sqrt((double)rms_aux / n) + .499);
- if (rms_mapq < 0) rms_mapq = rms_aux;
- }
- putchar('\t');
- // print quality
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pu + i;
- int c = bam1_qual(p->b)[p->qpos] + 33;
- if (c > 126) c = 126;
- putchar(c);
- }
- if (d->format & BAM_PLF_2ND) { // print 2nd calls and qualities
- const unsigned char *q;
- putchar('\t');
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pu + i;
- q = bam_aux_get(p->b, "E2");
- putchar(q? q[p->qpos + 1] : 'N');
- }
- putchar('\t');
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pu + i;
- q = bam_aux_get(p->b, "U2");
- putchar(q? q[p->qpos + 1] : '!');
- }
- }
- // print mapping quality if -s is flagged on the command line
- if (d->format & BAM_PLF_SIMPLE) {
- putchar('\t');
- for (i = 0; i < n; ++i) {
- int c = pu[i].b->core.qual + 33;
- if (c > 126) c = 126;
- putchar(c);
- }
- }
- // print read position
- if (d->format & BAM_PLF_READPOS) {
- putchar('\t');
- for (i = 0; i < n; ++i) {
- int x = pu[i].qpos;
- int l = pu[i].b->core.l_qseq;
- printw(x < l/2? x+1 : -((l-1)-x+1), stdout); putchar(',');
- }
- }
- putchar('\n');
- // print the indel line if r has been calculated. This only happens if:
- // a) -c or -i are flagged, AND b) the reference sequence is available
- if (r) {
- printf("%s\t%d\t*\t", d->h->target_name[tid], pos + 1);
- if (r->gt < 2) printf("%s/%s\t", r->s[r->gt], r->s[r->gt]);
- else printf("%s/%s\t", r->s[0], r->s[1]);
- printf("%d\t%d\t", r->q_cns, r->q_ref);
- printf("%d\t%d\t", rms_mapq, n);
- printf("%s\t%s\t", r->s[0], r->s[1]);
- //printf("%d\t%d\t", r->gl[0], r->gl[1]);
- printf("%d\t%d\t%d\t", r->cnt1, r->cnt2, r->cnt_anti);
- printf("%d\t%d\n", r->cnt_ref, r->cnt_ambi);
- bam_maqindel_ret_destroy(r);
- }
- return 0;
-}
-
-int bam_pileup(int argc, char *argv[])
-{
- int c, is_SAM = 0;
- char *fn_list = 0, *fn_fa = 0, *fn_pos = 0;
- pu_data_t *d = (pu_data_t*)calloc(1, sizeof(pu_data_t));
- d->max_depth = 1024; d->tid = -1; d->mask = BAM_DEF_MASK; d->min_baseQ = 13;
- d->c = bam_maqcns_init();
- d->c->errmod = BAM_ERRMOD_MAQ2; // change the default model
- d->ido = bam_maqindel_opt_init();
- while ((c = getopt(argc, argv, "st:f:cT:N:r:l:d:im:gI:G:vM:S2aR:PAQ:C:B")) >= 0) {
- switch (c) {
- case 'Q': d->c->min_baseQ = atoi(optarg); break;
- case 'C': d->capQ_thres = atoi(optarg); break;
- case 'B': d->format |= BAM_PLF_NOBAQ; break;
- case 'a': d->c->errmod = BAM_ERRMOD_SOAP; break;
- case 'A': d->c->errmod = BAM_ERRMOD_MAQ; break;
- case 's': d->format |= BAM_PLF_SIMPLE; break;
- case 't': fn_list = strdup(optarg); break;
- case 'l': fn_pos = strdup(optarg); break;
- case 'f': fn_fa = strdup(optarg); break;
- case 'T': d->c->theta = atof(optarg); break;
- case 'N': d->c->n_hap = atoi(optarg); break;
- case 'r': d->c->het_rate = atof(optarg); d->ido->r_snp = d->c->het_rate; break;
- case 'M': d->c->cap_mapQ = atoi(optarg); break;
- case 'd': d->max_depth = atoi(optarg); break;
- case 'c': d->format |= BAM_PLF_CNS; break;
- case 'i': d->format |= BAM_PLF_INDEL_ONLY; break;
- case 'v': d->format |= BAM_PLF_VAR_ONLY; break;
- case 'm': d->mask = strtol(optarg, 0, 0); break;
- case 'g': d->format |= BAM_PLF_GLF; break;
- case '2': d->format |= BAM_PLF_2ND; break;
- case 'P': d->format |= BAM_PLF_READPOS; break;
- case 'I': d->ido->q_indel = atoi(optarg); break;
- case 'G': d->ido->r_indel = atof(optarg); break;
- case 'S': is_SAM = 1; break;
- case 'R':
- if (strcmp(optarg, "random") == 0) d->format |= BAM_PLF_RANBASE;
- else if (strcmp(optarg, "first") == 0) d->format |= BAM_PLF_1STBASE;
- else if (strcmp(optarg, "all") == 0) d->format |= BAM_PLF_ALLBASE;
- else fprintf(stderr, "[bam_pileup] unrecognized -R\n");
- break;
- default: fprintf(stderr, "Unrecognizd option '-%c'.\n", c); return 1;
- }
- }
- if (d->c->errmod != BAM_ERRMOD_MAQ2) d->c->theta += 0.02;
- if (d->c->theta > 1.0) d->c->theta = 1.0;
- if (fn_list) is_SAM = 1;
- if (optind == argc) {
- fprintf(stderr, "\n");
- 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, " -B disable BAQ computation\n");
- fprintf(stderr, " -A use the original MAQ model for SNP calling (DEPRECATED)\n");
- fprintf(stderr, " -2 output the 2nd best call and quality\n");
- fprintf(stderr, " -i only show lines/consensus with indels\n");
- fprintf(stderr, " -Q INT min base quality (possibly capped by BAQ) [%d]\n", d->c->min_baseQ);
- fprintf(stderr, " -C INT coefficient for adjusting mapQ of poor mappings [%d]\n", d->capQ_thres);
- fprintf(stderr, " -m INT filtering reads with bits in INT [0x%x]\n", d->mask);
- fprintf(stderr, " -M INT cap mapping quality at INT [%d]\n", d->c->cap_mapQ);
- fprintf(stderr, " -d INT limit maximum depth for indels [%d]\n", d->max_depth);
- fprintf(stderr, " -t FILE list of reference sequences (force -S)\n");
- fprintf(stderr, " -l FILE list of sites at which pileup is output\n");
- fprintf(stderr, " -f FILE reference sequence in the FASTA format\n\n");
- fprintf(stderr, " -c compute the consensus sequence\n");
- fprintf(stderr, " -v print variants only (for -c)\n");
- fprintf(stderr, " -g output in the GLFv3 format (DEPRECATED)\n");
- fprintf(stderr, " -T FLOAT theta in maq consensus calling model (for -c) [%.4g]\n", d->c->theta);
- fprintf(stderr, " -N INT number of haplotypes in the sample (for -c) [%d]\n", d->c->n_hap);
- fprintf(stderr, " -r FLOAT prior of a difference between two haplotypes (for -c) [%.4g]\n", d->c->het_rate);
- fprintf(stderr, " -G FLOAT prior of an indel between two haplotypes (for -c) [%.4g]\n", d->ido->r_indel);
- fprintf(stderr, " -I INT phred prob. of an indel in sequencing/prep. (for -c) [%d]\n", d->ido->q_indel);
- fprintf(stderr, "\n");
- fprintf(stderr, "Warning: Please use the `mpileup' command instead `pileup'. `Pileup' is deprecated!\n\n");
- free(fn_list); free(fn_fa); free(d);
- return 1;
- }
- if (d->format & (BAM_PLF_RANBASE|BAM_PLF_1STBASE|BAM_PLF_ALLBASE)) d->format |= BAM_PLF_CNS;
- if (fn_fa) d->fai = fai_load(fn_fa);
- if (d->format & (BAM_PLF_CNS|BAM_PLF_GLF)) bam_maqcns_prepare(d->c); // consensus calling
- if (d->format & BAM_PLF_GLF) { // for glf output
- glf3_header_t *h;
- h = glf3_header_init();
- d->fp_glf = bgzf_fdopen(fileno(stdout), "w");
- glf3_header_write(d->fp_glf, h);
- glf3_header_destroy(h);
- }
- if (d->fai == 0 && (d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)))
- fprintf(stderr, "[bam_pileup] indels will not be called when -f is absent.\n");
- if (fn_fa && is_SAM && fn_list == 0) fn_list = samfaipath(fn_fa);
-
- {
- samfile_t *fp;
- fp = is_SAM? samopen(argv[optind], "r", fn_list) : samopen(argv[optind], "rb", 0);
- if (fp == 0 || fp->header == 0) {
- fprintf(stderr, "[bam_pileup] fail to read the header: non-exisiting file or wrong format.\n");
- return 1;
- }
- d->h = fp->header;
- if (fn_pos) d->hash = load_pos(fn_pos, d->h);
- { // run pileup
- extern int bam_prob_realn(bam1_t *b, const char *ref);
- extern int bam_cap_mapQ(bam1_t *b, char *ref, int thres);
- bam1_t *b;
- int ret, tid, pos, n_plp;
- bam_plp_t iter;
- const bam_pileup1_t *plp;
- b = bam_init1();
- iter = bam_plp_init(0, 0);
- bam_plp_set_mask(iter, d->mask);
- while ((ret = samread(fp, b)) >= 0) {
- int skip = 0;
- if ((int)b->core.tid < 0) break;
- // update d->ref if necessary
- if (d->fai && (int)b->core.tid != d->tid) {
- free(d->ref);
- d->ref = faidx_fetch_seq(d->fai, d->h->target_name[b->core.tid], 0, 0x7fffffff, &d->len);
- d->tid = b->core.tid;
- }
- if (d->ref && (d->format&BAM_PLF_CNS) && !(d->format&BAM_PLF_NOBAQ)) bam_prob_realn(b, d->ref);
- if (d->ref && (d->format&BAM_PLF_CNS) && d->capQ_thres > 10) {
- int q = bam_cap_mapQ(b, d->ref, d->capQ_thres);
- if (q < 0) skip = 1;
- else if (b->core.qual > q) b->core.qual = q;
- } else if (b->core.flag&BAM_FUNMAP) skip = 1;
- else if ((d->format&BAM_PLF_CNS) && (b->core.flag&1) && !(b->core.flag&2)) skip = 1;
- if (skip) continue;
- bam_plp_push(iter, b);
- while ((plp = bam_plp_next(iter, &tid, &pos, &n_plp)) != 0)
- pileup_func(tid, pos, n_plp, plp, d);
- }
- bam_plp_push(iter, 0);
- while ((plp = bam_plp_next(iter, &tid, &pos, &n_plp)) != 0)
- pileup_func(tid, pos, n_plp, plp, d);
- bam_plp_destroy(iter);
- bam_destroy1(b);
- }
- samclose(fp); // d->h will be destroyed here
- }
-
- // free
- if (d->format & BAM_PLF_GLF) bgzf_close(d->fp_glf);
- if (fn_pos) { // free the hash table
- khint_t k;
- for (k = kh_begin(d->hash); k < kh_end(d->hash); ++k)
- if (kh_exist(d->hash, k)) free(kh_val(d->hash, k));
- kh_destroy(64, d->hash);
- }
- free(fn_pos); free(fn_list); free(fn_fa);
- if (d->fai) fai_destroy(d->fai);
- bam_maqcns_destroy(d->c);
- free(d->ido); free(d->ref); free(d);
- return 0;
-}
-
-/***********
- * mpileup *
- ***********/
-
#include <assert.h>
#include "bam2bcf.h"
#include "sample.h"
#define MPLP_NO_INDEL 0x400
#define MPLP_EXT_BAQ 0x800
#define MPLP_ILLUMINA13 0x1000
+#define MPLP_IGNORE_RG 0x2000
+#define MPLP_PRINT_POS 0x4000
+#define MPLP_PRINT_MAPQ 0x8000
void *bed_read(const char *fn);
void bed_destroy(void *_h);
}
static void group_smpl(mplp_pileup_t *m, bam_sample_t *sm, kstring_t *buf,
- int n, char *const*fn, int *n_plp, const bam_pileup1_t **plp)
+ int n, char *const*fn, int *n_plp, const bam_pileup1_t **plp, int ignore_rg)
{
int i, j;
memset(m->n_plp, 0, m->n * sizeof(int));
const bam_pileup1_t *p = plp[i] + j;
uint8_t *q;
int id = -1;
- q = bam_aux_get(p->b, "RG");
+ q = ignore_rg? 0 : bam_aux_get(p->b, "RG");
if (q) id = bam_smpl_rg2smid(sm, fn[i], (char*)q+1, buf);
if (id < 0) id = bam_smpl_rg2smid(sm, fn[i], 0, buf);
if (id < 0 || id >= m->n) {
extern void *bcf_call_add_rg(void *rghash, const char *hdtext, const char *list);
extern void bcf_call_del_rghash(void *rghash);
mplp_aux_t **data;
- int i, tid, pos, *n_plp, beg0 = 0, end0 = 1u<<29, ref_len, ref_tid, max_depth, max_indel_depth;
+ int i, tid, pos, *n_plp, tid0 = -1, beg0 = 0, end0 = 1u<<29, ref_len, ref_tid = -1, max_depth, max_indel_depth;
const bam_pileup1_t **plp;
bam_mplp_t iter;
bam_header_t *h = 0;
data[i]->conf = conf;
h_tmp = bam_header_read(data[i]->fp);
data[i]->h = i? h : h_tmp; // for i==0, "h" has not been set yet
- bam_smpl_add(sm, fn[i], h_tmp->text);
+ bam_smpl_add(sm, fn[i], (conf->flag&MPLP_IGNORE_RG)? 0 : h_tmp->text);
rghash = bcf_call_add_rg(rghash, h_tmp->text, conf->pl_list);
if (conf->reg) {
int beg, end;
fprintf(stderr, "[%s] malformatted region or wrong seqname for %d-th input.\n", __func__, i+1);
exit(1);
}
- if (i == 0) beg0 = beg, end0 = end;
+ if (i == 0) tid0 = tid, beg0 = beg, end0 = end;
data[i]->iter = bam_iter_query(idx, tid, beg, end);
bam_index_destroy(idx);
}
bca->min_frac = conf->min_frac;
bca->min_support = conf->min_support;
}
- ref_tid = -1; ref = 0;
+ if (tid0 >= 0 && conf->fai) { // region is set
+ ref = faidx_fetch_seq(conf->fai, h->target_name[tid0], 0, 0x7fffffff, &ref_len);
+ ref_tid = tid0;
+ for (i = 0; i < n; ++i) data[i]->ref = ref, data[i]->ref_id = tid0;
+ } else ref_tid = -1, ref = 0;
iter = bam_mplp_init(n, mplp_func, (void**)data);
max_depth = conf->max_depth;
if (max_depth * sm->n > 1<<20)
if (conf->bed && tid >= 0 && !bed_overlap(conf->bed, h->target_name[tid], pos, pos+1)) continue;
if (tid != ref_tid) {
free(ref); ref = 0;
- if (conf->fai) ref = fai_fetch(conf->fai, h->target_name[tid], &ref_len);
+ if (conf->fai) ref = faidx_fetch_seq(conf->fai, h->target_name[tid], 0, 0x7fffffff, &ref_len);
for (i = 0; i < n; ++i) data[i]->ref = ref, data[i]->ref_id = tid;
ref_tid = tid;
}
int total_depth, _ref0, ref16;
bcf1_t *b = calloc(1, sizeof(bcf1_t));
for (i = total_depth = 0; i < n; ++i) total_depth += n_plp[i];
- group_smpl(&gplp, sm, &buf, n, fn, n_plp, plp);
+ group_smpl(&gplp, sm, &buf, n, fn, n_plp, plp, conf->flag & MPLP_IGNORE_RG);
_ref0 = (ref && pos < ref_len)? ref[pos] : 'N';
ref16 = bam_nt16_table[_ref0];
for (i = 0; i < gplp.n; ++i)
for (i = 0; i < n; ++i) {
int j;
printf("\t%d\t", n_plp[i]);
- if (n_plp[i] == 0) printf("*\t*");
- else {
+ if (n_plp[i] == 0) {
+ printf("*\t*"); // FIXME: printf() is very slow...
+ if (conf->flag & MPLP_PRINT_POS) printf("\t*");
+ } else {
for (j = 0; j < n_plp[i]; ++j)
pileup_seq(plp[i] + j, pos, ref_len, ref);
putchar('\t');
if (c > 126) c = 126;
putchar(c);
}
+ if (conf->flag & MPLP_PRINT_MAPQ) {
+ putchar('\t');
+ for (j = 0; j < n_plp[i]; ++j) {
+ int c = plp[i][j].b->core.qual + 33;
+ if (c > 126) c = 126;
+ putchar(c);
+ }
+ }
+ if (conf->flag & MPLP_PRINT_POS) {
+ putchar('\t');
+ for (j = 0; j < n_plp[i]; ++j) {
+ if (j > 0) putchar(',');
+ printf("%d", plp[i][j].qpos + 1); // FIXME: printf() is very slow...
+ }
+ }
}
}
putchar('\n');
int nfiles = 0, use_orphan = 0;
mplp_conf_t mplp;
memset(&mplp, 0, sizeof(mplp_conf_t));
+ #define MPLP_PRINT_POS 0x4000
mplp.max_mq = 60;
mplp.min_baseQ = 13;
mplp.capQ_thres = 0;
mplp.openQ = 40; mplp.extQ = 20; mplp.tandemQ = 100;
mplp.min_frac = 0.002; mplp.min_support = 1;
mplp.flag = MPLP_NO_ORPHAN | MPLP_REALN;
- while ((c = getopt(argc, argv, "Agf:r:l:M:q:Q:uaRC:BDSd:L:b:P:o:e:h:Im:F:EG:6")) >= 0) {
+ while ((c = getopt(argc, argv, "Agf:r:l:M:q:Q:uaRC:BDSd:L:b:P:o:e:h:Im:F:EG:6Os")) >= 0) {
switch (c) {
case 'f':
mplp.fai = fai_load(optarg);
case 'u': mplp.flag |= MPLP_NO_COMP | MPLP_GLF; break;
case 'a': mplp.flag |= MPLP_NO_ORPHAN | MPLP_REALN; break;
case 'B': mplp.flag &= ~MPLP_REALN; break;
- case 'R': mplp.flag |= MPLP_REALN; break;
case 'D': mplp.flag |= MPLP_FMT_DP; break;
case 'S': mplp.flag |= MPLP_FMT_SP; break;
case 'I': mplp.flag |= MPLP_NO_INDEL; break;
case 'E': mplp.flag |= MPLP_EXT_BAQ; break;
case '6': mplp.flag |= MPLP_ILLUMINA13; break;
+ case 'R': mplp.flag |= MPLP_IGNORE_RG; break;
+ case 's': mplp.flag |= MPLP_PRINT_MAPQ; break;
+ case 'O': mplp.flag |= MPLP_PRINT_POS; break;
case 'C': mplp.capQ_thres = atoi(optarg); break;
case 'M': mplp.max_mq = atoi(optarg); break;
case 'q': mplp.min_mq = atoi(optarg); break;
fprintf(stderr, " -A count anomalous read pairs\n");
fprintf(stderr, " -B disable BAQ computation\n");
fprintf(stderr, " -b FILE list of input BAM files [null]\n");
+ fprintf(stderr, " -C INT parameter for adjusting mapQ; 0 to disable [0]\n");
fprintf(stderr, " -d INT max per-BAM depth to avoid excessive memory usage [%d]\n", mplp.max_depth);
fprintf(stderr, " -E extended BAQ for higher sensitivity but lower specificity\n");
fprintf(stderr, " -f FILE faidx indexed reference sequence file [null]\n");
fprintf(stderr, " -l FILE list of positions (chr pos) or regions (BED) [null]\n");
fprintf(stderr, " -M INT cap mapping quality at INT [%d]\n", mplp.max_mq);
fprintf(stderr, " -r STR region in which pileup is generated [null]\n");
+ fprintf(stderr, " -R ignore RG tags\n");
fprintf(stderr, " -q INT skip alignments with mapQ smaller than INT [%d]\n", mplp.min_mq);
fprintf(stderr, " -Q INT skip bases with baseQ/BAQ smaller than INT [%d]\n", mplp.min_baseQ);
fprintf(stderr, "\nOutput options:\n\n");
fprintf(stderr, " -D output per-sample DP in BCF (require -g/-u)\n");
fprintf(stderr, " -g generate BCF output (genotype likelihoods)\n");
+ fprintf(stderr, " -O output base positions on reads (disabled by -g/-u)\n");
+ fprintf(stderr, " -s output mapping quality (disabled by -g/-u)\n");
fprintf(stderr, " -S output per-sample strand bias P-value in BCF (require -g/-u)\n");
fprintf(stderr, " -u generate uncompress BCF output\n");
fprintf(stderr, "\nSNP/INDEL genotype likelihoods options (effective with `-g' or `-u'):\n\n");
#include <ctype.h>
#include <assert.h>
#include <string.h>
+#include <math.h>
#include "bam.h"
#include "faidx.h"
-#include "bam_maqcns.h"
+#include "bam2bcf.h"
char bam_aux_getCEi(bam1_t *b, int i);
char bam_aux_getCSi(bam1_t *b, int i);
bamFile fp;
int curr_tid, left_pos;
faidx_t *fai;
- bam_maqcns_t *bmc;
+ bcf_callaux_t *bca;
int ccol, last_pos, row_shift, base_for, color_for, is_dot, l_ref, ins, no_skip, show_name;
char *ref;
int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data)
{
+ extern unsigned char bam_nt16_table[256];
tview_t *tv = (tview_t*)data;
int i, j, c, rb, attr, max_ins = 0;
uint32_t call = 0;
mvaddch(1, tv->ccol++, c);
}
if (pos%10 == 0 && tv->mcol - tv->ccol >= 10) mvprintw(0, tv->ccol, "%-d", pos+1);
- // print consensus
- call = bam_maqcns_call(n, pl, tv->bmc);
+ { // call consensus
+ bcf_callret1_t bcr;
+ int qsum[4], a1, a2, tmp;
+ double p[3], prior = 30;
+ bcf_call_glfgen(n, pl, bam_nt16_table[rb], tv->bca, &bcr);
+ for (i = 0; i < 4; ++i) qsum[i] = bcr.qsum[i]<<2 | i;
+ for (i = 1; i < 4; ++i) // insertion sort
+ for (j = i; j > 0 && qsum[j] > qsum[j-1]; --j)
+ tmp = qsum[j], qsum[j] = qsum[j-1], qsum[j-1] = tmp;
+ a1 = qsum[0]&3; a2 = qsum[1]&3;
+ p[0] = bcr.p[a1*5+a1]; p[1] = bcr.p[a1*5+a2] + prior; p[2] = bcr.p[a2*5+a2];
+ if ("ACGT"[a1] != toupper(rb)) p[0] += prior + 3;
+ if ("ACGT"[a2] != toupper(rb)) p[2] += prior + 3;
+ if (p[0] < p[1] && p[0] < p[2]) call = (1<<a1)<<16 | (int)((p[1]<p[2]?p[1]:p[2]) - p[0] + .499);
+ else if (p[2] < p[1] && p[2] < p[0]) call = (1<<a2)<<16 | (int)((p[0]<p[1]?p[0]:p[1]) - p[2] + .499);
+ else call = (1<<a1|1<<a2)<<16 | (int)((p[0]<p[2]?p[0]:p[2]) - p[1] + .499);
+ }
attr = A_UNDERLINE;
- c = ",ACMGRSVTWYHKDBN"[call>>28&0xf];
- i = (call>>8&0xff)/10+1;
+ c = ",ACMGRSVTWYHKDBN"[call>>16&0xf];
+ i = (call&0xffff)/10+1;
if (i > 4) i = 4;
attr |= COLOR_PAIR(i);
if (c == toupper(rb)) c = '.';
if (tv->idx == 0) exit(1);
tv->lplbuf = bam_lplbuf_init(tv_pl_func, tv);
if (fn_fa) tv->fai = fai_load(fn_fa);
- tv->bmc = bam_maqcns_init();
+ tv->bca = bcf_call_init(0.83, 13);
tv->ins = 1;
- bam_maqcns_prepare(tv->bmc);
initscr();
keypad(stdscr, TRUE);
endwin();
bam_lplbuf_destroy(tv->lplbuf);
- bam_maqcns_destroy(tv->bmc);
+ bcf_call_destroy(tv->bca);
bam_index_destroy(tv->idx);
if (tv->fai) fai_destroy(tv->fai);
free(tv->ref);
#endif
int bam_taf2baf(int argc, char *argv[]);
-int bam_pileup(int argc, char *argv[]);
int bam_mpileup(int argc, char *argv[]);
int bam_merge(int argc, char *argv[]);
int bam_index(int argc, char *argv[]);
int main_phase(int argc, char *argv[]);
int main_cat(int argc, char *argv[]);
int main_depth(int argc, char *argv[]);
+int main_bam2fq(int argc, char *argv[]);
int faidx_main(int argc, char *argv[]);
-int glf3_view_main(int argc, char *argv[]);
static int usage()
{
fprintf(stderr, "Usage: samtools <command> [options]\n\n");
fprintf(stderr, "Command: view SAM<->BAM conversion\n");
fprintf(stderr, " sort sort alignment file\n");
- fprintf(stderr, " pileup generate pileup output\n");
fprintf(stderr, " mpileup multi-way pileup\n");
fprintf(stderr, " depth compute the depth\n");
fprintf(stderr, " faidx index/extract FASTA\n");
fprintf(stderr, " index index alignment\n");
fprintf(stderr, " idxstats BAM index stats (r595 or later)\n");
fprintf(stderr, " fixmate fix mate information\n");
- fprintf(stderr, " glfview print GLFv3 file\n");
fprintf(stderr, " flagstat simple stats\n");
fprintf(stderr, " calmd recalculate MD/NM tags and '=' bases\n");
fprintf(stderr, " merge merge sorted alignments\n");
if (argc < 2) return usage();
if (strcmp(argv[1], "view") == 0) return main_samview(argc-1, argv+1);
else if (strcmp(argv[1], "import") == 0) return main_import(argc-1, argv+1);
- else if (strcmp(argv[1], "pileup") == 0) return bam_pileup(argc-1, argv+1);
else if (strcmp(argv[1], "mpileup") == 0) return bam_mpileup(argc-1, argv+1);
else if (strcmp(argv[1], "merge") == 0) return bam_merge(argc-1, argv+1);
else if (strcmp(argv[1], "sort") == 0) return bam_sort(argc-1, argv+1);
else if (strcmp(argv[1], "faidx") == 0) return faidx_main(argc-1, argv+1);
else if (strcmp(argv[1], "fixmate") == 0) return bam_mating(argc-1, argv+1);
else if (strcmp(argv[1], "rmdup") == 0) return bam_rmdup(argc-1, argv+1);
- else if (strcmp(argv[1], "glfview") == 0) return glf3_view_main(argc-1, argv+1);
else if (strcmp(argv[1], "flagstat") == 0) return bam_flagstat(argc-1, argv+1);
else if (strcmp(argv[1], "calmd") == 0) return bam_fillmd(argc-1, argv+1);
else if (strcmp(argv[1], "fillmd") == 0) return bam_fillmd(argc-1, argv+1);
else if (strcmp(argv[1], "targetcut") == 0) return main_cut_target(argc-1, argv+1);
else if (strcmp(argv[1], "phase") == 0) return main_phase(argc-1, argv+1);
else if (strcmp(argv[1], "depth") == 0) return main_depth(argc-1, argv+1);
+ else if (strcmp(argv[1], "bam2fq") == 0) return main_bam2fq(argc-1, argv+1);
+ else if (strcmp(argv[1], "pileup") == 0) {
+ fprintf(stderr, "[main] The `pileup' command has been removed. Please use `mpileup' instead.\n");
+ return 1;
+ }
#if _CURSES_LIB != 0
else if (strcmp(argv[1], "tview") == 0) return bam_tview_main(argc-1, argv+1);
#endif
CC= gcc
CFLAGS= -g -Wall -O2 #-m64 #-arch ppc
DFLAGS= -D_FILE_OFFSET_BITS=64 -D_USE_KNETFILE
-LOBJS= bcf.o vcf.o bcfutils.o prob1.o em.o kfunc.o kmin.o index.o fet.o bcf2qcall.o
+LOBJS= bcf.o vcf.o bcfutils.o prob1.o em.o kfunc.o kmin.o index.o fet.o mut.o bcf2qcall.o
OMISC= ..
AOBJS= call1.o main.o $(OMISC)/kstring.o $(OMISC)/bgzf.o $(OMISC)/knetfile.o $(OMISC)/bedidx.o
PROG= bcftools
lib:libbcf.a
libbcf.a:$(LOBJS)
- $(AR) -cru $@ $(LOBJS)
+ $(AR) -csru $@ $(LOBJS)
bcftools:lib $(AOBJS)
$(CC) $(CFLAGS) -o $@ $(AOBJS) -L. $(LIBPATH) -lbcf -lm -lz
#ifndef BCF_H
#define BCF_H
+#define BCF_VERSION "0.1.17 (r973:277)"
+
#include <stdint.h>
#include <zlib.h>
int bcf_fix_gt(bcf1_t *b);
// update PL generated by old samtools
int bcf_fix_pl(bcf1_t *b);
+ // convert PL to GLF-like 10-likelihood GL
+ int bcf_gl10(const bcf1_t *b, uint8_t *gl);
// string hash table
void *bcf_build_refhash(bcf_hdr_t *h);
+++ /dev/null
-.TH bcftools 1 "16 March 2011" "bcftools" "Bioinformatics tools"
-.SH NAME
-.PP
-bcftools - Utilities for the Binary Call Format (BCF) and VCF.
-.SH SYNOPSIS
-.PP
-bcftools index in.bcf
-.PP
-bcftools view in.bcf chr2:100-200 > out.vcf
-.PP
-bcftools view -vc in.bcf > out.vcf 2> out.afs
-
-.SH DESCRIPTION
-.PP
-Bcftools is a toolkit for processing VCF/BCF files, calling variants and
-estimating site allele frequencies and allele frequency spectrums.
-
-.SH COMMANDS AND OPTIONS
-
-.TP 10
-.B view
-.B bcftools view
-.RB [ \-AbFGNQSucgv ]
-.RB [ \-D
-.IR seqDict ]
-.RB [ \-l
-.IR listLoci ]
-.RB [ \-s
-.IR listSample ]
-.RB [ \-i
-.IR gapSNPratio ]
-.RB [ \-t
-.IR mutRate ]
-.RB [ \-p
-.IR varThres ]
-.RB [ \-P
-.IR prior ]
-.RB [ \-1
-.IR nGroup1 ]
-.RB [ \-d
-.IR minFrac ]
-.RB [ \-U
-.IR nPerm ]
-.RB [ \-X
-.IR permThres ]
-.I in.bcf
-.RI [ region ]
-
-Convert between BCF and VCF, call variant candidates and estimate allele
-frequencies.
-
-.RS
-.TP
-.B Input/Output Options:
-.TP 10
-.B -A
-Retain all possible alternate alleles at variant sites. By default, the view
-command discards unlikely alleles.
-.TP 10
-.B -b
-Output in the BCF format. The default is VCF.
-.TP
-.BI -D \ FILE
-Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null]
-.TP
-.B -F
-Indicate PL is generated by r921 or before (ordering is different).
-.TP
-.B -G
-Suppress all individual genotype information.
-.TP
-.BI -l \ FILE
-List of sites at which information are outputted [all sites]
-.TP
-.B -N
-Skip sites where the REF field is not A/C/G/T
-.TP
-.B -Q
-Output the QCALL likelihood format
-.TP
-.BI -s \ FILE
-List of samples to use. The first column in the input gives the sample names
-and the second gives the ploidy, which can only be 1 or 2. When the 2nd column
-is absent, the sample ploidy is assumed to be 2. In the output, the ordering of
-samples will be identical to the one in
-.IR FILE .
-[null]
-.TP
-.B -S
-The input is VCF instead of BCF.
-.TP
-.B -u
-Uncompressed BCF output (force -b).
-.TP
-.B Consensus/Variant Calling Options:
-.TP 10
-.B -c
-Call variants using Bayesian inference. This option automatically invokes option
-.BR -e .
-.TP
-.BI -d \ FLOAT
-When
-.B -v
-is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0]
-.TP
-.B -e
-Perform max-likelihood inference only, including estimating the site allele frequency,
-testing Hardy-Weinberg equlibrium and testing associations with LRT.
-.TP
-.B -g
-Call per-sample genotypes at variant sites (force -c)
-.TP
-.BI -i \ FLOAT
-Ratio of INDEL-to-SNP mutation rate [0.15]
-.TP
-.BI -p \ FLOAT
-A site is considered to be a variant if P(ref|D)<FLOAT [0.5]
-.TP
-.BI -P \ STR
-Prior or initial allele frequency spectrum. If STR can be
-.IR full ,
-.IR cond2 ,
-.I flat
-or the file consisting of error output from a previous variant calling
-run.
-.TP
-.BI -t \ FLOAT
-Scaled muttion rate for variant calling [0.001]
-.TP
-.B -v
-Output variant sites only (force -c)
-.TP
-.B Contrast Calling and Association Test Options:
-.TP
-.BI -1 \ INT
-Number of group-1 samples. This option is used for dividing the samples into
-two groups for contrast SNP calling or association test.
-When this option is in use, the following VCF INFO will be outputted:
-PC2, PCHI2 and QCHI2. [0]
-.TP
-.BI -U \ INT
-Number of permutations for association test (effective only with
-.BR -1 )
-[0]
-.TP
-.BI -X \ FLOAT
-Only perform permutations for P(chi^2)<FLOAT (effective only with
-.BR -U )
-[0.01]
-.RE
-
-.TP
-.B index
-.B bcftools index
-.I in.bcf
-
-Index sorted BCF for random access.
-.RE
-
-.TP
-.B cat
-.B bcftools cat
-.I in1.bcf
-.RI [ "in2.bcf " [ ... "]]]"
-
-Concatenate BCF files. The input files are required to be sorted and
-have identical samples appearing in the same order.
-.RE
-
-.SH BCFTOOLS SPECIFIC VCF TAGS
-
-.TS
-center box;
-cb | cb | cb
-l | l | l .
-Tag Format Description
-_
-AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
-CI95 double[2] Equal-tail Bayesian credible interval of AF at the 95% level
-DP int Raw read depth (without quality filtering)
-DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
-FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise
-MQ int Root-Mean-Square mapping quality of covering reads
-PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
-PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
-PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
-QCHI2 int Phred-scaled PCHI2
-RP int # permutations yielding a smaller PCHI2
-.TE
#include "khash.h"
KHASH_MAP_INIT_STR(str2id, int)
+// FIXME: valgrind report a memory leak in this function. Probably it does not get deallocated...
void *bcf_build_refhash(bcf_hdr_t *h)
{
khash_t(str2id) *hash;
b->n_smpl = n_smpl;
return 0;
}
+
+static int8_t nt4_table[128] = {
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 /*'-'*/, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 3, 4, 4, 4, -1, 4, 4, 4, 4, 4, 4, 4,
+ 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 3, 4, 4, 4, -1, 4, 4, 4, 4, 4, 4, 4
+};
+
+int bcf_gl10(const bcf1_t *b, uint8_t *gl)
+{
+ int a[4], k, l, map[4], k1, j, i;
+ const bcf_ginfo_t *PL;
+ char *s;
+ if (b->ref[1] != 0 || b->n_alleles > 4) return -1; // ref is not a single base or >4 alleles
+ for (i = 0; i < b->n_gi; ++i)
+ if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
+ if (i == b->n_gi) return -1; // no PL
+ PL = b->gi + i;
+ a[0] = nt4_table[(int)b->ref[0]];
+ if (a[0] > 3 || a[0] < 0) return -1; // ref is not A/C/G/T
+ a[1] = a[2] = a[3] = -2; // -1 has a special meaning
+ if (b->alt[0] == 0) return -1; // no alternate allele
+ map[0] = map[1] = map[2] = map[3] = -2;
+ map[a[0]] = 0;
+ for (k = 0, s = b->alt, k1 = -1; k < 3 && *s; ++k, s += 2) {
+ if (s[1] != ',' && s[1] != 0) return -1; // ALT is not single base
+ a[k+1] = nt4_table[(int)*s];
+ if (a[k+1] >= 0) map[a[k+1]] = k+1;
+ else k1 = k + 1;
+ if (s[1] == 0) break; // the end of the ALT string
+ }
+ for (k = 0; k < 4; ++k)
+ if (map[k] < 0) map[k] = k1;
+ for (i = 0; i < b->n_smpl; ++i) {
+ const uint8_t *p = PL->data + i * PL->len; // the PL for the i-th individual
+ uint8_t *g = gl + 10 * i;
+ for (j = 0; j < PL->len; ++j)
+ if (p[j]) break;
+ for (k = j = 0; k < 4; ++k) {
+ for (l = k; l < 4; ++l) {
+ int t, x = map[k], y = map[l];
+ if (x > y) t = x, x = y, y = t; // make sure x is the smaller
+ g[j++] = p[y * (y+1) / 2 + x];
+ }
+ }
+ }
+ return 0;
+}
#define VC_ANNO_MAX 16384
#define VC_FIX_PL 32768
#define VC_EM 0x10000
+#define VC_PAIRCALL 0x20000
+#define VC_QCNT 0x40000
typedef struct {
int flag, prior_type, n1, n_sub, *sublist, n_perm;
+ uint32_t *trio_aux;
char *prior_file, **subsam, *fn_dict;
uint8_t *ploidy;
double theta, pref, indel_frac, min_perm_p, min_smpl_frac, min_lrt;
bcf_sync(b);
}
-static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag, double em[9])
+static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag, double em[10], int cons_llr, int64_t cons_gt)
{
kstring_t s;
int has_I16, is_var;
if (em[4] >= 0 && em[4] <= 0.05) ksprintf(&s, ";G3=%.4g,%.4g,%.4g;HWE=%.3g", em[3], em[2], em[1], em[4]);
if (em[5] >= 0 && em[6] >= 0) ksprintf(&s, ";AF2=%.4g,%.4g", 1 - em[5], 1 - em[6]);
if (em[7] >= 0) ksprintf(&s, ";LRT=%.3g", em[7]);
+ if (em[8] >= 0) ksprintf(&s, ";LRT2=%.3g", em[8]);
+ //if (em[9] >= 0) ksprintf(&s, ";LRT1=%.3g", em[9]);
+ }
+ if (cons_llr > 0) {
+ ksprintf(&s, ";CLR=%d", cons_llr);
+ if (cons_gt > 0)
+ ksprintf(&s, ";UGT=%c%c%c;CGT=%c%c%c", cons_gt&0xff, cons_gt>>8&0xff, cons_gt>>16&0xff,
+ cons_gt>>32&0xff, cons_gt>>40&0xff, cons_gt>>48&0xff);
}
if (pr == 0) { // if pr is unset, return
kputc('\0', &s); kputs(b->fmt, &s); kputc('\0', &s);
is_var = (pr->p_ref < pref);
r = is_var? pr->p_ref : pr->p_var;
- ksprintf(&s, ";CI95=%.4g,%.4g", pr->cil, pr->cih); // FIXME: when EM is not used, ";" should be omitted!
+// ksprintf(&s, ";CI95=%.4g,%.4g", pr->cil, pr->cih); // FIXME: when EM is not used, ";" should be omitted!
+ ksprintf(&s, ";AC1=%d", pr->ac);
if (has_I16) ksprintf(&s, ";DP4=%d,%d,%d,%d;MQ=%d", a.d[0], a.d[1], a.d[2], a.d[3], a.mq);
fq = pr->p_ref_folded < 0.5? -4.343 * log(pr->p_ref_folded) : 4.343 * log(pr->p_var_folded);
if (fq < -999) fq = -999;
if (q[i] > 255) q[i] = 255;
}
if (pr->perm_rank >= 0) ksprintf(&s, ";PR=%d", pr->perm_rank);
+ // ksprintf(&s, ";LRT3=%.3g", pr->lrt);
ksprintf(&s, ";PCHI2=%.3g;PC2=%d,%d", q[1], q[2], pr->p_chi2);
}
if (has_I16 && a.is_tested) ksprintf(&s, ";PV4=%.2g,%.2g,%.2g,%.2g", a.p[0], a.p[1], a.p[2], a.p[3]);
if (!strstr(str.s, "##INFO=<ID=FQ,"))
kputs("##INFO=<ID=FQ,Number=1,Type=Float,Description=\"Phred probability of all samples being the same\">\n", &str);
if (!strstr(str.s, "##INFO=<ID=AF1,"))
- kputs("##INFO=<ID=AF1,Number=1,Type=Float,Description=\"Max-likelihood estimate of the site allele frequency of the first ALT allele\">\n", &str);
+ kputs("##INFO=<ID=AF1,Number=1,Type=Float,Description=\"Max-likelihood estimate of the first ALT allele frequency (assuming HWE)\">\n", &str);
+ if (!strstr(str.s, "##INFO=<ID=AC1,"))
+ kputs("##INFO=<ID=AC1,Number=1,Type=Float,Description=\"Max-likelihood estimate of the first ALT allele count (no HWE assumption)\">\n", &str);
if (!strstr(str.s, "##INFO=<ID=G3,"))
kputs("##INFO=<ID=G3,Number=3,Type=Float,Description=\"ML estimate of genotype frequencies\">\n", &str);
if (!strstr(str.s, "##INFO=<ID=HWE,"))
kputs("##INFO=<ID=HWE,Number=1,Type=Float,Description=\"Chi^2 based HWE test P-value based on G3\">\n", &str);
- if (!strstr(str.s, "##INFO=<ID=CI95,"))
- kputs("##INFO=<ID=CI95,Number=2,Type=Float,Description=\"Equal-tail Bayesian credible interval of the site allele frequency at the 95% level\">\n", &str);
+ if (!strstr(str.s, "##INFO=<ID=CLR,"))
+ kputs("##INFO=<ID=CLR,Number=1,Type=Integer,Description=\"Log ratio of genotype likelihoods with and without the constraint\">\n", &str);
+ if (!strstr(str.s, "##INFO=<ID=UGT,"))
+ kputs("##INFO=<ID=UGT,Number=1,Type=String,Description=\"The most probable unconstrained genotype configuration in the trio\">\n", &str);
+ if (!strstr(str.s, "##INFO=<ID=CGT,"))
+ kputs("##INFO=<ID=CGT,Number=1,Type=String,Description=\"The most probable constrained genotype configuration in the trio\">\n", &str);
+// if (!strstr(str.s, "##INFO=<ID=CI95,"))
+// kputs("##INFO=<ID=CI95,Number=2,Type=Float,Description=\"Equal-tail Bayesian credible interval of the site allele frequency at the 95% level\">\n", &str);
if (!strstr(str.s, "##INFO=<ID=PV4,"))
kputs("##INFO=<ID=PV4,Number=4,Type=Float,Description=\"P-values for strand bias, baseQ bias, mapQ bias and tail distance bias\">\n", &str);
if (!strstr(str.s, "##INFO=<ID=INDEL,"))
extern int bcf_fix_gt(bcf1_t *b);
extern int bcf_anno_max(bcf1_t *b);
extern int bcf_shuffle(bcf1_t *b, int seed);
+ extern uint32_t *bcf_trio_prep(int is_x, int is_son);
+ extern int bcf_trio_call(uint32_t *prep, const bcf1_t *b, int *llr, int64_t *gt);
+ extern int bcf_pair_call(const bcf1_t *b);
+ extern int bcf_min_diff(const bcf1_t *b);
+
bcf_t *bp, *bout = 0;
bcf1_t *b, *blast;
int c, *seeds = 0;
- uint64_t n_processed = 0;
+ uint64_t n_processed = 0, qcnt[256];
viewconf_t vc;
bcf_p1aux_t *p1 = 0;
bcf_hdr_t *hin, *hout;
tid = begin = end = -1;
memset(&vc, 0, sizeof(viewconf_t));
vc.prior_type = vc.n1 = -1; vc.theta = 1e-3; vc.pref = 0.5; vc.indel_frac = -1.; vc.n_perm = 0; vc.min_perm_p = 0.01; vc.min_smpl_frac = 0; vc.min_lrt = 1;
- while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:")) >= 0) {
+ memset(qcnt, 0, 8 * 256);
+ while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:T:Y")) >= 0) {
switch (c) {
case '1': vc.n1 = atoi(optarg); break;
case 'l': vc.bed = bed_read(optarg); break;
case 'g': vc.flag |= VC_CALL_GT | VC_CALL; break;
case 'I': vc.flag |= VC_NO_INDEL; break;
case 'M': vc.flag |= VC_ANNO_MAX; break;
+ case 'Y': vc.flag |= VC_QCNT; break;
case 't': vc.theta = atof(optarg); break;
case 'p': vc.pref = atof(optarg); break;
case 'i': vc.indel_frac = atof(optarg); break;
for (tid = 0; tid < vc.n_sub; ++tid) vc.ploidy[tid] = vc.subsam[tid][strlen(vc.subsam[tid]) + 1];
tid = -1;
break;
+ case 'T':
+ if (strcmp(optarg, "trioauto") == 0) vc.trio_aux = bcf_trio_prep(0, 0);
+ else if (strcmp(optarg, "trioxd") == 0) vc.trio_aux = bcf_trio_prep(1, 0);
+ else if (strcmp(optarg, "trioxs") == 0) vc.trio_aux = bcf_trio_prep(1, 1);
+ else if (strcmp(optarg, "pair") == 0) vc.flag |= VC_PAIRCALL;
+ else {
+ fprintf(stderr, "[%s] Option '-T' can only take value trioauto, trioxd or trioxs.\n", __func__);
+ return 1;
+ }
+ break;
case 'P':
if (strcmp(optarg, "full") == 0) vc.prior_type = MC_PTYPE_FULL;
else if (strcmp(optarg, "cond2") == 0) vc.prior_type = MC_PTYPE_COND2;
fprintf(stderr, " -p FLOAT variant if P(ref|D)<FLOAT [%.3g]\n", vc.pref);
fprintf(stderr, " -P STR type of prior: full, cond2, flat [full]\n");
fprintf(stderr, " -t FLOAT scaled substitution mutation rate [%.4g]\n", vc.theta);
+ fprintf(stderr, " -T STR constrained calling; STR can be: pair, trioauto, trioxd and trioxs (see manual) [null]\n");
fprintf(stderr, " -v output potential variant sites only (force -c)\n");
fprintf(stderr, "\nContrast calling and association test options:\n\n");
fprintf(stderr, " -1 INT number of group-1 samples [0]\n");
}
}
while (vcf_read(bp, hin, b) > 0) {
- int is_indel;
- double em[9];
+ int is_indel, cons_llr = -1;
+ int64_t cons_gt = -1;
+ double em[10];
if ((vc.flag & VC_VARONLY) && strcmp(b->alt, "X") == 0) continue;
if ((vc.flag & VC_VARONLY) && vc.min_smpl_frac > 0.) {
extern int bcf_smpl_covered(const bcf1_t *b);
if (!(l > begin && end > b->pos)) continue;
}
++n_processed;
+ if (vc.flag & VC_QCNT) { // summarize the difference
+ int x = bcf_min_diff(b);
+ if (x > 255) x = 255;
+ if (x >= 0) ++qcnt[x];
+ }
if (vc.flag & VC_QCALL) { // output QCALL format; STOP here
bcf_2qcall(hout, b);
continue;
}
- if (vc.flag & VC_EM) bcf_em1(b, vc.n1, 0xff, em);
+ if (vc.trio_aux) // do trio calling
+ bcf_trio_call(vc.trio_aux, b, &cons_llr, &cons_gt);
+ else if (vc.flag & VC_PAIRCALL)
+ cons_llr = bcf_pair_call(b);
+ if (vc.flag & (VC_CALL|VC_ADJLD|VC_EM)) bcf_gl2pl(b);
+ if (vc.flag & VC_EM) bcf_em1(b, vc.n1, 0x1ff, em);
else {
int i;
for (i = 0; i < 9; ++i) em[i] = -1.;
}
- if (vc.flag & (VC_CALL|VC_ADJLD)) bcf_gl2pl(b);
if (vc.flag & VC_CALL) { // call variants
bcf_p1rst_t pr;
int calret = bcf_p1_cal(b, (em[7] >= 0 && em[7] < vc.min_lrt), p1, &pr);
int i, n = 0;
for (i = 0; i < vc.n_perm; ++i) {
#ifdef BCF_PERM_LRT // LRT based permutation is much faster but less robust to artifacts
- double x[9];
+ double x[10];
bcf_shuffle(b, seeds[i]);
bcf_em1(b, vc.n1, 1<<7, x);
if (x[7] < em[7]) ++n;
}
pr.perm_rank = n;
}
- if (calret >= 0) update_bcf1(b, p1, &pr, vc.pref, vc.flag, em);
- } else if (vc.flag & VC_EM) update_bcf1(b, 0, 0, 0, vc.flag, em);
+ if (calret >= 0) update_bcf1(b, p1, &pr, vc.pref, vc.flag, em, cons_llr, cons_gt);
+ } else if (vc.flag & VC_EM) update_bcf1(b, 0, 0, 0, vc.flag, em, cons_llr, cons_gt);
if (vc.flag & VC_ADJLD) { // compute LD
double f[4], r2;
if ((r2 = bcf_pair_freq(blast, b, f)) >= 0) {
vcf_close(bp); vcf_close(bout);
if (vc.fn_dict) free(vc.fn_dict);
if (vc.ploidy) free(vc.ploidy);
+ if (vc.trio_aux) free(vc.trio_aux);
if (vc.n_sub) {
int i;
for (i = 0; i < vc.n_sub; ++i) free(vc.subsam[i]);
free(vc.subsam); free(vc.sublist);
}
if (vc.bed) bed_destroy(vc.bed);
+ if (vc.flag & VC_QCNT)
+ for (c = 0; c < 256; ++c)
+ fprintf(stderr, "QT\t%d\t%lld\n", c, (long long)qcnt[c]);
if (seeds) free(seeds);
if (p1) bcf_p1_destroy(p1);
return 0;
// x[5..6]: group1 freq, group2 freq
// x[7]: 1-degree P-value
// x[8]: 2-degree P-value
-int bcf_em1(const bcf1_t *b, int n1, int flag, double x[9])
+// x[9]: 1-degree P-value with freq estimated from genotypes
+int bcf_em1(const bcf1_t *b, int n1, int flag, double x[10])
{
double *pdg;
int i, n, n2;
n = b->n_smpl; n2 = n - n1;
pdg = get_pdg3(b);
if (pdg == 0) return -1;
- for (i = 0; i < 9; ++i) x[i] = -1.;
+ for (i = 0; i < 10; ++i) x[i] = -1.; // set to negative
{
if ((x[0] = est_freq(n, pdg)) < 0.) {
free(pdg);
}
x[0] = freqml(x[0], 0, n, pdg);
}
- if (flag & (0xf<<1|1<<8)) { // estimate the genotype frequency and test HWE
+ if (flag & (0xf<<1|3<<8)) { // estimate the genotype frequency and test HWE
double *g = x + 1, f3[3], r;
f3[0] = g[0] = (1 - x[0]) * (1 - x[0]);
f3[1] = g[1] = 2 * x[0] * (1 - x[0]);
f3[i][0] = (1-f[i])*(1-f[i]), f3[i][1] = 2*f[i]*(1-f[i]), f3[i][2] = f[i]*f[i];
x[7] = kf_gammaq(.5, log(lk_ratio_test(n, n1, pdg, f3)));
}
- if ((flag & 1<<8) && n1 > 0 && n1 < n) { // 2-degree P-value
- double g[3][3];
+ if ((flag & 3<<8) && n1 > 0 && n1 < n) { // 2-degree P-value
+ double g[3][3], tmp;
for (i = 0; i < 3; ++i) memcpy(g[i], x + 1, 3 * sizeof(double));
for (i = 0; i < ITER_MAX; ++i)
if (g3_iter(g[1], pdg, 0, n1) < EPS) break;
for (i = 0; i < ITER_MAX; ++i)
if (g3_iter(g[2], pdg, n1, n) < EPS) break;
- x[8] = kf_gammaq(1., log(lk_ratio_test(n, n1, pdg, g)));
+ tmp = log(lk_ratio_test(n, n1, pdg, g));
+ x[8] = kf_gammaq(1., tmp);
+ x[9] = kf_gammaq(.5, tmp);
}
// free
free(pdg);
{
if (argc == 1) {
fprintf(stderr, "\n");
+ fprintf(stderr, "Program: bcftools (Tools for data in the VCF/BCF formats)\n");
+ fprintf(stderr, "Version: %s\n\n", BCF_VERSION);
fprintf(stderr, "Usage: bcftools <command> <arguments>\n\n");
fprintf(stderr, "Command: view print, extract, convert and call SNPs from BCF\n");
fprintf(stderr, " index index BCF\n");
--- /dev/null
+#include <stdlib.h>
+#include <stdint.h>
+#include "bcf.h"
+
+#define MAX_GENO 359
+
+int8_t seq_bitcnt[] = { 4, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4 };
+char *seq_nt16rev = "XACMGRSVTWYHKDBN";
+
+uint32_t *bcf_trio_prep(int is_x, int is_son)
+{
+ int i, j, k, n, map[10];
+ uint32_t *ret;
+ ret = calloc(MAX_GENO, 4);
+ for (i = 0, k = 0; i < 4; ++i)
+ for (j = i; j < 4; ++j)
+ map[k++] = 1<<i|1<<j;
+ for (i = 0, n = 1; i < 10; ++i) { // father
+ if (is_x && seq_bitcnt[map[i]] != 1) continue;
+ if (is_x && is_son) {
+ for (j = 0; j < 10; ++j) // mother
+ for (k = 0; k < 10; ++k) // child
+ if (seq_bitcnt[map[k]] == 1 && (map[j]&map[k]))
+ ret[n++] = j<<16 | i<<8 | k;
+ } else {
+ for (j = 0; j < 10; ++j) // mother
+ for (k = 0; k < 10; ++k) // child
+ if ((map[i]&map[k]) && (map[j]&map[k]) && ((map[i]|map[j])&map[k]) == map[k])
+ ret[n++] = j<<16 | i<<8 | k;
+ }
+ }
+ ret[0] = n - 1;
+ return ret;
+}
+
+int bcf_trio_call(const uint32_t *prep, const bcf1_t *b, int *llr, int64_t *gt)
+{
+ int i, j, k;
+ const bcf_ginfo_t *PL;
+ uint8_t *gl10;
+ int map[10];
+ if (b->n_smpl != 3) return -1; // not a trio
+ for (i = 0; i < b->n_gi; ++i)
+ if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
+ if (i == b->n_gi) return -1; // no PL
+ gl10 = alloca(10 * b->n_smpl);
+ if (bcf_gl10(b, gl10) < 0) return -1;
+ PL = b->gi + i;
+ for (i = 0, k = 0; i < 4; ++i)
+ for (j = i; j < 4; ++j)
+ map[k++] = seq_nt16rev[1<<i|1<<j];
+ for (j = 0; j < 3; ++j) // check if ref hom is the most probable in all members
+ if (((uint8_t*)PL->data)[j * PL->len] != 0) break;
+ if (j < 3) { // we need to go through the complex procedure
+ uint8_t *g[3];
+ int minc = 1<<30, minc_j = -1, minf = 0, gtf = 0, gtc = 0;
+ g[0] = gl10;
+ g[1] = gl10 + 10;
+ g[2] = gl10 + 20;
+ for (j = 1; j <= (int)prep[0]; ++j) { // compute LK with constraint
+ int sum = g[0][prep[j]&0xff] + g[1][prep[j]>>8&0xff] + g[2][prep[j]>>16&0xff];
+ if (sum < minc) minc = sum, minc_j = j;
+ }
+ gtc |= map[prep[minc_j]&0xff]; gtc |= map[prep[minc_j]>>8&0xff]<<8; gtc |= map[prep[minc_j]>>16]<<16;
+ for (j = 0; j < 3; ++j) { // compute LK without constraint
+ int min = 1<<30, min_k = -1;
+ for (k = 0; k < 10; ++k)
+ if (g[j][k] < min) min = g[j][k], min_k = k;
+ gtf |= map[min_k]<<(j*8);
+ minf += min;
+ }
+ *llr = minc - minf; *gt = (int64_t)gtc<<32 | gtf;
+ } else *llr = 0, *gt = -1;
+ return 0;
+}
+
+int bcf_pair_call(const bcf1_t *b)
+{
+ int i, j, k;
+ const bcf_ginfo_t *PL;
+ if (b->n_smpl != 2) return -1; // not a pair
+ for (i = 0; i < b->n_gi; ++i)
+ if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
+ if (i == b->n_gi) return -1; // no PL
+ PL = b->gi + i;
+ for (j = 0; j < 2; ++j) // check if ref hom is the most probable in all members
+ if (((uint8_t*)PL->data)[j * PL->len] != 0) break;
+ if (j < 2) { // we need to go through the complex procedure
+ uint8_t *g[2];
+ int minc = 1<<30, minf = 0;
+ g[0] = PL->data;
+ g[1] = (uint8_t*)PL->data + PL->len;
+ for (j = 0; j < PL->len; ++j) // compute LK with constraint
+ minc = minc < g[0][j] + g[1][j]? minc : g[0][j] + g[1][j];
+ for (j = 0; j < 2; ++j) { // compute LK without constraint
+ int min = 1<<30;
+ for (k = 0; k < PL->len; ++k)
+ min = min < g[j][k]? min : g[j][k];
+ minf += min;
+ }
+ return minc - minf;
+ } else return 0;
+}
+
+int bcf_min_diff(const bcf1_t *b)
+{
+ int i, min = 1<<30;
+ const bcf_ginfo_t *PL;
+ for (i = 0; i < b->n_gi; ++i)
+ if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
+ if (i == b->n_gi) return -1; // no PL
+ PL = b->gi + i;
+ for (i = 0; i < b->n_smpl; ++i) {
+ int m1, m2, j;
+ const uint8_t *p = (uint8_t*)PL->data;
+ m1 = m2 = 1<<30;
+ for (j = 0; j < PL->len; ++j) {
+ if ((int)p[j] < m1) m2 = m1, m1 = p[j];
+ else if ((int)p[j] < m2) m2 = p[j];
+ }
+ min = min < m2 - m1? min : m2 - m1;
+ }
+ return min;
+}
for (k = 0, sum = 0.; k < ma->M; ++k)
sum += ma->afs1[k];
rst->p_var = (double)sum;
+ { // compute the allele count
+ double max = -1;
+ rst->ac = -1;
+ for (k = 0; k <= ma->M; ++k)
+ if (max < ma->z[k]) max = ma->z[k], rst->ac = k;
+ rst->ac = ma->M - rst->ac;
+ }
// calculate f_flat and f_em
for (k = 0, sum = 0.; k <= ma->M; ++k)
sum += (long double)ma->z[k];
{ // estimate equal-tail credible interval (95% level)
int l, h;
double p;
- for (i = 0, p = 0.; i < ma->M; ++i)
+ for (i = 0, p = 0.; i <= ma->M; ++i)
if (p + ma->afs1[i] > 0.025) break;
else p += ma->afs1[i];
l = i;
- for (i = ma->M-1, p = 0.; i >= 0; --i)
+ for (i = ma->M, p = 0.; i >= 0; --i)
if (p + ma->afs1[i] > 0.025) break;
else p += ma->afs1[i];
h = i;
rst->cil = (double)(ma->M - h) / ma->M; rst->cih = (double)(ma->M - l) / ma->M;
}
+ if (ma->n1 > 0) { // compute LRT
+ double max0, max1, max2;
+ for (k = 0, max0 = -1; k <= ma->M; ++k)
+ if (max0 < ma->z[k]) max0 = ma->z[k];
+ for (k = 0, max1 = -1; k <= ma->n1 * 2; ++k)
+ if (max1 < ma->z1[k]) max1 = ma->z1[k];
+ for (k = 0, max2 = -1; k <= ma->M - ma->n1 * 2; ++k)
+ if (max2 < ma->z2[k]) max2 = ma->z2[k];
+ rst->lrt = log(max1 * max2 / max0);
+ rst->lrt = rst->lrt < 0? 1 : kf_gammaq(.5, rst->lrt);
+ } else rst->lrt = -1.0;
rst->cmp[0] = rst->cmp[1] = rst->cmp[2] = rst->p_chi2 = -1.0;
if (do_contrast && rst->p_var > 0.5) // skip contrast2() if the locus is a strong non-variant
rst->p_chi2 = contrast2(ma, rst->cmp);
typedef struct {
int rank0, perm_rank; // NB: perm_rank is always set to -1 by bcf_p1_cal()
+ int ac; // ML alternative allele count
double f_exp, f_flat, p_ref_folded, p_ref, p_var_folded, p_var;
double cil, cih;
- double cmp[3], p_chi2; // used by contrast2()
+ double cmp[3], p_chi2, lrt; // used by contrast2()
} bcf_p1rst_t;
#define MC_PTYPE_FULL 1
int bcf_p1_set_n1(bcf_p1aux_t *b, int n1);
void bcf_p1_set_folded(bcf_p1aux_t *p1a); // only effective when set_n1() is not called
- int bcf_em1(const bcf1_t *b, int n1, int flag, double x[9]);
+ int bcf_em1(const bcf1_t *b, int n1, int flag, double x[10]);
#ifdef __cplusplus
}
print;
} else {
my @t = split;
- my @c = (0);
+ my @c = (0, 0);
my $n = 0;
my $s = -1;
@_ = split(":", $t[8]);
}
sub varFilter {
- my %opts = (d=>2, D=>10000000, a=>2, W=>10, Q=>10, w=>10, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, G=>0, S=>1000);
- getopts('pd:D:W:Q:w:a:1:2:3:4:G:S:', \%opts);
+ my %opts = (d=>2, D=>10000000, a=>2, W=>10, Q=>10, w=>3, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, G=>0, S=>1000, e=>1e-4);
+ getopts('pd:D:W:Q:w:a:1:2:3:4:G:S:e:', \%opts);
die(qq/
Usage: vcfutils.pl varFilter [options] <in.vcf>
-2 FLOAT min P-value for baseQ bias [$opts{2}]
-3 FLOAT min P-value for mapQ bias [$opts{3}]
-4 FLOAT min P-value for end distance bias [$opts{4}]
+ -e FLOAT min P-value for HWE (plus F<0) [$opts{e}]
-p print filtered variants
Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
$flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
&& ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
$flt = 8 if ($flt == 0 && ((/MXGQ=(\d+)/ && $1 < $opts{G}) || (/MXSP=(\d+)/ && $1 >= $opts{S})));
+ # HWE filter
+ if ($t[7] =~ /G3=([^;,]+),([^;,]+),([^;,]+).*HWE=([^;,]+)/ && $4 < $opts{e}) {
+ my $p = 2*$1 + $2;
+ my $f = ($p > 0 && $p < 1)? 1 - $2 / ($p * (1-$p)) : 0;
+ $flt = 9 if ($f < 0);
+ }
my $score = $t[5] * 100 + $dp_alt;
my $rlen = length($t[3]) - 1; # $indel_score<0 for SNPs
my ($b, $q);
$q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
if ($q < 0) {
- $_ = $1 if ($t[7] =~ /AF1=([\d\.]+)/);
+ $_ = ($t[7] =~ /AF1=([\d\.]+)/)? $1 : 0;
$b = ($_ < .5 || $alt eq '.')? $ref : $alt;
$q = -$q;
} else {
if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0])) {
beg = atoi(str->s); // begin
if (dret != '\n') {
- if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0]))
+ if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0])) {
end = atoi(str->s); // end
+ if (end < beg) end = -1;
+ }
}
}
}
#include "khash.h"
typedef struct {
- uint64_t len:32, line_len:16, line_blen:16;
+ int32_t line_len, line_blen;
+ int64_t len;
uint64_t offset;
} faidx1_t;
KHASH_MAP_INIT_STR(s, faidx1_t)
{
char c, *name;
int l_name, m_name, ret;
- int len, line_len, line_blen, state;
+ int line_len, line_blen, state;
int l1, l2;
faidx_t *idx;
uint64_t offset;
+ int64_t len;
idx = (faidx_t*)calloc(1, sizeof(faidx_t));
idx->hash = kh_init(s);
return 0;
}
++l1; len += l2;
- if (l2 >= 0x10000) {
- fprintf(stderr, "[fai_build_core] line length exceeds 65535 in sequence '%s'.\n", name);
- free(name); fai_destroy(idx);
- return 0;
- }
if (state == 1) line_len = l1, line_blen = l2, state = 0;
else if (state == 0) {
if (l1 != line_len || l2 != line_blen) state = 2;
char *fai_fetch(const faidx_t *fai, const char *str, int *len)
{
- char *s, *p, c;
- int i, l, k;
+ char *s, c;
+ int i, l, k, name_end;
khiter_t iter;
faidx1_t val;
khash_t(s) *h;
beg = end = -1;
h = fai->hash;
- l = strlen(str);
- p = s = (char*)malloc(l+1);
- /* squeeze out "," */
- for (i = k = 0; i != l; ++i)
- if (str[i] != ',' && !isspace(str[i])) s[k++] = str[i];
- s[k] = 0;
- for (i = 0; i != k; ++i) if (s[i] == ':') break;
- s[i] = 0;
- iter = kh_get(s, h, s); /* get the ref_id */
- if (iter == kh_end(h)) {
- *len = 0;
- free(s); return 0;
- }
+ name_end = l = strlen(str);
+ s = (char*)malloc(l+1);
+ // remove space
+ for (i = k = 0; i < l; ++i)
+ if (!isspace(str[i])) s[k++] = str[i];
+ s[k] = 0; l = k;
+ // determine the sequence name
+ for (i = l - 1; i >= 0; --i) if (s[i] == ':') break; // look for colon from the end
+ if (i >= 0) name_end = i;
+ if (name_end < l) { // check if this is really the end
+ int n_hyphen = 0;
+ for (i = name_end + 1; i < l; ++i) {
+ if (s[i] == '-') ++n_hyphen;
+ else if (!isdigit(s[i]) && s[i] != ',') break;
+ }
+ if (i < l || n_hyphen > 1) name_end = l; // malformated region string; then take str as the name
+ s[name_end] = 0;
+ iter = kh_get(s, h, s);
+ if (iter == kh_end(h)) { // cannot find the sequence name
+ iter = kh_get(s, h, str); // try str as the name
+ if (iter == kh_end(h)) {
+ *len = 0;
+ free(s); return 0;
+ } else s[name_end] = ':', name_end = l;
+ }
+ } else iter = kh_get(s, h, str);
val = kh_value(h, iter);
- if (i == k) { /* dump the whole sequence */
- beg = 0; end = val.len;
- } else {
- for (p = s + i + 1; i != k; ++i) if (s[i] == '-') break;
- beg = atoi(p);
- if (i < k) {
- p = s + i + 1;
- end = atoi(p);
- } else end = val.len;
- }
- if (beg > 0) --beg;
+ // parse the interval
+ if (name_end < l) {
+ for (i = k = name_end + 1; i < l; ++i)
+ if (s[i] != ',') s[k++] = s[i];
+ s[k] = 0;
+ beg = atoi(s + name_end + 1);
+ for (i = name_end + 1; i != k; ++i) if (s[i] == '-') break;
+ end = i < k? atoi(s + i + 1) : val.len;
+ if (beg > 0) --beg;
+ } else beg = 0, end = val.len;
if (beg >= val.len) beg = val.len;
if (end >= val.len) end = val.len;
if (beg > end) beg = end;
+++ /dev/null
-#include <string.h>
-#include <stdlib.h>
-#include "glf.h"
-
-#ifdef _NO_BGZF
-// then alias bgzf_*() functions
-#endif
-
-static int glf3_is_BE = 0;
-
-static inline uint32_t bam_swap_endian_4(uint32_t v)
-{
- v = ((v & 0x0000FFFFU) << 16) | (v >> 16);
- return ((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8);
-}
-
-static inline uint16_t bam_swap_endian_2(uint16_t v)
-{
- return (uint16_t)(((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8));
-}
-
-static inline int bam_is_big_endian()
-{
- long one= 1;
- return !(*((char *)(&one)));
-}
-
-glf3_header_t *glf3_header_init()
-{
- glf3_is_BE = bam_is_big_endian();
- return (glf3_header_t*)calloc(1, sizeof(glf3_header_t));
-}
-
-glf3_header_t *glf3_header_read(glfFile fp)
-{
- glf3_header_t *h;
- char magic[4];
- h = glf3_header_init();
- bgzf_read(fp, magic, 4);
- if (strncmp(magic, "GLF\3", 4)) {
- fprintf(stderr, "[glf3_header_read] invalid magic.\n");
- glf3_header_destroy(h);
- return 0;
- }
- bgzf_read(fp, &h->l_text, 4);
- if (glf3_is_BE) h->l_text = bam_swap_endian_4(h->l_text);
- if (h->l_text) {
- h->text = (uint8_t*)calloc(h->l_text + 1, 1);
- bgzf_read(fp, h->text, h->l_text);
- }
- return h;
-}
-
-void glf3_header_write(glfFile fp, const glf3_header_t *h)
-{
- int32_t x;
- bgzf_write(fp, "GLF\3", 4);
- x = glf3_is_BE? bam_swap_endian_4(h->l_text) : h->l_text;
- bgzf_write(fp, &x, 4);
- if (h->l_text) bgzf_write(fp, h->text, h->l_text);
-}
-
-void glf3_header_destroy(glf3_header_t *h)
-{
- free(h->text);
- free(h);
-}
-
-char *glf3_ref_read(glfFile fp, int *len)
-{
- int32_t n, x;
- char *str;
- *len = 0;
- if (bgzf_read(fp, &n, 4) != 4) return 0;
- if (glf3_is_BE) n = bam_swap_endian_4(n);
- if (n < 0) {
- fprintf(stderr, "[glf3_ref_read] invalid reference name length: %d.\n", n);
- return 0;
- }
- str = (char*)calloc(n + 1, 1); // not necesarily n+1 in fact
- x = bgzf_read(fp, str, n);
- x += bgzf_read(fp, len, 4);
- if (x != n + 4) {
- free(str); *len = -1; return 0; // truncated
- }
- if (glf3_is_BE) *len = bam_swap_endian_4(*len);
- return str;
-}
-
-void glf3_ref_write(glfFile fp, const char *str, int len)
-{
- int32_t m, n = strlen(str) + 1;
- m = glf3_is_BE? bam_swap_endian_4(n) : n;
- bgzf_write(fp, &m, 4);
- bgzf_write(fp, str, n);
- if (glf3_is_BE) len = bam_swap_endian_4(len);
- bgzf_write(fp, &len, 4);
-}
-
-void glf3_view1(const char *ref_name, const glf3_t *g3, int pos)
-{
- int j;
- if (g3->rtype == GLF3_RTYPE_END) return;
- printf("%s\t%d\t%c\t%d\t%d\t%d", ref_name, pos + 1,
- g3->rtype == GLF3_RTYPE_INDEL? '*' : "XACMGRSVTWYHKDBN"[g3->ref_base],
- g3->depth, g3->rms_mapQ, g3->min_lk);
- if (g3->rtype == GLF3_RTYPE_SUB)
- for (j = 0; j != 10; ++j) printf("\t%d", g3->lk[j]);
- else {
- printf("\t%d\t%d\t%d\t%d\t%d\t%s\t%s\t", g3->lk[0], g3->lk[1], g3->lk[2], g3->indel_len[0], g3->indel_len[1],
- g3->indel_len[0]? g3->indel_seq[0] : "*", g3->indel_len[1]? g3->indel_seq[1] : "*");
- }
- printf("\n");
-}
-
-int glf3_write1(glfFile fp, const glf3_t *g3)
-{
- int r;
- uint8_t c;
- uint32_t y[2];
- c = g3->rtype<<4 | g3->ref_base;
- r = bgzf_write(fp, &c, 1);
- if (g3->rtype == GLF3_RTYPE_END) return r;
- y[0] = g3->offset;
- y[1] = g3->min_lk<<24 | g3->depth;
- if (glf3_is_BE) {
- y[0] = bam_swap_endian_4(y[0]);
- y[1] = bam_swap_endian_4(y[1]);
- }
- r += bgzf_write(fp, y, 8);
- r += bgzf_write(fp, &g3->rms_mapQ, 1);
- if (g3->rtype == GLF3_RTYPE_SUB) r += bgzf_write(fp, g3->lk, 10);
- else {
- int16_t x[2];
- r += bgzf_write(fp, g3->lk, 3);
- x[0] = glf3_is_BE? bam_swap_endian_2(g3->indel_len[0]) : g3->indel_len[0];
- x[1] = glf3_is_BE? bam_swap_endian_2(g3->indel_len[1]) : g3->indel_len[1];
- r += bgzf_write(fp, x, 4);
- if (g3->indel_len[0]) r += bgzf_write(fp, g3->indel_seq[0], abs(g3->indel_len[0]));
- if (g3->indel_len[1]) r += bgzf_write(fp, g3->indel_seq[1], abs(g3->indel_len[1]));
- }
- return r;
-}
-
-#ifndef kv_roundup32
-#define kv_roundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
-#endif
-
-int glf3_read1(glfFile fp, glf3_t *g3)
-{
- int r;
- uint8_t c;
- uint32_t y[2];
- r = bgzf_read(fp, &c, 1);
- if (r == 0) return 0;
- g3->ref_base = c & 0xf;
- g3->rtype = c>>4;
- if (g3->rtype == GLF3_RTYPE_END) return r;
- r += bgzf_read(fp, y, 8);
- if (glf3_is_BE) {
- y[0] = bam_swap_endian_4(y[0]);
- y[1] = bam_swap_endian_4(y[1]);
- }
- g3->offset = y[0];
- g3->min_lk = y[1]>>24;
- g3->depth = y[1]<<8>>8;
- r += bgzf_read(fp, &g3->rms_mapQ, 1);
- if (g3->rtype == GLF3_RTYPE_SUB) r += bgzf_read(fp, g3->lk, 10);
- else {
- int16_t x[2], max;
- r += bgzf_read(fp, g3->lk, 3);
- r += bgzf_read(fp, x, 4);
- if (glf3_is_BE) {
- x[0] = bam_swap_endian_2(x[0]);
- x[1] = bam_swap_endian_2(x[1]);
- }
- g3->indel_len[0] = x[0];
- g3->indel_len[1] = x[1];
- x[0] = abs(x[0]); x[1] = abs(x[1]);
- max = (x[0] > x[1]? x[0] : x[1]) + 1;
- if (g3->max_len < max) {
- g3->max_len = max;
- kv_roundup32(g3->max_len);
- g3->indel_seq[0] = (char*)realloc(g3->indel_seq[0], g3->max_len);
- g3->indel_seq[1] = (char*)realloc(g3->indel_seq[1], g3->max_len);
- }
- r += bgzf_read(fp, g3->indel_seq[0], x[0]);
- r += bgzf_read(fp, g3->indel_seq[1], x[1]);
- g3->indel_seq[0][x[0]] = g3->indel_seq[1][x[1]] = 0;
- }
- return r;
-}
-
-void glf3_view(glfFile fp)
-{
- glf3_header_t *h;
- char *name;
- glf3_t *g3;
- int len;
- h = glf3_header_read(fp);
- g3 = glf3_init1();
- while ((name = glf3_ref_read(fp, &len)) != 0) {
- int pos = 0;
- while (glf3_read1(fp, g3) && g3->rtype != GLF3_RTYPE_END) {
- pos += g3->offset;
- glf3_view1(name, g3, pos);
- }
- free(name);
- }
- glf3_header_destroy(h);
- glf3_destroy1(g3);
-}
-
-int glf3_view_main(int argc, char *argv[])
-{
- glfFile fp;
- if (argc == 1) {
- fprintf(stderr, "Usage: glfview <in.glf>\n");
- return 1;
- }
- fp = (strcmp(argv[1], "-") == 0)? bgzf_fdopen(fileno(stdin), "r") : bgzf_open(argv[1], "r");
- if (fp == 0) {
- fprintf(stderr, "Fail to open file '%s'\n", argv[1]);
- return 1;
- }
- glf3_view(fp);
- bgzf_close(fp);
- return 0;
-}
-
-#ifdef GLFVIEW_MAIN
-int main(int argc, char *argv[])
-{
- return glf3_view_main(argc, argv);
-}
-#endif
+++ /dev/null
-#ifndef GLF_H_
-#define GLF_H_
-
-typedef struct {
- unsigned char ref_base:4, dummy:4; /** "XACMGRSVTWYHKDBN"[ref_base] gives the reference base */
- unsigned char max_mapQ; /** maximum mapping quality */
- unsigned char lk[10]; /** log likelihood ratio, capped at 255 */
- unsigned min_lk:8, depth:24; /** minimum lk capped at 255, and the number of mapped reads */
-} glf1_t;
-
-#include <stdint.h>
-#include "bgzf.h"
-typedef BGZF *glfFile;
-
-#define GLF3_RTYPE_END 0
-#define GLF3_RTYPE_SUB 1
-#define GLF3_RTYPE_INDEL 2
-
-typedef struct {
- uint8_t ref_base:4, rtype:4; /** "XACMGRSVTWYHKDBN"[ref_base] gives the reference base */
- uint8_t rms_mapQ; /** RMS mapping quality */
- uint8_t lk[10]; /** log likelihood ratio, capped at 255 */
- uint32_t min_lk:8, depth:24; /** minimum lk capped at 255, and the number of mapped reads */
- int32_t offset; /** the first base in a chromosome has offset zero. */
- // for indel (lkHom1, lkHom2 and lkHet are the first three elements in lk[10])
- int16_t indel_len[2];
- int32_t max_len; // maximum indel len; will be modified by glf3_read1()
- char *indel_seq[2];
-} glf3_t;
-
-typedef struct {
- int32_t l_text;
- uint8_t *text;
-} glf3_header_t;
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-#define glf3_init1() ((glf3_t*)calloc(1, sizeof(glf3_t)))
-#define glf3_destroy1(g3) do { free((g3)->indel_seq[0]); free((g3)->indel_seq[1]); free(g3); } while (0)
-
- glf3_header_t *glf3_header_init();
- glf3_header_t *glf3_header_read(glfFile fp);
- void glf3_header_write(glfFile fp, const glf3_header_t *h);
- void glf3_header_destroy(glf3_header_t *h);
- char *glf3_ref_read(glfFile fp, int *len);
- void glf3_ref_write(glfFile fp, const char *name, int len);
- int glf3_write1(glfFile fp, const glf3_t *g3);
- int glf3_read1(glfFile fp, glf3_t *g3);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
static rghash_t g_rghash = 0;
static int g_min_mapQ = 0, g_flag_on = 0, g_flag_off = 0;
static char *g_library, *g_rg;
-static int g_sol2sanger_tbl[128];
static void *g_bed;
void *bed_read(const char *fn);
void bed_destroy(void *_h);
int bed_overlap(const void *_h, const char *chr, int beg, int end);
-static void sol2sanger(bam1_t *b)
-{
- int l;
- uint8_t *qual = bam1_qual(b);
- if (g_sol2sanger_tbl[30] == 0) {
- for (l = 0; l != 128; ++l) {
- g_sol2sanger_tbl[l] = (int)(10.0 * log(1.0 + pow(10.0, (l - 64 + 33) / 10.0)) / log(10.0) + .499);
- if (g_sol2sanger_tbl[l] >= 93) g_sol2sanger_tbl[l] = 93;
- }
- }
- for (l = 0; l < b->core.l_qseq; ++l) {
- int q = qual[l];
- if (q > 127) q = 127;
- qual[l] = g_sol2sanger_tbl[q];
- }
-}
-
static inline int __g_skip_aln(const bam_header_t *h, const bam1_t *b)
{
if (b->core.qual < g_min_mapQ || ((b->core.flag & g_flag_on) != g_flag_on) || (b->core.flag & g_flag_off))
int main_samview(int argc, char *argv[])
{
- int c, is_header = 0, is_header_only = 0, is_bamin = 1, ret = 0, compress_level = -1, is_bamout = 0, slx2sngr = 0, is_count = 0;
+ int c, is_header = 0, is_header_only = 0, is_bamin = 1, ret = 0, compress_level = -1, is_bamout = 0, is_count = 0;
int of_type = BAM_OFDEC, is_long_help = 0;
int count = 0;
samfile_t *in = 0, *out = 0;
/* parse command-line options */
strcpy(in_mode, "r"); strcpy(out_mode, "w");
- while ((c = getopt(argc, argv, "Sbct:h1Ho:q:f:F:ul:r:xX?T:CR:L:")) >= 0) {
+ while ((c = getopt(argc, argv, "Sbct:h1Ho:q:f:F:ul:r:xX?T:R:L:")) >= 0) {
switch (c) {
case 'c': is_count = 1; break;
- case 'C': slx2sngr = 1; break;
case 'S': is_bamin = 0; break;
case 'b': is_bamout = 1; break;
case 't': fn_list = strdup(optarg); is_bamin = 0; break;
int r;
while ((r = samread(in, b)) >= 0) { // read one alignment from `in'
if (!__g_skip_aln(in->header, b)) {
- if (slx2sngr) sol2sanger(b);
if (!is_count) samwrite(out, b); // write the alignment to `out'
count++;
}
free(argv2);
return ret;
}
+
+int8_t seq_comp_table[16] = { 0, 8, 4, 12, 2, 10, 9, 14, 1, 6, 5, 13, 3, 11, 7, 15 };
+
+int main_bam2fq(int argc, char *argv[])
+{
+ bamFile fp;
+ bam_header_t *h;
+ bam1_t *b;
+ int8_t *buf;
+ int max_buf;
+ if (argc == 1) {
+ fprintf(stderr, "Usage: samtools bam2fq <in.bam>\n");
+ return 1;
+ }
+ fp = strcmp(argv[1], "-")? bam_open(argv[1], "r") : bam_dopen(fileno(stdin), "r");
+ if (fp == 0) return 1;
+ h = bam_header_read(fp);
+ b = bam_init1();
+ buf = 0;
+ max_buf = 0;
+ while (bam_read1(fp, b) >= 0) {
+ int i, qlen = b->core.l_qseq;
+ uint8_t *seq;
+ putchar('@'); fputs(bam1_qname(b), stdout);
+ if ((b->core.flag & 0x40) && !(b->core.flag & 0x80)) puts("/1");
+ else if ((b->core.flag & 0x80) && !(b->core.flag & 0x40)) puts("/2");
+ else putchar('\n');
+ if (max_buf < qlen + 1) {
+ max_buf = qlen + 1;
+ kroundup32(max_buf);
+ buf = realloc(buf, max_buf);
+ }
+ buf[qlen] = 0;
+ seq = bam1_seq(b);
+ for (i = 0; i < qlen; ++i)
+ buf[i] = bam1_seqi(seq, i);
+ if (b->core.flag & 16) { // reverse complement
+ for (i = 0; i < qlen>>1; ++i) {
+ int8_t t = seq_comp_table[buf[qlen - 1 - i]];
+ buf[qlen - 1 - i] = seq_comp_table[buf[i]];
+ buf[i] = t;
+ }
+ if (qlen&1) buf[i] = seq_comp_table[buf[i]];
+ }
+ for (i = 0; i < qlen; ++i)
+ buf[i] = bam_nt16_rev_table[buf[i]];
+ puts((char*)buf);
+ puts("+");
+ seq = bam1_qual(b);
+ for (i = 0; i < qlen; ++i)
+ buf[i] = 33 + seq[i];
+ if (b->core.flag & 16) { // reverse
+ for (i = 0; i < qlen>>1; ++i) {
+ int8_t t = buf[qlen - 1 - i];
+ buf[qlen - 1 - i] = buf[i];
+ buf[i] = t;
+ }
+ }
+ puts((char*)buf);
+ }
+ free(buf);
+ bam_destroy1(b);
+ bam_header_destroy(h);
+ bam_close(fp);
+ return 0;
+}
kstring_t buf;
int n = 0;
khash_t(sm) *sm2id = (khash_t(sm)*)sm->sm2id;
+ if (txt == 0) {
+ add_pair(sm, sm2id, fn, fn);
+ return 0;
+ }
memset(&buf, 0, sizeof(kstring_t));
while ((q = strstr(p, "@RG")) != 0) {
p = q + 3;
-.TH samtools 1 "21 April 2011" "samtools-0.1.16" "Bioinformatics tools"
+.TH samtools 1 "05 July 2011" "samtools-0.1.17" "Bioinformatics tools"
.SH NAME
.PP
samtools - Utilities for the Sequence Alignment/Map (SAM) format
+
+bcftools - Utilities for the Binary Call Format (BCF) and VCF
.SH SYNOPSIS
.PP
samtools view -bt ref_list.txt -o aln.bam aln.sam.gz
samtools mpileup -C50 -gf ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam
.PP
samtools tview aln.sorted.bam ref.fasta
+.PP
+bcftools index in.bcf
+.PP
+bcftools view in.bcf chr2:100-200 > out.vcf
+.PP
+bcftools view -vc in.bcf > out.vcf 2> out.afs
.SH DESCRIPTION
.PP
will download the index upon absence. Samtools does not retrieve the
entire alignment file unless it is asked to do so.
-.SH COMMANDS AND OPTIONS
+.SH SAMTOOLS COMMANDS AND OPTIONS
.TP 10
.B view
.TP
.B mpileup
-samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]
+.B samtools mpileup
+.RB [ \-EBug ]
+.RB [ \-C
+.IR capQcoef ]
+.RB [ \-r
+.IR reg ]
+.RB [ \-f
+.IR in.fa ]
+.RB [ \-l
+.IR list ]
+.RB [ \-M
+.IR capMapQ ]
+.RB [ \-Q
+.IR minBaseQ ]
+.RB [ \-q
+.IR minMapQ ]
+.I in.bam
+.RI [ in2.bam
+.RI [ ... ]]
Generate BCF or pileup for one or multiple BAM files. Alignment records
are grouped by sample identifiers in @RG header lines. If sample
identifiers are absent, each input file is regarded as one sample.
-.B OPTIONS:
+In the pileup format (without
+.BR -u or -g ),
+each
+line represents a genomic position, consisting of chromosome name,
+coordinate, reference base, read bases, read qualities and alignment
+mapping qualities. Information on match, mismatch, indel, strand,
+mapping quality and start and end of a read are all encoded at the read
+base column. At this column, a dot stands for a match to the reference
+base on the forward strand, a comma for a match on the reverse strand,
+a '>' or '<' for a reference skip, `ACGTN' for a mismatch on the forward
+strand and `acgtn' for a mismatch on the reverse strand. A pattern
+`\\+[0-9]+[ACGTNacgtn]+' indicates there is an insertion between this
+reference position and the next reference position. The length of the
+insertion is given by the integer in the pattern, followed by the
+inserted sequence. Similarly, a pattern `-[0-9]+[ACGTNacgtn]+'
+represents a deletion from the reference. The deleted bases will be
+presented as `*' in the following lines. Also at the read base column, a
+symbol `^' marks the start of a read. The ASCII of the character
+following `^' minus 33 gives the mapping quality. A symbol `$' marks the
+end of a read segment.
+
+.B Input Options:
.RS
.TP 10
+.B -6
+Assume the quality is in the Illumina 1.3+ encoding.
.B -A
Do not skip anomalous read pairs in variant calling.
.TP
misaligned. Applying this option greatly helps to reduce false SNPs
caused by misalignments.
.TP
+.BI -b \ FILE
+List of input BAM files, one file per line [null]
+.TP
.BI -C \ INT
Coefficient for downgrading mapping quality for reads containing
excessive mismatches. Given a read with a phred-scaled probability q of
.I INT
reads per input BAM. [250]
.TP
-.B -D
-Output per-sample read depth
-.TP
-.BI -e \ INT
-Phred-scaled gap extension sequencing error probability. Reducing
-.I INT
-leads to longer indels. [20]
-.TP
.B -E
Extended BAQ computation. This option helps sensitivity especially for MNPs, but may hurt
specificity a little bit.
.TP
.BI -f \ FILE
-The reference file [null]
+The
+.BR faidx -indexed
+reference file in the FASTA format. The file can be optionally compressed by
+.BR razip .
+[null]
+.TP
+.BI -l \ FILE
+BED or position list file containing a list of regions or sites where pileup or BCF should be generated [null]
+.TP
+.BI -q \ INT
+Minimum mapping quality for an alignment to be used [0]
+.TP
+.BI -Q \ INT
+Minimum base quality for a base to be considered [13]
+.TP
+.BI -r \ STR
+Only generate pileup in region
+.I STR
+[all sites]
+.TP
+.B Output Options:
+
+.TP
+.B -D
+Output per-sample read depth
.TP
.B -g
Compute genotype likelihoods and output them in the binary call format (BCF).
.TP
+.B -S
+Output per-sample Phred-scaled strand bias P-value
+.TP
+.B -u
+Similar to
+.B -g
+except that the output is uncompressed BCF, which is preferred for piping.
+
+.TP
+.B Options for Genotype Likelihood Computation (for -g or -u):
+
+.TP
+.BI -e \ INT
+Phred-scaled gap extension sequencing error probability. Reducing
+.I INT
+leads to longer indels. [20]
+.TP
.BI -h \ INT
Coefficient for modeling homopolymer errors. Given an
.IR l -long
.B -I
Do not perform INDEL calling
.TP
-.BI -l \ FILE
-File containing a list of sites where pileup or BCF is outputted [null]
-.TP
.BI -L \ INT
Skip INDEL calling if the average per-sample depth is above
.IR INT .
from which indel candidates are obtained. It is recommended to collect
indel candidates from sequencing technologies that have low indel error
rate such as ILLUMINA. [all]
-.TP
-.BI -q \ INT
-Minimum mapping quality for an alignment to be used [0]
-.TP
-.BI -Q \ INT
-Minimum base quality for a base to be considered [13]
-.TP
-.BI -r \ STR
-Only generate pileup in region
-.I STR
-[all sites]
-.TP
-.B -S
-Output per-sample Phred-scaled strand bias P-value
-.TP
-.B -u
-Similar to
-.B -g
-except that the output is uncompressed BCF, which is preferred for piping.
.RE
.TP
Minimum base quality to be used in het calling. [13]
.RE
-.TP
-.B pileup
-samtools pileup [-2sSBicv] [-f in.ref.fasta] [-t in.ref_list] [-l
-in.site_list] [-C capMapQ] [-M maxMapQ] [-T theta] [-N nHap] [-r
-pairDiffRate] [-m mask] [-d maxIndelDepth] [-G indelPrior]
-<in.bam>|<in.sam>
+.SH BCFTOOLS COMMANDS AND OPTIONS
-Print the alignment in the pileup format. In the pileup format, each
-line represents a genomic position, consisting of chromosome name,
-coordinate, reference base, read bases, read qualities and alignment
-mapping qualities. Information on match, mismatch, indel, strand,
-mapping quality and start and end of a read are all encoded at the read
-base column. At this column, a dot stands for a match to the reference
-base on the forward strand, a comma for a match on the reverse strand,
-a '>' or '<' for a reference skip, `ACGTN' for a mismatch on the forward
-strand and `acgtn' for a mismatch on the reverse strand. A pattern
-`\\+[0-9]+[ACGTNacgtn]+' indicates there is an insertion between this
-reference position and the next reference position. The length of the
-insertion is given by the integer in the pattern, followed by the
-inserted sequence. Similarly, a pattern `-[0-9]+[ACGTNacgtn]+'
-represents a deletion from the reference. The deleted bases will be
-presented as `*' in the following lines. Also at the read base column, a
-symbol `^' marks the start of a read. The ASCII of the character
-following `^' minus 33 gives the mapping quality. A symbol `$' marks the
-end of a read segment.
-
-If option
-.B -c
-is applied, the consensus base, Phred-scaled consensus quality, SNP
-quality (i.e. the Phred-scaled probability of the consensus being
-identical to the reference) and root mean square (RMS) mapping quality
-of the reads covering the site will be inserted between the `reference
-base' and the `read bases' columns. An indel occupies an additional
-line. Each indel line consists of chromosome name, coordinate, a star,
-the genotype, consensus quality, SNP quality, RMS mapping quality, #
-covering reads, the first alllele, the second allele, # reads supporting
-the first allele, # reads supporting the second allele and # reads
-containing indels different from the top two alleles.
-
-.B NOTE:
-Since 0.1.10, the `pileup' command is deprecated by `mpileup'.
+.TP 10
+.B view
+.B bcftools view
+.RB [ \-AbFGNQSucgv ]
+.RB [ \-D
+.IR seqDict ]
+.RB [ \-l
+.IR listLoci ]
+.RB [ \-s
+.IR listSample ]
+.RB [ \-i
+.IR gapSNPratio ]
+.RB [ \-t
+.IR mutRate ]
+.RB [ \-p
+.IR varThres ]
+.RB [ \-P
+.IR prior ]
+.RB [ \-1
+.IR nGroup1 ]
+.RB [ \-d
+.IR minFrac ]
+.RB [ \-U
+.IR nPerm ]
+.RB [ \-X
+.IR permThres ]
+.RB [ \-T
+.IR trioType ]
+.I in.bcf
+.RI [ region ]
+
+Convert between BCF and VCF, call variant candidates and estimate allele
+frequencies.
-.B OPTIONS:
.RS
+.TP
+.B Input/Output Options:
.TP 10
-.B -B
-Disable the BAQ computation. See the
-.B mpileup
-command for details.
+.B -A
+Retain all possible alternate alleles at variant sites. By default, the view
+command discards unlikely alleles.
+.TP 10
+.B -b
+Output in the BCF format. The default is VCF.
.TP
-.B -c
-Call the consensus sequence. Options
-.BR -T ", " -N ", " -I " and " -r
-are only effective when
-.BR -c " or " -g
-is in use.
+.BI -D \ FILE
+Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null]
.TP
-.BI -C \ INT
-Coefficient for downgrading the mapping quality of poorly mapped
-reads. See the
-.B mpileup
-command for details. [0]
+.B -F
+Indicate PL is generated by r921 or before (ordering is different).
.TP
-.BI -d \ INT
-Use the first
-.I NUM
-reads in the pileup for indel calling for speed up. Zero for unlimited. [1024]
+.B -G
+Suppress all individual genotype information.
.TP
-.BI -f \ FILE
-The reference sequence in the FASTA format. Index file
-.I FILE.fai
-will be created if
-absent.
+.BI -l \ FILE
+List of sites at which information are outputted [all sites]
.TP
-.B -g
-Generate genotype likelihood in the binary GLFv3 format. This option
-suppresses -c, -i and -s. This option is deprecated by the
-.B mpileup
-command.
+.B -N
+Skip sites where the REF field is not A/C/G/T
.TP
-.B -i
-Only output pileup lines containing indels.
+.B -Q
+Output the QCALL likelihood format
.TP
-.BI -I \ INT
-Phred probability of an indel in sequencing/prep. [40]
+.BI -s \ FILE
+List of samples to use. The first column in the input gives the sample names
+and the second gives the ploidy, which can only be 1 or 2. When the 2nd column
+is absent, the sample ploidy is assumed to be 2. In the output, the ordering of
+samples will be identical to the one in
+.IR FILE .
+[null]
.TP
-.BI -l \ FILE
-List of sites at which pileup is output. This file is space
-delimited. The first two columns are required to be chromosome and
-1-based coordinate. Additional columns are ignored. It is
-recommended to use option
+.B -S
+The input is VCF instead of BCF.
.TP
-.BI -m \ INT
-Filter reads with flag containing bits in
-.I INT
-[1796]
+.B -u
+Uncompressed BCF output (force -b).
.TP
-.BI -M \ INT
-Cap mapping quality at INT [60]
+.B Consensus/Variant Calling Options:
+.TP 10
+.B -c
+Call variants using Bayesian inference. This option automatically invokes option
+.BR -e .
.TP
-.BI -N \ INT
-Number of haplotypes in the sample (>=2) [2]
+.BI -d \ FLOAT
+When
+.B -v
+is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0]
.TP
-.BI -r \ FLOAT
-Expected fraction of differences between a pair of haplotypes [0.001]
+.B -e
+Perform max-likelihood inference only, including estimating the site allele frequency,
+testing Hardy-Weinberg equlibrium and testing associations with LRT.
.TP
-.B -s
-Print the mapping quality as the last column. This option makes the
-output easier to parse, although this format is not space efficient.
+.B -g
+Call per-sample genotypes at variant sites (force -c)
.TP
-.B -S
-The input file is in SAM.
+.BI -i \ FLOAT
+Ratio of INDEL-to-SNP mutation rate [0.15]
.TP
-.BI -t \ FILE
-List of reference names ane sequence lengths, in the format described
-for the
-.B import
-command. If this option is present, samtools assumes the input
-.I <in.alignment>
-is in SAM format; otherwise it assumes in BAM format.
+.BI -p \ FLOAT
+A site is considered to be a variant if P(ref|D)<FLOAT [0.5]
+.TP
+.BI -P \ STR
+Prior or initial allele frequency spectrum. If STR can be
+.IR full ,
+.IR cond2 ,
+.I flat
+or the file consisting of error output from a previous variant calling
+run.
+.TP
+.BI -t \ FLOAT
+Scaled muttion rate for variant calling [0.001]
+.TP
+.BI -T \ STR
+Enable pair/trio calling. For trio calling, option
.B -s
-together with
-.B -l
-as in the default format we may not know the mapping quality.
+is usually needed to be applied to configure the trio members and their ordering.
+In the file supplied to the option
+.BR -s ,
+the first sample must be the child, the second the father and the third the mother.
+The valid values of
+.I STR
+are `pair', `trioauto', `trioxd' and `trioxs', where `pair' calls differences between two input samples, and `trioxd' (`trioxs') specifies that the input
+is from the X chromosome non-PAR regions and the child is a female (male). [null]
.TP
-.BI -T \ FLOAT
-The theta parameter (error dependency coefficient) in the maq consensus
-calling model [0.85]
+.B -v
+Output variant sites only (force -c)
+.TP
+.B Contrast Calling and Association Test Options:
+.TP
+.BI -1 \ INT
+Number of group-1 samples. This option is used for dividing the samples into
+two groups for contrast SNP calling or association test.
+When this option is in use, the following VCF INFO will be outputted:
+PC2, PCHI2 and QCHI2. [0]
+.TP
+.BI -U \ INT
+Number of permutations for association test (effective only with
+.BR -1 )
+[0]
+.TP
+.BI -X \ FLOAT
+Only perform permutations for P(chi^2)<FLOAT (effective only with
+.BR -U )
+[0.01]
.RE
+.TP
+.B index
+.B bcftools index
+.I in.bcf
+
+Index sorted BCF for random access.
+.RE
+
+.TP
+.B cat
+.B bcftools cat
+.I in1.bcf
+.RI [ "in2.bcf " [ ... "]]]"
+
+Concatenate BCF files. The input files are required to be sorted and
+have identical samples appearing in the same order.
+.RE
.SH SAM FORMAT
-SAM is TAB-delimited. Apart from the header lines, which are started
+Sequence Alignment/Map (SAM) format is TAB-delimited. Apart from the header lines, which are started
with the `@' symbol, each alignment line consists of:
.TS
n | l | l .
Col Field Description
_
-1 QNAME Query (pair) NAME
+1 QNAME Query template/pair NAME
2 FLAG bitwise FLAG
3 RNAME Reference sequence NAME
4 POS 1-based leftmost POSition/coordinate of clipped sequence
6 CIAGR extended CIGAR string
7 MRNM Mate Reference sequence NaMe (`=' if same as RNAME)
8 MPOS 1-based Mate POSistion
-9 ISIZE Inferred insert SIZE
+9 TLEN inferred Template LENgth (insert size)
10 SEQ query SEQuence on the same strand as the reference
11 QUAL query QUALity (ASCII-33 gives the Phred base quality)
-12 OPT variable OPTional fields in the format TAG:VTYPE:VALUE
+12+ OPT variable OPTional fields in the format TAG:VTYPE:VALUE
.TE
.PP
0x0400 d the read is either a PCR or an optical duplicate
.TE
+where the second column gives the string representation of the FLAG field.
+
+.SH VCF FORMAT
+
+The Variant Call Format (VCF) is a TAB-delimited format with each data line consists of the following fields:
+.TS
+center box;
+cb | cb | cb
+n | l | l .
+Col Field Description
+_
+1 CHROM CHROMosome name
+2 POS the left-most POSition of the variant
+3 ID unique variant IDentifier
+4 REF the REFerence allele
+5 ALT the ALTernate allele(s), separated by comma
+6 QUAL variant/reference QUALity
+7 FILTER FILTers applied
+8 INFO INFOrmation related to the variant, separated by semi-colon
+9 FORMAT FORMAT of the genotype fields, separated by colon (optional)
+10+ SAMPLE SAMPLE genotypes and per-sample information (optional)
+.TE
+
+.PP
+The following table gives the
+.B INFO
+tags used by samtools and bcftools.
+
+.TS
+center box;
+cb | cb | cb
+l | l | l .
+Tag Format Description
+_
+AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
+DP int Raw read depth (without quality filtering)
+DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
+FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise
+MQ int Root-Mean-Square mapping quality of covering reads
+PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
+PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
+PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
+QCHI2 int Phred-scaled PCHI2
+RP int # permutations yielding a smaller PCHI2
+CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint
+UGT string Most probable genotype configuration without the trio constraint
+CGT string Most probable configuration with the trio constraint
+.TE
+
.SH EXAMPLES
.IP o 2
Import SAM to BAM when
.IR RG:Z:454 .
.IP o 2
-Call SNPs and short indels for one diploid individual:
+Call SNPs and short INDELs for one diploid individual:
samtools mpileup -ugf ref.fa aln.bam | bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -D 100 > var.flt.vcf
samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq
+.IP o 2
+Call somatic mutations from a pair of samples:
+
+ samtools mpileup -DSuf ref.fa aln.bam | bcftools view -bvcgT pair - > var.bcf
+
+In the output INFO field,
+.I CLR
+gives the Phred-log ratio between the likelihood by treating the
+two samples independently, and the likelihood by requiring the genotype to be identical.
+This
+.I CLR
+is effectively a score measuring the confidence of somatic calls. The higher the better.
+
+.IP o 2
+Call de novo and somatic mutations from a family trio:
+
+ samtools mpileup -DSuf ref.fa aln.bam | bcftools view -bvcgT pair -s samples.txt - > var.bcf
+
+File
+.I samples.txt
+should consist of three lines specifying the member and order of samples (in the order of child-father-mother).
+Similarly,
+.I CLR
+gives the Phred-log likelihood ratio with and without the trio constraint.
+.I UGT
+shows the most likely genotype configuration without the trio constraint, and
+.I CGT
+gives the most likely genotype configuration satisfying the trio constraint.
+
.IP o 2
Phase one individual:
.IP o 2
Unaligned words used in bam_import.c, bam_endian.h, bam.c and bam_aux.c.
.IP o 2
-In merging, the input files are required to have the same number of
-reference sequences. The requirement can be relaxed. In addition,
-merging does not reconstruct the header dictionaries
-automatically. Endusers have to provide the correct header. Picard is
-better at merging.
-.IP o 2
Samtools paired-end rmdup does not work for unpaired reads (e.g. orphan
reads or ends mapped to different chromosomes). If this is a concern,
please use Picard's MarkDuplicate which correctly handles these cases,
projects / samtools.git / commitdiff