From: Heng Li Date: Thu, 7 Jul 2011 02:59:05 +0000 (+0000) Subject: Release samtools-0.1.17 X-Git-Url: https://git.donarmstrong.com/?p=samtools.git;a=commitdiff_plain;h=3384f6c6c1642ada4edae9204ca1202672de7d5a Release samtools-0.1.17 --- 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 34b2795..246b020 100644 --- a/bam.h +++ b/bam.h @@ -40,7 +40,7 @@ @copyright Genome Research Ltd. */ -#define BAM_VERSION "0.1.16-dev (r969:252)" +#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 f37a783..2fdee9f 100644 --- a/bam2bcf_indel.c +++ b/bam2bcf_indel.c @@ -168,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; @@ -178,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_index.c b/bam_index.c index 9b71bdb..66d8eb8 100644 --- a/bam_index.c +++ b/bam_index.c @@ -176,7 +176,7 @@ bam_index_t *bam_index_core(bamFile fp) 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) { diff --git a/bam_plcmd.c b/bam_plcmd.c index 4870805..17e5806 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -71,6 +71,9 @@ static inline void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, cons #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); @@ -145,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)); @@ -154,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) { @@ -176,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, max_depth, max_indel_depth; const bam_pileup1_t **plp; bam_mplp_t iter; bam_header_t *h = 0; @@ -209,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; @@ -223,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); } @@ -271,7 +274,10 @@ 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); + 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) @@ -295,7 +301,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) @@ -322,8 +328,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'); @@ -333,6 +341,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'); @@ -413,6 +436,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; @@ -420,7 +444,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); @@ -434,12 +458,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; @@ -474,6 +500,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"); @@ -481,11 +508,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/bamtk.c b/bamtk.c index 79b0dbc..8ba2581 100644 --- a/bamtk.c +++ b/bamtk.c @@ -26,6 +26,7 @@ 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[]); @@ -92,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 dff439d..aeae6ca 100644 --- a/bcftools/bcf.h +++ b/bcftools/bcf.h @@ -28,7 +28,7 @@ #ifndef BCF_H #define BCF_H -#define BCF_VERSION "0.1.16-dev (r963:240)" +#define BCF_VERSION "0.1.17 (r973:277)" #include #include @@ -152,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 1c35ee8..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[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; @@ -125,6 +128,12 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, 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); @@ -240,6 +249,12 @@ static void write_header(bcf_hdr_t *h) 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== 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; @@ -309,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; @@ -323,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; @@ -357,6 +389,7 @@ int bcfview(int argc, char *argv[]) fprintf(stderr, " -p FLOAT variant if P(ref|D) 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.) { @@ -456,10 +490,19 @@ 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.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 { @@ -491,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) { @@ -521,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/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/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 fd03893..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; 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)