From b301e959d73eee0955c57004f344f17af00703f4 Mon Sep 17 00:00:00 2001 From: Charles Plessy Date: Mon, 11 Jul 2011 20:18:58 +0900 Subject: [PATCH] Imported Upstream version 0.1.17 --- Makefile | 12 +- NEWS | 44 +++ bam.h | 2 +- bam2bcf.c | 2 +- bam2bcf_indel.c | 12 +- bam_aux.c | 75 +++--- bam_import.c | 2 +- bam_index.c | 20 +- bam_maqcns.c | 626 ------------------------------------------- bam_maqcns.h | 61 ----- bam_plcmd.c | 520 +++-------------------------------- bam_tview.c | 34 ++- bamtk.c | 12 +- bcftools/Makefile | 4 +- bcftools/bcf.h | 4 + bcftools/bcftools.1 | 189 ------------- bcftools/bcfutils.c | 53 ++++ bcftools/call1.c | 81 +++++- bcftools/em.c | 15 +- bcftools/main.c | 2 + bcftools/mut.c | 124 +++++++++ bcftools/prob1.c | 22 +- bcftools/prob1.h | 5 +- bcftools/vcfutils.pl | 15 +- bedidx.c | 4 +- faidx.c | 75 +++--- glf.c | 236 ---------------- glf.h | 56 ---- sam_view.c | 90 +++++-- sample.c | 4 + samtools.1 | 476 +++++++++++++++++++++----------- 31 files changed, 926 insertions(+), 1951 deletions(-) delete mode 100644 bam_maqcns.c delete mode 100644 bam_maqcns.h delete mode 100644 bcftools/bcftools.1 create mode 100644 bcftools/mut.c delete mode 100644 glf.c delete mode 100644 glf.h diff --git a/Makefile b/Makefile index 4fc03d5..db18333 100644 --- a/Makefile +++ b/Makefile @@ -3,9 +3,9 @@ CFLAGS= -g -Wall -O2 #-m64 #-arch ppc 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 @@ -38,7 +38,7 @@ all:$(PROG) 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 @@ -54,14 +54,12 @@ bam.o:bam.h razf.h bam_endian.h kstring.h sam_header.h 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 diff --git a/NEWS b/NEWS index a600bb1..1aa2359 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,47 @@ +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) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/bam.h b/bam.h index e7360bb..246b020 100644 --- a/bam.h +++ b/bam.h @@ -40,7 +40,7 @@ @copyright Genome Research Ltd. */ -#define BAM_VERSION "0.1.16 (r963:234)" +#define BAM_VERSION "0.1.17 (r973:277)" #include #include diff --git a/bam2bcf.c b/bam2bcf.c index bd5503d..f263fad 100644 --- a/bam2bcf.c +++ b/bam2bcf.c @@ -236,7 +236,7 @@ int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bc 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; diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c index 899428f..2fdee9f 100644 --- a/bam2bcf_indel.c +++ b/bam2bcf_indel.c @@ -3,12 +3,14 @@ #include #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 @@ -110,7 +112,6 @@ static inline int est_indelreg(int pos, const char *ref, int l, char *ins4) 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; @@ -167,6 +168,12 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla 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; @@ -177,7 +184,6 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla 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; diff --git a/bam_aux.c b/bam_aux.c index b6a8d88..2247bdf 100644 --- a/bam_aux.c +++ b/bam_aux.c @@ -87,47 +87,56 @@ int32_t bam_get_tid(const bam_header_t *header, const char *seq_name) 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) diff --git a/bam_import.c b/bam_import.c index 8c24692..29516c9 100644 --- a/bam_import.c +++ b/bam_import.c @@ -445,7 +445,7 @@ int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b) 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"); diff --git a/bam_index.c b/bam_index.c index 51c2701..66d8eb8 100644 --- a/bam_index.c +++ b/bam_index.c @@ -172,17 +172,21 @@ bam_index_t *bam_index_core(bamFile fp) 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 @@ -203,7 +207,7 @@ bam_index_t *bam_index_core(bamFile fp) 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; @@ -222,7 +226,7 @@ bam_index_t *bam_index_core(bamFile fp) ++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; } } } @@ -473,6 +477,10 @@ int bam_index_build2(const char *fn, const char *_fnidx) } 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"); diff --git a/bam_maqcns.c b/bam_maqcns.c deleted file mode 100644 index 9931036..0000000 --- a/bam_maqcns.c +++ /dev/null @@ -1,626 +0,0 @@ -#include -#include -#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() = \theta \sum_{i=1}^{N-1} 1/i - P(D|) = \sum_{k=1}^{N-1} p_k 1/2 [(k/N)^n_2(1-k/N)^n_1 + (k/N)^n1(1-k/N)^n_2] - p_k = 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<= 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; -} diff --git a/bam_maqcns.h b/bam_maqcns.h deleted file mode 100644 index 291ae53..0000000 --- a/bam_maqcns.h +++ /dev/null @@ -1,61 +0,0 @@ -#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 diff --git a/bam_plcmd.c b/bam_plcmd.c index f6381c3..cbf6ae8 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -6,93 +6,8 @@ #include #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]; @@ -108,75 +23,6 @@ static inline int printw(int c, FILE *fp) 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; @@ -212,317 +58,6 @@ static inline void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, cons 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] |\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 #include "bam2bcf.h" #include "sample.h" @@ -536,6 +71,9 @@ int bam_pileup(int argc, char *argv[]) #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); @@ -610,7 +148,7 @@ static int mplp_func(void *data, bam1_t *b) } 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)); @@ -619,7 +157,7 @@ static void group_smpl(mplp_pileup_t *m, bam_sample_t *sm, kstring_t *buf, 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) { @@ -641,7 +179,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) 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; @@ -674,7 +212,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) 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; @@ -688,7 +226,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) 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); } @@ -736,7 +274,11 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) 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) @@ -752,7 +294,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) 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; } @@ -760,7 +302,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) 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) @@ -787,8 +329,10 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) 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'); @@ -798,6 +342,21 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) 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'); @@ -878,6 +437,7 @@ int bam_mpileup(int argc, char *argv[]) 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; @@ -885,7 +445,7 @@ int bam_mpileup(int argc, char *argv[]) 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); @@ -899,12 +459,14 @@ int bam_mpileup(int argc, char *argv[]) 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; @@ -939,6 +501,7 @@ int bam_mpileup(int argc, char *argv[]) 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"); @@ -946,11 +509,14 @@ int bam_mpileup(int argc, char *argv[]) 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"); diff --git a/bam_tview.c b/bam_tview.c index bf01e15..4eea955 100644 --- a/bam_tview.c +++ b/bam_tview.c @@ -19,9 +19,10 @@ #include #include #include +#include #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); @@ -50,7 +51,7 @@ typedef struct { 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; @@ -58,6 +59,7 @@ typedef struct { 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; @@ -70,11 +72,26 @@ int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void 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<>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 = '.'; @@ -191,9 +208,8 @@ tview_t *tv_init(const char *fn, const char *fn_fa) 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); @@ -224,7 +240,7 @@ void tv_destroy(tview_t *tv) 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); diff --git a/bamtk.c b/bamtk.c index 0941f6f..8ba2581 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,6 @@ #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[]); @@ -27,9 +26,9 @@ int main_cut_target(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() { @@ -39,7 +38,6 @@ static int usage() fprintf(stderr, "Usage: samtools [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"); @@ -49,7 +47,6 @@ static int usage() 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"); @@ -80,7 +77,6 @@ int main(int argc, char *argv[]) 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); @@ -89,7 +85,6 @@ int main(int argc, char *argv[]) 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); @@ -98,6 +93,11 @@ int main(int argc, char *argv[]) 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 diff --git a/bcftools/Makefile b/bcftools/Makefile index fd2e2df..9b6f863 100644 --- a/bcftools/Makefile +++ b/bcftools/Makefile @@ -1,7 +1,7 @@ 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 @@ -28,7 +28,7 @@ all:$(PROG) lib:libbcf.a libbcf.a:$(LOBJS) - $(AR) -cru $@ $(LOBJS) + $(AR) -csru $@ $(LOBJS) bcftools:lib $(AOBJS) $(CC) $(CFLAGS) -o $@ $(AOBJS) -L. $(LIBPATH) -lbcf -lm -lz diff --git a/bcftools/bcf.h b/bcftools/bcf.h index ba6dbe9..aeae6ca 100644 --- a/bcftools/bcf.h +++ b/bcftools/bcf.h @@ -28,6 +28,8 @@ #ifndef BCF_H #define BCF_H +#define BCF_VERSION "0.1.17 (r973:277)" + #include #include @@ -150,6 +152,8 @@ extern "C" { 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); diff --git a/bcftools/bcftools.1 b/bcftools/bcftools.1 deleted file mode 100644 index c6b4968..0000000 --- a/bcftools/bcftools.1 +++ /dev/null @@ -1,189 +0,0 @@ -.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)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; +} diff --git a/bcftools/call1.c b/bcftools/call1.c index 8e57aa1..b2d7d4a 100644 --- a/bcftools/call1.c +++ b/bcftools/call1.c @@ -26,9 +26,12 @@ KSTREAM_INIT(gzFile, gzread, 16384) #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; @@ -102,7 +105,7 @@ static void rm_info(bcf1_t *b, const char *key) 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; @@ -122,6 +125,14 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, 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); @@ -134,7 +145,8 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, 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; @@ -148,6 +160,7 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, 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]); @@ -229,13 +242,21 @@ static void write_header(bcf_hdr_t *h) if (!strstr(str.s, "##INFO=\n", &str); if (!strstr(str.s, "##INFO=\n", &str); + kputs("##INFO=\n", &str); + if (!strstr(str.s, "##INFO=\n", &str); if (!strstr(str.s, "##INFO=\n", &str); if (!strstr(str.s, "##INFO=\n", &str); - if (!strstr(str.s, "##INFO=\n", &str); + if (!strstr(str.s, "##INFO=\n", &str); + if (!strstr(str.s, "##INFO=\n", &str); + if (!strstr(str.s, "##INFO=\n", &str); +// if (!strstr(str.s, "##INFO=\n", &str); if (!strstr(str.s, "##INFO=\n", &str); if (!strstr(str.s, "##INFO== 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; @@ -303,6 +330,7 @@ int bcfview(int argc, char *argv[]) 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; @@ -317,6 +345,16 @@ int bcfview(int argc, char *argv[]) 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; @@ -351,6 +389,7 @@ int bcfview(int argc, char *argv[]) fprintf(stderr, " -p FLOAT variant if P(ref|D) 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); @@ -450,16 +490,25 @@ int bcfview(int argc, char *argv[]) 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); @@ -473,7 +522,7 @@ int bcfview(int argc, char *argv[]) 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; @@ -485,8 +534,8 @@ int bcfview(int argc, char *argv[]) } 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) { @@ -515,12 +564,16 @@ int bcfview(int argc, char *argv[]) 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; diff --git a/bcftools/em.c b/bcftools/em.c index f5d2fd9..32b400f 100644 --- a/bcftools/em.c +++ b/bcftools/em.c @@ -168,7 +168,8 @@ static double lk_ratio_test(int n, int n1, const double *pdg, double f3[3][3]) // 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; @@ -180,7 +181,7 @@ int bcf_em1(const bcf1_t *b, int n1, int flag, double x[9]) 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); @@ -188,7 +189,7 @@ int bcf_em1(const bcf1_t *b, int n1, int flag, double x[9]) } 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]); @@ -213,14 +214,16 @@ int bcf_em1(const bcf1_t *b, int n1, int flag, double x[9]) 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); diff --git a/bcftools/main.c b/bcftools/main.c index 7293374..fcd94b8 100644 --- a/bcftools/main.c +++ b/bcftools/main.c @@ -166,6 +166,8 @@ int main(int argc, char *argv[]) { 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 \n\n"); fprintf(stderr, "Command: view print, extract, convert and call SNPs from BCF\n"); fprintf(stderr, " index index BCF\n"); diff --git a/bcftools/mut.c b/bcftools/mut.c new file mode 100644 index 0000000..4b7cd10 --- /dev/null +++ b/bcftools/mut.c @@ -0,0 +1,124 @@ +#include +#include +#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<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<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; +} diff --git a/bcftools/prob1.c b/bcftools/prob1.c index fc9cb29..a380484 100644 --- a/bcftools/prob1.c +++ b/bcftools/prob1.c @@ -497,6 +497,13 @@ int bcf_p1_cal(const bcf1_t *b, int do_contrast, bcf_p1aux_t *ma, bcf_p1rst_t *r 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]; @@ -509,16 +516,27 @@ int bcf_p1_cal(const bcf1_t *b, int do_contrast, bcf_p1aux_t *ma, bcf_p1rst_t *r { // 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); diff --git a/bcftools/prob1.h b/bcftools/prob1.h index 571f42f..0a51a0a 100644 --- a/bcftools/prob1.h +++ b/bcftools/prob1.h @@ -8,9 +8,10 @@ typedef struct __bcf_p1aux_t bcf_p1aux_t; 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 @@ -32,7 +33,7 @@ extern "C" { 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 } diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl index 6ced66a..2b7ba0b 100755 --- a/bcftools/vcfutils.pl +++ b/bcftools/vcfutils.pl @@ -86,7 +86,7 @@ sub fillac { print; } else { my @t = split; - my @c = (0); + my @c = (0, 0); my $n = 0; my $s = -1; @_ = split(":", $t[8]); @@ -215,8 +215,8 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions # } 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] @@ -230,6 +230,7 @@ Options: -Q INT minimum RMS mapping quality for SNPs [$opts{Q}] -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. @@ -291,6 +292,12 @@ 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 @@ -499,7 +506,7 @@ Options: -d INT minimum depth [$opts{d}] 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 { diff --git a/bedidx.c b/bedidx.c index 722877d..34f0f2f 100644 --- a/bedidx.c +++ b/bedidx.c @@ -119,8 +119,10 @@ void *bed_read(const char *fn) 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; + } } } } diff --git a/faidx.c b/faidx.c index 38292a1..f0798fc 100644 --- a/faidx.c +++ b/faidx.c @@ -7,7 +7,8 @@ #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) @@ -64,10 +65,11 @@ faidx_t *fai_build_core(RAZF *rz) { 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); @@ -119,11 +121,6 @@ faidx_t *fai_build_core(RAZF *rz) 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; @@ -305,8 +302,8 @@ faidx_t *fai_load(const char *fn) 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; @@ -314,31 +311,43 @@ char *fai_fetch(const faidx_t *fai, const char *str, int *len) 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; diff --git a/glf.c b/glf.c deleted file mode 100644 index 8d5346a..0000000 --- a/glf.c +++ /dev/null @@ -1,236 +0,0 @@ -#include -#include -#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 \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 diff --git a/glf.h b/glf.h deleted file mode 100644 index 12e5400..0000000 --- a/glf.h +++ /dev/null @@ -1,56 +0,0 @@ -#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 -#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 diff --git a/sam_view.c b/sam_view.c index a579614..0821a61 100644 --- a/sam_view.c +++ b/sam_view.c @@ -22,30 +22,12 @@ typedef khash_t(rg) *rghash_t; 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)) @@ -121,7 +103,7 @@ static int usage(int is_long_help); 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; @@ -129,10 +111,9 @@ int main_samview(int argc, char *argv[]) /* 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; @@ -216,7 +197,6 @@ int main_samview(int argc, char *argv[]) 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++; } @@ -349,3 +329,69 @@ int main_import(int argc, char *argv[]) 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 \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; +} diff --git a/sample.c b/sample.c index 3efc559..4e4b8a4 100644 --- a/sample.c +++ b/sample.c @@ -55,6 +55,10 @@ int bam_smpl_add(bam_sample_t *sm, const char *fn, const char *txt) 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; diff --git a/samtools.1 b/samtools.1 index c71fc87..98ce9d0 100644 --- a/samtools.1 +++ b/samtools.1 @@ -1,7 +1,9 @@ -.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 @@ -23,6 +25,12 @@ samtools pileup -vcf ref.fasta aln.sorted.bam 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 @@ -43,7 +51,7 @@ Samtools checks the current working directory for the index file and 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 @@ -137,15 +145,56 @@ viewing the same reference sequence. .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 @@ -155,6 +204,9 @@ quality (BAQ). BAQ is the Phred-scaled probability of a read base being 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 @@ -167,24 +219,57 @@ At a position, read maximally .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 @@ -198,9 +283,6 @@ is modeled as .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 . @@ -217,25 +299,6 @@ Comma dilimited list of platforms (determined by 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 @@ -485,139 +548,174 @@ Minimum Phred-scaled LOD to call a heterozygote. [40] 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] -| +.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 -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)