From 63af098c78ef4163356c281e87bfb467adc8824d Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 23 Aug 2010 14:58:53 +0000 Subject: [PATCH] remove VCF output in mpileup --- Makefile | 7 +- bam_mcns.c | 278 ---------------------------------------------------- bam_mcns.h | 36 ------- bam_plcmd.c | 110 +-------------------- 4 files changed, 8 insertions(+), 423 deletions(-) delete mode 100644 bam_mcns.c delete mode 100644 bam_mcns.h diff --git a/Makefile b/Makefile index c14d663..6a15c28 100644 --- a/Makefile +++ b/Makefile @@ -7,7 +7,7 @@ LOBJS= bgzf.o kstring.o bam_aux.o bam.o bam_import.o sam.o bam_index.o \ $(KNETFILE_O) bam_sort.o sam_header.o bam_reheader.o AOBJS= bam_tview.o bam_maqcns.o bam_plcmd.o sam_view.o \ bam_rmdup.o bam_rmdupse.o bam_mate.o bam_stat.o bam_color.o \ - bamtk.o kaln.o bam_mcns.o bam2bcf.o errmod.o + bamtk.o kaln.o bam2bcf.o errmod.o PROG= samtools INCLUDES= -I. SUBDIRS= . bcftools misc @@ -53,7 +53,7 @@ bam.o:bam.h razf.h bam_endian.h kstring.h sam_header.h sam.o:sam.h bam.h bam_import.o:bam.h kseq.h khash.h razf.h bam_pileup.o:bam.h razf.h ksort.h -bam_plcmd.o:bam.h faidx.h bam_maqcns.h glf.h bam_mcns.h bcftools/bcf.h bam2bcf.h +bam_plcmd.o:bam.h faidx.h bam_maqcns.h glf.h bcftools/bcf.h bam2bcf.h bam_index.o:bam.h khash.h ksort.h razf.h bam_endian.h bam_lpileup.o:bam.h ksort.h bam_tview.o:bam.h faidx.h bam_maqcns.h @@ -63,7 +63,8 @@ bam_md.o:bam.h faidx.h glf.o:glf.h sam_header.o:sam_header.h khash.h bcf.o:bcftools/bcf.h -bam2bcf.o:bam2bcf.h bcftools/bcf.h +bam2bcf.o:bam2bcf.h errmod.h bcftools/bcf.h +errmod.o:errmod.h faidx.o:faidx.h razf.h khash.h faidx_main.o:faidx.h razf.h diff --git a/bam_mcns.c b/bam_mcns.c deleted file mode 100644 index 893d2ef..0000000 --- a/bam_mcns.c +++ /dev/null @@ -1,278 +0,0 @@ -#include -#include -#include -#include "bam_mcns.h" - -#define MC_MIN_QUAL 13 -#define MC_AVG_ERR 0.007 -#define MC_MAX_SUMQ 3000 -#define MC_MAX_SUMQP 1e-300 -#define MC_MAX_EM_ITER 16 -#define MC_EM_EPS 1e-4 - -struct __mc_aux_t { - int n, M; - int ref, alt, alt2; - double *q2p, *pdg; // pdg -> P(D|g) - double *phi, *CMk; // CMk=\binom{M}{k} - double *z, *zswap; // aux for afs - double *afs, *afs1; // afs: accumulative AFS; afs1: site posterior distribution - int *qsum, *bcnt; -}; - -void mc_init_prior(mc_aux_t *ma, int type, double theta) -{ - int i; - if (type == MC_PTYPE_COND2) { - for (i = 0; i <= 2 * ma->n; ++i) - ma->phi[i] = 2. * (i + 1) / (2 * ma->n + 1) / (2 * ma->n + 2); - } else if (type == MC_PTYPE_FLAT) { - for (i = 0; i <= ma->M; ++i) - ma->phi[i] = 1. / (ma->M + 1); - } else { - double sum; - for (i = 0, sum = 0.; i < 2 * ma->n; ++i) - sum += (ma->phi[i] = theta / (2 * ma->n - i)); - ma->phi[2 * ma->n] = 1. - sum; - } -} - -mc_aux_t *mc_init(int n) // FIXME: assuming diploid -{ - mc_aux_t *ma; - int i; - ma = calloc(1, sizeof(mc_aux_t)); - ma->n = n; ma->M = 2 * n; - ma->q2p = calloc(MC_MAX_SUMQ + 1, sizeof(double)); - ma->qsum = calloc(4 * ma->n, sizeof(int)); - ma->bcnt = calloc(4 * ma->n, sizeof(int)); - ma->pdg = calloc(3 * ma->n, sizeof(double)); - ma->phi = calloc(ma->M + 1, sizeof(double)); - ma->CMk = calloc(ma->M + 1, sizeof(double)); - ma->z = calloc(2 * ma->n + 1, sizeof(double)); - ma->zswap = calloc(2 * ma->n + 1, sizeof(double)); - ma->afs = calloc(2 * ma->n + 1, sizeof(double)); - ma->afs1 = calloc(2 * ma->n + 1, sizeof(double)); - for (i = 0; i <= MC_MAX_SUMQ; ++i) - ma->q2p[i] = pow(10., -i / 10.); - for (i = 0; i <= ma->M; ++i) - ma->CMk[i] = exp(lgamma(ma->M + 1) - lgamma(i + 1) - lgamma(ma->M - i + 1)); - mc_init_prior(ma, MC_PTYPE_FULL, 1e-3); // the simplest prior - return ma; -} - -void mc_destroy(mc_aux_t *ma) -{ - if (ma) { - free(ma->qsum); free(ma->bcnt); - free(ma->q2p); free(ma->pdg); - free(ma->phi); free(ma->CMk); - free(ma->z); free(ma->zswap); - free(ma->afs); free(ma->afs1); - free(ma); - } -} - -static int sum_err(int *n, const bam_pileup1_t **plp, mc_aux_t *ma) -{ - int i, j, tot = 0; - memset(ma->qsum, 0, sizeof(int) * 4 * ma->n); - memset(ma->bcnt, 0, sizeof(int) * 4 * ma->n); - for (j = 0; j < ma->n; ++j) { - int *qsum = ma->qsum + j * 4; - int *bcnt = ma->bcnt + j * 4; - for (i = 0; i < n[j]; ++i) { - const bam_pileup1_t *p = plp[j] + i; - int q, b; - if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue; - q = bam1_qual(p->b)[p->qpos]; - if (p->b->core.qual < q) q = p->b->core.qual; - if (q < MC_MIN_QUAL) continue; // small qual - b = bam_nt16_nt4_table[(int)bam1_seqi(bam1_seq(p->b), p->qpos)]; - if (b > 3) continue; // N - qsum[b] += q; - ++bcnt[b]; - ++tot; - } - } - return tot; -} - -static void set_allele(int ref, mc_aux_t *ma) -{ - int i, j, sum[4], tmp; - sum[0] = sum[1] = sum[2] = sum[3] = 0; - for (i = 0; i < ma->n; ++i) - for (j = 0; j < 4; ++j) - sum[j] += ma->qsum[i * 4 + j]; - for (j = 0; j < 4; ++j) sum[j] = sum[j]<<2 | j; - for (i = 1; i < 4; ++i) // insertion sort - for (j = i; j > 0 && sum[j] < sum[j-1]; --j) - tmp = sum[j], sum[j] = sum[j-1], sum[j-1] = tmp; - ma->ref = sum[3]&3; ma->alt = sum[2]&3; ma->alt2 = -1; - if (ma->ref != ref) { // the best base is not ref - if (ref >= 0 && ref <= 3) { // ref is not N - if (ma->alt == ref) tmp = ma->ref, ma->ref = ma->alt, ma->alt = tmp; // then switch alt and ref - else ma->alt2 = ma->alt, ma->alt = ma->ref, ma->ref = ref; // then set ref as ref - } else ma->alt2 = ma->alt, ma->alt = ma->ref, ma->ref = sum[0]&3; // then set the weakest as ref - } -} - -static void cal_pdg(mc_aux_t *ma) -{ - int i, j; - for (j = 0; j < ma->n; ++j) { - int pi[3], *qsum, *bcnt; - double *pdg = ma->pdg + j * 3; - qsum = ma->qsum + j * 4; - bcnt = ma->bcnt + j * 4; - pi[1] = 3 * (bcnt[ma->ref] + bcnt[ma->alt]); - pi[0] = qsum[ma->ref]; - pi[2] = qsum[ma->alt]; - for (i = 0; i < 3; ++i) - pdg[i] = pi[i] > MC_MAX_SUMQ? MC_MAX_SUMQP : ma->q2p[pi[i]]; - } -} -// this calculates the naive allele frequency and Nielsen's frequency -static double mc_freq0(const mc_aux_t *ma, double *_f) -{ - int i, cnt; - double f, f_nielsen, w_sum; - *_f = -1.; - for (i = cnt = 0, f = f_nielsen = w_sum = 0.; i < ma->n; ++i) { - int *bcnt = ma->bcnt + i * 4; - int x = bcnt[ma->ref] + bcnt[ma->alt]; - if (x) { - double w, p; - ++cnt; - f += (double)bcnt[ma->ref] / x; - p = (bcnt[ma->ref] - MC_AVG_ERR * x) / (1. - 2. * MC_AVG_ERR) / x; - w = 2. * x / (1. + x); - w_sum += w; - f_nielsen += p * w; - } - } - if (cnt) { - f_nielsen /= w_sum; - if (f_nielsen < 0.) f_nielsen = 0.; - if (f_nielsen > 1.) f_nielsen = 1.; - *_f = f_nielsen; - return f / cnt; - } else return -1.; -} -// f0 is the reference allele frequency -static double mc_freq_iter(double f0, const mc_aux_t *ma) -{ - double f, f3[3]; - int i; - f3[0] = (1.-f0)*(1.-f0); f3[1] = 2.*f0*(1.-f0); f3[2] = f0*f0; - for (i = 0, f = 0.; i < ma->n; ++i) { - double *pdg; - pdg = ma->pdg + i * 3; - f += (pdg[1] * f3[1] + 2. * pdg[2] * f3[2]) - / (pdg[0] * f3[0] + pdg[1] * f3[1] + pdg[2] * f3[2]); - } - f /= ma->n * 2.; - return f; -} - -int mc_call_gt(const mc_aux_t *ma, double f0, int k) -{ - double sum, g[3]; - double max, f3[3], *pdg = ma->pdg + k * 3; - int q, i, max_i; - f3[0] = (1.-f0)*(1.-f0); f3[1] = 2.*f0*(1.-f0); f3[2] = f0*f0; - for (i = 0, sum = 0.; i < 3; ++i) - sum += (g[i] = pdg[i] * f3[i]); - for (i = 0, max = -1., max_i = 0; i < 3; ++i) { - g[i] /= sum; - if (g[i] > max) max = g[i], max_i = i; - } - max = 1. - max; - if (max < 1e-308) max = 1e-308; - q = (int)(-3.434 * log(max) + .499); - if (q > 99) q = 99; - return q<<2|max_i; -} - -static void mc_cal_z(mc_aux_t *ma) -{ - double *z[2], *tmp, *pdg; - int i, j; - z[0] = ma->z; - z[1] = ma->zswap; - pdg = ma->pdg; - z[0][0] = 1.; z[0][1] = z[0][2] = 0.; - for (j = 0; j < ma->n; ++j) { - int max = (j + 1) * 2; - double p[3]; - pdg = ma->pdg + j * 3; - p[0] = pdg[0]; p[1] = 2. * pdg[1]; p[2] = pdg[2]; - z[1][0] = p[0] * z[0][0]; - z[1][1] = p[0] * z[0][1] + p[1] * z[0][0]; - for (i = 2; i <= max; ++i) - z[1][i] = p[0] * z[0][i] + p[1] * z[0][i-1] + p[2] * z[0][i-2]; - if (j < ma->n - 1) z[1][max+1] = z[1][max+2] = 0.; -// int k; for (k = 0; k <= max; ++k) printf("%d:%.3lg ", k, z[1][k]); putchar('\n'); - tmp = z[0]; z[0] = z[1]; z[1] = tmp; - } - if (z[0] != ma->z) memcpy(ma->z, z[0], sizeof(double) * (2 * ma->n + 1)); -} - -static double mc_add_afs(mc_aux_t *ma) -{ - int k; - long double sum = 0.; - memset(ma->afs1, 0, sizeof(double) * (ma->M + 1)); - mc_cal_z(ma); - for (k = 0, sum = 0.; k <= ma->M; ++k) - sum += (long double)ma->phi[k] * ma->z[k] / ma->CMk[k]; - for (k = 0; k <= ma->M; ++k) { - ma->afs1[k] = ma->phi[k] * ma->z[k] / ma->CMk[k] / sum; - if (isnan(ma->afs1[k]) || isinf(ma->afs1[k])) return -1.; - } - for (k = 0, sum = 0.; k <= ma->M; ++k) { - ma->afs[k] += ma->afs1[k]; - sum += k * ma->afs1[k]; - } - return sum / ma->M; -} - -int mc_cal(int ref, int *n, const bam_pileup1_t **plp, mc_aux_t *ma, mc_rst_t *rst, int level) -{ - int i, tot; - memset(rst, 0, sizeof(mc_rst_t)); - rst->f_em = rst->f_exp = -1.; rst->ref = rst->alt = -1; - // precalculation - tot = sum_err(n, plp, ma); - if (tot == 0) return 0; // no good bases - set_allele(ref, ma); - cal_pdg(ma); - // set ref/major allele - rst->ref = ma->ref; rst->alt = ma->alt; rst->alt2 = ma->alt2; - // calculate naive and Nielsen's freq - rst->f_naive = mc_freq0(ma, &rst->f_nielsen); - { // calculate f_em - double flast = rst->f_naive; - for (i = 0; i < MC_MAX_EM_ITER; ++i) { - rst->f_em = mc_freq_iter(flast, ma); - if (fabs(rst->f_em - flast) < MC_EM_EPS) break; - flast = rst->f_em; - } - } - if (level >= 2) { - rst->f_exp = mc_add_afs(ma); - rst->p_ref = ma->afs1[ma->M]; - } - return tot; -} - -void mc_dump_afs(mc_aux_t *ma) -{ - int k; - fprintf(stderr, "[afs]"); - for (k = 0; k <= ma->M; ++k) - fprintf(stderr, " %d:%.3lf", k, ma->afs[ma->M - k]); - fprintf(stderr, "\n"); - memset(ma->afs, 0, sizeof(double) * (ma->M + 1)); -} diff --git a/bam_mcns.h b/bam_mcns.h deleted file mode 100644 index 6f1a310..0000000 --- a/bam_mcns.h +++ /dev/null @@ -1,36 +0,0 @@ -#ifndef BAM_MCNS_H -#define BAM_MCNS_H - -#include "bam.h" - -struct __mc_aux_t; -typedef struct __mc_aux_t mc_aux_t; - -typedef struct { - // O(n) - int ref, alt, alt2; - double f_em, f_naive, f_nielsen; - // O(n^2) - double p_ref, f_exp; -} mc_rst_t; - -#define MC_PTYPE_FULL 1 -#define MC_PTYPE_COND2 2 -#define MC_PTYPE_FLAT 3 - -#ifdef __cplusplus -extern "C" { -#endif - - mc_aux_t *mc_init(int n); - void mc_init_prior(mc_aux_t *ma, int type, double theta); - void mc_destroy(mc_aux_t *ma); - int mc_cal(int ref, int *n, const bam_pileup1_t **plp, mc_aux_t *ma, mc_rst_t *rst, int level); - int mc_call_gt(const mc_aux_t *ma, double f0, int k); - void mc_dump_afs(mc_aux_t *ma); - -#ifdef __cplusplus -} -#endif - -#endif diff --git a/bam_plcmd.c b/bam_plcmd.c index aef10d3..80799e3 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -5,7 +5,6 @@ #include "sam.h" #include "faidx.h" #include "bam_maqcns.h" -#include "bam_mcns.h" #include "bam2bcf.h" #include "khash.h" #include "glf.h" @@ -450,16 +449,11 @@ int bam_pileup(int argc, char *argv[]) * mpileup * ***********/ -#define MPLP_VCF 0x1 -#define MPLP_VAR 0x2 -#define MPLP_AFALL 0x8 #define MPLP_GLF 0x10 #define MPLP_NO_COMP 0x20 -#define MPLP_AFS_BLOCK 0x10000 - typedef struct { - int max_mq, min_mq, prior_type, flag, min_baseQ; + int max_mq, min_mq, flag, min_baseQ; double theta; char *reg, *fn_pos; faidx_t *fai; @@ -489,12 +483,9 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) const bam_pileup1_t **plp; bam_mplp_t iter; bam_header_t *h = 0; - uint64_t N = 0; char *ref; khash_t(64) *hash = 0; - mc_aux_t *ma = 0; - bcf_callaux_t *bca = 0; bcf_callret1_t *bcr = 0; bcf_call_t bc; @@ -564,33 +555,11 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) free(s.s); bcf_hdr_sync(bh); bcf_hdr_write(bp, bh); - } else if (conf->flag & MPLP_VCF) { - kstring_t s; - s.l = s.m = 0; s.s = 0; - puts("##fileformat=VCFv4.0"); - puts("##INFO="); - puts("##INFO="); - puts("##INFO="); - puts("##FILTER="); - puts("##FILTER="); - kputs("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT", &s); - for (i = 0; i < n; ++i) { - const char *p; - kputc('\t', &s); - if ((p = strstr(fn[i], ".bam")) != 0) - kputsn(fn[i], p - fn[i], &s); - else kputs(fn[i], &s); - } - puts(s.s); - free(s.s); } // mpileup if (conf->flag & MPLP_GLF) { bca = bcf_call_init(-1., conf->min_baseQ); bcr = calloc(n, sizeof(bcf_callret1_t)); - } else if (conf->flag & MPLP_VCF) { - ma = mc_init(n); - mc_init_prior(ma, conf->prior_type, conf->theta); } ref_tid = -1; ref = 0; iter = bam_mplp_init(n, mplp_func, (void**)data); @@ -618,57 +587,6 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) bcf_write(bp, bh, b); //fprintf(stderr, "%d,%d,%d\n", b->tid, b->pos, b->l_str); bcf_destroy(b); - } else if (conf->flag & MPLP_VCF) { - mc_rst_t r; - int j, _ref0, depth, rms_q, _ref0b, is_var = 0, qref = 0, level = 2, tot; - uint64_t sqr_sum; - _ref0 = _ref0b = (ref && pos < ref_len)? ref[pos] : 'N'; - _ref0 = bam_nt16_nt4_table[bam_nt16_table[_ref0]]; - tot = mc_cal(_ref0, n_plp, plp, ma, &r, level); - if (tot) { // has good bases - double q; - is_var = (r.p_ref < .5); - q = is_var? r.p_ref : 1. - r.p_ref; - if (q < 1e-308) q = 1e-308; - qref = (int)(-3.434 * log(q) + .499); - if (qref > 99) qref = 99; - } - if ((conf->flag & MPLP_VAR) && !is_var) continue; - ++N; // number of processed lines - printf("%s\t%d\t.\t%c\t", h->target_name[tid], pos + 1, _ref0b); - if (is_var) { - putchar("ACGTN"[r.alt]); - if (r.alt2 >= 0 && r.alt2 < 4) printf(",%c", "ACGT"[r.alt2]); - } else putchar('.'); - printf("\t%d\t", qref); - if (!tot) printf("Q13\t"); - else if (r.f_exp < 0.) printf("FPE\t"); - else printf(".\t"); - for (i = depth = 0, sqr_sum = 0; i < n; ++i) { - depth += n_plp[i]; - for (j = 0; j < n_plp[i]; ++j) { - int q = plp[i][j].b->core.qual; - if (q > conf->max_mq) q = conf->max_mq; - sqr_sum += q * q; - } - } - rms_q = (int)(sqrt((double)sqr_sum / depth) + .499); - printf("DP=%d;MQ=%d", depth, rms_q); - if (tot) { - printf(";AF=%.3lf", 1. - r.f_em); - if (level >= 2) printf(";AFE=%.3lf", 1-r.f_exp); - if (conf->flag & MPLP_AFALL) - printf(";AF0=%.3lf;AFN=%.3lf", 1-r.f_naive, 1-r.f_nielsen); - } - printf("\tGT:GQ:DP"); - if (tot) { - for (i = 0; i < n; ++i) { - int x = mc_call_gt(ma, r.f_exp, i); - printf("\t%c/%c:%d:%d", "10"[((x&3)==2)], "10"[((x&3)>0)], x>>2, n_plp[i]); - } - } else for (i = 0; i < n; ++i) printf("\t./.:0:0"); - putchar('\n'); - if (N % MPLP_AFS_BLOCK == 0) mc_dump_afs(ma); } else { printf("%s\t%d\t%c", h->target_name[tid], pos + 1, (ref && pos < ref_len)? ref[pos] : 'N'); for (i = 0; i < n; ++i) { @@ -691,7 +609,6 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) } } bcf_close(bp); - if (conf->flag&MPLP_VCF) mc_dump_afs(ma); if (hash) { // free the hash table khint_t k; for (k = kh_begin(hash); k < kh_end(hash); ++k) @@ -699,7 +616,6 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) kh_destroy(64, hash); } bcf_hdr_destroy(bh); bcf_call_destroy(bca); free(bc.PL); free(bcr); - mc_destroy(ma); bam_mplp_destroy(iter); bam_header_destroy(h); for (i = 0; i < n; ++i) { @@ -717,21 +633,11 @@ int bam_mpileup(int argc, char *argv[]) mplp_conf_t mplp; memset(&mplp, 0, sizeof(mplp_conf_t)); mplp.max_mq = 60; - mplp.prior_type = MC_PTYPE_FULL; mplp.theta = 1e-3; mplp.min_baseQ = 13; - while ((c = getopt(argc, argv, "gvVcFSP:f:r:l:VM:q:t:Q:u")) >= 0) { + while ((c = getopt(argc, argv, "gf:r:l:M:q:t:Q:u")) >= 0) { switch (c) { case 't': mplp.theta = atof(optarg); break; - case 'P': - if (strcmp(optarg, "full") == 0) mplp.prior_type = MC_PTYPE_FULL; - else if (strcmp(optarg, "cond2") == 0) mplp.prior_type = MC_PTYPE_COND2; - else if (strcmp(optarg, "flat") == 0) mplp.prior_type = MC_PTYPE_FLAT; - else { - fprintf(stderr, "[%s] unrecognized prior type.\n", __func__); - return 1; - } - break; case 'f': mplp.fai = fai_load(optarg); if (mplp.fai == 0) return 1; @@ -739,17 +645,12 @@ int bam_mpileup(int argc, char *argv[]) case 'r': mplp.reg = strdup(optarg); break; case 'l': mplp.fn_pos = strdup(optarg); break; case 'g': mplp.flag |= MPLP_GLF; break; - case 'V': - case 'c': mplp.flag |= MPLP_VCF; break; - case 'F': mplp.flag |= MPLP_AFALL; break; - case 'v': mplp.flag |= MPLP_VAR; break; case 'u': mplp.flag |= MPLP_NO_COMP; break; case 'M': mplp.max_mq = atoi(optarg); break; case 'q': mplp.min_mq = atoi(optarg); break; case 'Q': mplp.min_baseQ = atoi(optarg); break; } } - if (mplp.flag&MPLP_GLF) mplp.flag &= ~MPLP_VCF; if (argc == 1) { fprintf(stderr, "\n"); fprintf(stderr, "Usage: samtools mpileup [options] in1.bam [in2.bam [...]]\n\n"); @@ -757,15 +658,12 @@ int bam_mpileup(int argc, char *argv[]) fprintf(stderr, " -r STR region in which pileup is generated [null]\n"); fprintf(stderr, " -l FILE list of positions (format: chr pos) [null]\n"); fprintf(stderr, " -M INT cap mapping quality at INT [%d]\n", mplp.max_mq); + fprintf(stderr, " -Q INT min base quality [%d]\n", mplp.min_baseQ); fprintf(stderr, " -q INT filter out alignment with MQ smaller than INT [%d]\n", mplp.min_mq); fprintf(stderr, " -t FLOAT scaled mutation rate [%lg]\n", mplp.theta); - fprintf(stderr, " -P STR prior: full, flat, cond2 [full]\n"); - fprintf(stderr, " -Q INT min base quality [%d]\n", mplp.min_baseQ); - fprintf(stderr, " -c generate VCF output (consensus calling)\n"); fprintf(stderr, " -g generate GLF output\n"); - fprintf(stderr, " -v show variant sites only\n"); fprintf(stderr, "\n"); - fprintf(stderr, "Notes: Assuming error independency and diploid individuals.\n\n"); + fprintf(stderr, "Notes: Assuming diploid individuals.\n\n"); return 1; } mpileup(&mplp, argc - optind, argv + optind); -- 2.39.2