+Beta Release 0.1.17 (6 July, 2011)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+With the maturity of `mpileup' and the lack of update in the `pileup' command,
+the `pileup' command is now formally dropped. Most of the pileup functionality,
+such as outputting mapping quality and read positions, have been added
+`mpileup'.
+
+Since this release, `bcftools view' is able to perform contrast SNP calling
+(option -T) for discovering de novo and/or somatic mutations between a pair of
+samples or in a family trio. Potential mutations are scored by a log likelihood
+ratio, which is very simple in math, but should be comparable to more
+sophisticated methods. Note that getting the score is only the very first step.
+A lot more need to be done to reduce systematical errors due to mapping and
+reference errors and structural variations.
+
+Other notable changes in samtools:
+
+ * Improved sorting order checking during indexing.
+
+ * Improved region parsing. Colons in reference sequence names are parsed
+ properly.
+
+ * Fixed an issue where mpileup does not apply BAQ for the first few reads when
+ a region is specified.
+
+ * Fixed an issue where `faidx' does not work with FASTA files with long lines.
+
+ * Bugfix: wrong SP genotype information in the BCF output.
+
+Other notable changes in bcftools:
+
+ * Output the ML esitmate of the allele count.
+
+ * Added the HWE plus F<0 filter to varFilter. For multiple samples, it
+ effectively filters false heterozygous calls around centromeres.
+
+ * For association mapping, perform both 1-degree and 2-degree test. The
+ 2-degree test is conservative but more robust to HWE violation.
+
+(0.1.17: 6 July 2011, r973:277)
+
+
+
Beta Release 0.1.16 (21 April, 2011)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@copyright Genome Research Ltd.
*/
-#define BAM_VERSION "0.1.16-dev (r969:252)"
+#define BAM_VERSION "0.1.17 (r973:277)"
#include <stdint.h>
#include <stdlib.h>
memcpy(b->gi[0].data, bc->PL, b->gi[0].len * bc->n);
if (bcr) {
uint16_t *dp = (uint16_t*)b->gi[1].data;
- uint8_t *sp = is_SP? b->gi[2].data : 0;
+ int32_t *sp = is_SP? b->gi[2].data : 0;
for (i = 0; i < bc->n; ++i) {
bcf_callret1_t *p = bcr + i;
dp[i] = p->depth < 0xffff? p->depth : 0xffff;
if (n_types == 1 || (double)n_alt / n_tot < bca->min_frac || n_alt < bca->min_support) { // then skip
free(aux); return -1;
}
+ if (n_types >= 64) {
+ free(aux);
+ if (bam_verbose >= 2)
+ fprintf(stderr, "[%s] excessive INDEL alleles at position %d. Skip the position.\n", __func__, pos + 1);
+ return -1;
+ }
types = (int*)calloc(n_types, sizeof(int));
t = 0;
types[t++] = aux[0] - MINUS_CONST;
for (t = 0; t < n_types; ++t)
if (types[t] == 0) break;
ref_type = t; // the index of the reference type (0)
- assert(n_types < 64);
}
{ // calculate left and right boundary
left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;
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 ((uint32_t)last_tid > (uint32_t)c->tid) {
#define MPLP_NO_INDEL 0x400
#define MPLP_EXT_BAQ 0x800
#define MPLP_ILLUMINA13 0x1000
+#define MPLP_IGNORE_RG 0x2000
+#define MPLP_PRINT_POS 0x4000
+#define MPLP_PRINT_MAPQ 0x8000
void *bed_read(const char *fn);
void bed_destroy(void *_h);
}
static void group_smpl(mplp_pileup_t *m, bam_sample_t *sm, kstring_t *buf,
- int n, char *const*fn, int *n_plp, const bam_pileup1_t **plp)
+ int n, char *const*fn, int *n_plp, const bam_pileup1_t **plp, int ignore_rg)
{
int i, j;
memset(m->n_plp, 0, m->n * sizeof(int));
const bam_pileup1_t *p = plp[i] + j;
uint8_t *q;
int id = -1;
- q = bam_aux_get(p->b, "RG");
+ q = ignore_rg? 0 : bam_aux_get(p->b, "RG");
if (q) id = bam_smpl_rg2smid(sm, fn[i], (char*)q+1, buf);
if (id < 0) id = bam_smpl_rg2smid(sm, fn[i], 0, buf);
if (id < 0 || id >= m->n) {
extern void *bcf_call_add_rg(void *rghash, const char *hdtext, const char *list);
extern void bcf_call_del_rghash(void *rghash);
mplp_aux_t **data;
- int i, tid, pos, *n_plp, beg0 = 0, end0 = 1u<<29, ref_len, ref_tid, max_depth, max_indel_depth;
+ int i, tid, pos, *n_plp, tid0 = -1, beg0 = 0, end0 = 1u<<29, ref_len, ref_tid, max_depth, max_indel_depth;
const bam_pileup1_t **plp;
bam_mplp_t iter;
bam_header_t *h = 0;
data[i]->conf = conf;
h_tmp = bam_header_read(data[i]->fp);
data[i]->h = i? h : h_tmp; // for i==0, "h" has not been set yet
- bam_smpl_add(sm, fn[i], h_tmp->text);
+ bam_smpl_add(sm, fn[i], (conf->flag&MPLP_IGNORE_RG)? 0 : h_tmp->text);
rghash = bcf_call_add_rg(rghash, h_tmp->text, conf->pl_list);
if (conf->reg) {
int beg, end;
fprintf(stderr, "[%s] malformatted region or wrong seqname for %d-th input.\n", __func__, i+1);
exit(1);
}
- if (i == 0) beg0 = beg, end0 = end;
+ if (i == 0) tid0 = tid, beg0 = beg, end0 = end;
data[i]->iter = bam_iter_query(idx, tid, beg, end);
bam_index_destroy(idx);
}
bca->min_frac = conf->min_frac;
bca->min_support = conf->min_support;
}
- ref_tid = -1; ref = 0;
+ if (tid0 >= 0 && conf->fai) { // region is set
+ ref = faidx_fetch_seq(conf->fai, h->target_name[tid0], 0, 0x7fffffff, &ref_len);
+ 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)
int total_depth, _ref0, ref16;
bcf1_t *b = calloc(1, sizeof(bcf1_t));
for (i = total_depth = 0; i < n; ++i) total_depth += n_plp[i];
- group_smpl(&gplp, sm, &buf, n, fn, n_plp, plp);
+ group_smpl(&gplp, sm, &buf, n, fn, n_plp, plp, conf->flag & MPLP_IGNORE_RG);
_ref0 = (ref && pos < ref_len)? ref[pos] : 'N';
ref16 = bam_nt16_table[_ref0];
for (i = 0; i < gplp.n; ++i)
for (i = 0; i < n; ++i) {
int j;
printf("\t%d\t", n_plp[i]);
- if (n_plp[i] == 0) printf("*\t*");
- else {
+ if (n_plp[i] == 0) {
+ printf("*\t*"); // FIXME: printf() is very slow...
+ if (conf->flag & MPLP_PRINT_POS) printf("\t*");
+ } else {
for (j = 0; j < n_plp[i]; ++j)
pileup_seq(plp[i] + j, pos, ref_len, ref);
putchar('\t');
if (c > 126) c = 126;
putchar(c);
}
+ if (conf->flag & MPLP_PRINT_MAPQ) {
+ putchar('\t');
+ for (j = 0; j < n_plp[i]; ++j) {
+ int c = plp[i][j].b->core.qual + 33;
+ if (c > 126) c = 126;
+ putchar(c);
+ }
+ }
+ if (conf->flag & MPLP_PRINT_POS) {
+ putchar('\t');
+ for (j = 0; j < n_plp[i]; ++j) {
+ if (j > 0) putchar(',');
+ printf("%d", plp[i][j].qpos + 1); // FIXME: printf() is very slow...
+ }
+ }
}
}
putchar('\n');
int nfiles = 0, use_orphan = 0;
mplp_conf_t mplp;
memset(&mplp, 0, sizeof(mplp_conf_t));
+ #define MPLP_PRINT_POS 0x4000
mplp.max_mq = 60;
mplp.min_baseQ = 13;
mplp.capQ_thres = 0;
mplp.openQ = 40; mplp.extQ = 20; mplp.tandemQ = 100;
mplp.min_frac = 0.002; mplp.min_support = 1;
mplp.flag = MPLP_NO_ORPHAN | MPLP_REALN;
- while ((c = getopt(argc, argv, "Agf:r:l:M:q:Q:uaRC:BDSd:L:b:P:o:e:h:Im:F:EG:6")) >= 0) {
+ while ((c = getopt(argc, argv, "Agf:r:l:M:q:Q:uaRC:BDSd:L:b:P:o:e:h:Im:F:EG:6Os")) >= 0) {
switch (c) {
case 'f':
mplp.fai = fai_load(optarg);
case 'u': mplp.flag |= MPLP_NO_COMP | MPLP_GLF; break;
case 'a': mplp.flag |= MPLP_NO_ORPHAN | MPLP_REALN; break;
case 'B': mplp.flag &= ~MPLP_REALN; break;
- case 'R': mplp.flag |= MPLP_REALN; break;
case 'D': mplp.flag |= MPLP_FMT_DP; break;
case 'S': mplp.flag |= MPLP_FMT_SP; break;
case 'I': mplp.flag |= MPLP_NO_INDEL; break;
case 'E': mplp.flag |= MPLP_EXT_BAQ; break;
case '6': mplp.flag |= MPLP_ILLUMINA13; break;
+ case 'R': mplp.flag |= MPLP_IGNORE_RG; break;
+ case 's': mplp.flag |= MPLP_PRINT_MAPQ; break;
+ case 'O': mplp.flag |= MPLP_PRINT_POS; break;
case 'C': mplp.capQ_thres = atoi(optarg); break;
case 'M': mplp.max_mq = atoi(optarg); break;
case 'q': mplp.min_mq = atoi(optarg); break;
fprintf(stderr, " -A count anomalous read pairs\n");
fprintf(stderr, " -B disable BAQ computation\n");
fprintf(stderr, " -b FILE list of input BAM files [null]\n");
+ fprintf(stderr, " -C INT parameter for adjusting mapQ; 0 to disable [0]\n");
fprintf(stderr, " -d INT max per-BAM depth to avoid excessive memory usage [%d]\n", mplp.max_depth);
fprintf(stderr, " -E extended BAQ for higher sensitivity but lower specificity\n");
fprintf(stderr, " -f FILE faidx indexed reference sequence file [null]\n");
fprintf(stderr, " -l FILE list of positions (chr pos) or regions (BED) [null]\n");
fprintf(stderr, " -M INT cap mapping quality at INT [%d]\n", mplp.max_mq);
fprintf(stderr, " -r STR region in which pileup is generated [null]\n");
+ fprintf(stderr, " -R ignore RG tags\n");
fprintf(stderr, " -q INT skip alignments with mapQ smaller than INT [%d]\n", mplp.min_mq);
fprintf(stderr, " -Q INT skip bases with baseQ/BAQ smaller than INT [%d]\n", mplp.min_baseQ);
fprintf(stderr, "\nOutput options:\n\n");
fprintf(stderr, " -D output per-sample DP in BCF (require -g/-u)\n");
fprintf(stderr, " -g generate BCF output (genotype likelihoods)\n");
+ fprintf(stderr, " -O output base positions on reads (disabled by -g/-u)\n");
+ fprintf(stderr, " -s output mapping quality (disabled by -g/-u)\n");
fprintf(stderr, " -S output per-sample strand bias P-value in BCF (require -g/-u)\n");
fprintf(stderr, " -u generate uncompress BCF output\n");
fprintf(stderr, "\nSNP/INDEL genotype likelihoods options (effective with `-g' or `-u'):\n\n");
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[]);
else if (strcmp(argv[1], "targetcut") == 0) return main_cut_target(argc-1, argv+1);
else if (strcmp(argv[1], "phase") == 0) return main_phase(argc-1, argv+1);
else if (strcmp(argv[1], "depth") == 0) return main_depth(argc-1, argv+1);
+ else if (strcmp(argv[1], "bam2fq") == 0) return main_bam2fq(argc-1, argv+1);
+ else if (strcmp(argv[1], "pileup") == 0) {
+ fprintf(stderr, "[main] The `pileup' command has been removed. Please use `mpileup' instead.\n");
+ return 1;
+ }
#if _CURSES_LIB != 0
else if (strcmp(argv[1], "tview") == 0) return bam_tview_main(argc-1, argv+1);
#endif
CC= gcc
CFLAGS= -g -Wall -O2 #-m64 #-arch ppc
DFLAGS= -D_FILE_OFFSET_BITS=64 -D_USE_KNETFILE
-LOBJS= bcf.o vcf.o bcfutils.o prob1.o em.o kfunc.o kmin.o index.o fet.o bcf2qcall.o
+LOBJS= bcf.o vcf.o bcfutils.o prob1.o em.o kfunc.o kmin.o index.o fet.o mut.o bcf2qcall.o
OMISC= ..
AOBJS= call1.o main.o $(OMISC)/kstring.o $(OMISC)/bgzf.o $(OMISC)/knetfile.o $(OMISC)/bedidx.o
PROG= bcftools
lib:libbcf.a
libbcf.a:$(LOBJS)
- $(AR) -cru $@ $(LOBJS)
+ $(AR) -csru $@ $(LOBJS)
bcftools:lib $(AOBJS)
$(CC) $(CFLAGS) -o $@ $(AOBJS) -L. $(LIBPATH) -lbcf -lm -lz
#ifndef BCF_H
#define BCF_H
-#define BCF_VERSION "0.1.16-dev (r963:240)"
+#define BCF_VERSION "0.1.17 (r973:277)"
#include <stdint.h>
#include <zlib.h>
int bcf_fix_gt(bcf1_t *b);
// update PL generated by old samtools
int bcf_fix_pl(bcf1_t *b);
+ // convert PL to GLF-like 10-likelihood GL
+ int bcf_gl10(const bcf1_t *b, uint8_t *gl);
// string hash table
void *bcf_build_refhash(bcf_hdr_t *h);
+++ /dev/null
-.TH bcftools 1 "16 March 2011" "bcftools" "Bioinformatics tools"
-.SH NAME
-.PP
-bcftools - Utilities for the Binary Call Format (BCF) and VCF.
-.SH SYNOPSIS
-.PP
-bcftools index in.bcf
-.PP
-bcftools view in.bcf chr2:100-200 > out.vcf
-.PP
-bcftools view -vc in.bcf > out.vcf 2> out.afs
-
-.SH DESCRIPTION
-.PP
-Bcftools is a toolkit for processing VCF/BCF files, calling variants and
-estimating site allele frequencies and allele frequency spectrums.
-
-.SH COMMANDS AND OPTIONS
-
-.TP 10
-.B view
-.B bcftools view
-.RB [ \-AbFGNQSucgv ]
-.RB [ \-D
-.IR seqDict ]
-.RB [ \-l
-.IR listLoci ]
-.RB [ \-s
-.IR listSample ]
-.RB [ \-i
-.IR gapSNPratio ]
-.RB [ \-t
-.IR mutRate ]
-.RB [ \-p
-.IR varThres ]
-.RB [ \-P
-.IR prior ]
-.RB [ \-1
-.IR nGroup1 ]
-.RB [ \-d
-.IR minFrac ]
-.RB [ \-U
-.IR nPerm ]
-.RB [ \-X
-.IR permThres ]
-.I in.bcf
-.RI [ region ]
-
-Convert between BCF and VCF, call variant candidates and estimate allele
-frequencies.
-
-.RS
-.TP
-.B Input/Output Options:
-.TP 10
-.B -A
-Retain all possible alternate alleles at variant sites. By default, the view
-command discards unlikely alleles.
-.TP 10
-.B -b
-Output in the BCF format. The default is VCF.
-.TP
-.BI -D \ FILE
-Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null]
-.TP
-.B -F
-Indicate PL is generated by r921 or before (ordering is different).
-.TP
-.B -G
-Suppress all individual genotype information.
-.TP
-.BI -l \ FILE
-List of sites at which information are outputted [all sites]
-.TP
-.B -N
-Skip sites where the REF field is not A/C/G/T
-.TP
-.B -Q
-Output the QCALL likelihood format
-.TP
-.BI -s \ FILE
-List of samples to use. The first column in the input gives the sample names
-and the second gives the ploidy, which can only be 1 or 2. When the 2nd column
-is absent, the sample ploidy is assumed to be 2. In the output, the ordering of
-samples will be identical to the one in
-.IR FILE .
-[null]
-.TP
-.B -S
-The input is VCF instead of BCF.
-.TP
-.B -u
-Uncompressed BCF output (force -b).
-.TP
-.B Consensus/Variant Calling Options:
-.TP 10
-.B -c
-Call variants using Bayesian inference. This option automatically invokes option
-.BR -e .
-.TP
-.BI -d \ FLOAT
-When
-.B -v
-is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0]
-.TP
-.B -e
-Perform max-likelihood inference only, including estimating the site allele frequency,
-testing Hardy-Weinberg equlibrium and testing associations with LRT.
-.TP
-.B -g
-Call per-sample genotypes at variant sites (force -c)
-.TP
-.BI -i \ FLOAT
-Ratio of INDEL-to-SNP mutation rate [0.15]
-.TP
-.BI -p \ FLOAT
-A site is considered to be a variant if P(ref|D)<FLOAT [0.5]
-.TP
-.BI -P \ STR
-Prior or initial allele frequency spectrum. If STR can be
-.IR full ,
-.IR cond2 ,
-.I flat
-or the file consisting of error output from a previous variant calling
-run.
-.TP
-.BI -t \ FLOAT
-Scaled muttion rate for variant calling [0.001]
-.TP
-.B -v
-Output variant sites only (force -c)
-.TP
-.B Contrast Calling and Association Test Options:
-.TP
-.BI -1 \ INT
-Number of group-1 samples. This option is used for dividing the samples into
-two groups for contrast SNP calling or association test.
-When this option is in use, the following VCF INFO will be outputted:
-PC2, PCHI2 and QCHI2. [0]
-.TP
-.BI -U \ INT
-Number of permutations for association test (effective only with
-.BR -1 )
-[0]
-.TP
-.BI -X \ FLOAT
-Only perform permutations for P(chi^2)<FLOAT (effective only with
-.BR -U )
-[0.01]
-.RE
-
-.TP
-.B index
-.B bcftools index
-.I in.bcf
-
-Index sorted BCF for random access.
-.RE
-
-.TP
-.B cat
-.B bcftools cat
-.I in1.bcf
-.RI [ "in2.bcf " [ ... "]]]"
-
-Concatenate BCF files. The input files are required to be sorted and
-have identical samples appearing in the same order.
-.RE
-
-.SH BCFTOOLS SPECIFIC VCF TAGS
-
-.TS
-center box;
-cb | cb | cb
-l | l | l .
-Tag Format Description
-_
-AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
-CI95 double[2] Equal-tail Bayesian credible interval of AF at the 95% level
-DP int Raw read depth (without quality filtering)
-DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
-FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise
-MQ int Root-Mean-Square mapping quality of covering reads
-PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
-PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
-PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
-QCHI2 int Phred-scaled PCHI2
-RP int # permutations yielding a smaller PCHI2
-.TE
b->n_smpl = n_smpl;
return 0;
}
+
+static int8_t nt4_table[128] = {
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 /*'-'*/, 4, 4,
+ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 3, 4, 4, 4, -1, 4, 4, 4, 4, 4, 4, 4,
+ 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 3, 4, 4, 4, -1, 4, 4, 4, 4, 4, 4, 4
+};
+
+int bcf_gl10(const bcf1_t *b, uint8_t *gl)
+{
+ int a[4], k, l, map[4], k1, j, i;
+ const bcf_ginfo_t *PL;
+ char *s;
+ if (b->ref[1] != 0 || b->n_alleles > 4) return -1; // ref is not a single base or >4 alleles
+ for (i = 0; i < b->n_gi; ++i)
+ if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
+ if (i == b->n_gi) return -1; // no PL
+ PL = b->gi + i;
+ a[0] = nt4_table[(int)b->ref[0]];
+ if (a[0] > 3 || a[0] < 0) return -1; // ref is not A/C/G/T
+ a[1] = a[2] = a[3] = -2; // -1 has a special meaning
+ if (b->alt[0] == 0) return -1; // no alternate allele
+ map[0] = map[1] = map[2] = map[3] = -2;
+ map[a[0]] = 0;
+ for (k = 0, s = b->alt, k1 = -1; k < 3 && *s; ++k, s += 2) {
+ if (s[1] != ',' && s[1] != 0) return -1; // ALT is not single base
+ a[k+1] = nt4_table[(int)*s];
+ if (a[k+1] >= 0) map[a[k+1]] = k+1;
+ else k1 = k + 1;
+ if (s[1] == 0) break; // the end of the ALT string
+ }
+ for (k = 0; k < 4; ++k)
+ if (map[k] < 0) map[k] = k1;
+ for (i = 0; i < b->n_smpl; ++i) {
+ const uint8_t *p = PL->data + i * PL->len; // the PL for the i-th individual
+ uint8_t *g = gl + 10 * i;
+ for (j = 0; j < PL->len; ++j)
+ if (p[j]) break;
+ for (k = j = 0; k < 4; ++k) {
+ for (l = k; l < 4; ++l) {
+ int t, x = map[k], y = map[l];
+ if (x > y) t = x, x = y, y = t; // make sure x is the smaller
+ g[j++] = p[y * (y+1) / 2 + x];
+ }
+ }
+ }
+ return 0;
+}
#define VC_ANNO_MAX 16384
#define VC_FIX_PL 32768
#define VC_EM 0x10000
+#define VC_PAIRCALL 0x20000
+#define VC_QCNT 0x40000
typedef struct {
int flag, prior_type, n1, n_sub, *sublist, n_perm;
+ uint32_t *trio_aux;
char *prior_file, **subsam, *fn_dict;
uint8_t *ploidy;
double theta, pref, indel_frac, min_perm_p, min_smpl_frac, min_lrt;
bcf_sync(b);
}
-static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag, double em[10])
+static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag, double em[10], int cons_llr, int64_t cons_gt)
{
kstring_t s;
int has_I16, is_var;
if (em[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);
free(b->str);
kputs("##INFO=<ID=G3,Number=3,Type=Float,Description=\"ML estimate of genotype frequencies\">\n", &str);
if (!strstr(str.s, "##INFO=<ID=HWE,"))
kputs("##INFO=<ID=HWE,Number=1,Type=Float,Description=\"Chi^2 based HWE test P-value based on G3\">\n", &str);
+ if (!strstr(str.s, "##INFO=<ID=CLR,"))
+ kputs("##INFO=<ID=CLR,Number=1,Type=Integer,Description=\"Log ratio of genotype likelihoods with and without the constraint\">\n", &str);
+ if (!strstr(str.s, "##INFO=<ID=UGT,"))
+ kputs("##INFO=<ID=UGT,Number=1,Type=String,Description=\"The most probable unconstrained genotype configuration in the trio\">\n", &str);
+ if (!strstr(str.s, "##INFO=<ID=CGT,"))
+ kputs("##INFO=<ID=CGT,Number=1,Type=String,Description=\"The most probable constrained genotype configuration in the trio\">\n", &str);
// if (!strstr(str.s, "##INFO=<ID=CI95,"))
// kputs("##INFO=<ID=CI95,Number=2,Type=Float,Description=\"Equal-tail Bayesian credible interval of the site allele frequency at the 95% level\">\n", &str);
if (!strstr(str.s, "##INFO=<ID=PV4,"))
extern int bcf_fix_gt(bcf1_t *b);
extern int bcf_anno_max(bcf1_t *b);
extern int bcf_shuffle(bcf1_t *b, int seed);
+ extern uint32_t *bcf_trio_prep(int is_x, int is_son);
+ extern int bcf_trio_call(uint32_t *prep, const bcf1_t *b, int *llr, int64_t *gt);
+ extern int bcf_pair_call(const bcf1_t *b);
+ extern int bcf_min_diff(const bcf1_t *b);
+
bcf_t *bp, *bout = 0;
bcf1_t *b, *blast;
int c, *seeds = 0;
- uint64_t n_processed = 0;
+ uint64_t n_processed = 0, qcnt[256];
viewconf_t vc;
bcf_p1aux_t *p1 = 0;
bcf_hdr_t *hin, *hout;
tid = begin = end = -1;
memset(&vc, 0, sizeof(viewconf_t));
vc.prior_type = vc.n1 = -1; vc.theta = 1e-3; vc.pref = 0.5; vc.indel_frac = -1.; vc.n_perm = 0; vc.min_perm_p = 0.01; vc.min_smpl_frac = 0; vc.min_lrt = 1;
- while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:")) >= 0) {
+ memset(qcnt, 0, 8 * 256);
+ while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:T:Y")) >= 0) {
switch (c) {
case '1': vc.n1 = atoi(optarg); break;
case 'l': vc.bed = bed_read(optarg); break;
case 'g': vc.flag |= VC_CALL_GT | VC_CALL; break;
case 'I': vc.flag |= VC_NO_INDEL; break;
case 'M': vc.flag |= VC_ANNO_MAX; break;
+ case 'Y': vc.flag |= VC_QCNT; break;
case 't': vc.theta = atof(optarg); break;
case 'p': vc.pref = atof(optarg); break;
case 'i': vc.indel_frac = atof(optarg); break;
for (tid = 0; tid < vc.n_sub; ++tid) vc.ploidy[tid] = vc.subsam[tid][strlen(vc.subsam[tid]) + 1];
tid = -1;
break;
+ case 'T':
+ if (strcmp(optarg, "trioauto") == 0) vc.trio_aux = bcf_trio_prep(0, 0);
+ else if (strcmp(optarg, "trioxd") == 0) vc.trio_aux = bcf_trio_prep(1, 0);
+ else if (strcmp(optarg, "trioxs") == 0) vc.trio_aux = bcf_trio_prep(1, 1);
+ else if (strcmp(optarg, "pair") == 0) vc.flag |= VC_PAIRCALL;
+ else {
+ fprintf(stderr, "[%s] Option '-T' can only take value trioauto, trioxd or trioxs.\n", __func__);
+ return 1;
+ }
+ break;
case 'P':
if (strcmp(optarg, "full") == 0) vc.prior_type = MC_PTYPE_FULL;
else if (strcmp(optarg, "cond2") == 0) vc.prior_type = MC_PTYPE_COND2;
fprintf(stderr, " -p FLOAT variant if P(ref|D)<FLOAT [%.3g]\n", vc.pref);
fprintf(stderr, " -P STR type of prior: full, cond2, flat [full]\n");
fprintf(stderr, " -t FLOAT scaled substitution mutation rate [%.4g]\n", vc.theta);
+ fprintf(stderr, " -T STR constrained calling; STR can be: pair, trioauto, trioxd and trioxs (see manual) [null]\n");
fprintf(stderr, " -v output potential variant sites only (force -c)\n");
fprintf(stderr, "\nContrast calling and association test options:\n\n");
fprintf(stderr, " -1 INT number of group-1 samples [0]\n");
}
}
while (vcf_read(bp, hin, b) > 0) {
- int is_indel;
+ 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.) {
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.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 {
}
pr.perm_rank = n;
}
- if (calret >= 0) update_bcf1(b, p1, &pr, vc.pref, vc.flag, em);
- } else if (vc.flag & VC_EM) update_bcf1(b, 0, 0, 0, vc.flag, em);
+ if (calret >= 0) update_bcf1(b, p1, &pr, vc.pref, vc.flag, em, cons_llr, cons_gt);
+ } else if (vc.flag & VC_EM) update_bcf1(b, 0, 0, 0, vc.flag, em, cons_llr, cons_gt);
if (vc.flag & VC_ADJLD) { // compute LD
double f[4], r2;
if ((r2 = bcf_pair_freq(blast, b, f)) >= 0) {
vcf_close(bp); vcf_close(bout);
if (vc.fn_dict) free(vc.fn_dict);
if (vc.ploidy) free(vc.ploidy);
+ if (vc.trio_aux) free(vc.trio_aux);
if (vc.n_sub) {
int i;
for (i = 0; i < vc.n_sub; ++i) free(vc.subsam[i]);
free(vc.subsam); free(vc.sublist);
}
if (vc.bed) bed_destroy(vc.bed);
+ if (vc.flag & VC_QCNT)
+ for (c = 0; c < 256; ++c)
+ fprintf(stderr, "QT\t%d\t%lld\n", c, (long long)qcnt[c]);
if (seeds) free(seeds);
if (p1) bcf_p1_destroy(p1);
return 0;
--- /dev/null
+#include <stdlib.h>
+#include <stdint.h>
+#include "bcf.h"
+
+#define MAX_GENO 359
+
+int8_t seq_bitcnt[] = { 4, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4 };
+char *seq_nt16rev = "XACMGRSVTWYHKDBN";
+
+uint32_t *bcf_trio_prep(int is_x, int is_son)
+{
+ int i, j, k, n, map[10];
+ uint32_t *ret;
+ ret = calloc(MAX_GENO, 4);
+ for (i = 0, k = 0; i < 4; ++i)
+ for (j = i; j < 4; ++j)
+ map[k++] = 1<<i|1<<j;
+ for (i = 0, n = 1; i < 10; ++i) { // father
+ if (is_x && seq_bitcnt[map[i]] != 1) continue;
+ if (is_x && is_son) {
+ for (j = 0; j < 10; ++j) // mother
+ for (k = 0; k < 10; ++k) // child
+ if (seq_bitcnt[map[k]] == 1 && (map[j]&map[k]))
+ ret[n++] = j<<16 | i<<8 | k;
+ } else {
+ for (j = 0; j < 10; ++j) // mother
+ for (k = 0; k < 10; ++k) // child
+ if ((map[i]&map[k]) && (map[j]&map[k]) && ((map[i]|map[j])&map[k]) == map[k])
+ ret[n++] = j<<16 | i<<8 | k;
+ }
+ }
+ ret[0] = n - 1;
+ return ret;
+}
+
+int bcf_trio_call(const uint32_t *prep, const bcf1_t *b, int *llr, int64_t *gt)
+{
+ int i, j, k;
+ const bcf_ginfo_t *PL;
+ uint8_t *gl10;
+ int map[10];
+ if (b->n_smpl != 3) return -1; // not a trio
+ for (i = 0; i < b->n_gi; ++i)
+ if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
+ if (i == b->n_gi) return -1; // no PL
+ gl10 = alloca(10 * b->n_smpl);
+ if (bcf_gl10(b, gl10) < 0) return -1;
+ PL = b->gi + i;
+ for (i = 0, k = 0; i < 4; ++i)
+ for (j = i; j < 4; ++j)
+ map[k++] = seq_nt16rev[1<<i|1<<j];
+ for (j = 0; j < 3; ++j) // check if ref hom is the most probable in all members
+ if (((uint8_t*)PL->data)[j * PL->len] != 0) break;
+ if (j < 3) { // we need to go through the complex procedure
+ uint8_t *g[3];
+ int minc = 1<<30, minc_j = -1, minf = 0, gtf = 0, gtc = 0;
+ g[0] = gl10;
+ g[1] = gl10 + 10;
+ g[2] = gl10 + 20;
+ for (j = 1; j <= (int)prep[0]; ++j) { // compute LK with constraint
+ int sum = g[0][prep[j]&0xff] + g[1][prep[j]>>8&0xff] + g[2][prep[j]>>16&0xff];
+ if (sum < minc) minc = sum, minc_j = j;
+ }
+ gtc |= map[prep[minc_j]&0xff]; gtc |= map[prep[minc_j]>>8&0xff]<<8; gtc |= map[prep[minc_j]>>16]<<16;
+ for (j = 0; j < 3; ++j) { // compute LK without constraint
+ int min = 1<<30, min_k = -1;
+ for (k = 0; k < 10; ++k)
+ if (g[j][k] < min) min = g[j][k], min_k = k;
+ gtf |= map[min_k]<<(j*8);
+ minf += min;
+ }
+ *llr = minc - minf; *gt = (int64_t)gtc<<32 | gtf;
+ } else *llr = 0, *gt = -1;
+ return 0;
+}
+
+int bcf_pair_call(const bcf1_t *b)
+{
+ int i, j, k;
+ const bcf_ginfo_t *PL;
+ if (b->n_smpl != 2) return -1; // not a pair
+ for (i = 0; i < b->n_gi; ++i)
+ if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
+ if (i == b->n_gi) return -1; // no PL
+ PL = b->gi + i;
+ for (j = 0; j < 2; ++j) // check if ref hom is the most probable in all members
+ if (((uint8_t*)PL->data)[j * PL->len] != 0) break;
+ if (j < 2) { // we need to go through the complex procedure
+ uint8_t *g[2];
+ int minc = 1<<30, minf = 0;
+ g[0] = PL->data;
+ g[1] = (uint8_t*)PL->data + PL->len;
+ for (j = 0; j < PL->len; ++j) // compute LK with constraint
+ minc = minc < g[0][j] + g[1][j]? minc : g[0][j] + g[1][j];
+ for (j = 0; j < 2; ++j) { // compute LK without constraint
+ int min = 1<<30;
+ for (k = 0; k < PL->len; ++k)
+ min = min < g[j][k]? min : g[j][k];
+ minf += min;
+ }
+ return minc - minf;
+ } else return 0;
+}
+
+int bcf_min_diff(const bcf1_t *b)
+{
+ int i, min = 1<<30;
+ const bcf_ginfo_t *PL;
+ for (i = 0; i < b->n_gi; ++i)
+ if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
+ if (i == b->n_gi) return -1; // no PL
+ PL = b->gi + i;
+ for (i = 0; i < b->n_smpl; ++i) {
+ int m1, m2, j;
+ const uint8_t *p = (uint8_t*)PL->data;
+ m1 = m2 = 1<<30;
+ for (j = 0; j < PL->len; ++j) {
+ if ((int)p[j] < m1) m2 = m1, m1 = p[j];
+ else if ((int)p[j] < m2) m2 = p[j];
+ }
+ min = min < m2 - m1? min : m2 - m1;
+ }
+ return min;
+}
print;
} else {
my @t = split;
- my @c = (0);
+ my @c = (0, 0);
my $n = 0;
my $s = -1;
@_ = split(":", $t[8]);
}
sub varFilter {
- my %opts = (d=>2, D=>10000000, a=>2, W=>10, Q=>10, w=>10, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, G=>0, S=>1000);
- getopts('pd:D:W:Q:w:a:1:2:3:4:G:S:', \%opts);
+ my %opts = (d=>2, D=>10000000, a=>2, W=>10, Q=>10, w=>3, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, G=>0, S=>1000, e=>1e-4);
+ getopts('pd:D:W:Q:w:a:1:2:3:4:G:S:e:', \%opts);
die(qq/
Usage: vcfutils.pl varFilter [options] <in.vcf>
-2 FLOAT min P-value for baseQ bias [$opts{2}]
-3 FLOAT min P-value for mapQ bias [$opts{3}]
-4 FLOAT min P-value for end distance bias [$opts{4}]
+ -e FLOAT min P-value for HWE (plus F<0) [$opts{e}]
-p print filtered variants
Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
$flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
&& ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
$flt = 8 if ($flt == 0 && ((/MXGQ=(\d+)/ && $1 < $opts{G}) || (/MXSP=(\d+)/ && $1 >= $opts{S})));
+ # HWE filter
+ if ($t[7] =~ /G3=([^;,]+),([^;,]+),([^;,]+).*HWE=([^;,]+)/ && $4 < $opts{e}) {
+ my $p = 2*$1 + $2;
+ my $f = ($p > 0 && $p < 1)? 1 - $2 / ($p * (1-$p)) : 0;
+ $flt = 9 if ($f < 0);
+ }
my $score = $t[5] * 100 + $dp_alt;
my $rlen = length($t[3]) - 1; # $indel_score<0 for SNPs
my ($b, $q);
$q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
if ($q < 0) {
- $_ = $1 if ($t[7] =~ /AF1=([\d\.]+)/);
+ $_ = ($t[7] =~ /AF1=([\d\.]+)/)? $1 : 0;
$b = ($_ < .5 || $alt eq '.')? $ref : $alt;
$q = -$q;
} else {
if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0])) {
beg = atoi(str->s); // begin
if (dret != '\n') {
- if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0]))
+ if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0])) {
end = atoi(str->s); // end
+ if (end < beg) end = -1;
+ }
}
}
}
#include "khash.h"
typedef struct {
- uint64_t len:32, line_len:16, line_blen:16;
+ int32_t line_len, line_blen;
+ int64_t len;
uint64_t offset;
} faidx1_t;
KHASH_MAP_INIT_STR(s, faidx1_t)
{
char c, *name;
int l_name, m_name, ret;
- int len, line_len, line_blen, state;
+ int line_len, line_blen, state;
int l1, l2;
faidx_t *idx;
uint64_t offset;
+ int64_t len;
idx = (faidx_t*)calloc(1, sizeof(faidx_t));
idx->hash = kh_init(s);
return 0;
}
++l1; len += l2;
- if (l2 >= 0x10000) {
- fprintf(stderr, "[fai_build_core] line length exceeds 65535 in sequence '%s'.\n", name);
- free(name); fai_destroy(idx);
- return 0;
- }
if (state == 1) line_len = l1, line_blen = l2, state = 0;
else if (state == 0) {
if (l1 != line_len || l2 != line_blen) state = 2;
static rghash_t g_rghash = 0;
static int g_min_mapQ = 0, g_flag_on = 0, g_flag_off = 0;
static char *g_library, *g_rg;
-static int g_sol2sanger_tbl[128];
static void *g_bed;
void *bed_read(const char *fn);
void bed_destroy(void *_h);
int bed_overlap(const void *_h, const char *chr, int beg, int end);
-static void sol2sanger(bam1_t *b)
-{
- int l;
- uint8_t *qual = bam1_qual(b);
- if (g_sol2sanger_tbl[30] == 0) {
- for (l = 0; l != 128; ++l) {
- g_sol2sanger_tbl[l] = (int)(10.0 * log(1.0 + pow(10.0, (l - 64 + 33) / 10.0)) / log(10.0) + .499);
- if (g_sol2sanger_tbl[l] >= 93) g_sol2sanger_tbl[l] = 93;
- }
- }
- for (l = 0; l < b->core.l_qseq; ++l) {
- int q = qual[l];
- if (q > 127) q = 127;
- qual[l] = g_sol2sanger_tbl[q];
- }
-}
-
static inline int __g_skip_aln(const bam_header_t *h, const bam1_t *b)
{
if (b->core.qual < g_min_mapQ || ((b->core.flag & g_flag_on) != g_flag_on) || (b->core.flag & g_flag_off))
int main_samview(int argc, char *argv[])
{
- int c, is_header = 0, is_header_only = 0, is_bamin = 1, ret = 0, compress_level = -1, is_bamout = 0, slx2sngr = 0, is_count = 0;
+ int c, is_header = 0, is_header_only = 0, is_bamin = 1, ret = 0, compress_level = -1, is_bamout = 0, is_count = 0;
int of_type = BAM_OFDEC, is_long_help = 0;
int count = 0;
samfile_t *in = 0, *out = 0;
/* parse command-line options */
strcpy(in_mode, "r"); strcpy(out_mode, "w");
- while ((c = getopt(argc, argv, "Sbct:h1Ho:q:f:F:ul:r:xX?T:CR:L:")) >= 0) {
+ while ((c = getopt(argc, argv, "Sbct:h1Ho:q:f:F:ul:r:xX?T:R:L:")) >= 0) {
switch (c) {
case 'c': is_count = 1; break;
- case 'C': slx2sngr = 1; break;
case 'S': is_bamin = 0; break;
case 'b': is_bamout = 1; break;
case 't': fn_list = strdup(optarg); is_bamin = 0; break;
int r;
while ((r = samread(in, b)) >= 0) { // read one alignment from `in'
if (!__g_skip_aln(in->header, b)) {
- if (slx2sngr) sol2sanger(b);
if (!is_count) samwrite(out, b); // write the alignment to `out'
count++;
}
free(argv2);
return ret;
}
+
+int8_t seq_comp_table[16] = { 0, 8, 4, 12, 2, 10, 9, 14, 1, 6, 5, 13, 3, 11, 7, 15 };
+
+int main_bam2fq(int argc, char *argv[])
+{
+ bamFile fp;
+ bam_header_t *h;
+ bam1_t *b;
+ int8_t *buf;
+ int max_buf;
+ if (argc == 1) {
+ fprintf(stderr, "Usage: samtools bam2fq <in.bam>\n");
+ return 1;
+ }
+ fp = strcmp(argv[1], "-")? bam_open(argv[1], "r") : bam_dopen(fileno(stdin), "r");
+ if (fp == 0) return 1;
+ h = bam_header_read(fp);
+ b = bam_init1();
+ buf = 0;
+ max_buf = 0;
+ while (bam_read1(fp, b) >= 0) {
+ int i, qlen = b->core.l_qseq;
+ uint8_t *seq;
+ putchar('@'); fputs(bam1_qname(b), stdout);
+ if ((b->core.flag & 0x40) && !(b->core.flag & 0x80)) puts("/1");
+ else if ((b->core.flag & 0x80) && !(b->core.flag & 0x40)) puts("/2");
+ else putchar('\n');
+ if (max_buf < qlen + 1) {
+ max_buf = qlen + 1;
+ kroundup32(max_buf);
+ buf = realloc(buf, max_buf);
+ }
+ buf[qlen] = 0;
+ seq = bam1_seq(b);
+ for (i = 0; i < qlen; ++i)
+ buf[i] = bam1_seqi(seq, i);
+ if (b->core.flag & 16) { // reverse complement
+ for (i = 0; i < qlen>>1; ++i) {
+ int8_t t = seq_comp_table[buf[qlen - 1 - i]];
+ buf[qlen - 1 - i] = seq_comp_table[buf[i]];
+ buf[i] = t;
+ }
+ if (qlen&1) buf[i] = seq_comp_table[buf[i]];
+ }
+ for (i = 0; i < qlen; ++i)
+ buf[i] = bam_nt16_rev_table[buf[i]];
+ puts((char*)buf);
+ puts("+");
+ seq = bam1_qual(b);
+ for (i = 0; i < qlen; ++i)
+ buf[i] = 33 + seq[i];
+ if (b->core.flag & 16) { // reverse
+ for (i = 0; i < qlen>>1; ++i) {
+ int8_t t = buf[qlen - 1 - i];
+ buf[qlen - 1 - i] = buf[i];
+ buf[i] = t;
+ }
+ }
+ puts((char*)buf);
+ }
+ free(buf);
+ bam_destroy1(b);
+ bam_header_destroy(h);
+ bam_close(fp);
+ return 0;
+}
kstring_t buf;
int n = 0;
khash_t(sm) *sm2id = (khash_t(sm)*)sm->sm2id;
+ if (txt == 0) {
+ add_pair(sm, sm2id, fn, fn);
+ return 0;
+ }
memset(&buf, 0, sizeof(kstring_t));
while ((q = strstr(p, "@RG")) != 0) {
p = q + 3;
-.TH samtools 1 "21 April 2011" "samtools-0.1.16" "Bioinformatics tools"
+.TH samtools 1 "05 July 2011" "samtools-0.1.17" "Bioinformatics tools"
.SH NAME
.PP
samtools - Utilities for the Sequence Alignment/Map (SAM) format
+
+bcftools - Utilities for the Binary Call Format (BCF) and VCF
.SH SYNOPSIS
.PP
samtools view -bt ref_list.txt -o aln.bam aln.sam.gz
samtools mpileup -C50 -gf ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam
.PP
samtools tview aln.sorted.bam ref.fasta
+.PP
+bcftools index in.bcf
+.PP
+bcftools view in.bcf chr2:100-200 > out.vcf
+.PP
+bcftools view -vc in.bcf > out.vcf 2> out.afs
.SH DESCRIPTION
.PP
will download the index upon absence. Samtools does not retrieve the
entire alignment file unless it is asked to do so.
-.SH COMMANDS AND OPTIONS
+.SH SAMTOOLS COMMANDS AND OPTIONS
.TP 10
.B view
.TP
.B mpileup
-samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]
+.B samtools mpileup
+.RB [ \-EBug ]
+.RB [ \-C
+.IR capQcoef ]
+.RB [ \-r
+.IR reg ]
+.RB [ \-f
+.IR in.fa ]
+.RB [ \-l
+.IR list ]
+.RB [ \-M
+.IR capMapQ ]
+.RB [ \-Q
+.IR minBaseQ ]
+.RB [ \-q
+.IR minMapQ ]
+.I in.bam
+.RI [ in2.bam
+.RI [ ... ]]
Generate BCF or pileup for one or multiple BAM files. Alignment records
are grouped by sample identifiers in @RG header lines. If sample
identifiers are absent, each input file is regarded as one sample.
-.B OPTIONS:
+In the pileup format (without
+.BR -u or -g ),
+each
+line represents a genomic position, consisting of chromosome name,
+coordinate, reference base, read bases, read qualities and alignment
+mapping qualities. Information on match, mismatch, indel, strand,
+mapping quality and start and end of a read are all encoded at the read
+base column. At this column, a dot stands for a match to the reference
+base on the forward strand, a comma for a match on the reverse strand,
+a '>' or '<' for a reference skip, `ACGTN' for a mismatch on the forward
+strand and `acgtn' for a mismatch on the reverse strand. A pattern
+`\\+[0-9]+[ACGTNacgtn]+' indicates there is an insertion between this
+reference position and the next reference position. The length of the
+insertion is given by the integer in the pattern, followed by the
+inserted sequence. Similarly, a pattern `-[0-9]+[ACGTNacgtn]+'
+represents a deletion from the reference. The deleted bases will be
+presented as `*' in the following lines. Also at the read base column, a
+symbol `^' marks the start of a read. The ASCII of the character
+following `^' minus 33 gives the mapping quality. A symbol `$' marks the
+end of a read segment.
+
+.B Input Options:
.RS
.TP 10
+.B -6
+Assume the quality is in the Illumina 1.3+ encoding.
.B -A
Do not skip anomalous read pairs in variant calling.
.TP
misaligned. Applying this option greatly helps to reduce false SNPs
caused by misalignments.
.TP
+.BI -b \ FILE
+List of input BAM files, one file per line [null]
+.TP
.BI -C \ INT
Coefficient for downgrading mapping quality for reads containing
excessive mismatches. Given a read with a phred-scaled probability q of
.I INT
reads per input BAM. [250]
.TP
-.B -D
-Output per-sample read depth
-.TP
-.BI -e \ INT
-Phred-scaled gap extension sequencing error probability. Reducing
-.I INT
-leads to longer indels. [20]
-.TP
.B -E
Extended BAQ computation. This option helps sensitivity especially for MNPs, but may hurt
specificity a little bit.
.TP
.BI -f \ FILE
-The reference file [null]
+The
+.BR faidx -indexed
+reference file in the FASTA format. The file can be optionally compressed by
+.BR razip .
+[null]
+.TP
+.BI -l \ FILE
+BED or position list file containing a list of regions or sites where pileup or BCF should be generated [null]
+.TP
+.BI -q \ INT
+Minimum mapping quality for an alignment to be used [0]
+.TP
+.BI -Q \ INT
+Minimum base quality for a base to be considered [13]
+.TP
+.BI -r \ STR
+Only generate pileup in region
+.I STR
+[all sites]
+.TP
+.B Output Options:
+
+.TP
+.B -D
+Output per-sample read depth
.TP
.B -g
Compute genotype likelihoods and output them in the binary call format (BCF).
.TP
+.B -S
+Output per-sample Phred-scaled strand bias P-value
+.TP
+.B -u
+Similar to
+.B -g
+except that the output is uncompressed BCF, which is preferred for piping.
+
+.TP
+.B Options for Genotype Likelihood Computation (for -g or -u):
+
+.TP
+.BI -e \ INT
+Phred-scaled gap extension sequencing error probability. Reducing
+.I INT
+leads to longer indels. [20]
+.TP
.BI -h \ INT
Coefficient for modeling homopolymer errors. Given an
.IR l -long
.B -I
Do not perform INDEL calling
.TP
-.BI -l \ FILE
-File containing a list of sites where pileup or BCF is outputted [null]
-.TP
.BI -L \ INT
Skip INDEL calling if the average per-sample depth is above
.IR INT .
from which indel candidates are obtained. It is recommended to collect
indel candidates from sequencing technologies that have low indel error
rate such as ILLUMINA. [all]
-.TP
-.BI -q \ INT
-Minimum mapping quality for an alignment to be used [0]
-.TP
-.BI -Q \ INT
-Minimum base quality for a base to be considered [13]
-.TP
-.BI -r \ STR
-Only generate pileup in region
-.I STR
-[all sites]
-.TP
-.B -S
-Output per-sample Phred-scaled strand bias P-value
-.TP
-.B -u
-Similar to
-.B -g
-except that the output is uncompressed BCF, which is preferred for piping.
.RE
.TP
Minimum base quality to be used in het calling. [13]
.RE
-.TP
-.B pileup
-samtools pileup [-2sSBicv] [-f in.ref.fasta] [-t in.ref_list] [-l
-in.site_list] [-C capMapQ] [-M maxMapQ] [-T theta] [-N nHap] [-r
-pairDiffRate] [-m mask] [-d maxIndelDepth] [-G indelPrior]
-<in.bam>|<in.sam>
+.SH BCFTOOLS COMMANDS AND OPTIONS
-Print the alignment in the pileup format. In the pileup format, each
-line represents a genomic position, consisting of chromosome name,
-coordinate, reference base, read bases, read qualities and alignment
-mapping qualities. Information on match, mismatch, indel, strand,
-mapping quality and start and end of a read are all encoded at the read
-base column. At this column, a dot stands for a match to the reference
-base on the forward strand, a comma for a match on the reverse strand,
-a '>' or '<' for a reference skip, `ACGTN' for a mismatch on the forward
-strand and `acgtn' for a mismatch on the reverse strand. A pattern
-`\\+[0-9]+[ACGTNacgtn]+' indicates there is an insertion between this
-reference position and the next reference position. The length of the
-insertion is given by the integer in the pattern, followed by the
-inserted sequence. Similarly, a pattern `-[0-9]+[ACGTNacgtn]+'
-represents a deletion from the reference. The deleted bases will be
-presented as `*' in the following lines. Also at the read base column, a
-symbol `^' marks the start of a read. The ASCII of the character
-following `^' minus 33 gives the mapping quality. A symbol `$' marks the
-end of a read segment.
-
-If option
-.B -c
-is applied, the consensus base, Phred-scaled consensus quality, SNP
-quality (i.e. the Phred-scaled probability of the consensus being
-identical to the reference) and root mean square (RMS) mapping quality
-of the reads covering the site will be inserted between the `reference
-base' and the `read bases' columns. An indel occupies an additional
-line. Each indel line consists of chromosome name, coordinate, a star,
-the genotype, consensus quality, SNP quality, RMS mapping quality, #
-covering reads, the first alllele, the second allele, # reads supporting
-the first allele, # reads supporting the second allele and # reads
-containing indels different from the top two alleles.
-
-.B NOTE:
-Since 0.1.10, the `pileup' command is deprecated by `mpileup'.
+.TP 10
+.B view
+.B bcftools view
+.RB [ \-AbFGNQSucgv ]
+.RB [ \-D
+.IR seqDict ]
+.RB [ \-l
+.IR listLoci ]
+.RB [ \-s
+.IR listSample ]
+.RB [ \-i
+.IR gapSNPratio ]
+.RB [ \-t
+.IR mutRate ]
+.RB [ \-p
+.IR varThres ]
+.RB [ \-P
+.IR prior ]
+.RB [ \-1
+.IR nGroup1 ]
+.RB [ \-d
+.IR minFrac ]
+.RB [ \-U
+.IR nPerm ]
+.RB [ \-X
+.IR permThres ]
+.RB [ \-T
+.IR trioType ]
+.I in.bcf
+.RI [ region ]
+
+Convert between BCF and VCF, call variant candidates and estimate allele
+frequencies.
-.B OPTIONS:
.RS
+.TP
+.B Input/Output Options:
.TP 10
-.B -B
-Disable the BAQ computation. See the
-.B mpileup
-command for details.
+.B -A
+Retain all possible alternate alleles at variant sites. By default, the view
+command discards unlikely alleles.
+.TP 10
+.B -b
+Output in the BCF format. The default is VCF.
.TP
-.B -c
-Call the consensus sequence. Options
-.BR -T ", " -N ", " -I " and " -r
-are only effective when
-.BR -c " or " -g
-is in use.
+.BI -D \ FILE
+Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null]
.TP
-.BI -C \ INT
-Coefficient for downgrading the mapping quality of poorly mapped
-reads. See the
-.B mpileup
-command for details. [0]
+.B -F
+Indicate PL is generated by r921 or before (ordering is different).
.TP
-.BI -d \ INT
-Use the first
-.I NUM
-reads in the pileup for indel calling for speed up. Zero for unlimited. [1024]
+.B -G
+Suppress all individual genotype information.
.TP
-.BI -f \ FILE
-The reference sequence in the FASTA format. Index file
-.I FILE.fai
-will be created if
-absent.
+.BI -l \ FILE
+List of sites at which information are outputted [all sites]
.TP
-.B -g
-Generate genotype likelihood in the binary GLFv3 format. This option
-suppresses -c, -i and -s. This option is deprecated by the
-.B mpileup
-command.
+.B -N
+Skip sites where the REF field is not A/C/G/T
.TP
-.B -i
-Only output pileup lines containing indels.
+.B -Q
+Output the QCALL likelihood format
.TP
-.BI -I \ INT
-Phred probability of an indel in sequencing/prep. [40]
+.BI -s \ FILE
+List of samples to use. The first column in the input gives the sample names
+and the second gives the ploidy, which can only be 1 or 2. When the 2nd column
+is absent, the sample ploidy is assumed to be 2. In the output, the ordering of
+samples will be identical to the one in
+.IR FILE .
+[null]
.TP
-.BI -l \ FILE
-List of sites at which pileup is output. This file is space
-delimited. The first two columns are required to be chromosome and
-1-based coordinate. Additional columns are ignored. It is
-recommended to use option
+.B -S
+The input is VCF instead of BCF.
.TP
-.BI -m \ INT
-Filter reads with flag containing bits in
-.I INT
-[1796]
+.B -u
+Uncompressed BCF output (force -b).
.TP
-.BI -M \ INT
-Cap mapping quality at INT [60]
+.B Consensus/Variant Calling Options:
+.TP 10
+.B -c
+Call variants using Bayesian inference. This option automatically invokes option
+.BR -e .
.TP
-.BI -N \ INT
-Number of haplotypes in the sample (>=2) [2]
+.BI -d \ FLOAT
+When
+.B -v
+is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0]
.TP
-.BI -r \ FLOAT
-Expected fraction of differences between a pair of haplotypes [0.001]
+.B -e
+Perform max-likelihood inference only, including estimating the site allele frequency,
+testing Hardy-Weinberg equlibrium and testing associations with LRT.
.TP
-.B -s
-Print the mapping quality as the last column. This option makes the
-output easier to parse, although this format is not space efficient.
+.B -g
+Call per-sample genotypes at variant sites (force -c)
.TP
-.B -S
-The input file is in SAM.
+.BI -i \ FLOAT
+Ratio of INDEL-to-SNP mutation rate [0.15]
.TP
-.BI -t \ FILE
-List of reference names ane sequence lengths, in the format described
-for the
-.B import
-command. If this option is present, samtools assumes the input
-.I <in.alignment>
-is in SAM format; otherwise it assumes in BAM format.
+.BI -p \ FLOAT
+A site is considered to be a variant if P(ref|D)<FLOAT [0.5]
+.TP
+.BI -P \ STR
+Prior or initial allele frequency spectrum. If STR can be
+.IR full ,
+.IR cond2 ,
+.I flat
+or the file consisting of error output from a previous variant calling
+run.
+.TP
+.BI -t \ FLOAT
+Scaled muttion rate for variant calling [0.001]
+.TP
+.BI -T \ STR
+Enable pair/trio calling. For trio calling, option
.B -s
-together with
-.B -l
-as in the default format we may not know the mapping quality.
+is usually needed to be applied to configure the trio members and their ordering.
+In the file supplied to the option
+.BR -s ,
+the first sample must be the child, the second the father and the third the mother.
+The valid values of
+.I STR
+are `pair', `trioauto', `trioxd' and `trioxs', where `pair' calls differences between two input samples, and `trioxd' (`trioxs') specifies that the input
+is from the X chromosome non-PAR regions and the child is a female (male). [null]
.TP
-.BI -T \ FLOAT
-The theta parameter (error dependency coefficient) in the maq consensus
-calling model [0.85]
+.B -v
+Output variant sites only (force -c)
+.TP
+.B Contrast Calling and Association Test Options:
+.TP
+.BI -1 \ INT
+Number of group-1 samples. This option is used for dividing the samples into
+two groups for contrast SNP calling or association test.
+When this option is in use, the following VCF INFO will be outputted:
+PC2, PCHI2 and QCHI2. [0]
+.TP
+.BI -U \ INT
+Number of permutations for association test (effective only with
+.BR -1 )
+[0]
+.TP
+.BI -X \ FLOAT
+Only perform permutations for P(chi^2)<FLOAT (effective only with
+.BR -U )
+[0.01]
.RE
+.TP
+.B index
+.B bcftools index
+.I in.bcf
+
+Index sorted BCF for random access.
+.RE
+
+.TP
+.B cat
+.B bcftools cat
+.I in1.bcf
+.RI [ "in2.bcf " [ ... "]]]"
+
+Concatenate BCF files. The input files are required to be sorted and
+have identical samples appearing in the same order.
+.RE
.SH SAM FORMAT
-SAM is TAB-delimited. Apart from the header lines, which are started
+Sequence Alignment/Map (SAM) format is TAB-delimited. Apart from the header lines, which are started
with the `@' symbol, each alignment line consists of:
.TS
n | l | l .
Col Field Description
_
-1 QNAME Query (pair) NAME
+1 QNAME Query template/pair NAME
2 FLAG bitwise FLAG
3 RNAME Reference sequence NAME
4 POS 1-based leftmost POSition/coordinate of clipped sequence
6 CIAGR extended CIGAR string
7 MRNM Mate Reference sequence NaMe (`=' if same as RNAME)
8 MPOS 1-based Mate POSistion
-9 ISIZE Inferred insert SIZE
+9 TLEN inferred Template LENgth (insert size)
10 SEQ query SEQuence on the same strand as the reference
11 QUAL query QUALity (ASCII-33 gives the Phred base quality)
-12 OPT variable OPTional fields in the format TAG:VTYPE:VALUE
+12+ OPT variable OPTional fields in the format TAG:VTYPE:VALUE
.TE
.PP
0x0400 d the read is either a PCR or an optical duplicate
.TE
+where the second column gives the string representation of the FLAG field.
+
+.SH VCF FORMAT
+
+The Variant Call Format (VCF) is a TAB-delimited format with each data line consists of the following fields:
+.TS
+center box;
+cb | cb | cb
+n | l | l .
+Col Field Description
+_
+1 CHROM CHROMosome name
+2 POS the left-most POSition of the variant
+3 ID unique variant IDentifier
+4 REF the REFerence allele
+5 ALT the ALTernate allele(s), separated by comma
+6 QUAL variant/reference QUALity
+7 FILTER FILTers applied
+8 INFO INFOrmation related to the variant, separated by semi-colon
+9 FORMAT FORMAT of the genotype fields, separated by colon (optional)
+10+ SAMPLE SAMPLE genotypes and per-sample information (optional)
+.TE
+
+.PP
+The following table gives the
+.B INFO
+tags used by samtools and bcftools.
+
+.TS
+center box;
+cb | cb | cb
+l | l | l .
+Tag Format Description
+_
+AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
+DP int Raw read depth (without quality filtering)
+DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
+FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise
+MQ int Root-Mean-Square mapping quality of covering reads
+PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
+PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
+PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
+QCHI2 int Phred-scaled PCHI2
+RP int # permutations yielding a smaller PCHI2
+CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint
+UGT string Most probable genotype configuration without the trio constraint
+CGT string Most probable configuration with the trio constraint
+.TE
+
.SH EXAMPLES
.IP o 2
Import SAM to BAM when
.IR RG:Z:454 .
.IP o 2
-Call SNPs and short indels for one diploid individual:
+Call SNPs and short INDELs for one diploid individual:
samtools mpileup -ugf ref.fa aln.bam | bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -D 100 > var.flt.vcf
samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq
+.IP o 2
+Call somatic mutations from a pair of samples:
+
+ samtools mpileup -DSuf ref.fa aln.bam | bcftools view -bvcgT pair - > var.bcf
+
+In the output INFO field,
+.I CLR
+gives the Phred-log ratio between the likelihood by treating the
+two samples independently, and the likelihood by requiring the genotype to be identical.
+This
+.I CLR
+is effectively a score measuring the confidence of somatic calls. The higher the better.
+
+.IP o 2
+Call de novo and somatic mutations from a family trio:
+
+ samtools mpileup -DSuf ref.fa aln.bam | bcftools view -bvcgT pair -s samples.txt - > var.bcf
+
+File
+.I samples.txt
+should consist of three lines specifying the member and order of samples (in the order of child-father-mother).
+Similarly,
+.I CLR
+gives the Phred-log likelihood ratio with and without the trio constraint.
+.I UGT
+shows the most likely genotype configuration without the trio constraint, and
+.I CGT
+gives the most likely genotype configuration satisfying the trio constraint.
+
.IP o 2
Phase one individual:
.IP o 2
Unaligned words used in bam_import.c, bam_endian.h, bam.c and bam_aux.c.
.IP o 2
-In merging, the input files are required to have the same number of
-reference sequences. The requirement can be relaxed. In addition,
-merging does not reconstruct the header dictionaries
-automatically. Endusers have to provide the correct header. Picard is
-better at merging.
-.IP o 2
Samtools paired-end rmdup does not work for unpaired reads (e.g. orphan
reads or ends mapped to different chromosomes). If this is a concern,
please use Picard's MarkDuplicate which correctly handles these cases,
projects / samtools.git / commitdiff