From 0f4aed2936fdfc9df7f81e94ea8f22071acf236d Mon Sep 17 00:00:00 2001 From: Petr Danecek Date: Mon, 7 Dec 2020 14:00:13 +0000 Subject: [PATCH] [NEWS] Revamp of `bcftools call -G` Sample grouping by population was not truly independent and could still be influenced by the presence of other sample groups. With this commit, -G should work as expected. Other changes are: - Optional addition of INFO/PV4 annotation with `call -a INFO/PV4` - Remove generation of useless HOB and ICB annotation; use +fill-tags -- -t HWE,ExcHet` instead - The `call -f` option was renamed to `-a` to (1) make it consistent with `mpileup` and (2) to indicate that it includes both INFO and FORMAT annotations, not just FORMAT as previously - Any sensible Number=R,Type=Integer annotation can be used with -G, such as AD or QS - Don't trim QUAL; although usefuleness of this change is questionable for true probabilistic interpretation (such high precision is unrealistic), using QUAL as a score rather than probability *is* helpful and permits more fine-grained filtering - Fix a suspected bug in `call -F` in the worst case, for certain improve readability - Note that `call -C trio` is temporarily disabled by this commit Fixes #1332 --- call.h | 34 +- ccall.c | 4 +- doc/bcftools.1 | 38 +- doc/bcftools.html | 34 +- doc/bcftools.txt | 10 +- mcall.c | 777 +++++++++++++++++------------------- test/call-G.1.out | 17 + test/call-G.2.1.out | 27 ++ test/call-G.2.out | 17 + test/call-G.2.vcf | 24 ++ test/call-G.vcf | 13 + test/call.af-fixation.1.out | 13 + test/call.af-fixation.2.out | 13 + test/call.af-fixation.3.out | 15 + test/call.af-fixation.txt | 35 ++ test/call.af-fixation.vcf | 9 + test/check.chk | 22 +- test/check_merge.chk | 20 +- test/mpileup.1.out | 24 +- test/mpileup.2.out | 24 +- test/mpileup.3.out | 24 +- test/mpileup.4.out | 24 +- test/mpileup.5.out | 24 +- test/mpileup.X.2.out | 24 +- test/mpileup.X.out | 24 +- test/mpileup.cAls.2.out | 12 +- test/mpileup.cAls.3.out | 4 +- test/mpileup.cAls.4.out | 6 +- test/mpileup.cAls.5.out | 4 +- test/mpileup.cAls.6.out | 6 +- test/mpileup.cAls.7.out | 8 +- test/mpileup.cAls.out | 20 +- test/mpileup.cals.8.out | 4 +- test/mpileup.cals.9.out | 2 - test/mpileup.hwe.1.out | 4 +- test/mpileup.hwe.1b.out | 25 ++ test/mpileup.hwe.2.out | 4 +- test/mpileup.hwe.3.out | 6 +- test/mpileup.hwe.4.out | 4 +- test/test.pl | 12 +- test/trio-dnm.1.vcf | 31 -- test/trio-dnm.2.vcf | 31 -- vcfcall.c | 90 +++-- 43 files changed, 852 insertions(+), 711 deletions(-) create mode 100644 test/call-G.1.out create mode 100644 test/call-G.2.1.out create mode 100644 test/call-G.2.out create mode 100644 test/call-G.2.vcf create mode 100644 test/call-G.vcf create mode 100644 test/call.af-fixation.1.out create mode 100644 test/call.af-fixation.2.out create mode 100644 test/call.af-fixation.3.out create mode 100644 test/call.af-fixation.txt create mode 100644 test/call.af-fixation.vcf create mode 100644 test/mpileup.hwe.1b.out delete mode 100644 test/trio-dnm.1.vcf delete mode 100644 test/trio-dnm.2.vcf diff --git a/call.h b/call.h index 50e4815ab..4ea5443c3 100644 --- a/call.h +++ b/call.h @@ -34,7 +34,7 @@ THE SOFTWARE. */ #define CALL_CONSTR_TRIO (1<<2) #define CALL_CONSTR_ALLELES (1<<3) // -// +#define CALL_FMT_PV4 (1<<5) #define CALL_FMT_GQ (1<<6) #define CALL_FMT_GP (1<<7) @@ -52,18 +52,13 @@ family_t; // For the single-sample and grouped -G calling typedef struct { + double ref_lk, max_lk, lk_sum; float *qsum; // QS(quality sum) values - int nqsum, dp; - double fa,fb,fc,fa2,fb2,fc2,fab,fac,fbc; -} -grp1_t; -typedef struct -{ - grp1_t *grp; - int ngrp; - int *smpl2grp; + int nqsum; + uint32_t *smpl, nsmpl; + uint32_t nals, als; } -grp_t; +smpl_grp_t; // For the `-C alleles -i` constrained calling typedef struct @@ -82,6 +77,7 @@ typedef struct int *pl_map, npl_map; // same as above for PLs, but reverse (new -> old) char **als; // array to hold the trimmed set of alleles to appear on output int nals; // size of the als array + int als_new, nals_new; // bitmask with final alleles and their number family_t *fams; // list of families and samples for trio calling int nfams, mfams; int ntrio[5][5]; // possible trio genotype combinations and their counts; first idx: @@ -96,18 +92,16 @@ typedef struct int32_t *ugts, *cgts; // unconstraind and constrained GTs uint32_t output_tags; char *prior_AN, *prior_AC; // reference panel AF tags (AF=AC/AN) - tgt_als_t *tgt_als; // for CALL_CONSTR_ALLELES - char *sample_groups; // for single-sample or grouped calling with -G - grp_t smpl_grp; - float *qsum; - int nqsum; + tgt_als_t *tgt_als; // for CALL_CONSTR_ALLELES + char *sample_groups; // for single-sample or grouped calling with -G + char *sample_groups_tag; // for -G [AD|QS:] + smpl_grp_t *smpl_grp; + int nsmpl_grp; // ccall only double indel_frac, min_perm_p, min_lrt; double prior_type, pref; - double ref_lk, lk_sum; int ngrp1_samples, n_perm; - int nhets, ndiploid; char *prior_file; ccall_t *cdat; @@ -149,7 +143,7 @@ void qcall_destroy(call_t *call); void call_init_pl2p(call_t *call); uint32_t *call_trio_prep(int is_x, int is_son); -void init_allele_trimming_maps(call_t *call, int als, int nals); -void mcall_trim_numberR(call_t *call, bcf1_t *rec, int nals, int nout_als, int out_als); +void init_allele_trimming_maps(call_t *call, int nals_ori, int als_out); +void mcall_trim_and_update_numberR(call_t *call, bcf1_t *rec, int nals_ori, int nals_new); #endif diff --git a/ccall.c b/ccall.c index 126ed92d0..6bf987b69 100644 --- a/ccall.c +++ b/ccall.c @@ -303,8 +303,8 @@ static int update_bcf1(call_t *call, bcf1_t *rec, const bcf_p1rst_t *pr, double // trim Number=R tags int out_als = 0; for (i=0; i -.\" Date: 2020-09-22 +.\" Date: 2020-11-25 16:08 GMT .\" Manual: \ \& .\" Source: \ \& .\" Language: English .\" -.TH "BCFTOOLS" "1" "2020\-09\-22" "\ \&" "\ \&" +.TH "BCFTOOLS" "1" "2020\-11\-25 16:08 GMT" "\ \&" "\ \&" .\" ----------------------------------------------------------------- .\" * Define some portability stuff .\" ----------------------------------------------------------------- @@ -41,7 +41,7 @@ Most commands accept VCF, bgzipped VCF and BCF with filetype detected automatica BCFtools is designed to work on a stream\&. It regards an input file "\-" as the standard input (stdin) and outputs to the standard output (stdout)\&. Several commands can thus be combined with Unix pipes\&. .SS "VERSION" .sp -This manual page was last updated \fB2020\-09\-22\fR and refers to bcftools git version \fB1\&.11\fR\&. +This manual page was last updated \fB2020\-11\-25 16:08 GMT\fR and refers to bcftools git version \fB1\&.11\-24\-g9718479+\fR\&. .SS "BCF1" .sp The BCF1 format output by versions of samtools <= 0\&.1\&.19 is \fBnot\fR compatible with this version of bcftools\&. To read BCF1 files one can use the view command from old versions of bcftools packaged with samtools versions <= 0\&.1\&.19 to convert to VCF, which can then be read by this version of bcftools\&. @@ -664,6 +664,15 @@ See also the option\&. .RE .PP +\fB\-C, \-\-columns\-file\fR \fIfile\fR +.RS 4 +Read the list of columns from a file (normally given via the +\fB\-c, \-\-columns\fR +option)\&. "\-" to skip a column of the annotation file\&. One column name per row, an additional space\- or tab\-separated field can be present to indicate the merge logic (normally given via the +\fB\-l, \-\-merge\-logic\fR +option)\&. This is useful when many annotations are added at once\&. +.RE +.PP \fB\-e, \-\-exclude\fR \fIEXPRESSION\fR .RS 4 exclude sites for which @@ -730,11 +739,11 @@ and expressions instead of discarding them .RE .PP -\fB\-l, \-\-merge\-logic\fR \fItag\fR:\*(Aqfirst\*(Aq|\fIappend\fR|\fIunique\fR|\fIsum\fR|\fIavg\fR|\fImin\fR|\fImax\fR[,\&...] +\fB\-l, \-\-merge\-logic\fR \fItag\fR:\*(Aqfirst\*(Aq|\fIappend\fR|\fIappend\-missing\fR|\fIunique\fR|\fIsum\fR|\fIavg\fR|\fImin\fR|\fImax\fR[,\&...] .RS 4 if multiple regions overlap a single record, the option defines how to treat multiple annotation values when setting \fItag\fR -in the destination file: use the first encountered value ignoring the rest (\fIfirst\fR); append allowing duplicates (\fIappend\fR); append discarding duplicate values (\fIunique\fR); sum the values (\fIsum\fR, numeric fields only); average the values (\fIavg\fR); use the minimum value (\fImin\fR) or the maximum (\fImax\fR)\&. +in the destination file: use the first encountered value ignoring the rest (\fIfirst\fR); append allowing duplicates (\fIappend\fR); append even if the appended value is missing, i\&.e\&. is a dot (\fIappend\-missing\fR); append discarding duplicate values (\fIunique\fR); sum the values (\fIsum\fR, numeric fields only); average the values (\fIavg\fR); use the minimum value (\fImin\fR) or the maximum (\fImax\fR)\&. Note that this option is intended for use with BED or TAB\-delimited annotation files only\&. Moreover, it is effective only when either \fIREF\fR @@ -787,6 +796,12 @@ see \fBCommon Options\fR .RE .PP +\fB\-\-rename\-annots\fR \fIfile\fR +.RS 4 +rename annotations according to the map in +\fIfile\fR, with "old_name new_name\en" pairs separated by whitespaces, each on a separate line\&. The old name must be prefixed with the annotation type: INFO, FORMAT, or FILTER\&. +.RE +.PP \fB\-\-rename\-chrs\fR \fIfile\fR .RS 4 rename chromosomes according to the map in @@ -976,7 +991,7 @@ output all alternate alleles present in the alignments even if they do not appea .PP \fB\-f, \-\-format\-fields\fR \fIlist\fR .RS 4 -comma\-separated list of FORMAT fields to output for each sample\&. Currently GQ and GP fields are supported\&. For convenience, the fields can be given as lower case letters\&. +comma\-separated list of FORMAT fields to output for each sample\&. Currently GQ and GP fields are supported\&. For convenience, the fields can be given as lower case letters\&. Prefixed with "^" indicates a request for tag removal of auxiliary tags useful only for calling\&. .RE .PP \fB\-F, \-\-prior\-freqs\fR \fIAN\fR,\fIAC\fR @@ -1015,11 +1030,10 @@ is a tab\-delimited text file with sample names in the first column and group na \fI\-\fR is given instead, no HWE assumption is made at all and single\-sample calling is performed\&. (Note that in low coverage data this inflates the rate of false positives\&.) The \fB\-G\fR -option requires the presence of FORMAT/AD generated at the -\fBbcftools mpileup\fR -step by providing the -\fB\-a FMT/AD\fR -option\&. +option requires the presence of per\-sample FORMAT/QS or FORMAT/AD tag generated with +\fBbcftools mpileup \-a QS\fR +(or +\fB\-a AD\fR)\&. .RE .PP \fB\-g, \-\-gvcf\fR \fIINT\fR @@ -2279,7 +2293,7 @@ Homozygous genotypes only, useful with low coverage data (requires .PP \fB\-\-n\-matches\fR \fIINT\fR .RS 4 -Print only top INT matches for each sample, 0 for unlimited\&. Use negative value to sort by HWE probability rather than the number of discordant sites\&. +Print only top INT matches for each sample, 0 for unlimited\&. Use negative value to sort by HWE probability rather than the number of discordant sites\&. Note that average score is used to determine the top matches, not absolute values\&. .RE .PP \fB\-\-no\-HWE\-prob\fR diff --git a/doc/bcftools.html b/doc/bcftools.html index fd918301a..95b873b91 100644 --- a/doc/bcftools.html +++ b/doc/bcftools.html @@ -1,6 +1,6 @@ -bcftools

Name

bcftools — utilities for variant calling and manipulating VCFs and BCFs.

Synopsis

bcftools [--version|--version-only] [--help] [COMMAND] [OPTIONS]

DESCRIPTION

BCFtools is a set of utilities that manipulate variant calls in the Variant +bcftools

Name

bcftools — utilities for variant calling and manipulating VCFs and BCFs.

Synopsis

bcftools [--version|--version-only] [--help] [COMMAND] [OPTIONS]

DESCRIPTION

BCFtools is a set of utilities that manipulate variant calls in the Variant Call Format (VCF) and its binary counterpart BCF. All commands work transparently with both VCFs and BCFs, both uncompressed and BGZF-compressed.

Most commands accept VCF, bgzipped VCF and BCF with filetype detected automatically even when streaming from a pipe. Indexed VCF and BCF @@ -10,7 +10,7 @@ (Note that files with non-standard index names can be accessed as e.g. "bcftools view -r X:2928329 file.vcf.gz##idx##non-standard-index-name".)

BCFtools is designed to work on a stream. It regards an input file "-" as the standard input (stdin) and outputs to the standard output (stdout). Several -commands can thus be combined with Unix pipes.

VERSION

This manual page was last updated 2020-09-22 and refers to bcftools git version 1.11.

BCF1

The BCF1 format output by versions of samtools <= 0.1.19 is not +commands can thus be combined with Unix pipes.

VERSION

This manual page was last updated 2020-11-25 16:08 GMT and refers to bcftools git version 1.11-24-g9718479+.

BCF1

The BCF1 format output by versions of samtools <= 0.1.19 is not compatible with this version of bcftools. To read BCF1 files one can use the view command from old versions of bcftools packaged with samtools versions <= 0.1.19 to convert to VCF, which can then be read by @@ -286,6 +286,14 @@ See also the -l, --merge-logic option.

+-C, --columns-file file +
+ Read the list of columns from a file (normally given via the -c, --columns option). + "-" to skip a column of the annotation file. + One column name per row, an additional space- or tab-separated field can + be present to indicate the merge logic (normally given via the -l, --merge-logic option). + This is useful when many annotations are added at once. +
-e, --exclude EXPRESSION
exclude sites for which EXPRESSION is true. For valid expressions see @@ -318,11 +326,12 @@
keep sites which do not pass -i and -e expressions instead of discarding them
--l, --merge-logic tag:'first'|append|unique|sum|avg|min|max[,…] +-l, --merge-logic tag:'first'|append|append-missing|unique|sum|avg|min|max[,…]
if multiple regions overlap a single record, the option defines how to treat multiple annotation values when setting tag in the destination file: use the first encountered value ignoring - the rest (first); append allowing duplicates (append); append discarding duplicate values (unique); + the rest (first); append allowing duplicates (append); append even if the appended value is missing, + i.e. is a dot (append-missing); append discarding duplicate values (unique); sum the values (sum, numeric fields only); average the values (avg); use the minimum value (min) or the maximum (max). @@ -356,6 +365,13 @@
see Common Options
+--rename-annots file +
+ rename annotations according to the map in file, with + "old_name new_name\n" pairs separated by whitespaces, each on a separate + line. The old name must be prefixed with the annotation type: + INFO, FORMAT, or FILTER. +
--rename-chrs file
rename chromosomes according to the map in file, with @@ -490,7 +506,8 @@
comma-separated list of FORMAT fields to output for each sample. Currently GQ and GP fields are supported. For convenience, the fields can be given - as lower case letters. + as lower case letters. Prefixed with "^" indicates a request for tag + removal of auxiliary tags useful only for calling.
-F, --prior-freqs AN,AC
@@ -510,14 +527,14 @@ # Now before calling, stream the raw mpileup output through `bcftools annotate` to add the frequencies bcftools mpileup [...] -Ou | bcftools annotate -a AFs.tab.gz -h AFs.hdr -c CHROM,POS,REF,ALT,REF_AN,REF_AC -Ou | bcftools call -mv -F REF_AN,REF_AC [...]
--G, --group-samples FILE|- +-G, --group-samples FILE|-
by default, all samples are assumed to come from a single population. This option allows to group samples into populations and apply the HWE assumption within but not across the populations. FILE is a tab-delimited text file with sample names in the first column and group names in the second column. If - is given instead, no HWE assumption is made at all and single-sample calling is performed. (Note that in low coverage data this inflates the rate of false positives.) The -G option requires the presence of - FORMAT/AD generated at the bcftools mpileup step by providing the -a FMT/AD option. + per-sample FORMAT/QS or FORMAT/AD tag generated with bcftools mpileup -a QS (or -a AD).
-g, --gvcf INT
@@ -1361,7 +1378,8 @@ --n-matches INT
Print only top INT matches for each sample, 0 for unlimited. Use negative value - to sort by HWE probability rather than the number of discordant sites. + to sort by HWE probability rather than the number of discordant sites. Note + that average score is used to determine the top matches, not absolute values.
--no-HWE-prob
diff --git a/doc/bcftools.txt b/doc/bcftools.txt index a0a0d3ac2..5ece9eb8f 100644 --- a/doc/bcftools.txt +++ b/doc/bcftools.txt @@ -552,7 +552,8 @@ demand. The original calling model can be invoked with the *-c* option. *-f, --format-fields* 'list':: comma-separated list of FORMAT fields to output for each sample. Currently GQ and GP fields are supported. For convenience, the fields can be given - as lower case letters. + as lower case letters. Prefixed with "^" indicates a request for tag + removal of auxiliary tags useful only for calling. *-F, --prior-freqs* 'AN','AC':: take advantage of prior knowledge of population allele frequencies. The @@ -574,13 +575,13 @@ demand. The original calling model can be invoked with the *-c* option. bcftools mpileup [...] -Ou | bcftools annotate -a AFs.tab.gz -h AFs.hdr -c CHROM,POS,REF,ALT,REF_AN,REF_AC -Ou | bcftools call -mv -F REF_AN,REF_AC [...] ---- -*-G, --group-samples* 'FILE'|'-':: +*-G, --group-samples* [TAG:]'FILE'|'-':: by default, all samples are assumed to come from a single population. This option allows to group samples into populations and apply the HWE assumption within but not across the populations. 'FILE' is a tab-delimited text file with sample names in the first column and group names in the second column. If '-' is given instead, no HWE assumption is made at all and single-sample calling is performed. (Note that in low coverage data this inflates the rate of false positives.) The *-G* option requires the presence of - FORMAT/AD generated at the *bcftools mpileup* step by providing the *-a FMT/AD* option. + per-sample FORMAT/QS or FORMAT/AD tag generated with *bcftools mpileup -a QS* (or *-a AD*). *-g, --gvcf* 'INT':: output also gVCF blocks of homozygous REF calls. The parameter 'INT' is the @@ -1398,7 +1399,8 @@ Without the *-g* option, multi-sample cross-check of samples in 'query.vcf.gz' i *--n-matches* 'INT':: Print only top INT matches for each sample, 0 for unlimited. Use negative value - to sort by HWE probability rather than the number of discordant sites. + to sort by HWE probability rather than the number of discordant sites. Note + that average score is used to determine the top matches, not absolute values. *--no-HWE-prob*:: Disable calculation of HWE probability to reduce memory requirements with diff --git a/mcall.c b/mcall.c index c24a5a751..e02157a85 100644 --- a/mcall.c +++ b/mcall.c @@ -1,6 +1,6 @@ /* mcall.c -- multiallelic and rare variant calling. - Copyright (C) 2012-2016 Genome Research Ltd. + Copyright (C) 2012-2020 Genome Research Ltd. Author: Petr Danecek @@ -25,9 +25,11 @@ THE SOFTWARE. */ #include #include #include +#include #include #include #include "call.h" +#include "prob1.h" // Using priors for GTs does not seem to be mathematically justified. Although // it seems effective in removing false calls, it also flips a significant @@ -39,6 +41,7 @@ THE SOFTWARE. */ // genotypes is reported instead. #define FLAT_PDG_FOR_MISSING 0 +int test16(float *anno16, anno16_t *a); void qcall_init(call_t *call) { return; } void qcall_destroy(call_t *call) { return; } @@ -250,60 +253,104 @@ static void init_sample_groups(call_t *call) if ( !call->sample_groups ) { // standard pooled calling, all samples in the same group - grp_t *grps = &call->smpl_grp; - grps->ngrp = 1; - grps->grp = (grp1_t*)calloc(grps->ngrp, sizeof(grp1_t)); - grps->smpl2grp = (int*)calloc(nsmpl,sizeof(int)); + call->nsmpl_grp = 1; + call->smpl_grp = (smpl_grp_t*)calloc(1,sizeof(*call->smpl_grp)); + call->smpl_grp[0].nsmpl = nsmpl; + call->smpl_grp[0].smpl = (uint32_t*)calloc(call->smpl_grp[0].nsmpl,sizeof(uint32_t)); + for (i=0; ismpl_grp[0].smpl[i] = i; + return; + } + + // Parse tag (optional) and file name + char *fname = call->sample_groups; + while ( *fname && *fname!=':' ) fname++; + if ( *fname ) + { + call->sample_groups_tag = call->sample_groups; + *fname = 0; + fname++; + + // Is the tag defined in the header? + int tag_id = bcf_hdr_id2int(call->hdr,BCF_DT_ID,call->sample_groups_tag); + if ( tag_id==-1 ) error("No such tag \"%s\"\n",call->sample_groups_tag); + if ( !bcf_hdr_idinfo_exists(call->hdr,BCF_HL_FMT,tag_id) ) error("No such FORMAT tag \"%s\"\n", call->sample_groups_tag); + } + else + { + int tag_id = bcf_hdr_id2int(call->hdr,BCF_DT_ID,"QS"); + if ( tag_id >= 0 && bcf_hdr_idinfo_exists(call->hdr,BCF_HL_FMT,tag_id) ) call->sample_groups_tag = "QS"; + else + { + tag_id = bcf_hdr_id2int(call->hdr,BCF_DT_ID,"AD"); + if ( tag_id >= 0 && bcf_hdr_idinfo_exists(call->hdr,BCF_HL_FMT,tag_id) ) call->sample_groups_tag = "AD"; + else error("Error: neither \"AD\" nor \"QS\" FORMAT tag exists and no alternative given with -G\n"); + } + fname = call->sample_groups; } - else if ( !strcmp("-",call->sample_groups) ) + + // Read samples/groups + if ( !strcmp("-",fname) ) { // single-sample calling, each sample creates its own group - grp_t *grps = &call->smpl_grp; - grps->ngrp = nsmpl; - grps->grp = (grp1_t*)calloc(grps->ngrp, sizeof(grp1_t)); - grps->smpl2grp = (int*)malloc(nsmpl*sizeof(int)); - for (i=0; ismpl2grp[i] = i; + call->nsmpl_grp = nsmpl; + call->smpl_grp = (smpl_grp_t*)calloc(nsmpl,sizeof(*call->smpl_grp)); + for (i=0; ismpl_grp[i].nsmpl = 1; + call->smpl_grp[i].smpl = (uint32_t*)calloc(call->smpl_grp[i].nsmpl,sizeof(uint32_t)); + call->smpl_grp[i].smpl[0] = i; + } } else { int nlines; - char **lines = hts_readlist(call->sample_groups, 1, &nlines); - if ( !lines ) error("Could not read the file: %s\n", call->sample_groups); + char **lines = hts_readlist(fname, 1, &nlines); + if ( !lines ) error("Could not read the file: %s\n", fname); - uint32_t *smpl2grp1 = (uint32_t*)calloc(nsmpl,sizeof(uint32_t)); + uint32_t *smpl2grp = (uint32_t*)calloc(nsmpl,sizeof(uint32_t)); + uint32_t *grp2n = (uint32_t*)calloc(nsmpl,sizeof(uint32_t)); void *grp2idx = khash_str2int_init(); - grp_t *grps = &call->smpl_grp; + call->nsmpl_grp = 0; for (i=0; isample_groups,lines[i]); - *ptr = 0; + while ( *ptr && !isspace(*ptr) ) ptr++; + if ( !*ptr ) error("Could not parse the line in %s, expected a sample name followed by tab and a population name: %s\n",fname,lines[i]); + char *tmp = ptr; + while ( *ptr && isspace(*ptr) ) ptr++; + if ( !*ptr ) error("Could not parse the line in %s, expected a sample name followed by tab and a population name: %s\n",fname,lines[i]); + *tmp = 0; int ismpl = bcf_hdr_id2int(call->hdr, BCF_DT_SAMPLE, lines[i]); if ( ismpl<0 ) continue; - if ( smpl2grp1[ismpl] ) error("Error: the sample \"%s\" is listed twice in %s\n", lines[i],call->sample_groups); + if ( smpl2grp[ismpl] ) error("Error: the sample \"%s\" is listed twice in %s\n", lines[i],fname); if ( !khash_str2int_has_key(grp2idx,ptr+1) ) { - khash_str2int_inc(grp2idx, ptr+1); - grps->ngrp++; + khash_str2int_set(grp2idx, ptr+1, call->nsmpl_grp); + call->nsmpl_grp++; } - int igrp; - if ( khash_str2int_get(grp2idx, ptr+1, &igrp)==0 ) - smpl2grp1[ismpl] = igrp+1; - else + int igrp = -1; + if ( khash_str2int_get(grp2idx, ptr+1, &igrp)!=0 ) error("This should not happen, fixme: %s\n",ptr+1); + grp2n[igrp]++; + smpl2grp[ismpl] = igrp+1; // +1 to distinguish unlisted samples } khash_str2int_destroy(grp2idx); + if ( !call->nsmpl_grp ) error("Could not parse the file, no matching samples found: %s\n", fname); - grps->grp = (grp1_t*)calloc(grps->ngrp, sizeof(grp1_t)); - grps->smpl2grp = (int*)malloc(nsmpl*sizeof(int)); + call->smpl_grp = (smpl_grp_t*)calloc(call->nsmpl_grp,sizeof(*call->smpl_grp)); for (i=0; ihdr->samples[i],call->sample_groups); - grps->smpl2grp[i] = smpl2grp1[i] - 1; + if ( !smpl2grp[i] ) error("Error: The sample \"%s\" is not listed in %s\n",call->hdr->samples[i],fname); + int igrp = smpl2grp[i] - 1; + if ( !call->smpl_grp[igrp].nsmpl ) + call->smpl_grp[igrp].smpl = (uint32_t*)calloc(grp2n[igrp],sizeof(uint32_t)); + call->smpl_grp[igrp].smpl[call->smpl_grp[igrp].nsmpl] = i; + call->smpl_grp[igrp].nsmpl++; } - free(smpl2grp1); + free(smpl2grp); + free(grp2n); for (i=0; ismpl_grp; - for (i=0; ingrp; i++) - free(grps->grp[i].qsum); - free(grps->grp); - free(grps->smpl2grp); + for (i=0; insmpl_grp; i++) + { + free(call->smpl_grp[i].qsum); + free(call->smpl_grp[i].smpl); + } + free(call->smpl_grp); } void mcall_init(call_t *call) { + init_sample_groups(call); call_init_pl2p(call); call->nals_map = 5; @@ -342,15 +391,15 @@ void mcall_init(call_t *call) if ( call->output_tags & CALL_FMT_GQ ) bcf_hdr_append(call->hdr,"##FORMAT="); if ( call->output_tags & CALL_FMT_GP ) - bcf_hdr_append(call->hdr,"##FORMAT="); + bcf_hdr_append(call->hdr,"##FORMAT="); if ( call->output_tags & (CALL_FMT_GQ|CALL_FMT_GP) ) call->GQs = (int32_t*) malloc(sizeof(int32_t)*bcf_hdr_nsamples(call->hdr)); - bcf_hdr_append(call->hdr,"##INFO="); - bcf_hdr_append(call->hdr,"##INFO="); bcf_hdr_append(call->hdr,"##INFO="); bcf_hdr_append(call->hdr,"##INFO="); bcf_hdr_append(call->hdr,"##INFO="); bcf_hdr_append(call->hdr,"##INFO="); + if ( call->output_tags & CALL_FMT_PV4 ) + bcf_hdr_append(call->hdr,"##INFO=\n"); // init the prior if ( call->theta>0 ) @@ -373,8 +422,6 @@ void mcall_init(call_t *call) } call->theta = log(call->theta); } - - init_sample_groups(call); } void mcall_destroy(call_t *call) @@ -395,7 +442,6 @@ void mcall_destroy(call_t *call) free(call->pdg); free(call->als); free(call->ac); - free(call->qsum); return; } @@ -506,14 +552,14 @@ void set_pdg(double *pl2p, int *PLs, double *pdg, int n_smpl, int n_gt, int unse } // Create mapping between old and new (trimmed) alleles -void init_allele_trimming_maps(call_t *call, int als, int nals) +void init_allele_trimming_maps(call_t *call, int nals_ori, int als_out) { - int i, j; + int i, j, nout = 0; // als_map: old(i) -> new(j) - for (i=0, j=0; ials_map[i] = j++; + if ( als_out & (1<als_map[i] = nout++; else call->als_map[i] = -1; } @@ -521,85 +567,16 @@ void init_allele_trimming_maps(call_t *call, int als, int nals) // pl_map: new(k) -> old(l) int k = 0, l = 0; - for (i=0; ipl_map[k++] = l; + if ( (als_out & (1<pl_map[k++] = l; l++; } } } -double binom_dist(int N, double p, int k) -{ - int mean = (int) (N*p); - if ( mean==k ) return 1.0; - - double log_p = (k-mean)*log(p) + (mean-k)*log(1.0-p); - if ( k > N - k ) k = N - k; - if ( mean > N - mean ) mean = N - mean; - - if ( k < mean ) { int tmp = k; k = mean; mean = tmp; } - double diff = k - mean; - - double val = 1.0; - int i; - for (i=0; i10 && (1-q)*ndiploid>10 ) || ndiploid>200 ) - { - //fprintf(stderr,"out: mean=%e p=%e\n", mean,exp(-0.5*(nhets-mean)*(nhets-mean)/(mean*(1-q)))); - return exp(-0.5*(nhets-mean)*(nhets-mean)/(mean*(1-q))); - } - - return binom_dist(ndiploid, q, nhets); -} - -float calc_HOB(int nref, int nalt, int nhets, int ndiploid) -{ - if ( !nref || !nalt || !ndiploid ) return HUGE_VAL; - - double fref = (double)nref/(nref+nalt); // fraction of reference allelels - double falt = (double)nalt/(nref+nalt); // non-ref als - return fabs((double)nhets/ndiploid - 2*fref*falt); -} - -/** - * log(sum_i exp(a_i)) - */ -// static inline double logsumexp(double *vals, int nvals) -// { -// int i; -// double max_exp = vals[0]; -// for (i=1; ihdr); + int nsmpl = grp->nsmpl; int ngts = nals*(nals+1)/2; // Single allele @@ -635,60 +611,45 @@ static int mcall_find_best_alleles(call_t *call, int nals, int *out_als) double lk_tot = 0; int lk_tot_set = 0; int iaa = (ia+1)*(ia+2)/2-1; // index in PL which corresponds to the homozygous "ia/ia" genotype - int isample; - double *pdg = call->pdg + iaa; - for (isample=0; isamplepdg + grp->smpl[ismpl]*ngts + iaa; if ( *pdg ) { lk_tot += log(*pdg); lk_tot_set = 1; } - pdg += ngts; } if ( ia==0 ) ref_lk = lk_tot; // likelihood of 0/0 for all samples else lk_tot += call->theta; // the prior UPDATE_MAX_LKs(1<0 && lk_tot_set); } - grp_t *grps = &call->smpl_grp; - // Two alleles if ( nals>1 ) { for (ia=0; iangrp==1 && grps->grp[0].qsum[ia]==0 ) continue; + if ( grp->qsum[ia]==0 ) continue; int iaa = (ia+1)*(ia+2)/2-1; for (ib=0; ibngrp==1 && grps->grp[0].qsum[ib]==0 ) continue; + if ( grp->qsum[ib]==0 ) continue; double lk_tot = 0; int lk_tot_set = 0; - int ia_cov = 0, ib_cov = 0; - for (j=0; jngrp; j++) + double fa = grp->qsum[ia]/(grp->qsum[ia] + grp->qsum[ib]); + double fb = grp->qsum[ib]/(grp->qsum[ia] + grp->qsum[ib]); + double fa2 = fa*fa; + double fb2 = fb*fb; + double fab = 2*fa*fb; + int is, ibb = (ib+1)*(ib+2)/2-1, iab = iaa - ia + ib; + for (is=0; isgrp[j]; - if ( grp->qsum[ia] ) ia_cov = 1; - if ( grp->qsum[ib] ) ib_cov = 1; - if ( !grp->qsum[ia] && !grp->qsum[ib] ) { grp->dp = 0; continue; } - grp->dp = 1; - grp->fa = grp->qsum[ia]/(grp->qsum[ia]+grp->qsum[ib]); - grp->fb = grp->qsum[ib]/(grp->qsum[ia]+grp->qsum[ib]); - grp->fa2 = grp->fa*grp->fa; - grp->fb2 = grp->fb*grp->fb; - grp->fab = 2*grp->fa*grp->fb; - } - if ( !ia_cov || !ib_cov ) continue; - int isample, ibb = (ib+1)*(ib+2)/2-1, iab = iaa - ia + ib; - double *pdg = call->pdg; - for (isample=0; isamplegrp[grps->smpl2grp[isample]]; - if ( !grp->dp ) continue; + int ismpl = grp->smpl[is]; + double *pdg = call->pdg + ismpl*ngts; double val = 0; - if ( !call->ploidy || call->ploidy[isample]==2 ) - val = grp->fa2*pdg[iaa] + grp->fb2*pdg[ibb] + grp->fab*pdg[iab]; - else if ( call->ploidy && call->ploidy[isample]==1 ) - val = grp->fa*pdg[iaa] + grp->fb*pdg[ibb]; + if ( !call->ploidy || call->ploidy[ismpl]==2 ) + val = fa2*pdg[iaa] + fb2*pdg[ibb] + fab*pdg[iab]; + else if ( call->ploidy && call->ploidy[ismpl]==1 ) + val = fa*pdg[iaa] + fb*pdg[ibb]; if ( val ) { lk_tot += log(val); lk_tot_set = 1; } - pdg += ngts; } if ( ia!=0 ) lk_tot += call->theta; // the prior if ( ib!=0 ) lk_tot += call->theta; @@ -702,50 +663,38 @@ static int mcall_find_best_alleles(call_t *call, int nals, int *out_als) { for (ia=0; iangrp==1 && grps->grp[0].qsum[ia]==0 ) continue; + if ( grp->qsum[ia]==0 ) continue; int iaa = (ia+1)*(ia+2)/2-1; for (ib=0; ibngrp==1 && grps->grp[0].qsum[ib]==0 ) continue; + if ( grp->qsum[ib]==0 ) continue; int ibb = (ib+1)*(ib+2)/2-1; int iab = iaa - ia + ib; for (ic=0; icngrp==1 && grps->grp[0].qsum[ic]==0 ) continue; + if ( grp->qsum[ic]==0 ) continue; double lk_tot = 0; - int lk_tot_set = 1; - int ia_cov = 0, ib_cov = 0, ic_cov = 0; - for (j=0; jngrp; j++) - { - grp1_t *grp = &grps->grp[j]; - if ( grp->qsum[ia] ) ia_cov = 1; - if ( grp->qsum[ib] ) ib_cov = 1; - if ( grp->qsum[ic] ) ic_cov = 1; - if ( !grp->qsum[ia] && !grp->qsum[ib] && !grp->qsum[ic] ) { grp->dp = 0; continue; } - grp->dp = 1; - grp->fa = grp->qsum[ia]/(grp->qsum[ia]+grp->qsum[ib]+grp->qsum[ic]); - grp->fb = grp->qsum[ib]/(grp->qsum[ia]+grp->qsum[ib]+grp->qsum[ic]); - grp->fc = grp->qsum[ic]/(grp->qsum[ia]+grp->qsum[ib]+grp->qsum[ic]); - grp->fa2 = grp->fa*grp->fa; - grp->fb2 = grp->fb*grp->fb; - grp->fc2 = grp->fc*grp->fc; - grp->fab = 2*grp->fa*grp->fb, grp->fac = 2*grp->fa*grp->fc, grp->fbc = 2*grp->fb*grp->fc; - } - if ( !ia_cov || !ib_cov || !ic_cov ) continue; - int isample, icc = (ic+1)*(ic+2)/2-1; + int lk_tot_set = 0; + + double fa = grp->qsum[ia]/(grp->qsum[ia] + grp->qsum[ib] + grp->qsum[ic]); + double fb = grp->qsum[ib]/(grp->qsum[ia] + grp->qsum[ib] + grp->qsum[ic]); + double fc = grp->qsum[ic]/(grp->qsum[ia] + grp->qsum[ib] + grp->qsum[ic]); + double fa2 = fa*fa; + double fb2 = fb*fb; + double fc2 = fc*fc; + double fab = 2*fa*fb, fac = 2*fa*fc, fbc = 2*fb*fc; + int is, icc = (ic+1)*(ic+2)/2-1; int iac = iaa - ia + ic, ibc = ibb - ib + ic; - double *pdg = call->pdg; - for (isample=0; isamplegrp[grps->smpl2grp[isample]]; - if ( !grp->dp ) continue; + int ismpl = grp->smpl[is]; + double *pdg = call->pdg + ismpl*ngts; double val = 0; - if ( !call->ploidy || call->ploidy[isample]==2 ) - val = grp->fa2*pdg[iaa] + grp->fb2*pdg[ibb] + grp->fc2*pdg[icc] + grp->fab*pdg[iab] + grp->fac*pdg[iac] + grp->fbc*pdg[ibc]; - else if ( call->ploidy && call->ploidy[isample]==1 ) - val = grp->fa*pdg[iaa] + grp->fb*pdg[ibb] + grp->fc*pdg[icc]; + if ( !call->ploidy || call->ploidy[ismpl]==2 ) + val = fa2*pdg[iaa] + fb2*pdg[ibb] + fc2*pdg[icc] + fab*pdg[iab] + fac*pdg[iac] + fbc*pdg[ibc]; + else if ( call->ploidy && call->ploidy[ismpl]==1 ) + val = fa*pdg[iaa] + fb*pdg[ibb] + fc*pdg[icc]; if ( val ) { lk_tot += log(val); lk_tot_set = 1; } - pdg += ngts; } if ( ia!=0 ) lk_tot += call->theta; // the prior if ( ib!=0 ) lk_tot += call->theta; // the prior @@ -756,25 +705,26 @@ static int mcall_find_best_alleles(call_t *call, int nals, int *out_als) } } - call->ref_lk = ref_lk; - call->lk_sum = lk_sum; - *out_als = max_als; - int i, n = 0; for (i=0; imax_lk = max_lk; + grp->ref_lk = ref_lk; + grp->lk_sum = lk_sum; + grp->als = max_als; + grp->nals = n; + return n; } -static void mcall_set_ref_genotypes(call_t *call, int nals) +// Sets GT=0/0 or GT=. if PL=0,0,0 +static void mcall_set_ref_genotypes(call_t *call, int nals_ori) { int i; - int ngts = nals*(nals+1)/2; + int ngts = nals_ori*(nals_ori+1)/2; // need this to distinguish between GT=0/0 vs GT=. int nsmpl = bcf_hdr_nsamples(call->hdr); - for (i=0; iac[i] = 0; - call->nhets = 0; - call->ndiploid = 0; + for (i=0; iac[i] = 0; // nals_new<=nals_ori, never mind setting extra 0's // Set all genotypes to 0/0 or 0 int *gts = call->gts; @@ -800,34 +750,27 @@ static void mcall_set_ref_genotypes(call_t *call, int nals) } } -static void mcall_call_genotypes(call_t *call, bcf1_t *rec, int nals, int nout_als, int out_als) +static void mcall_call_genotypes(call_t *call, int nals_ori, smpl_grp_t *grp) { int ia, ib, i; - int ngts = nals*(nals+1)/2; - int nsmpl = bcf_hdr_nsamples(call->hdr); - int nout_gts = nout_als*(nout_als+1)/2; - hts_expand(float,nout_gts*nsmpl,call->nGPs,call->GPs); - - for (i=0; iac[i] = 0; - call->nhets = 0; - call->ndiploid = 0; + int ngts_ori = nals_ori*(nals_ori+1)/2; + int ngts_new = call->nals_new*(call->nals_new+1)/2; + int nsmpl = grp->nsmpl; #if USE_PRIOR_FOR_GTS float prior = exp(call->theta); #endif - float *gps = call->GPs - nout_gts; - double *pdg = call->pdg - ngts; - int *gts = call->gts - 2; - int isample; - for (isample = 0; isample < nsmpl; isample++) + int is; + for (is = 0; is < nsmpl; is++) { - int ploidy = call->ploidy ? call->ploidy[isample] : 2; - assert( ploidy>=0 && ploidy<=2 ); + int ismpl = grp->smpl[is]; + double *pdg = call->pdg + ismpl*ngts_ori; + float *gps = call->GPs + ismpl*ngts_new; + int *gts = call->gts + ismpl*2; - pdg += ngts; - gts += 2; - gps += nout_gts; + int ploidy = call->ploidy ? call->ploidy[ismpl] : 2; + assert( ploidy>=0 && ploidy<=2 ); if ( !ploidy ) { @@ -839,8 +782,8 @@ static void mcall_call_genotypes(call_t *call, bcf1_t *rec, int nals, int nout_a #if !FLAT_PDG_FOR_MISSING // Skip samples with zero depth, they have all pdg's equal to 0 - for (i=0; indiploid++; - // Default fallback for the case all LKs are the same gts[0] = bcf_gt_unphased(0); gts[1] = ploidy==2 ? bcf_gt_unphased(0) : bcf_int32_vector_end; // Non-zero depth, determine the most likely genotype - grp1_t *grp = &call->smpl_grp.grp[call->smpl_grp.smpl2grp[isample]]; double best_lk = 0; - for (ia=0; iaals & 1<qsum[ia]*grp->qsum[ia] : pdg[iaa]*grp->qsum[ia]; #if USE_PRIOR_FOR_GTS if ( ia!=0 ) lk *= prior; @@ -877,13 +817,13 @@ static void mcall_call_genotypes(call_t *call, bcf1_t *rec, int nals, int nout_a if ( ploidy==2 ) { gts[1] = gts[0]; - for (ia=0; iaals & 1<als & 1<qsum[ia]*grp->qsum[ib]; #if USE_PRIOR_FOR_GTS @@ -900,7 +840,6 @@ static void mcall_call_genotypes(call_t *call, bcf1_t *rec, int nals, int nout_a } } } - if ( gts[0] != gts[1] ) call->nhets++; } else gts[1] = bcf_int32_vector_end; @@ -908,55 +847,50 @@ static void mcall_call_genotypes(call_t *call, bcf1_t *rec, int nals, int nout_a call->ac[ bcf_gt_allele(gts[0]) ]++; if ( gts[1]!=bcf_int32_vector_end ) call->ac[ bcf_gt_allele(gts[1]) ]++; } - if ( call->output_tags & (CALL_FMT_GQ|CALL_FMT_GP) ) + if ( !(call->output_tags & (CALL_FMT_GQ|CALL_FMT_GP)) ) return; + double max, sum; + for (is=0; isGPs + isample*nout_gts; + int ismpl = grp->smpl[is]; + float *gps = call->GPs + ismpl*ngts_new; - int nmax; - if ( call->ploidy ) - { - if ( call->ploidy[isample]==2 ) nmax = nout_gts; - else if ( call->ploidy[isample]==1 ) nmax = nout_als; - else nmax = 0; - } - else nmax = nout_gts; + int nmax; + if ( call->ploidy ) + { + if ( call->ploidy[ismpl]==2 ) nmax = ngts_new; + else if ( call->ploidy[ismpl]==1 ) nmax = grp->nals; + else nmax = 0; + } + else nmax = ngts_new; - max = gps[0]; - if ( max<0 || nmax==0 ) - { - // no call - if ( call->output_tags & CALL_FMT_GP ) - { - for (i=0; iGQs[isample] = 0; - continue; - } - sum = gps[0]; - for (i=1; iGQs[isample] = max<=INT8_MAX ? max : INT8_MAX; + max = gps[0]; + if ( max<0 || nmax==0 ) + { + // no call if ( call->output_tags & CALL_FMT_GP ) { - assert( max ); - for (i=0; iGQs[ismpl] = 0; + continue; + } + sum = gps[0]; + for (i=1; iGQs[ismpl] = max<=INT8_MAX ? max : INT8_MAX; + if ( call->output_tags & CALL_FMT_GP ) + { + assert( max ); + for (i=0; ioutput_tags & CALL_FMT_GP ) - bcf_update_format_float(call->hdr, rec, "GP", call->GPs, nsmpl*nout_gts); - if ( call->output_tags & CALL_FMT_GQ ) - bcf_update_format_int32(call->hdr, rec, "GQ", call->GQs, nsmpl); } @@ -979,12 +913,13 @@ static void mcall_call_genotypes(call_t *call, bcf1_t *rec, int nals, int nout_a Individual qualities are calculated as GQ(F=i,M=j,K=k) = P(F=i,M=j,K=k) / \sum_{x,y} P(F=i,M=x,K=y) */ -static void mcall_call_trio_genotypes(call_t *call, bcf1_t *rec, int nals, int nout_als, int out_als) +#if 0 +static void mcall_call_trio_genotypes(call_t *call, bcf1_t *rec, int nals, int nals_new, int als_new) { int ia, ib, i; int nsmpl = bcf_hdr_nsamples(call->hdr); int ngts = nals*(nals+1)/2; - int nout_gts = nout_als*(nout_als+1)/2; + int nout_gts = nals_new*(nals_new+1)/2; double *gls = call->GLs - nout_gts; double *pdg = call->pdg - ngts; @@ -1014,7 +949,7 @@ static void mcall_call_trio_genotypes(call_t *call, bcf1_t *rec, int nals, int n double best_lk = 0; for (ia=0; iaals_map[ia],call->als_map[ia]); double lk = ploidy==2 ? pdg[iaa]*grp->qsum[ia]*grp->qsum[ia] : pdg[iaa]*grp->qsum[ia]; @@ -1030,10 +965,10 @@ static void mcall_call_trio_genotypes(call_t *call, bcf1_t *rec, int nals, int n { for (ia=0; iaals_map[ia],call->als_map[ib]); double lk = 2*pdg[iab]*grp->qsum[ia]*grp->qsum[ib]; @@ -1077,8 +1012,8 @@ static void mcall_call_trio_genotypes(call_t *call, bcf1_t *rec, int nals, int n for (ifm=0; ifmnfams; ifm++) { family_t *fam = &call->fams[ifm]; - int ntrio = call->ntrio[fam->type][nout_als]; - uint16_t *trio = call->trio[fam->type][nout_als]; + int ntrio = call->ntrio[fam->type][nals_new]; + uint16_t *trio = call->trio[fam->type][nals_new]; // Unconstrained likelihood int uc_itr = 0; @@ -1226,11 +1161,12 @@ static void mcall_call_trio_genotypes(call_t *call, bcf1_t *rec, int nals, int n bcf_update_format_int32(call->hdr,rec,"CGT",call->cgts,nsmpl); } } +#endif -static void mcall_trim_PLs(call_t *call, bcf1_t *rec, int nals, int nout_als, int out_als) +static void mcall_trim_and_update_PLs(call_t *call, bcf1_t *rec, int nals_ori, int nals_new) { - int ngts = nals*(nals+1)/2; - int npls_src = ngts, npls_dst = nout_als*(nout_als+1)/2; // number of PL values in diploid samples, ori and new + int npls_src = nals_ori*(nals_ori+1)/2; + int npls_dst = nals_new*(nals_new+1)/2; // number of PL values in diploid samples, ori and new if ( call->all_diploid && npls_src == npls_dst ) return; int *pls_src = call->PLs, *pls_dst = call->PLs; @@ -1247,7 +1183,7 @@ static void mcall_trim_PLs(call_t *call, bcf1_t *rec, int nals, int nout_als, in } else if ( ploidy==1 ) { - for (ia=0; iapl_map[isrc] ]; @@ -1257,7 +1193,7 @@ static void mcall_trim_PLs(call_t *call, bcf1_t *rec, int nals, int nout_als, in else { pls_dst[0] = bcf_int32_missing; - pls_dst[1] = bcf_int32_vector_end; // relying on nout_als>1 in mcall() + pls_dst[1] = bcf_int32_vector_end; // relying on nals_new>1 in mcall() } pls_src += npls_src; pls_dst += npls_dst; @@ -1265,9 +1201,9 @@ static void mcall_trim_PLs(call_t *call, bcf1_t *rec, int nals, int nout_als, in bcf_update_format_int32(call->hdr, rec, "PL", call->PLs, npls_dst*nsmpl); } -void mcall_trim_numberR(call_t *call, bcf1_t *rec, int nals, int nout_als, int out_als) +void mcall_trim_and_update_numberR(call_t *call, bcf1_t *rec, int nals_ori, int nals_new) { - if ( nals==nout_als ) return; + if ( nals_ori==nals_new ) return; int i,j, nret, size = sizeof(float); @@ -1286,17 +1222,17 @@ void mcall_trim_numberR(call_t *call, bcf1_t *rec, int nals, int nout_als, int o nret = bcf_get_info_values(call->hdr, rec, key, &tmp_ori, &ntmp_ori, type); if ( nret<=0 ) continue; - if ( nout_als==1 ) + if ( nals_new==1 ) bcf_update_info_int32(call->hdr, rec, key, tmp_ori, 1); // has to be the REF, the order could not change else { - for (j=0; jals_map[j]; if ( k==-1 ) continue; // to be dropped memcpy((char *)tmp_new+size*k, (char *)tmp_ori+size*j, size); } - bcf_update_info_int32(call->hdr, rec, key, tmp_new, nout_als); + bcf_update_info_int32(call->hdr, rec, key, tmp_new, nals_new); } } @@ -1313,21 +1249,21 @@ void mcall_trim_numberR(call_t *call, bcf1_t *rec, int nals, int nout_als, int o if (nret<=0) continue; int nsmpl = bcf_hdr_nsamples(call->hdr); - assert( nret==nals*nsmpl ); + assert( nret==nals_ori*nsmpl ); for (j=0; jals_map[k]; if ( l==-1 ) continue; // to be dropped memcpy(ptr_dst+size*l, ptr_src+size*k, size); } } - bcf_update_format_int32(call->hdr, rec, key, tmp_new, nout_als*nsmpl); + bcf_update_format_int32(call->hdr, rec, key, tmp_new, nals_new*nsmpl); } call->PLs = (int32_t*) tmp_new; @@ -1442,12 +1378,12 @@ static int mcall_constrain_alleles(call_t *call, bcf1_t *rec, int *unseen) } bcf_update_format_int32(call->hdr, rec, "PL", call->itmp, npls_new*nsmpl); - // update QS - int nqs = bcf_get_info_float(call->hdr, rec, "QS", &call->smpl_grp.grp[0].qsum, &call->smpl_grp.grp[0].nqsum); - hts_expand(float,nals,call->nqsum,call->qsum); + // update QS, use temporarily call->GPs to store the values + int nqs = bcf_get_info_float(call->hdr, rec, "QS", &call->smpl_grp[0].qsum, &call->smpl_grp[0].nqsum); + hts_expand(float,nals,call->nGPs,call->GPs); for (i=0; iqsum[i] = call->als_map[i]smpl_grp.grp[0].qsum[call->als_map[i]] : 0; - bcf_update_info_float(call->hdr, rec, "QS", call->qsum, nals); + call->GPs[i] = call->als_map[i]smpl_grp[0].qsum[call->als_map[i]] : 0; + bcf_update_info_float(call->hdr, rec, "QS", call->GPs, nals); // update any Number=R tags void *tmp_ori = call->itmp, *tmp_new = call->PLs; // reusing PLs storage which is not used at this point @@ -1488,7 +1424,6 @@ static int mcall_constrain_alleles(call_t *call, bcf1_t *rec, int *unseen) call->itmp = (int32_t*) tmp_ori; call->n_itmp = ntmp_ori; - if ( *unseen ) *unseen = nals-1; return 0; } @@ -1507,203 +1442,229 @@ int mcall(call_t *call, bcf1_t *rec) // Force alleles when calling genotypes given alleles was requested if ( call->flag & CALL_CONSTR_ALLELES && mcall_constrain_alleles(call, rec, &unseen)!=0 ) return -2; - int nsmpl = bcf_hdr_nsamples(call->hdr); - int nals = rec->n_allele; - hts_expand(int,nals,call->nac,call->ac); - hts_expand(int,nals,call->nals_map,call->als_map); - hts_expand(int,nals*(nals+1)/2,call->npl_map,call->pl_map); + int nsmpl = bcf_hdr_nsamples(call->hdr); + int nals_ori = rec->n_allele; + hts_expand(int,nals_ori,call->nac,call->ac); + hts_expand(int,nals_ori,call->nals_map,call->als_map); + hts_expand(int,nals_ori*(nals_ori+1)/2,call->npl_map,call->pl_map); // Get the genotype likelihoods call->nPLs = bcf_get_format_int32(call->hdr, rec, "PL", &call->PLs, &call->mPLs); - if ( call->nPLs!=nsmpl*nals*(nals+1)/2 && call->nPLs!=nsmpl*nals ) // a mixture of diploid and haploid or haploid only - error("Wrong number of PL fields? nals=%d npl=%d\n", nals,call->nPLs); + if ( call->nPLs!=nsmpl*nals_ori*(nals_ori+1)/2 && call->nPLs!=nsmpl*nals_ori ) // a mixture of diploid and haploid or haploid only + error("Wrong number of PL fields? nals=%d npl=%d\n", nals_ori,call->nPLs); // Convert PLs to probabilities - int ngts = nals*(nals+1)/2; + int ngts_ori = nals_ori*(nals_ori+1)/2; hts_expand(double, call->nPLs, call->npdg, call->pdg); - set_pdg(call->pl2p, call->PLs, call->pdg, nsmpl, ngts, unseen); + set_pdg(call->pl2p, call->PLs, call->pdg, nsmpl, ngts_ori, unseen); // Get sum of qualities, serves as an AF estimate, f_x = QS/N in Eq. 1 in call-m math notes. - if ( call->smpl_grp.ngrp == 1 ) + if ( call->nsmpl_grp == 1 ) { - int nqs = bcf_get_info_float(call->hdr, rec, "QS", &call->smpl_grp.grp[0].qsum, &call->smpl_grp.grp[0].nqsum); + int nqs = bcf_get_info_float(call->hdr, rec, "QS", &call->smpl_grp[0].qsum, &call->smpl_grp[0].nqsum); if ( nqs<=0 ) error("The QS annotation not present at %s:%d\n", bcf_seqname(call->hdr,rec),rec->pos+1); - if ( nqs < nals ) + if ( nqs < nals_ori ) { // Some of the listed alleles do not have the corresponding QS field. This is // typically ref-only site with <*> in ALT. - hts_expand(float,nals,call->smpl_grp.grp[0].nqsum,call->smpl_grp.grp[0].qsum); - for (i=nqs; ismpl_grp.grp[0].qsum[i] = 0; + hts_expand(float,nals_ori,call->smpl_grp[0].nqsum,call->smpl_grp[0].qsum); + for (i=nqs; ismpl_grp[0].qsum[i] = 0; } } else { - for (j=0; jsmpl_grp.ngrp; j++) + for (j=0; jnsmpl_grp; j++) { - hts_expand(float,nals,call->smpl_grp.grp[j].nqsum,call->smpl_grp.grp[j].qsum); - memset(call->smpl_grp.grp[j].qsum, 0, sizeof(float)*nals); + hts_expand(float,nals_ori,call->smpl_grp[j].nqsum,call->smpl_grp[j].qsum); + memset(call->smpl_grp[j].qsum, 0, sizeof(float)*nals_ori); } - int nad = bcf_get_format_int32(call->hdr, rec, "AD", &call->ADs, &call->nADs); - if ( nad<1 ) error("Error: FORMAT/AD is required with the -G option, mpileup must be run with -a AD\n"); + // Use FORMAT/AD or FORMAT/QS + int nad = bcf_get_format_int32(call->hdr, rec, call->sample_groups_tag, &call->ADs, &call->nADs); + if ( nad<1 ) error("Error: FORMAT/%s is required with the -G option, mpileup must be run with \"-a AD\" or \"-a QS\"\n",call->sample_groups_tag); nad /= bcf_hdr_nsamples(call->hdr); - hts_expand(float,nals,call->nqsum,call->qsum); - float qsum = 0; - for (i=0; ihdr); i++) + for (i=0; insmpl_grp; i++) { - int32_t *ptr = call->ADs + i*nad; - for (j=0; jsmpl_grp[i]; + hts_expand(float,nals_ori,grp->nqsum,grp->qsum); + for (j=0; jqsum[j] = 0; + for (is=0; isnsmpl; is++) { - if ( ptr[j]==bcf_int32_vector_end ) break; - if ( ptr[j]==bcf_int32_missing ) call->qsum[j] = 0; - else { call->qsum[j] = ptr[j]; qsum += ptr[j]; } + int ismpl = grp->smpl[is]; + int32_t *ptr = call->ADs + ismpl*nad; + float sum = 0; + for (j=0; jqsum[j] += ptr[j]/sum; + } + } } - for (; jqsum[j] = 0; - if ( qsum ) - for (j=0; jqsum[j] /= qsum; - - grp1_t *grp = &call->smpl_grp.grp[call->smpl_grp.smpl2grp[i]]; - for (j=0; jqsum[j] += call->qsum[j]; } } // If available, take into account reference panel AFs if ( call->prior_AN && bcf_get_info_int32(call->hdr, rec, call->prior_AN ,&call->ac, &call->nac)==1 ) { - int an = call->ac[0]; - if ( bcf_get_info_int32(call->hdr, rec, call->prior_AC ,&call->ac, &call->nac)==nals-1 ) + int an = call->ac[0]; // number of alleles total, procede only if not zero; reuse call->ac + if ( an > 0 && bcf_get_info_int32(call->hdr, rec, call->prior_AC ,&call->ac, &call->nac)==nals_ori-1 ) // number of ALT alleles { - int ac0 = an; // number of alleles in the reference population - for (i=0; iac[i]==bcf_int32_vector_end ) break; if ( call->ac[i]==bcf_int32_missing ) continue; ac0 -= call->ac[i]; - for (j=0; jsmpl_grp.ngrp; j++) - call->smpl_grp.grp[j].qsum[i+1] += call->ac[i]*0.5; + + // here an*0.5 is the number of samples in the populatio and ac*0.5 is the AF weighted by the number of samples + for (j=0; jnsmpl_grp; j++) + call->smpl_grp[j].qsum[i+1] = (call->smpl_grp[j].qsum[i+1] + 0.5*call->ac[i]) / (call->smpl_grp[j].nsmpl + 0.5*an); } if ( ac0<0 ) error("Incorrect %s,%s values at %s:%d\n", call->prior_AN,call->prior_AC,bcf_seqname(call->hdr,rec),rec->pos+1); - for (j=0; jsmpl_grp.ngrp; j++) - call->smpl_grp.grp[j].qsum[0] += ac0*0.5; - for (i=0; ismpl_grp.ngrp; j++) - call->smpl_grp.grp[j].qsum[i] /= nsmpl + 0.5*an; - } + for (j=0; jnsmpl_grp; j++) + call->smpl_grp[j].qsum[0] = (call->smpl_grp[j].qsum[0] + 0.5*ac0) / (call->smpl_grp[j].nsmpl + 0.5*an); } } - for (j=0; jsmpl_grp.ngrp; j++) + // normalize so that QS sums to 1 for each group + for (j=0; jnsmpl_grp; j++) { - float qsum_tot = 0; - for (i=0; ismpl_grp.grp[j].qsum[i]; - if ( qsum_tot ) for (i=0; ismpl_grp.grp[j].qsum[i] /= qsum_tot; + float sum = 0; + for (i=0; ismpl_grp[j].qsum[i]; + if ( sum ) for (i=0; ismpl_grp[j].qsum[i] /= sum; } bcf_update_info_int32(call->hdr, rec, "QS", NULL, 0); // remove QS tag - // Find the best combination of alleles - int out_als, nout; - if ( nals > 8*sizeof(out_als) ) + if ( nals_ori > 8*sizeof(call->als_new) ) { fprintf(stderr,"Too many alleles at %s:%"PRId64", skipping.\n", bcf_seqname(call->hdr,rec),(int64_t) rec->pos+1); return 0; } - nout = mcall_find_best_alleles(call, nals, &out_als); - // Make sure the REF allele is always present - if ( !(out_als&1) ) + // For each group find the best combination of alleles + call->als_new = 0; + double ref_lk = -HUGE_VAL, lk_sum = -HUGE_VAL, max_qual = -HUGE_VAL; + for (j=0; jnsmpl_grp; j++) { - out_als |= 1; - nout++; + smpl_grp_t *grp = &call->smpl_grp[j]; + mcall_find_best_alleles(call, nals_ori, grp); + call->als_new |= grp->als; + if ( grp->max_lk==-HUGE_VAL ) continue; + double qual = -4.343*(grp->ref_lk - logsumexp2(grp->lk_sum,grp->ref_lk)); + if ( max_qual < qual ) + { + max_qual = qual; + lk_sum = grp->lk_sum; + ref_lk = grp->ref_lk; + } } - int is_variant = out_als==1 ? 0 : 1; + + // Make sure the REF allele is always present + if ( !(call->als_new&1) ) call->als_new |= 1; + + int is_variant = call->als_new==1 ? 0 : 1; if ( call->flag & CALL_VARONLY && !is_variant ) return 0; - // With -A, keep all ALTs except X - if ( call->flag & CALL_KEEPALT ) + call->nals_new = 0; + for (i=0; i0 && i==unseen ) continue; - out_als |= 1<0 && i==unseen ) continue; + if ( call->flag & CALL_KEEPALT ) call->als_new |= 1<als_new & (1<nals_new++; } + init_allele_trimming_maps(call,nals_ori,call->als_new); + int nAC = 0; - if ( out_als==1 ) // only REF allele on output + if ( call->als_new==1 ) // only REF allele on output { - init_allele_trimming_maps(call, 1, nals); - mcall_set_ref_genotypes(call,nals); + mcall_set_ref_genotypes(call,nals_ori); bcf_update_format_int32(call->hdr, rec, "PL", NULL, 0); // remove PL, useless now } + else if ( !is_variant ) + { + mcall_set_ref_genotypes(call,nals_ori); // running with -A, prevent mcall_call_genotypes from putting some ALT back + mcall_trim_and_update_PLs(call, rec, nals_ori, call->nals_new); + } else { // The most likely set of alleles includes non-reference allele (or was enforced), call genotypes. // Note that it is a valid outcome if the called genotypes exclude some of the ALTs. - init_allele_trimming_maps(call, out_als, nals); - if ( !is_variant ) - mcall_set_ref_genotypes(call,nals); // running with -A, prevent mcall_call_genotypes from putting some ALT back - else if ( call->flag & CALL_CONSTR_TRIO ) + int ngts_new = call->nals_new*(call->nals_new+1)/2; + hts_expand(float,ngts_new*nsmpl,call->nGPs,call->GPs); + for (i=0; inals_new; i++) call->ac[i] = 0; + + if ( call->flag & CALL_CONSTR_TRIO && call->nals_new>4 ) + { + fprintf(stderr,"Too many alleles at %s:%"PRId64", skipping.\n", bcf_seqname(call->hdr,rec),(int64_t) rec->pos+1); + return 0; + } + if ( call->output_tags & (CALL_FMT_GQ|CALL_FMT_GP) ) { - if ( nout>4 ) - { - fprintf(stderr,"Too many alleles at %s:%"PRId64", skipping.\n", bcf_seqname(call->hdr,rec),(int64_t) rec->pos+1); - return 0; - } - mcall_call_trio_genotypes(call, rec, nals,nout,out_als); + memset(call->GPs,0,nsmpl*ngts_new*sizeof(*call->GPs)); + memset(call->GQs,0,nsmpl*sizeof(*call->GQs)); + } + for (i=0; insmpl_grp; i++) + { + if ( call->flag & CALL_CONSTR_TRIO ) + error("todo: constrained trio calling temporarily disabled\n"); //mcall_call_trio_genotypes(call,rec,nals,&call->smpl_grp[i]); + else + mcall_call_genotypes(call,nals_ori,&call->smpl_grp[i]); } - else - mcall_call_genotypes(call,rec,nals,nout,out_als); // Skip the site if all samples are 0/0. This can happen occasionally. - nAC = 0; - for (i=1; iac[i]; + for (i=1; inals_new; i++) nAC += call->ac[i]; if ( !nAC && call->flag & CALL_VARONLY ) return 0; - mcall_trim_PLs(call, rec, nals, nout, out_als); + + if ( call->output_tags & CALL_FMT_GP ) + bcf_update_format_float(call->hdr, rec, "GP", call->GPs, nsmpl*ngts_new); + if ( call->output_tags & CALL_FMT_GQ ) + bcf_update_format_int32(call->hdr, rec, "GQ", call->GQs, nsmpl); + + mcall_trim_and_update_PLs(call,rec,nals_ori,call->nals_new); } - if ( nals!=nout ) mcall_trim_numberR(call, rec, nals, nout, out_als); + if ( nals_ori!=call->nals_new ) + mcall_trim_and_update_numberR(call,rec,nals_ori,call->nals_new); - // Set QUAL and calculate HWE-related annotations + // Set QUAL if ( nAC ) { - float icb = calc_ICB(call->ac[0],nAC, call->nhets, call->ndiploid); - if ( icb != HUGE_VAL ) bcf_update_info_float(call->hdr, rec, "ICB", &icb, 1); - - float hob = calc_HOB(call->ac[0],nAC, call->nhets, call->ndiploid); - if ( hob != HUGE_VAL ) bcf_update_info_float(call->hdr, rec, "HOB", &hob, 1); - // Quality of a variant site. fabs() to avoid negative zeros in VCF output when CALL_KEEPALT is set - rec->qual = -4.343*(call->ref_lk - logsumexp2(call->lk_sum,call->ref_lk)); + rec->qual = max_qual; } else { // Set the quality of a REF site - if ( call->lk_sum==-HUGE_VAL ) // no support from (high quality) reads, so QUAL=1-prior + if ( lk_sum!=-HUGE_VAL ) // no support from (high quality) reads, so QUAL=1-prior + rec->qual = -4.343*(lk_sum - logsumexp2(lk_sum,ref_lk)); + else if ( call->ac[0] ) rec->qual = call->theta ? -4.343*call->theta : 0; else - rec->qual = -4.343*(call->lk_sum - logsumexp2(call->lk_sum,call->ref_lk)); + bcf_float_set_missing(rec->qual); } - if ( rec->qual>999 ) rec->qual = 999; - if ( rec->qual>50 ) rec->qual = rint(rec->qual); - // AC, AN - if ( nout>1 ) bcf_update_info_int32(call->hdr, rec, "AC", call->ac+1, nout-1); + if ( call->nals_new>1 ) bcf_update_info_int32(call->hdr, rec, "AC", call->ac+1, call->nals_new-1); nAC += call->ac[0]; bcf_update_info_int32(call->hdr, rec, "AN", &nAC, 1); // Remove unused alleles - hts_expand(char*,nout,call->nals,call->als); - for (i=0; inals_new,call->nals,call->als); + for (i=0; ials_map[i]>=0 ) call->als[call->als_map[i]] = rec->d.allele[i]; - bcf_update_alleles(call->hdr, rec, (const char**)call->als, nout); + bcf_update_alleles(call->hdr, rec, (const char**)call->als, call->nals_new); bcf_update_genotypes(call->hdr, rec, call->gts, nsmpl*2); - // DP4 tag + // DP4 and PV4 tags if ( bcf_get_info_float(call->hdr, rec, "I16", &call->anno16, &call->n16)==16 ) { int32_t dp[4]; dp[0] = call->anno16[0]; dp[1] = call->anno16[1]; dp[2] = call->anno16[2]; dp[3] = call->anno16[3]; @@ -1711,10 +1672,22 @@ int mcall(call_t *call, bcf1_t *rec) int32_t mq = (call->anno16[8]+call->anno16[10])/(call->anno16[0]+call->anno16[1]+call->anno16[2]+call->anno16[3]); bcf_update_info_int32(call->hdr, rec, "MQ", &mq, 1); + + if ( call->output_tags & CALL_FMT_PV4 ) + { + anno16_t a; + float tmpf[4]; + int is_tested = test16(call->anno16, &a) >= 0 && a.is_tested ? 1 : 0; + if ( is_tested ) + { + for (i=0; i<4; i++) tmpf[i] = a.p[i]; + bcf_update_info_float(call->hdr, rec, "PV4", tmpf, 4); + } + } } bcf_update_info_int32(call->hdr, rec, "I16", NULL, 0); // remove I16 tag - return nout; + return call->nals_new; } diff --git a/test/call-G.1.out b/test/call-G.1.out new file mode 100644 index 000000000..171313184 --- /dev/null +++ b/test/call-G.1.out @@ -0,0 +1,17 @@ +##fileformat=VCFv4.2 +##FILTER= +##reference=test.fna +##contig= +##INFO= +##FORMAT= +##FORMAT= +##FORMAT= +##INFO= +##INFO= +##FORMAT= +##INFO= +##INFO= +##INFO= +##INFO= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT smpl1 smpl2 smpl3 smpl4 smpl5 smpl6 smpl7 smpl8 smpl9 smpl10 smpl11 smpl12 smpl13 smpl14 smpl15 smpl16 smpl17 smpl18 smpl19 smpl20 smpl21 smpl22 smpl23 smpl24 smpl25 smpl26 culprit smpl28 smpl29 smpl30 smpl31 smpl32 smpl33 +1 378 . G A,T 461.915 . DP=1418;AN_POP=1000;AC_POP=500,0,0;AC=5,5;AN=48;DP4=675,356,36,17;MQ=57 GT:PL:AD:QS 0/0:0,3,60,3,60,60:1,0,0:60,0,0 0/0:0,33,255,33,255,255:11,0,0:425,0,0 0/2:98,119,255,0,173,155:7,0,6:224,0,175 0/0:0,12,155,12,155,155:4,0,0:174,0,0 0/0:0,24,244,24,244,244:8,0,0:324,0,0 ./.:0,0,0,0,0,0:0,0,0:0,0,0 ./.:0,0,0,0,0,0:0,0,0:0,0,0 ./.:0,0,0,0,0,0:0,0,0:0,0,0 ./.:0,0,0,0,0,0:0,0,0:0,0,0 ./.:0,0,0,0,0,0:0,0,0:0,0,0 ./.:0,0,0,0,0,0:0,0,0:0,0,0 0/2:37,37,37,3,3,0:0,0,1:0,0,37 0/0:0,15,161,15,161,161:5,0,0:187,0,0 0/0:0,15,192,15,192,192:5,0,0:219,0,0 0/0:0,6,94,6,94,94:2,0,0:101,0,0 ./.:0,0,0,0,0,0:0,0,0:0,0,0 0/0:0,27,255,27,255,255:9,0,0:399,0,0 0/1:31,0,19,34,22,53:1,1,0:25,37,0 0/0:0,18,165,18,165,165:6,0,0:222,0,0 ./.:0,0,0,0,0,0:0,0,0:0,0,0 0/1:49,3,0,49,3,49:0,1,0:0,49,0 ./.:0,0,0,0,0,0:0,0,0:0,0,0 0/0:0,3,41,3,41,41:1,0,0:41,0,0 0/0:0,39,255,39,255,255:13,0,0:502,0,0 0/0:11,0,227,38,230,255:9,1,0:399,41,0 0/0:0,3,41,3,41,41:1,0,0:41,0,0 1/2:218,98,86,135,0,126:0,4,3:0,160,112 0/0:0,21,205,21,205,205:7,0,0:250,0,0 0/0:0,18,200,18,200,200:6,0,0:261,0,0 0/2:82,91,200,0,119,113:3,0,2:133,0,97 0/0:0,9,137,9,137,137:3,0,0:161,0,0 0/1:73,38,35,73,38,73:0,1,0:0,41,0 1/2:107,19,47,57,0,88:1,3,1:41,123,60 diff --git a/test/call-G.2.1.out b/test/call-G.2.1.out new file mode 100644 index 000000000..1c5a720f7 --- /dev/null +++ b/test/call-G.2.1.out @@ -0,0 +1,27 @@ +##fileformat=VCFv4.2 +##FILTER= +##reference=file://hwe.fa +##contig= +##ALT= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##FORMAT= +##FORMAT= +##FORMAT= +##INFO= +##INFO= +##INFO= +##INFO= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA07000 NA07056 NA12046 NA12144 NA12156 NA12234 NA12249 NA12282 NA12283 NA12286 NA12340 NA12341 NA12342 NA12399 NA12414 NA12546 NA12717 NA12718 NA12748 NA12749 NA12750 NA12776 NA12828 NA12843 NA18502 NA18867 NA18908 NA18933 NA18934 NA19108 NA19117 NA19119 NA19129 NA19131 NA19138 NA19147 NA19149 NA19159 NA19185 NA19189 NA19197 NA19198 NA19206 NA19213 NA19214 NA19223 NA19235 NA19236 NA19238 NA19239 NA19248 NA19256 +15 201 . A G 858.4 . DP=80;AN_POP=1000;AC_POP=1000,0;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;AC=79;AN=104;DP4=19,17,21,15;MQ=59 GT:PL:AD 0/1:0,3,36:1,0 0/1:0,6,93:2,0 0/1:0,3,33:1,0 0/1:0,3,32:1,0 0/1:0,3,36:1,0 0/1:0,3,32:1,0 0/1:0,3,40:1,0 0/1:0,6,79:2,0 0/1:0,3,42:1,0 0/1:0,6,64:2,0 0/1:0,6,67:2,0 0/1:0,6,66:2,0 0/1:0,9,96:3,0 0/1:0,3,32:1,0 0/1:0,3,36:1,0 0/1:0,3,31:1,0 0/1:0,3,44:1,0 0/1:0,3,44:1,0 0/1:0,6,61:2,0 0/1:0,3,60:1,0 0/1:0,6,60:2,0 0/1:0,6,57:2,0 0/1:0,3,31:1,0 0/1:0,6,68:2,0 0/1:0,3,60:1,0 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:37,3,0:0,1 1/1:60,3,0:0,1 1/1:38,3,0:0,1 1/1:36,3,0:0,1 1/1:37,3,0:0,1 1/1:85,9,0:0,3 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:63,6,0:0,2 1/1:29,3,0:0,1 1/1:32,3,0:0,1 1/1:43,3,0:0,1 1/1:27,3,0:0,1 1/1:43,3,0:0,1 1/1:34,3,0:0,1 1/1:46,3,0:0,1 1/1:144,12,0:0,4 1/1:59,3,0:0,1 1/1:16,3,0:0,1 1/1:37,3,0:0,1 1/1:57,6,0:0,2 1/1:91,6,0:0,2 1/1:21,3,0:0,1 1/1:52,6,0:0,2 1/1:26,3,0:0,1 diff --git a/test/call-G.2.out b/test/call-G.2.out new file mode 100644 index 000000000..8709e0515 --- /dev/null +++ b/test/call-G.2.out @@ -0,0 +1,17 @@ +##fileformat=VCFv4.2 +##FILTER= +##reference=test.fna +##contig= +##INFO= +##FORMAT= +##FORMAT= +##FORMAT= +##INFO= +##INFO= +##FORMAT= +##INFO= +##INFO= +##INFO= +##INFO= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT smpl1 smpl2 smpl3 smpl4 smpl5 smpl6 smpl7 smpl8 smpl9 smpl10 smpl11 smpl12 smpl13 smpl14 smpl15 smpl16 smpl17 smpl18 smpl19 smpl20 smpl21 smpl22 smpl23 smpl24 smpl25 smpl26 culprit smpl28 smpl29 smpl30 smpl31 smpl32 smpl33 +1 378 . G A,T,C 169.281 . DP=1418;AN_POP=1000;AC_POP=500,0,0;AC=6,5,1;AN=48;DP4=675,356,36,17;MQ=57 GT:PL:AD:QS 0/0:0,3,60,3,60,60,3,60,60,60:1,0,0,0:60,0,0,0 0/0:0,33,255,33,255,255,33,255,255,255:11,0,0,0:425,0,0,0 0/2:98,119,255,0,173,155,119,255,173,255:7,0,6,0:224,0,175,0 0/0:0,12,155,12,155,155,12,155,155,155:4,0,0,0:174,0,0,0 0/0:0,24,244,24,244,244,24,244,244,244:8,0,0,0:324,0,0,0 ./.:0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 ./.:0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 ./.:0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 ./.:0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 ./.:0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 ./.:0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 2/2:37,37,37,3,3,0,37,37,3,37:0,0,1,0:0,0,37,0 0/0:0,15,161,15,161,161,15,161,161,161:5,0,0,0:187,0,0,0 0/0:0,15,192,15,192,192,15,192,192,192:5,0,0,0:219,0,0,0 0/0:0,6,94,6,94,94,6,94,94,94:2,0,0,0:101,0,0,0 ./.:0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 0/0:0,27,255,27,255,255,27,255,255,255:9,0,0,0:399,0,0,0 0/1:31,0,19,34,22,53,34,22,53,53:1,1,0,0:25,37,0,0 0/0:0,18,165,18,165,165,18,165,165,165:6,0,0,0:222,0,0,0 ./.:0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 1/1:49,3,0,49,3,49,49,3,49,49:0,1,0,0:0,49,0,0 ./.:0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 0/0:0,3,41,3,41,41,3,41,41,41:1,0,0,0:41,0,0,0 0/0:0,39,255,39,255,255,39,255,255,255:13,0,0,0:502,0,0,0 0/0:11,0,227,38,230,255,38,230,255,255:9,1,0,0:399,41,0,0 0/0:0,3,41,3,41,41,3,41,41,41:1,0,0,0:41,0,0,0 1/2:218,98,86,135,0,126,218,98,135,218:0,4,3,0:0,160,112,0 0/0:0,21,205,21,205,205,21,205,205,205:7,0,0,0:250,0,0,0 0/0:0,18,200,18,200,200,18,200,200,200:6,0,0,0:261,0,0,0 0/2:82,91,200,0,119,113,91,200,119,200:3,0,2,0:133,0,97,0 0/0:0,9,137,9,137,137,9,137,137,137:3,0,0,0:161,0,0,0 1/3:73,38,35,73,38,73,38,0,38,35:0,1,0,1:0,41,0,41 0/1:107,19,47,57,0,88,110,56,91,141:1,3,1,0:41,123,60,0 diff --git a/test/call-G.2.vcf b/test/call-G.2.vcf new file mode 100644 index 000000000..531e3bd23 --- /dev/null +++ b/test/call-G.2.vcf @@ -0,0 +1,24 @@ +##fileformat=VCFv4.2 +##FILTER= +##reference=file://hwe.fa +##contig= +##ALT= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##FORMAT= +##FORMAT= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA07000 NA07056 NA12046 NA12144 NA12156 NA12234 NA12249 NA12282 NA12283 NA12286 NA12340 NA12341 NA12342 NA12399 NA12414 NA12546 NA12717 NA12718 NA12748 NA12749 NA12750 NA12776 NA12828 NA12843 NA18502 NA18867 NA18908 NA18933 NA18934 NA19108 NA19117 NA19119 NA19129 NA19131 NA19138 NA19147 NA19149 NA19159 NA19185 NA19189 NA19197 NA19198 NA19206 NA19213 NA19214 NA19223 NA19235 NA19236 NA19238 NA19239 NA19248 NA19256 +15 201 . A G,<*> 0 . DP=80;AN_POP=1000;AC_POP=1000,0;I16=19,17,21,15,1380,57328,1308,52148,2160,129600,2137,127369,613,12599,514,10028;QS=25,27,0;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0 PL:AD 0,3,36,3,36,36:1,0,0 0,6,93,6,93,93:2,0,0 0,3,33,3,33,33:1,0,0 0,3,32,3,32,32:1,0,0 0,3,36,3,36,36:1,0,0 0,3,32,3,32,32:1,0,0 0,3,40,3,40,40:1,0,0 0,6,79,6,79,79:2,0,0 0,3,42,3,42,42:1,0,0 0,6,64,6,64,64:2,0,0 0,6,67,6,67,67:2,0,0 0,6,66,6,66,66:2,0,0 0,9,96,9,96,96:3,0,0 0,3,32,3,32,32:1,0,0 0,3,36,3,36,36:1,0,0 0,3,31,3,31,31:1,0,0 0,3,44,3,44,44:1,0,0 0,3,44,3,44,44:1,0,0 0,6,61,6,61,61:2,0,0 0,3,60,3,60,60:1,0,0 0,6,60,6,60,60:2,0,0 0,6,57,6,57,57:2,0,0 0,3,31,3,31,31:1,0,0 0,6,68,6,68,68:2,0,0 0,3,60,3,60,60:1,0,0 39,3,0,39,3,39:0,1,0 31,3,0,31,3,31:0,1,0 37,3,0,37,3,37:0,1,0 60,3,0,60,3,60:0,1,0 38,3,0,38,3,38:0,1,0 36,3,0,36,3,36:0,1,0 37,3,0,37,3,37:0,1,0 85,9,0,85,9,85:0,3,0 39,3,0,39,3,39:0,1,0 31,3,0,31,3,31:0,1,0 63,6,0,63,6,63:0,2,0 29,3,0,29,3,29:0,1,0 32,3,0,32,3,32:0,1,0 43,3,0,43,3,43:0,1,0 27,3,0,27,3,27:0,1,0 43,3,0,43,3,43:0,1,0 34,3,0,34,3,34:0,1,0 46,3,0,46,3,46:0,1,0 144,12,0,144,12,144:0,4,0 59,3,0,59,3,59:0,1,0 16,3,0,16,3,16:0,1,0 37,3,0,37,3,37:0,1,0 57,6,0,57,6,57:0,2,0 91,6,0,91,6,91:0,2,0 21,3,0,21,3,21:0,1,0 52,6,0,52,6,52:0,2,0 26,3,0,26,3,26:0,1,0 diff --git a/test/call-G.vcf b/test/call-G.vcf new file mode 100644 index 000000000..d90044735 --- /dev/null +++ b/test/call-G.vcf @@ -0,0 +1,13 @@ +##fileformat=VCFv4.2 +##reference=test.fna +##contig= +##INFO= +##INFO= +##INFO= +##FORMAT= +##FORMAT= +##FORMAT= +##INFO= +##INFO= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT smpl1 smpl2 smpl3 smpl4 smpl5 smpl6 smpl7 smpl8 smpl9 smpl10 smpl11 smpl12 smpl13 smpl14 smpl15 smpl16 smpl17 smpl18 smpl19 smpl20 smpl21 smpl22 smpl23 smpl24 smpl25 smpl26 culprit smpl28 smpl29 smpl30 smpl31 smpl32 smpl33 +1 378 . G A,T,C 0 . DP=1418;AN_POP=1000;AC_POP=500,0,0;I16=675,356,36,17,49274,2.62616e+06,2313,113051,58808,3.52801e+06,3000,180000,21513,499753,1170,27336;QS=247.89,7.1438,4.46622,0.5 PL:AD:QS 0,3,60,3,60,60,3,60,60,60:1,0,0,0:60,0,0,0 0,33,255,33,255,255,33,255,255,255:11,0,0,0:425,0,0,0 98,119,255,0,173,155,119,255,173,255:7,0,6,0:224,0,175,0 0,12,155,12,155,155,12,155,155,155:4,0,0,0:174,0,0,0 0,24,244,24,244,244,24,244,244,244:8,0,0,0:324,0,0,0 0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 37,37,37,3,3,0,37,37,3,37:0,0,1,0:0,0,37,0 0,15,161,15,161,161,15,161,161,161:5,0,0,0:187,0,0,0 0,15,192,15,192,192,15,192,192,192:5,0,0,0:219,0,0,0 0,6,94,6,94,94,6,94,94,94:2,0,0,0:101,0,0,0 0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 0,27,255,27,255,255,27,255,255,255:9,0,0,0:399,0,0,0 31,0,19,34,22,53,34,22,53,53:1,1,0,0:25,37,0,0 0,18,165,18,165,165,18,165,165,165:6,0,0,0:222,0,0,0 0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 49,3,0,49,3,49,49,3,49,49:0,1,0,0:0,49,0,0 0,0,0,0,0,0,0,0,0,0:0,0,0,0:0,0,0,0 0,3,41,3,41,41,3,41,41,41:1,0,0,0:41,0,0,0 0,39,255,39,255,255,39,255,255,255:13,0,0,0:502,0,0,0 11,0,227,38,230,255,38,230,255,255:9,1,0,0:399,41,0,0 0,3,41,3,41,41,3,41,41,41:1,0,0,0:41,0,0,0 218,98,86,135,0,126,218,98,135,218:0,4,3,0:0,160,112,0 0,21,205,21,205,205,21,205,205,205:7,0,0,0:250,0,0,0 0,18,200,18,200,200,18,200,200,200:6,0,0,0:261,0,0,0 82,91,200,0,119,113,91,200,119,200:3,0,2,0:133,0,97,0 0,9,137,9,137,137,9,137,137,137:3,0,0,0:161,0,0,0 73,38,35,73,38,73,38,0,38,35:0,1,0,1:0,41,0,41 107,19,47,57,0,88,110,56,91,141:1,3,1,0:41,123,60,0 diff --git a/test/call.af-fixation.1.out b/test/call.af-fixation.1.out new file mode 100644 index 000000000..ede414aa5 --- /dev/null +++ b/test/call.af-fixation.1.out @@ -0,0 +1,13 @@ +##fileformat=VCFv4.2 +##FILTER= +##reference=file:///lustre/scratch116/vr/ref/human/GRCh37/hs37d5.fa +##contig= +##FORMAT= +##FORMAT= +##FORMAT= +##INFO= +##INFO= +##INFO= +##INFO= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA07051_CEU NA12843_CEU NA18488_YRI NA18489_YRI NA18498_YRI NA18504_YRI NA18510_YRI NA18519_YRI NA18858_YRI NA18861_YRI NA18868_YRI NA18871_YRI NA18874_YRI NA18881_YRI NA18907_YRI NA18909_YRI NA18917_YRI NA18933_YRI NA18934_YRI NA19108_YRI NA19119_YRI NA19129_YRI NA19131_YRI NA19138_YRI NA19147_YRI NA19189_YRI NA19197_YRI NA19198_YRI NA19206_YRI NA19214_YRI NA19223_YRI NA19238_YRI NA19239_YRI NA19248_YRI NA19256_YRI +2 136608646 . G . 24.5121 . AN=70;DP4=95,24,45,11;MQ=59 GT:AD 0/0:0 0/0:1 0/0:2 0/0:1 0/0:2 0/0:2 0/0:2 0/0:1 0/0:7 0/0:8 0/0:2 0/0:1 0/0:2 0/0:6 0/0:7 0/0:4 0/0:7 0/0:3 0/0:2 0/0:5 0/0:1 0/0:4 0/0:1 0/0:4 0/0:2 0/0:6 0/0:1 0/0:2 0/0:1 0/0:1 0/0:1 0/0:4 0/0:2 0/0:1 0/0:2 diff --git a/test/call.af-fixation.2.out b/test/call.af-fixation.2.out new file mode 100644 index 000000000..de69f88c9 --- /dev/null +++ b/test/call.af-fixation.2.out @@ -0,0 +1,13 @@ +##fileformat=VCFv4.2 +##FILTER= +##reference=file:///lustre/scratch116/vr/ref/human/GRCh37/hs37d5.fa +##contig= +##FORMAT= +##FORMAT= +##FORMAT= +##INFO= +##INFO= +##INFO= +##INFO= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA07051_CEU NA12843_CEU NA18488_YRI NA18489_YRI NA18498_YRI NA18504_YRI NA18510_YRI NA18519_YRI NA18858_YRI NA18861_YRI NA18868_YRI NA18871_YRI NA18874_YRI NA18881_YRI NA18907_YRI NA18909_YRI NA18917_YRI NA18933_YRI NA18934_YRI NA19108_YRI NA19119_YRI NA19129_YRI NA19131_YRI NA19138_YRI NA19147_YRI NA19189_YRI NA19197_YRI NA19198_YRI NA19206_YRI NA19214_YRI NA19223_YRI NA19238_YRI NA19239_YRI NA19248_YRI NA19256_YRI +2 136608646 . G A 4.92443 . AC=2;AN=70;DP4=95,24,45,11;MQ=59 GT:PL:AD 0/1:32,3,0:0,1 0/1:0,3,37:1,0 0/0:0,6,77:2,0 0/0:0,3,33:1,0 0/0:0,6,67:2,0 0/0:0,6,67:2,0 0/0:0,6,67:2,0 0/0:0,3,60:1,0 0/0:0,21,147:7,0 0/0:0,24,161:8,0 0/0:0,6,52:2,0 0/0:0,3,38:1,0 0/0:0,6,72:2,0 0/0:0,18,197:6,0 0/0:0,21,175:7,0 0/0:0,12,129:4,0 0/0:0,21,220:7,0 0/0:0,9,111:3,0 0/0:0,6,76:2,0 0/0:0,15,125:5,0 0/0:0,3,32:1,0 0/0:0,12,110:4,0 0/0:0,3,41:1,0 0/0:0,12,124:4,0 0/0:0,6,70:2,0 0/0:0,18,140:6,0 0/0:0,3,22:1,0 0/0:0,6,63:2,0 0/0:0,3,46:1,0 0/0:0,3,39:1,0 0/0:0,3,42:1,0 0/0:0,12,142:4,0 0/0:0,6,67:2,0 0/0:0,3,39:1,0 0/0:0,6,54:2,0 diff --git a/test/call.af-fixation.3.out b/test/call.af-fixation.3.out new file mode 100644 index 000000000..37644417b --- /dev/null +++ b/test/call.af-fixation.3.out @@ -0,0 +1,15 @@ +##fileformat=VCFv4.2 +##FILTER= +##reference=file:///lustre/scratch116/vr/ref/human/GRCh37/hs37d5.fa +##contig= +##FORMAT= +##FORMAT= +##FORMAT= +##FORMAT= +##FORMAT= +##INFO= +##INFO= +##INFO= +##INFO= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA07051_CEU NA12843_CEU NA18488_YRI NA18489_YRI NA18498_YRI NA18504_YRI NA18510_YRI NA18519_YRI NA18858_YRI NA18861_YRI NA18868_YRI NA18871_YRI NA18874_YRI NA18881_YRI NA18907_YRI NA18909_YRI NA18917_YRI NA18933_YRI NA18934_YRI NA19108_YRI NA19119_YRI NA19129_YRI NA19131_YRI NA19138_YRI NA19147_YRI NA19189_YRI NA19197_YRI NA19198_YRI NA19206_YRI NA19214_YRI NA19223_YRI NA19238_YRI NA19239_YRI NA19248_YRI NA19256_YRI +2 136608646 . G A 4.92443 . AC=2;AN=70;DP4=95,24,45,11;MQ=59 GT:PL:AD:GP:GQ 0/1:32,3,0:0,1:0.000315005,0.500435,0.49925:3 0/1:0,3,37:1,0:0.499357,0.500543,9.96349e-05:3 0/0:0,6,77:2,0:1,0,0:127 0/0:0,3,33:1,0:1,0,0:127 0/0:0,6,67:2,0:1,0,0:127 0/0:0,6,67:2,0:1,0,0:127 0/0:0,6,67:2,0:1,0,0:127 0/0:0,3,60:1,0:1,0,0:127 0/0:0,21,147:7,0:1,0,0:127 0/0:0,24,161:8,0:1,0,0:127 0/0:0,6,52:2,0:1,0,0:127 0/0:0,3,38:1,0:1,0,0:127 0/0:0,6,72:2,0:1,0,0:127 0/0:0,18,197:6,0:1,0,0:127 0/0:0,21,175:7,0:1,0,0:127 0/0:0,12,129:4,0:1,0,0:127 0/0:0,21,220:7,0:1,0,0:127 0/0:0,9,111:3,0:1,0,0:127 0/0:0,6,76:2,0:1,0,0:127 0/0:0,15,125:5,0:1,0,0:127 0/0:0,3,32:1,0:1,0,0:127 0/0:0,12,110:4,0:1,0,0:127 0/0:0,3,41:1,0:1,0,0:127 0/0:0,12,124:4,0:1,0,0:127 0/0:0,6,70:2,0:1,0,0:127 0/0:0,18,140:6,0:1,0,0:127 0/0:0,3,22:1,0:1,0,0:127 0/0:0,6,63:2,0:1,0,0:127 0/0:0,3,46:1,0:1,0,0:127 0/0:0,3,39:1,0:1,0,0:127 0/0:0,3,42:1,0:1,0,0:127 0/0:0,12,142:4,0:1,0,0:127 0/0:0,6,67:2,0:1,0,0:127 0/0:0,3,39:1,0:1,0,0:127 0/0:0,6,54:2,0:1,0,0:127 diff --git a/test/call.af-fixation.txt b/test/call.af-fixation.txt new file mode 100644 index 000000000..03892bb55 --- /dev/null +++ b/test/call.af-fixation.txt @@ -0,0 +1,35 @@ +NA07051_CEU CEU +NA12843_CEU CEU +NA18488_YRI YRI +NA18489_YRI YRI +NA18498_YRI YRI +NA18504_YRI YRI +NA18510_YRI YRI +NA18519_YRI YRI +NA18858_YRI YRI +NA18861_YRI YRI +NA18868_YRI YRI +NA18871_YRI YRI +NA18874_YRI YRI +NA18881_YRI YRI +NA18907_YRI YRI +NA18909_YRI YRI +NA18917_YRI YRI +NA18933_YRI YRI +NA18934_YRI YRI +NA19108_YRI YRI +NA19119_YRI YRI +NA19129_YRI YRI +NA19131_YRI YRI +NA19138_YRI YRI +NA19147_YRI YRI +NA19189_YRI YRI +NA19197_YRI YRI +NA19198_YRI YRI +NA19206_YRI YRI +NA19214_YRI YRI +NA19223_YRI YRI +NA19238_YRI YRI +NA19239_YRI YRI +NA19248_YRI YRI +NA19256_YRI YRI diff --git a/test/call.af-fixation.vcf b/test/call.af-fixation.vcf new file mode 100644 index 000000000..700178c7e --- /dev/null +++ b/test/call.af-fixation.vcf @@ -0,0 +1,9 @@ +##fileformat=VCFv4.2 +##reference=file:///lustre/scratch116/vr/ref/human/GRCh37/hs37d5.fa +##contig= +##INFO= +##INFO= +##FORMAT= +##FORMAT= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA07051_CEU NA12843_CEU NA18488_YRI NA18489_YRI NA18498_YRI NA18504_YRI NA18510_YRI NA18519_YRI NA18858_YRI NA18861_YRI NA18868_YRI NA18871_YRI NA18874_YRI NA18881_YRI NA18907_YRI NA18909_YRI NA18917_YRI NA18933_YRI NA18934_YRI NA19108_YRI NA19119_YRI NA19129_YRI NA19131_YRI NA19138_YRI NA19147_YRI NA19189_YRI NA19197_YRI NA19198_YRI NA19206_YRI NA19214_YRI NA19223_YRI NA19238_YRI NA19239_YRI NA19248_YRI NA19256_YRI +2 136608646 . G A,T,<*> 0 . I16=95,24,45,11,4466,172210,1917,71261,7086,423410,3283,194379,2130,47014,1038,22600;QS=40.8826,20.7918,0.325581,0 PL:AD 32,3,0,32,3,32,32,3,32,32:0,1,0,0 0,3,37,3,37,37,3,37,37,37:1,0,0,0 0,6,77,6,77,77,6,77,77,77:2,0,0,0 0,3,33,3,33,33,3,33,33,33:1,0,0,0 0,6,67,6,67,67,6,67,67,67:2,0,0,0 0,6,67,6,67,67,6,67,67,67:2,0,0,0 0,6,67,6,67,67,6,67,67,67:2,0,0,0 0,3,60,3,60,60,3,60,60,60:1,0,0,0 0,21,147,21,147,147,21,147,147,147:7,0,0,0 0,24,161,24,161,161,24,161,161,161:8,0,0,0 0,6,52,6,52,52,6,52,52,52:2,0,0,0 0,3,38,3,38,38,3,38,38,38:1,0,0,0 0,6,72,6,72,72,6,72,72,72:2,0,0,0 0,18,197,18,197,197,18,197,197,197:6,0,0,0 0,21,175,21,175,175,21,175,175,175:7,0,0,0 0,12,129,12,129,129,12,129,129,129:4,0,0,0 0,21,220,21,220,220,21,220,220,220:7,0,0,0 0,9,111,9,111,111,9,111,111,111:3,0,0,0 0,6,76,6,76,76,6,76,76,76:2,0,0,0 0,15,125,15,125,125,15,125,125,125:5,0,0,0 0,3,32,3,32,32,3,32,32,32:1,0,0,0 0,12,110,12,110,110,12,110,110,110:4,0,0,0 0,3,41,3,41,41,3,41,41,41:1,0,0,0 0,12,124,12,124,124,12,124,124,124:4,0,0,0 0,6,70,6,70,70,6,70,70,70:2,0,0,0 0,18,140,18,140,140,18,140,140,140:6,0,0,0 0,3,22,3,22,22,3,22,22,22:1,0,0,0 0,6,63,6,63,63,6,63,63,63:2,0,0,0 0,3,46,3,46,46,3,46,46,46:1,0,0,0 0,3,39,3,39,39,3,39,39,39:1,0,0,0 0,3,42,3,42,42,3,42,42,42:1,0,0,0 0,12,142,12,142,142,12,142,142,142:4,0,0,0 0,6,67,6,67,67,6,67,67,67:2,0,0,0 0,3,39,3,39,39,3,39,39,39:1,0,0,0 0,6,54,6,54,54,6,54,54,54:2,0,0,0 diff --git a/test/check.chk b/test/check.chk index 8a2ce22f6..31e4eb35b 100644 --- a/test/check.chk +++ b/test/check.chk @@ -38,17 +38,19 @@ AF 0 0.000000 3 1 2 2 0 0 2 AF 0 0.490000 0 0 0 12 0 0 12 AF 0 0.740000 1 1 0 0 0 0 0 AF 0 0.990000 1 1 0 0 0 0 0 -# QUAL, Stats by quality: +# QUAL, Stats by quality # QUAL [2]id [3]Quality [4]number of SNPs [5]number of transitions (1st ALT) [6]number of transversions (1st ALT) [7]number of indels -QUAL 0 12 1 0 1 1 -QUAL 0 45 0 0 0 1 -QUAL 0 59 2 1 0 3 -QUAL 0 60 1 1 0 0 -QUAL 0 61 0 0 0 1 -QUAL 0 79 0 0 0 1 -QUAL 0 82 0 0 0 1 -QUAL 0 90 1 1 0 0 -QUAL 0 342 0 0 0 1 +QUAL 0 12.6 1 0 1 0 +QUAL 0 12.9 0 0 0 1 +QUAL 0 45.0 0 0 0 1 +QUAL 0 59.2 1 1 0 0 +QUAL 0 59.9 0 0 0 3 +QUAL 0 60.2 1 1 0 0 +QUAL 0 61.5 0 0 0 1 +QUAL 0 79.0 0 0 0 1 +QUAL 0 82.7 0 0 0 1 +QUAL 0 90.6 1 1 0 0 +QUAL 0 342.0 0 0 0 1 # IDD, InDel distribution: # IDD [2]id [3]length (deletions negative) [4]number of sites [5]number of genotypes [6]mean VAF IDD 0 -10 1 0 . diff --git a/test/check_merge.chk b/test/check_merge.chk index 72fab0929..2f54579e4 100644 --- a/test/check_merge.chk +++ b/test/check_merge.chk @@ -54,12 +54,14 @@ DP 0 60 0 0.000000 1 33.333333 DP 0 62 0 0.000000 2 66.666667 # QUAL, Stats by quality # QUAL [2]id [3]Quality [4]number of SNPs [5]number of transitions (1st ALT) [6]number of transversions (1st ALT) [7]number of indels -QUAL 0 12 1 0 1 1 -QUAL 0 45 0 0 0 1 -QUAL 0 59 2 1 0 3 -QUAL 0 60 1 1 0 0 -QUAL 0 61 0 0 0 1 -QUAL 0 79 0 0 0 1 -QUAL 0 82 0 0 0 1 -QUAL 0 90 1 1 0 0 -QUAL 0 342 0 0 0 1 +QUAL 0 12.6 1 0 1 0 +QUAL 0 12.9 0 0 0 1 +QUAL 0 45.0 0 0 0 1 +QUAL 0 59.2 1 1 0 0 +QUAL 0 59.9 0 0 0 3 +QUAL 0 60.2 1 1 0 0 +QUAL 0 61.5 0 0 0 1 +QUAL 0 79.0 0 0 0 1 +QUAL 0 82.7 0 0 0 1 +QUAL 0 90.6 1 1 0 0 +QUAL 0 342.0 0 0 0 1 diff --git a/test/mpileup.1.out b/test/mpileup.1.out index 8cedcb346..c586fb4d6 100644 --- a/test/mpileup.1.out +++ b/test/mpileup.1.out @@ -20,21 +20,19 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100 HG00101 HG00102 -17 302 . T TA 488 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 0/1:157,0,9:7:6 1/1:201,21,0:7:7 -17 828 . T C 409 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 0/1:116,0,91:9:5 1/1:120,12,0:4:4 -17 834 . G A 364 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 0/1:128,0,59:8:5 1/1:89,9,0:3:3 -17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 -17 1869 . A T 138 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0/1:16,0,104:5:1 1/1:42,3,0:1:1 -17 2041 . G A 447 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 0/1:32,0,24:2:1 1/1:223,21,0:7:7 -17 2220 . G A 303 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 0/1:69,0,46:4:2 1/1:131,12,0:4:4 -17 2564 . A G 233 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 0/1:57,0,56:4:2 1/1:124,12,0:4:4 -17 3104 . C T 24.2837 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0/0:0,12,144:4:0 0/1:59,0,93:5:2 -17 3587 . G A 358 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0/1:22,0,118:5:1 1/1:212,24,0:8:8 -17 3936 . A G 469 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 0/1:77,0,58:6:4 1/1:196,24,0:8:8 +17 302 . T TA 487.586 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;AC=4;AN=6;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 0/1:157,0,9:7:6 1/1:201,21,0:7:7 +17 828 . T C 409.29 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 0/1:116,0,91:9:5 1/1:120,12,0:4:4 +17 834 . G A 363.72 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;AC=4;AN=6;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 0/1:128,0,59:8:5 1/1:89,9,0:3:3 +17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 +17 1869 . A T 138.104 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0/1:16,0,104:5:1 1/1:42,3,0:1:1 +17 2041 . G A 447.444 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 0/1:32,0,24:2:1 1/1:223,21,0:7:7 +17 2220 . G A 302.575 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 0/1:69,0,46:4:2 1/1:131,12,0:4:4 +17 2564 . A G 232.697 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 0/1:57,0,56:4:2 1/1:124,12,0:4:4 +17 3104 . C T 24.2837 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;AC=1;AN=6;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0/0:0,12,144:4:0 0/1:59,0,93:5:2 +17 3587 . G A 357.834 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;AC=4;AN=6;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0/1:22,0,118:5:1 1/1:212,24,0:8:8 +17 3936 . A G 469.356 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;AC=4;AN=6;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 0/1:77,0,58:6:4 1/1:196,24,0:8:8 diff --git a/test/mpileup.2.out b/test/mpileup.2.out index d6571607a..43e022fd4 100644 --- a/test/mpileup.2.out +++ b/test/mpileup.2.out @@ -22,33 +22,31 @@ ##INFO= ##INFO= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100 HG00101 HG00102 17 1 . A . . . END=301;MinDP=1 GT:DP ./.:5 ./.:1 ./.:3 -17 302 . T TA 488 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 0/1:157,0,9:7:6 1/1:201,21,0:7:7 +17 302 . T TA 487.586 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;AC=4;AN=6;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 0/1:157,0,9:7:6 1/1:201,21,0:7:7 17 303 . G . . . END=827;MinDP=2 GT:DP 0/0:9 0/0:2 0/0:3 -17 828 . T C 409 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 0/1:116,0,91:9:5 1/1:120,12,0:4:4 +17 828 . T C 409.29 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 0/1:116,0,91:9:5 1/1:120,12,0:4:4 17 829 . T . . . END=833;MinDP=4 GT:DP 0/0:11 0/0:8 0/0:4 -17 834 . G A 364 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 0/1:128,0,59:8:5 1/1:89,9,0:3:3 +17 834 . G A 363.72 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;AC=4;AN=6;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 0/1:128,0,59:8:5 1/1:89,9,0:3:3 17 835 . T . . . END=1664;MinDP=1 GT:DP 0/0:5 0/0:2 0/0:1 -17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 +17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 17 1666 . G . . . END=1868;MinDP=0 GT:DP 0/0:6 0/0:0 0/0:1 -17 1869 . A T 138 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0/1:16,0,104:5:1 1/1:42,3,0:1:1 +17 1869 . A T 138.104 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0/1:16,0,104:5:1 1/1:42,3,0:1:1 17 1870 . C . . . END=2040;MinDP=1 GT:DP 0/0:13 0/0:2 0/0:1 -17 2041 . G A 447 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 0/1:32,0,24:2:1 1/1:223,21,0:7:7 +17 2041 . G A 447.444 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 0/1:32,0,24:2:1 1/1:223,21,0:7:7 17 2042 . G . . . END=2219;MinDP=1 GT:DP 0/0:8 0/0:1 0/0:3 -17 2220 . G A 303 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 0/1:69,0,46:4:2 1/1:131,12,0:4:4 +17 2220 . G A 302.575 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 0/1:69,0,46:4:2 1/1:131,12,0:4:4 17 2221 . G . . . END=2563;MinDP=0 GT:DP 0/0:5 0/0:0 0/0:2 -17 2564 . A G 233 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 0/1:57,0,56:4:2 1/1:124,12,0:4:4 +17 2564 . A G 232.697 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 0/1:57,0,56:4:2 1/1:124,12,0:4:4 17 2565 . A . . . END=3103;MinDP=0 GT:DP 0/0:6 0/0:0 0/0:1 -17 3104 . C T 24.2837 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0/0:0,12,144:4:0 0/1:59,0,93:5:2 +17 3104 . C T 24.2837 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;AC=1;AN=6;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0/0:0,12,144:4:0 0/1:59,0,93:5:2 17 3105 . T . . . END=3586;MinDP=2 GT:DP 0/0:5 0/0:2 0/0:3 -17 3587 . G A 358 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0/1:22,0,118:5:1 1/1:212,24,0:8:8 +17 3587 . G A 357.834 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;AC=4;AN=6;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0/1:22,0,118:5:1 1/1:212,24,0:8:8 17 3588 . A . . . END=3935;MinDP=2 GT:DP 0/0:10 0/0:2 0/0:3 -17 3936 . A G 469 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 0/1:77,0,58:6:4 1/1:196,24,0:8:8 +17 3936 . A G 469.356 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;AC=4;AN=6;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 0/1:77,0,58:6:4 1/1:196,24,0:8:8 17 3937 . C . . . END=4101;MinDP=0 GT:DP 0/0:1 0/0:0 0/0:0 diff --git a/test/mpileup.3.out b/test/mpileup.3.out index 8cedcb346..c586fb4d6 100644 --- a/test/mpileup.3.out +++ b/test/mpileup.3.out @@ -20,21 +20,19 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100 HG00101 HG00102 -17 302 . T TA 488 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 0/1:157,0,9:7:6 1/1:201,21,0:7:7 -17 828 . T C 409 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 0/1:116,0,91:9:5 1/1:120,12,0:4:4 -17 834 . G A 364 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 0/1:128,0,59:8:5 1/1:89,9,0:3:3 -17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 -17 1869 . A T 138 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0/1:16,0,104:5:1 1/1:42,3,0:1:1 -17 2041 . G A 447 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 0/1:32,0,24:2:1 1/1:223,21,0:7:7 -17 2220 . G A 303 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 0/1:69,0,46:4:2 1/1:131,12,0:4:4 -17 2564 . A G 233 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 0/1:57,0,56:4:2 1/1:124,12,0:4:4 -17 3104 . C T 24.2837 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0/0:0,12,144:4:0 0/1:59,0,93:5:2 -17 3587 . G A 358 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0/1:22,0,118:5:1 1/1:212,24,0:8:8 -17 3936 . A G 469 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 0/1:77,0,58:6:4 1/1:196,24,0:8:8 +17 302 . T TA 487.586 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;AC=4;AN=6;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 0/1:157,0,9:7:6 1/1:201,21,0:7:7 +17 828 . T C 409.29 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 0/1:116,0,91:9:5 1/1:120,12,0:4:4 +17 834 . G A 363.72 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;AC=4;AN=6;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 0/1:128,0,59:8:5 1/1:89,9,0:3:3 +17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 +17 1869 . A T 138.104 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0/1:16,0,104:5:1 1/1:42,3,0:1:1 +17 2041 . G A 447.444 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 0/1:32,0,24:2:1 1/1:223,21,0:7:7 +17 2220 . G A 302.575 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 0/1:69,0,46:4:2 1/1:131,12,0:4:4 +17 2564 . A G 232.697 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 0/1:57,0,56:4:2 1/1:124,12,0:4:4 +17 3104 . C T 24.2837 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;AC=1;AN=6;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0/0:0,12,144:4:0 0/1:59,0,93:5:2 +17 3587 . G A 357.834 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;AC=4;AN=6;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0/1:22,0,118:5:1 1/1:212,24,0:8:8 +17 3936 . A G 469.356 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;AC=4;AN=6;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 0/1:77,0,58:6:4 1/1:196,24,0:8:8 diff --git a/test/mpileup.4.out b/test/mpileup.4.out index 9166dfa13..b8ded5132 100644 --- a/test/mpileup.4.out +++ b/test/mpileup.4.out @@ -20,21 +20,19 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00102 HG00101 HG00100 -17 302 . T TA 488 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 1/1:201,21,0:7:7 0/1:157,0,9:7:6 0/1:167,0,96:11:6 -17 828 . T C 409 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 1/1:120,12,0:4:4 0/1:116,0,91:9:5 0/1:211,0,35:12:10 -17 834 . G A 364 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 1/1:89,9,0:3:3 0/1:128,0,59:8:5 0/1:185,0,46:11:9 -17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/1:35,0,51:4:2 0/0:0,27,222:9:0 0/0:0,21,185:7:0 -17 1869 . A T 138 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 1/1:42,3,0:1:1 0/1:16,0,104:5:1 0/1:115,0,224:18:7 -17 2041 . G A 447 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 1/1:223,21,0:7:7 0/1:32,0,24:2:1 0/1:229,0,212:21:11 -17 2220 . G A 303 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 1/1:131,12,0:4:4 0/1:69,0,46:4:2 0/1:139,0,130:12:6 -17 2564 . A G 233 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 1/1:124,12,0:4:4 0/1:57,0,56:4:2 0/1:88,0,78:6:3 -17 3104 . C T 24.2837 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/1:59,0,93:5:2 0/0:0,12,144:4:0 0/0:0,48,255:16:0 -17 3587 . G A 358 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 1/1:212,24,0:8:8 0/1:22,0,118:5:1 0/1:161,0,184:14:7 -17 3936 . A G 469 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 1/1:196,24,0:8:8 0/1:77,0,58:6:4 0/1:233,0,206:20:11 +17 302 . T TA 487.586 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;AC=4;AN=6;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 1/1:201,21,0:7:7 0/1:157,0,9:7:6 0/1:167,0,96:11:6 +17 828 . T C 409.29 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 1/1:120,12,0:4:4 0/1:116,0,91:9:5 0/1:211,0,35:12:10 +17 834 . G A 363.72 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;AC=4;AN=6;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 1/1:89,9,0:3:3 0/1:128,0,59:8:5 0/1:185,0,46:11:9 +17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/1:35,0,51:4:2 0/0:0,27,222:9:0 0/0:0,21,185:7:0 +17 1869 . A T 138.104 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 1/1:42,3,0:1:1 0/1:16,0,104:5:1 0/1:115,0,224:18:7 +17 2041 . G A 447.444 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 1/1:223,21,0:7:7 0/1:32,0,24:2:1 0/1:229,0,212:21:11 +17 2220 . G A 302.575 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 1/1:131,12,0:4:4 0/1:69,0,46:4:2 0/1:139,0,130:12:6 +17 2564 . A G 232.697 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 1/1:124,12,0:4:4 0/1:57,0,56:4:2 0/1:88,0,78:6:3 +17 3104 . C T 24.2837 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;AC=1;AN=6;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/1:59,0,93:5:2 0/0:0,12,144:4:0 0/0:0,48,255:16:0 +17 3587 . G A 357.834 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;AC=4;AN=6;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 1/1:212,24,0:8:8 0/1:22,0,118:5:1 0/1:161,0,184:14:7 +17 3936 . A G 469.356 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;AC=4;AN=6;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 1/1:196,24,0:8:8 0/1:77,0,58:6:4 0/1:233,0,206:20:11 diff --git a/test/mpileup.5.out b/test/mpileup.5.out index 7606e647b..489408504 100644 --- a/test/mpileup.5.out +++ b/test/mpileup.5.out @@ -20,21 +20,19 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00102 HG00101 -17 302 . T TA 327 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;ICB=0.3;HOB=0.125;AC=3;AN=4;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 1/1:201,21,0:7:7 0/1:157,0,9:7:6 -17 828 . T C 202 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;ICB=0.3;HOB=0.125;AC=3;AN=4;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 1/1:120,12,0:4:4 0/1:116,0,91:9:5 -17 834 . G A 183 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;ICB=0.3;HOB=0.125;AC=3;AN=4;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 1/1:89,9,0:3:3 0/1:128,0,59:8:5 -17 1665 . T C 3.26671 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;ICB=0.3;HOB=0.125;AC=1;AN=4;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/1:35,0,51:4:2 0/0:0,27,222:9:0 -17 1869 . A T 25.1744 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;ICB=0.3;HOB=0.125;AC=3;AN=4;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 1/1:42,3,0:1:1 0/1:16,0,104:5:1 -17 2041 . G A 221 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;ICB=0.3;HOB=0.125;AC=3;AN=4;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 1/1:223,21,0:7:7 0/1:32,0,24:2:1 -17 2220 . G A 166 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;ICB=0.3;HOB=0.125;AC=3;AN=4;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 1/1:131,12,0:4:4 0/1:69,0,46:4:2 -17 2564 . A G 147 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;ICB=0.3;HOB=0.125;AC=3;AN=4;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 1/1:124,12,0:4:4 0/1:57,0,56:4:2 -17 3104 . C T 24.6136 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;ICB=0.3;HOB=0.125;AC=1;AN=4;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/1:59,0,93:5:2 0/0:0,12,144:4:0 -17 3587 . G A 199 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;ICB=0.3;HOB=0.125;AC=3;AN=4;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 1/1:212,24,0:8:8 0/1:22,0,118:5:1 -17 3936 . A G 239 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;ICB=0.3;HOB=0.125;AC=3;AN=4;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 1/1:196,24,0:8:8 0/1:77,0,58:6:4 +17 302 . T TA 326.859 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;AC=3;AN=4;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 1/1:201,21,0:7:7 0/1:157,0,9:7:6 +17 828 . T C 202.292 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;AC=3;AN=4;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 1/1:120,12,0:4:4 0/1:116,0,91:9:5 +17 834 . G A 183.191 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;AC=3;AN=4;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 1/1:89,9,0:3:3 0/1:128,0,59:8:5 +17 1665 . T C 3.26671 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;AC=1;AN=4;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/1:35,0,51:4:2 0/0:0,27,222:9:0 +17 1869 . A T 25.1744 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;AC=3;AN=4;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 1/1:42,3,0:1:1 0/1:16,0,104:5:1 +17 2041 . G A 221.292 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;AC=3;AN=4;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 1/1:223,21,0:7:7 0/1:32,0,24:2:1 +17 2220 . G A 166.425 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=3;AN=4;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 1/1:131,12,0:4:4 0/1:69,0,46:4:2 +17 2564 . A G 147.291 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=3;AN=4;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 1/1:124,12,0:4:4 0/1:57,0,56:4:2 +17 3104 . C T 24.6136 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;AC=1;AN=4;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/1:59,0,93:5:2 0/0:0,12,144:4:0 +17 3587 . G A 198.948 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;AC=3;AN=4;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 1/1:212,24,0:8:8 0/1:22,0,118:5:1 +17 3936 . A G 239.224 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;AC=3;AN=4;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 1/1:196,24,0:8:8 0/1:77,0,58:6:4 diff --git a/test/mpileup.X.2.out b/test/mpileup.X.2.out index ab82d4629..a727f707f 100644 --- a/test/mpileup.X.2.out +++ b/test/mpileup.X.2.out @@ -20,21 +20,19 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100 HG00101 HG00102 -X 302 . T TA 482 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 1:157,9:7:6 1/1:201,21,0:7:7 -X 828 . T C 322 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 1:116,91:9:5 1/1:120,12,0:4:4 -X 834 . G A 309 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 1:128,59:8:5 1/1:89,9,0:3:3 -X 1665 . T C 3.44176 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;ICB=0.235294;HOB=0.18;AC=1;AN=5;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0:0,222:9:0 0/1:35,0,51:4:2 -X 1869 . A T 122 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;ICB=0.461538;HOB=0.02;AC=3;AN=5;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0:16,104:5:1 1/1:42,3,0:1:1 -X 2041 . G A 426 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 1:32,24:2:1 1/1:223,21,0:7:7 -X 2220 . G A 259 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 1:69,46:4:2 1/1:131,12,0:4:4 -X 2564 . A G 180 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 1:57,56:4:2 1/1:124,12,0:4:4 -X 3104 . C T 24.8375 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;ICB=0.235294;HOB=0.18;AC=1;AN=5;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0:0,144:4:0 0/1:59,0,93:5:2 -X 3587 . G A 335 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;ICB=0.461538;HOB=0.02;AC=3;AN=5;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0:22,118:5:1 1/1:212,24,0:8:8 -X 3936 . A G 414 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 1:77,58:6:4 1/1:196,24,0:8:8 +X 302 . T TA 482.1 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;AC=4;AN=5;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 1:157,9:7:6 1/1:201,21,0:7:7 +X 828 . T C 322.296 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;AC=4;AN=5;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 1:116,91:9:5 1/1:120,12,0:4:4 +X 834 . G A 309.32 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;AC=4;AN=5;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 1:128,59:8:5 1/1:89,9,0:3:3 +X 1665 . T C 3.44176 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;AC=1;AN=5;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0:0,222:9:0 0/1:35,0,51:4:2 +X 1869 . A T 121.973 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;AC=3;AN=5;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0:16,104:5:1 1/1:42,3,0:1:1 +X 2041 . G A 425.853 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;AC=4;AN=5;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 1:32,24:2:1 1/1:223,21,0:7:7 +X 2220 . G A 258.867 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=4;AN=5;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 1:69,46:4:2 1/1:131,12,0:4:4 +X 2564 . A G 179.94 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=4;AN=5;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 1:57,56:4:2 1/1:124,12,0:4:4 +X 3104 . C T 24.8375 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;AC=1;AN=5;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0:0,144:4:0 0/1:59,0,93:5:2 +X 3587 . G A 335.348 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;AC=3;AN=5;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0:22,118:5:1 1/1:212,24,0:8:8 +X 3936 . A G 413.695 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;AC=4;AN=5;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 1:77,58:6:4 1/1:196,24,0:8:8 diff --git a/test/mpileup.X.out b/test/mpileup.X.out index 0e63b8efa..abce4a154 100644 --- a/test/mpileup.X.out +++ b/test/mpileup.X.out @@ -20,21 +20,19 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100 HG00101 HG00102 -X 302 . T TA 482 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 1:157,9:7:6 1/1:201,21,0:7:7 -X 828 . T C 322 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 1:116,91:9:5 1/1:120,12,0:4:4 -X 834 . G A 309 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 1:128,59:8:5 1/1:89,9,0:3:3 -X 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 -X 1869 . A T 138 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0/1:16,0,104:5:1 1/1:42,3,0:1:1 -X 2041 . G A 447 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 0/1:32,0,24:2:1 1/1:223,21,0:7:7 -X 2220 . G A 303 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 0/1:69,0,46:4:2 1/1:131,12,0:4:4 -X 2564 . A G 233 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 0/1:57,0,56:4:2 1/1:124,12,0:4:4 -X 3104 . C T 24.8375 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;ICB=0.235294;HOB=0.18;AC=1;AN=5;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0:0,144:4:0 0/1:59,0,93:5:2 -X 3587 . G A 335 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;ICB=0.461538;HOB=0.02;AC=3;AN=5;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0:22,118:5:1 1/1:212,24,0:8:8 -X 3936 . A G 414 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;ICB=0.235294;HOB=0.18;AC=4;AN=5;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 1:77,58:6:4 1/1:196,24,0:8:8 +X 302 . T TA 482.1 . INDEL;IDV=7;IMF=1;DP=25;VDB=0.27613;SGB=-4.22417;MQSB=0.0443614;MQ0F=0;AC=4;AN=5;DP4=2,4,8,11;MQ=49 GT:PL:DP:DV 0/1:167,0,96:11:6 1:157,9:7:6 1/1:201,21,0:7:7 +X 828 . T C 322.296 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;AC=4;AN=5;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 1:116,91:9:5 1/1:120,12,0:4:4 +X 834 . G A 309.32 . DP=25;VDB=0.788006;SGB=-4.01214;RPB=0.999233;MQB=1;MQSB=1;BQB=0.821668;MQ0F=0;AC=4;AN=5;DP4=2,3,7,10;MQ=60 GT:PL:DP:DV 0/1:185,0,46:11:9 1:128,59:8:5 1/1:89,9,0:3:3 +X 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 +X 1869 . A T 138.104 . DP=24;VDB=0.928022;SGB=-11.9537;RPB=0.984127;MQB=0.96464;MQSB=0.931547;BQB=0.359155;MQ0F=0;AC=4;AN=6;DP4=6,9,5,4;MQ=58 GT:PL:DP:DV 0/1:115,0,224:18:7 0/1:16,0,104:5:1 1/1:42,3,0:1:1 +X 2041 . G A 447.444 . DP=31;VDB=0.816435;SGB=-4.18892;RPB=0.88473;MQB=0.972375;MQSB=0.968257;BQB=0.311275;MQ0F=0;AC=4;AN=6;DP4=6,5,12,7;MQ=58 GT:PL:DP:DV 0/1:229,0,212:21:11 0/1:32,0,24:2:1 1/1:223,21,0:7:7 +X 2220 . G A 302.575 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=4;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/1:139,0,130:12:6 0/1:69,0,46:4:2 1/1:131,12,0:4:4 +X 2564 . A G 232.697 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=4;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/1:88,0,78:6:3 0/1:57,0,56:4:2 1/1:124,12,0:4:4 +X 3104 . C T 24.8375 . DP=25;VDB=0.8;SGB=0.346553;RPB=0.717391;MQB=0.956522;MQSB=0.962269;BQB=0.978261;MQ0F=0;AC=1;AN=5;DP4=8,15,2,0;MQ=58 GT:PL:DP:DV 0/0:0,48,255:16:0 0:0,144:4:0 0/1:59,0,93:5:2 +X 3587 . G A 335.348 . DP=29;VDB=0.902044;SGB=-3.91326;RPB=0.800999;MQB=1;MQSB=1;BQB=0.156944;MQ0F=0;AC=3;AN=5;DP4=4,7,10,6;MQ=60 GT:PL:DP:DV 0/1:161,0,184:14:7 0:22,118:5:1 1/1:212,24,0:8:8 +X 3936 . A G 413.695 . DP=37;VDB=0.0574114;SGB=-4.60123;RPB=0.741697;MQB=0.812605;MQSB=0.143788;BQB=0.883831;MQ0F=0;AC=4;AN=5;DP4=5,6,6,17;MQ=56 GT:PL:DP:DV 0/1:233,0,206:20:11 1:77,58:6:4 1/1:196,24,0:8:8 diff --git a/test/mpileup.cAls.2.out b/test/mpileup.cAls.2.out index 994bf4291..b27a92c8a 100644 --- a/test/mpileup.cAls.2.out +++ b/test/mpileup.cAls.2.out @@ -18,16 +18,14 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample1 sample2 -chr1 212740 . A G 483 . DP=73;VDB=0.520868;SGB=-1.38232;MQSB=1;MQ0F=0;AC=4;AN=4;DP4=0,0,39,4;MQ=60 GT:PL:DP:AD 1/1:255,72,0:24:0,24 1/1:255,57,0:19:0,19 -chr1 320055 . A G 534 . DP=101;MQSB=1;MQ0F=0;AC=0;AN=4;DP4=52,9,0,0;MQ=60 GT:PL:DP:AD 0/0:0,87,255:29:29,0 0/0:0,96,255:32:32,0 -chr1 486173 . A T 106 . DP=13;VDB=0.074936;SGB=0.620439;RPB=0.810265;MQB=1.01283;MQSB=1;BQB=0.810265;MQ0F=0;ICB=0.3;HOB=0.125;AC=1;AN=4;DP4=3,1,3,0;MQ=60 GT:PL:DP:AD 0/0:0,9,151:3:3,0 0/1:140,0,48:4:1,3 -chr1 511277 . A G 483 . DP=50;VDB=0.0722735;SGB=-1.26186;MQSB=1;MQ0F=0;AC=4;AN=4;DP4=0,0,25,4;MQ=60 GT:PL:DP:AD 1/1:255,30,0:10:0,10 1/1:255,57,0:19:0,19 +chr1 212740 . A G 483.052 . DP=73;VDB=0.520868;SGB=-1.38232;MQSB=1;MQ0F=0;AC=4;AN=4;DP4=0,0,39,4;MQ=60 GT:PL:DP:AD 1/1:255,72,0:24:0,24 1/1:255,57,0:19:0,19 +chr1 320055 . A G 533.95 . DP=101;MQSB=1;MQ0F=0;AC=0;AN=4;DP4=52,9,0,0;MQ=60 GT:PL:DP:AD 0/0:0,87,255:29:29,0 0/0:0,96,255:32:32,0 +chr1 486173 . A T 106.286 . DP=13;VDB=0.074936;SGB=0.620439;RPB=0.810265;MQB=1.01283;MQSB=1;BQB=0.810265;MQ0F=0;AC=1;AN=4;DP4=3,1,3,0;MQ=60 GT:PL:DP:AD 0/0:0,9,151:3:3,0 0/1:140,0,48:4:1,3 +chr1 511277 . A G 483.052 . DP=50;VDB=0.0722735;SGB=-1.26186;MQSB=1;MQ0F=0;AC=4;AN=4;DP4=0,0,25,4;MQ=60 GT:PL:DP:AD 1/1:255,30,0:10:0,10 1/1:255,57,0:19:0,19 chr1 602567 . A G 7.04355 . DP=9;SGB=-0.516033;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0;AC=0;AN=4;DP4=3,1,1,0;MQ=60 GT:PL:DP:AD 0/0:0,3,60:1:1,0 0/0:29,0,140:4:3,1 -chr1 639707 . T A 483 . DP=50;VDB=0.563111;SGB=-1.37269;MQSB=1;MQ0F=0;AC=4;AN=4;DP4=0,0,23,8;MQ=60 GT:PL:DP:AD 1/1:255,42,0:14:0,14 1/1:255,51,0:17:0,17 +chr1 639707 . T A 483.052 . DP=50;VDB=0.563111;SGB=-1.37269;MQSB=1;MQ0F=0;AC=4;AN=4;DP4=0,0,23,8;MQ=60 GT:PL:DP:AD 1/1:255,42,0:14:0,14 1/1:255,51,0:17:0,17 diff --git a/test/mpileup.cAls.3.out b/test/mpileup.cAls.3.out index 3550a5b88..66d58b6b2 100644 --- a/test/mpileup.cAls.3.out +++ b/test/mpileup.cAls.3.out @@ -7,14 +7,12 @@ ##INFO= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample -20 69066 . C G 60 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,33 +20 69066 . C G 59.5765 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,33 20 69093 . G A . . . GT . 20 69094 . G A . . . GT . 20 69408 . C T . . . GT . diff --git a/test/mpileup.cAls.4.out b/test/mpileup.cAls.4.out index 8a5124a63..ef884a104 100644 --- a/test/mpileup.cAls.4.out +++ b/test/mpileup.cAls.4.out @@ -7,15 +7,13 @@ ##INFO= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample -20 69066 . C G 60 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,33 -20 69076 . A G 55 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,28 +20 69066 . C G 59.5765 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,33 +20 69076 . A G 54.5765 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,28 20 69093 . G A . . . GT . 20 69094 . G A . . . GT . 20 69408 . C T . . . GT . diff --git a/test/mpileup.cAls.5.out b/test/mpileup.cAls.5.out index 084486ef1..3eeb85c51 100644 --- a/test/mpileup.cAls.5.out +++ b/test/mpileup.cAls.5.out @@ -7,8 +7,6 @@ ##INFO= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= @@ -17,6 +15,6 @@ 20 68799 . T C . . . GT . 20 68800 . A G . . . GT . 20 68810 . G A . . . GT . -20 69066 . C G 60 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,33 +20 69066 . C G 59.5765 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,33 20 69094 . G A . . . GT . 20 69408 . C T . . . GT . diff --git a/test/mpileup.cAls.6.out b/test/mpileup.cAls.6.out index 857d0663b..1a507325b 100644 --- a/test/mpileup.cAls.6.out +++ b/test/mpileup.cAls.6.out @@ -11,18 +11,16 @@ ##INFO= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample -1 1368828 . GT G 62 . DP=1;MQ0F=0;AC=0;AN=2;DP4=1,0,0,0;MQ=60 GT:PL 0/0:0,3,35 +1 1368828 . GT G 61.5766 . DP=1;MQ0F=0;AC=0;AN=2;DP4=1,0,0,0;MQ=60 GT:PL 0/0:0,3,35 1 1368833 . T C 39.5768 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,13 1 1368833 . TAAAAAAAAAAAAAAAA TAAAAAAAAAAAAAA 24.1741 . DP=2;AC=0;AN=2;DP4=0,1,1,0;MQ=60 GT:PL 0/0:6,0,6 1 1368833 . T G . . . GT . -1 1368834 . A T 60 . DP=1;MQ0F=0;AC=0;AN=2;DP4=1,0,0,0;MQ=60 GT:PL 0/0:0,3,33 +1 1368834 . A T 59.5765 . DP=1;MQ0F=0;AC=0;AN=2;DP4=1,0,0,0;MQ=60 GT:PL 0/0:0,3,33 16 60288 . C A . . . GT . 17 355 . G A . . . GT . 20 58799 . T C . . . GT . diff --git a/test/mpileup.cAls.7.out b/test/mpileup.cAls.7.out index 594df92b5..7168edf56 100644 --- a/test/mpileup.cAls.7.out +++ b/test/mpileup.cAls.7.out @@ -37,25 +37,23 @@ ##INFO= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT A 1 790667 . C T . . . GT . -1 116844268 . GA G 59 . DP=1;MQ0F=0;AC=0;AN=2;DP4=1,0,0,0;MQ=60 GT:PL 0/0:0,3,32 +1 116844268 . GA G 58.5765 . DP=1;MQ0F=0;AC=0;AN=2;DP4=1,0,0,0;MQ=60 GT:PL 0/0:0,3,32 1 116844279 . GAA G 39.5768 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,13 1 116844450 . T A . . . GT . 20 68799 . T C . . . GT . 20 68810 . G A . . . GT . -20 69066 . C G 60 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,33 +20 69066 . C G 59.5765 . DP=1;MQ0F=0;AC=0;AN=2;DP4=0,1,0,0;MQ=60 GT:PL 0/0:0,3,33 20 69094 . G A . . . GT . 20 69408 . C T . . . GT . 21 9411409 . T C . . . GT . 21 9411485 . C A . . . GT . 21 9411497 . A G . . . GT . -21 9412485 . C G 79 . DP=2;MQ0F=0;AC=0;AN=2;DP4=1,0,0,0;MQ=52 GT:PL 0/0:0,3,52 +21 9412485 . C G 78.5768 . DP=2;MQ0F=0;AC=0;AN=2;DP4=1,0,0,0;MQ=52 GT:PL 0/0:0,3,52 16 60288 . C A . . . GT . 17 355 . G A . . . GT . diff --git a/test/mpileup.cAls.out b/test/mpileup.cAls.out index 65979687e..b267c70f2 100644 --- a/test/mpileup.cAls.out +++ b/test/mpileup.cAls.out @@ -20,21 +20,19 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100 HG00101 HG00102 -17 1 . A G,T 52 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 -17 2 . A T,G 52 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 -17 3 . A C 26.0007 . DP=11;MQ0F=0;AC=0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 +17 1 . A G,T . . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 +17 2 . A T,G . . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 +17 3 . A C . . DP=11;MQ0F=0;AC=0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 17 4 . A G,T,C 21.815 . DP=11;MQ0F=0;AC=0,0,0;AN=2;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV 0/0:1,2,3,7,8,10,11,12,14,15:5:0 ./.:.:3:0 ./.:.:3:0 -17 5 . A G,T 26.0007 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 -17 6 . A T,G 26.0007 . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 +17 5 . A G,T . . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 +17 6 . A T,G . . DP=11;MQ0F=0;AC=0,0;AN=0;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV ./.:0,0,0,0,0,0:5:0 ./.:.:3:0 ./.:.:3:0 17 7 . A T,G,C 21.5769 . DP=11;MQ0F=0;AC=0,0,0;AN=2;DP4=11,0,0,0;MQ=29 GT:PL:DP:DV 0/0:1,2,3,4,5,6,2,3,5,3:5:0 ./.:.:3:0 ./.:.:3:0 -17 828 . T C 409 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;ICB=0.8;HOB=0.222222;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 0/1:116,0,91:9:5 1/1:120,12,0:4:4 -17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;ICB=0.128205;HOB=0.0555556;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 -17 2220 . G C 189 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=0;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/0:139,157,255:12:6 0/0:69,75,119:4:2 0/0:131,131,131:4:4 -17 2564 . A AG 166 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=0;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/0:88,98,171:6:3 0/0:57,63,117:4:2 0/0:124,124,124:4:4 +17 828 . T C 409.29 . DP=25;VDB=0.842082;SGB=-4.20907;RPB=0.950652;MQB=1;MQSB=1;BQB=0.929717;MQ0F=0;AC=4;AN=6;DP4=2,4,8,11;MQ=60 GT:PL:DP:DV 0/1:211,0,35:12:10 0/1:116,0,91:9:5 1/1:120,12,0:4:4 +17 1665 . T C 3.10665 . DP=20;VDB=0.1;SGB=0.346553;RPB=0.222222;MQB=0.611111;MQSB=0.988166;BQB=0.944444;MQ0F=0;AC=1;AN=6;DP4=7,11,1,1;MQ=55 GT:PL:DP:DV 0/0:0,21,185:7:0 0/0:0,27,222:9:0 0/1:35,0,51:4:2 +17 2220 . G C 188.992 . DP=21;VDB=0.532753;SGB=-3.51597;RPB=0.964198;MQB=0.898397;MQSB=0.875769;BQB=0.0354359;MQ0F=0;AC=0;AN=6;DP4=6,2,1,11;MQ=58 GT:PL:DP:DV 0/0:139,157,255:12:6 0/0:69,75,119:4:2 0/0:131,131,131:4:4 +17 2564 . A AG 165.992 . DP=15;VDB=0.690812;SGB=-3.20711;RPB=0.197899;MQB=1;MQSB=1;BQB=0.965069;MQ0F=0;AC=0;AN=6;DP4=1,4,4,5;MQ=60 GT:PL:DP:DV 0/0:88,98,171:6:3 0/0:57,63,117:4:2 0/0:124,124,124:4:4 diff --git a/test/mpileup.cals.8.out b/test/mpileup.cals.8.out index 169c8f78a..e07af78b9 100644 --- a/test/mpileup.cals.8.out +++ b/test/mpileup.cals.8.out @@ -21,11 +21,9 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA12878 -6 263708 . CTT C,CT 280 . DP=211;MQSB=0.901445;MQ0F=0;AC=0,0;AN=2;DP4=117,27,0,0;MQ=57 GT:PL:DP:SP:ADF:ADR:AD 0/0:0,255,255,255,255,255:144:0:117,0,0:27,0,0:144,0,0 +6 263708 . CTT C,CT 279.818 . DP=211;MQSB=0.901445;MQ0F=0;AC=0,0;AN=2;DP4=117,27,0,0;MQ=57 GT:PL:DP:SP:ADF:ADR:AD 0/0:0,255,255,255,255,255:144:0:117,0,0:27,0,0:144,0,0 diff --git a/test/mpileup.cals.9.out b/test/mpileup.cals.9.out index 35871e794..da7d723a9 100644 --- a/test/mpileup.cals.9.out +++ b/test/mpileup.cals.9.out @@ -11,8 +11,6 @@ ##INFO= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= diff --git a/test/mpileup.hwe.1.out b/test/mpileup.hwe.1.out index 684a18e65..2bc7803d3 100644 --- a/test/mpileup.hwe.1.out +++ b/test/mpileup.hwe.1.out @@ -19,9 +19,7 @@ ##INFO= ##INFO= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA19213 NA19129 -15 201 . A G 202 . DP=80;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;AC=4;AN=4;DP4=19,17,21,15;MQ=59 GT:PL:AD 1/1:144,12,0:0,4 1/1:85,9,0:0,3 +15 201 . A G 202.049 . DP=80;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;AC=4;AN=4;DP4=19,17,21,15;MQ=59 GT:PL:AD 1/1:144,12,0:0,4 1/1:85,9,0:0,3 diff --git a/test/mpileup.hwe.1b.out b/test/mpileup.hwe.1b.out new file mode 100644 index 000000000..56032579d --- /dev/null +++ b/test/mpileup.hwe.1b.out @@ -0,0 +1,25 @@ +##fileformat=VCFv4.2 +##FILTER= +##reference=file://hwe.fa +##contig= +##ALT= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##FORMAT= +##FORMAT= +##INFO= +##INFO= +##FORMAT= +##INFO= +##INFO= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA19213 NA19129 +15 201 . A G 117.048 . DP=80;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;AC=4;AN=4;DP4=19,17,21,15;MQ=59 GT:PL:AD 1/1:144,12,0:0,4 1/1:85,9,0:0,3 diff --git a/test/mpileup.hwe.2.out b/test/mpileup.hwe.2.out index d3cc90753..9b3d3e2ef 100644 --- a/test/mpileup.hwe.2.out +++ b/test/mpileup.hwe.2.out @@ -17,11 +17,9 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA07000 NA07056 NA12046 NA12144 NA12156 NA12234 NA12249 NA12282 NA12283 NA12286 NA12340 NA12341 NA12342 NA12399 NA12414 NA12546 NA12717 NA12718 NA12748 NA12749 NA12750 NA12776 NA12828 NA12843 NA18502 NA18867 NA18908 NA18933 NA18934 NA19108 NA19117 NA19119 NA19129 NA19131 NA19138 NA19147 NA19149 NA19159 NA19185 NA19189 NA19197 NA19198 NA19206 NA19213 NA19214 NA19223 NA19235 NA19236 NA19238 NA19239 NA19248 NA19256 -15 201 . A G 999 . DP=80;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;ICB=0.0721108;HOB=0.158099;AC=69;AN=104;DP4=19,17,21,15;MQ=59 GT:PL:AD 0/1:0,3,36:1,0 0/0:0,6,93:2,0 0/1:0,3,33:1,0 0/1:0,3,32:1,0 0/1:0,3,36:1,0 0/1:0,3,32:1,0 0/1:0,3,40:1,0 0/0:0,6,79:2,0 0/1:0,3,42:1,0 0/0:0,6,64:2,0 0/0:0,6,67:2,0 0/0:0,6,66:2,0 0/0:0,9,96:3,0 0/1:0,3,32:1,0 0/1:0,3,36:1,0 0/1:0,3,31:1,0 0/1:0,3,44:1,0 0/1:0,3,44:1,0 0/0:0,6,61:2,0 0/1:0,3,60:1,0 0/0:0,6,60:2,0 0/0:0,6,57:2,0 0/1:0,3,31:1,0 0/0:0,6,68:2,0 0/1:0,3,60:1,0 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:37,3,0:0,1 1/1:60,3,0:0,1 1/1:38,3,0:0,1 1/1:36,3,0:0,1 1/1:37,3,0:0,1 1/1:85,9,0:0,3 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:63,6,0:0,2 1/1:29,3,0:0,1 1/1:32,3,0:0,1 1/1:43,3,0:0,1 1/1:27,3,0:0,1 1/1:43,3,0:0,1 1/1:34,3,0:0,1 1/1:46,3,0:0,1 1/1:144,12,0:0,4 1/1:59,3,0:0,1 1/1:16,3,0:0,1 1/1:37,3,0:0,1 1/1:57,6,0:0,2 1/1:91,6,0:0,2 1/1:21,3,0:0,1 1/1:52,6,0:0,2 1/1:26,3,0:0,1 +15 201 . A G 1051.64 . DP=80;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;AC=69;AN=104;DP4=19,17,21,15;MQ=59 GT:PL:AD 0/1:0,3,36:1,0 0/0:0,6,93:2,0 0/1:0,3,33:1,0 0/1:0,3,32:1,0 0/1:0,3,36:1,0 0/1:0,3,32:1,0 0/1:0,3,40:1,0 0/0:0,6,79:2,0 0/1:0,3,42:1,0 0/0:0,6,64:2,0 0/0:0,6,67:2,0 0/0:0,6,66:2,0 0/0:0,9,96:3,0 0/1:0,3,32:1,0 0/1:0,3,36:1,0 0/1:0,3,31:1,0 0/1:0,3,44:1,0 0/1:0,3,44:1,0 0/0:0,6,61:2,0 0/1:0,3,60:1,0 0/0:0,6,60:2,0 0/0:0,6,57:2,0 0/1:0,3,31:1,0 0/0:0,6,68:2,0 0/1:0,3,60:1,0 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:37,3,0:0,1 1/1:60,3,0:0,1 1/1:38,3,0:0,1 1/1:36,3,0:0,1 1/1:37,3,0:0,1 1/1:85,9,0:0,3 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:63,6,0:0,2 1/1:29,3,0:0,1 1/1:32,3,0:0,1 1/1:43,3,0:0,1 1/1:27,3,0:0,1 1/1:43,3,0:0,1 1/1:34,3,0:0,1 1/1:46,3,0:0,1 1/1:144,12,0:0,4 1/1:59,3,0:0,1 1/1:16,3,0:0,1 1/1:37,3,0:0,1 1/1:57,6,0:0,2 1/1:91,6,0:0,2 1/1:21,3,0:0,1 1/1:52,6,0:0,2 1/1:26,3,0:0,1 diff --git a/test/mpileup.hwe.3.out b/test/mpileup.hwe.3.out index 5645f8479..76166969e 100644 --- a/test/mpileup.hwe.3.out +++ b/test/mpileup.hwe.3.out @@ -17,12 +17,10 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA07000 NA07056 NA12046 NA12144 NA12156 NA12234 NA12249 NA12282 NA12283 NA12286 NA12340 NA12341 NA12342 NA12399 NA12414 NA12546 NA12717 NA12718 NA12748 NA12749 NA12750 NA12776 NA12828 NA12843 NA18502 NA18867 NA18908 NA18933 NA18934 NA19108 NA19117 NA19119 NA19129 NA19131 NA19138 NA19147 NA19149 NA19159 NA19185 NA19189 NA19197 NA19198 NA19206 NA19213 NA19214 NA19223 NA19235 NA19236 NA19238 NA19239 NA19248 NA19256 -15 198 . C A 12.8211 . DP=79;SGB=-0.627953;RPB=1;MQB=1;MQSB=0.974642;BQB=1;MQ0F=0;ICB=0.490573;HOB=0.0377219;AC=2;AN=104;DP4=40,31,0,1;MQ=59 GT:PL:AD 0/0:0,3,44:1,0 0/0:0,6,95:2,0 0/0:0,3,39:1,0 0/0:0,3,43:1,0 0/0:0,3,34:1,0 0/0:0,3,37:1,0 0/0:0,3,41:1,0 0/0:0,6,87:2,0 0/0:0,3,46:1,0 0/0:0,6,76:2,0 0/0:0,6,80:2,0 0/0:0,6,78:2,0 0/0:0,9,112:3,0 0/0:0,3,41:1,0 1/1:35,3,0:0,1 0/0:0,3,33:1,0 0/0:0,3,44:1,0 0/0:0,3,45:1,0 0/0:0,6,69:2,0 0/0:0,3,60:1,0 0/0:0,6,70:2,0 0/0:0,6,57:2,0 0/0:0,3,37:1,0 0/0:0,6,79:2,0 0/0:0,3,60:1,0 0/0:0,3,42:1,0 0/0:0,3,39:1,0 0/0:0,3,44:1,0 0/0:0,3,60:1,0 0/0:0,3,39:1,0 0/0:0,3,41:1,0 0/0:0,3,44:1,0 0/0:0,9,98:3,0 0/0:0,3,40:1,0 0/0:0,3,38:1,0 0/0:0,6,72:2,0 0/0:0,3,29:1,0 0/0:0,3,34:1,0 0/0:0,3,43:1,0 0/0:0,3,34:1,0 0/0:0,3,44:1,0 0/0:0,3,36:1,0 0/0:0,3,47:1,0 0/0:0,12,163:4,0 0/0:0,3,60:1,0 0/0:0,3,41:1,0 0/0:0,3,37:1,0 0/0:0,6,73:2,0 0/0:0,6,96:2,0 0/0:0,3,40:1,0 0/0:0,6,67:2,0 0/0:0,3,33:1,0 -15 201 . A G 999 . DP=80;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;ICB=5.51698e-12;HOB=0.49926;AC=54;AN=104;DP4=19,17,21,15;MQ=59 GT:PL:AD 0/0:0,3,36:1,0 0/0:0,6,93:2,0 0/0:0,3,33:1,0 0/0:0,3,32:1,0 0/0:0,3,36:1,0 0/0:0,3,32:1,0 0/0:0,3,40:1,0 0/0:0,6,79:2,0 0/0:0,3,42:1,0 0/0:0,6,64:2,0 0/0:0,6,67:2,0 0/0:0,6,66:2,0 0/0:0,9,96:3,0 0/0:0,3,32:1,0 0/0:0,3,36:1,0 0/0:0,3,31:1,0 0/0:0,3,44:1,0 0/0:0,3,44:1,0 0/0:0,6,61:2,0 0/0:0,3,60:1,0 0/0:0,6,60:2,0 0/0:0,6,57:2,0 0/0:0,3,31:1,0 0/0:0,6,68:2,0 0/0:0,3,60:1,0 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:37,3,0:0,1 1/1:60,3,0:0,1 1/1:38,3,0:0,1 1/1:36,3,0:0,1 1/1:37,3,0:0,1 1/1:85,9,0:0,3 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:63,6,0:0,2 1/1:29,3,0:0,1 1/1:32,3,0:0,1 1/1:43,3,0:0,1 1/1:27,3,0:0,1 1/1:43,3,0:0,1 1/1:34,3,0:0,1 1/1:46,3,0:0,1 1/1:144,12,0:0,4 1/1:59,3,0:0,1 1/1:16,3,0:0,1 1/1:37,3,0:0,1 1/1:57,6,0:0,2 1/1:91,6,0:0,2 1/1:21,3,0:0,1 1/1:52,6,0:0,2 1/1:26,3,0:0,1 +15 198 . C A 12.8224 . DP=79;SGB=-0.627953;RPB=1;MQB=1;MQSB=0.974642;BQB=1;MQ0F=0;AC=2;AN=104;DP4=40,31,0,1;MQ=59 GT:PL:AD 0/0:0,3,44:1,0 0/0:0,6,95:2,0 0/0:0,3,39:1,0 0/0:0,3,43:1,0 0/0:0,3,34:1,0 0/0:0,3,37:1,0 0/0:0,3,41:1,0 0/0:0,6,87:2,0 0/0:0,3,46:1,0 0/0:0,6,76:2,0 0/0:0,6,80:2,0 0/0:0,6,78:2,0 0/0:0,9,112:3,0 0/0:0,3,41:1,0 1/1:35,3,0:0,1 0/0:0,3,33:1,0 0/0:0,3,44:1,0 0/0:0,3,45:1,0 0/0:0,6,69:2,0 0/0:0,3,60:1,0 0/0:0,6,70:2,0 0/0:0,6,57:2,0 0/0:0,3,37:1,0 0/0:0,6,79:2,0 0/0:0,3,60:1,0 0/0:0,3,42:1,0 0/0:0,3,39:1,0 0/0:0,3,44:1,0 0/0:0,3,60:1,0 0/0:0,3,39:1,0 0/0:0,3,41:1,0 0/0:0,3,44:1,0 0/0:0,9,98:3,0 0/0:0,3,40:1,0 0/0:0,3,38:1,0 0/0:0,6,72:2,0 0/0:0,3,29:1,0 0/0:0,3,34:1,0 0/0:0,3,43:1,0 0/0:0,3,34:1,0 0/0:0,3,44:1,0 0/0:0,3,36:1,0 0/0:0,3,47:1,0 0/0:0,12,163:4,0 0/0:0,3,60:1,0 0/0:0,3,41:1,0 0/0:0,3,37:1,0 0/0:0,6,73:2,0 0/0:0,6,96:2,0 0/0:0,3,40:1,0 0/0:0,6,67:2,0 0/0:0,3,33:1,0 +15 201 . A G 121.59 . DP=80;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;AC=50;AN=104;DP4=19,17,21,15;MQ=59 GT:PL:AD 0/0:0,3,36:1,0 0/0:0,6,93:2,0 0/0:0,3,33:1,0 0/0:0,3,32:1,0 0/0:0,3,36:1,0 0/0:0,3,32:1,0 0/0:0,3,40:1,0 0/0:0,6,79:2,0 0/0:0,3,42:1,0 0/0:0,6,64:2,0 0/0:0,6,67:2,0 0/0:0,6,66:2,0 0/0:0,9,96:3,0 0/0:0,3,32:1,0 0/0:0,3,36:1,0 0/0:0,3,31:1,0 0/0:0,3,44:1,0 0/0:0,3,44:1,0 0/0:0,6,61:2,0 0/0:0,3,60:1,0 0/0:0,6,60:2,0 0/0:0,6,57:2,0 0/0:0,3,31:1,0 0/0:0,6,68:2,0 0/0:0,3,60:1,0 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:37,3,0:0,1 1/1:60,3,0:0,1 1/1:38,3,0:0,1 1/1:36,3,0:0,1 1/1:37,3,0:0,1 1/1:85,9,0:0,3 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:63,6,0:0,2 1/1:29,3,0:0,1 1/1:32,3,0:0,1 1/1:43,3,0:0,1 1/1:27,3,0:0,1 1/1:43,3,0:0,1 1/1:34,3,0:0,1 1/1:46,3,0:0,1 1/1:144,12,0:0,4 1/1:59,3,0:0,1 0/0:16,3,0:0,1 1/1:37,3,0:0,1 1/1:57,6,0:0,2 1/1:91,6,0:0,2 0/0:21,3,0:0,1 1/1:52,6,0:0,2 1/1:26,3,0:0,1 diff --git a/test/mpileup.hwe.4.out b/test/mpileup.hwe.4.out index b1509eabc..b84853c81 100644 --- a/test/mpileup.hwe.4.out +++ b/test/mpileup.hwe.4.out @@ -17,11 +17,9 @@ ##FORMAT= ##FORMAT= ##FORMAT= -##INFO= -##INFO= ##INFO= ##INFO= ##INFO= ##INFO= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA07000 NA07056 NA12046 NA12144 NA12156 NA12234 NA12249 NA12282 NA12283 NA12286 NA12340 NA12341 NA12342 NA12399 NA12414 NA12546 NA12717 NA12718 NA12748 NA12749 NA12750 NA12776 NA12828 NA12843 NA18502 NA18867 NA18908 NA18933 NA18934 NA19108 NA19117 NA19119 NA19129 NA19131 NA19138 NA19147 NA19149 NA19159 NA19185 NA19189 NA19197 NA19198 NA19206 NA19213 NA19214 NA19223 NA19235 NA19236 NA19238 NA19239 NA19248 NA19256 -15 201 . A G 999 . DP=80;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;ICB=4.28809e-11;HOB=0.479105;AC=55;AN=104;DP4=19,17,21,15;MQ=59 GT:PL:AD 0/0:0,3,36:1,0 0/0:0,6,93:2,0 0/0:0,3,33:1,0 0/0:0,3,32:1,0 0/0:0,3,36:1,0 0/0:0,3,32:1,0 0/0:0,3,40:1,0 0/0:0,6,79:2,0 0/0:0,3,42:1,0 0/0:0,6,64:2,0 0/0:0,6,67:2,0 0/0:0,6,66:2,0 0/0:0,9,96:3,0 0/0:0,3,32:1,0 0/0:0,3,36:1,0 0/0:0,3,31:1,0 0/0:0,3,44:1,0 0/0:0,3,44:1,0 0/0:0,6,61:2,0 0/0:0,3,60:1,0 0/0:0,6,60:2,0 0/0:0,6,57:2,0 0/0:0,3,31:1,0 0/0:0,6,68:2,0 0/1:0,3,60:1,0 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:37,3,0:0,1 1/1:60,3,0:0,1 1/1:38,3,0:0,1 1/1:36,3,0:0,1 1/1:37,3,0:0,1 1/1:85,9,0:0,3 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:63,6,0:0,2 1/1:29,3,0:0,1 1/1:32,3,0:0,1 1/1:43,3,0:0,1 1/1:27,3,0:0,1 1/1:43,3,0:0,1 1/1:34,3,0:0,1 1/1:46,3,0:0,1 1/1:144,12,0:0,4 1/1:59,3,0:0,1 1/1:16,3,0:0,1 1/1:37,3,0:0,1 1/1:57,6,0:0,2 1/1:91,6,0:0,2 1/1:21,3,0:0,1 1/1:52,6,0:0,2 1/1:26,3,0:0,1 +15 201 . A G 1211.32 . DP=80;VDB=0.895839;SGB=-4.94122;RPB=0.999937;MQB=0.979662;MQSB=0.974642;BQB=0.739382;MQ0F=0;AC=55;AN=104;DP4=19,17,21,15;MQ=59 GT:PL:AD 0/0:0,3,36:1,0 0/0:0,6,93:2,0 0/0:0,3,33:1,0 0/0:0,3,32:1,0 0/0:0,3,36:1,0 0/0:0,3,32:1,0 0/0:0,3,40:1,0 0/0:0,6,79:2,0 0/0:0,3,42:1,0 0/0:0,6,64:2,0 0/0:0,6,67:2,0 0/0:0,6,66:2,0 0/0:0,9,96:3,0 0/0:0,3,32:1,0 0/0:0,3,36:1,0 0/0:0,3,31:1,0 0/0:0,3,44:1,0 0/0:0,3,44:1,0 0/0:0,6,61:2,0 0/0:0,3,60:1,0 0/0:0,6,60:2,0 0/0:0,6,57:2,0 0/0:0,3,31:1,0 0/0:0,6,68:2,0 0/1:0,3,60:1,0 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:37,3,0:0,1 1/1:60,3,0:0,1 1/1:38,3,0:0,1 1/1:36,3,0:0,1 1/1:37,3,0:0,1 1/1:85,9,0:0,3 1/1:39,3,0:0,1 1/1:31,3,0:0,1 1/1:63,6,0:0,2 1/1:29,3,0:0,1 1/1:32,3,0:0,1 1/1:43,3,0:0,1 1/1:27,3,0:0,1 1/1:43,3,0:0,1 1/1:34,3,0:0,1 1/1:46,3,0:0,1 1/1:144,12,0:0,4 1/1:59,3,0:0,1 1/1:16,3,0:0,1 1/1:37,3,0:0,1 1/1:57,6,0:0,2 1/1:91,6,0:0,2 1/1:21,3,0:0,1 1/1:52,6,0:0,2 1/1:26,3,0:0,1 diff --git a/test/test.pl b/test/test.pl index 8fade275b..633d1d5db 100755 --- a/test/test.pl +++ b/test/test.pl @@ -278,10 +278,10 @@ test_vcf_call($opts,in=>'mpileup.X',out=>'mpileup.X.out',args=>'-mv --ploidy-file {PATH}/mpileup.ploidy -S {PATH}/mpileup.ped'); test_vcf_call($opts,in=>'mpileup.X',out=>'mpileup.X.2.out',args=>'-mv --ploidy-file {PATH}/mpileup.ploidy -S {PATH}/mpileup.2.samples'); test_vcf_call($opts,in=>'mpileup.NA19213.NA19129',out=>'mpileup.hwe.1.out',args=>'-mv'); -test_vcf_call($opts,in=>'mpileup.NA19213.NA19129',out=>'mpileup.hwe.1.out',args=>'-mv -G -'); +test_vcf_call($opts,in=>'mpileup.NA19213.NA19129',out=>'mpileup.hwe.1b.out',args=>'-mv -G AD:-'); test_vcf_call($opts,in=>'mpileup.hwe',out=>'mpileup.hwe.2.out',args=>'-mv'); -test_vcf_call($opts,in=>'mpileup.hwe',out=>'mpileup.hwe.3.out',args=>'-mv -G -'); -test_vcf_call($opts,in=>'mpileup.hwe',out=>'mpileup.hwe.4.out',args=>'-mv -G {PATH}/mpileup.hwe.samples'); +test_vcf_call($opts,in=>'mpileup.hwe',out=>'mpileup.hwe.3.out',args=>'-mv -G AD:-'); # 21,3,0 becomes 0/0 because of the prior -P +test_vcf_call($opts,in=>'mpileup.hwe',out=>'mpileup.hwe.4.out',args=>'-mv -G AD:{PATH}/mpileup.hwe.samples'); test_vcf_call_cAls($opts,in=>'mpileup',out=>'mpileup.cAls.out',tab=>'mpileup'); test_vcf_call_cAls($opts,in=>'mpileup.2',out=>'mpileup.cAls.2.out',tab=>'mpileup.2'); test_vcf_call_cAls($opts,in=>'mpileup.3',out=>'mpileup.cAls.3.out',tab=>'mpileup.3',args=>'-i'); @@ -296,6 +296,12 @@ test_vcf_call($opts,in=>'mpileup.c.X',out=>'mpileup.c.X.out',args=>'-cv --ploidy-file {PATH}/mpileup.ploidy -S {PATH}/mpileup.samples'); test_vcf_call($opts,in=>'mpileup.c.X',out=>'mpileup.c.X.out',args=>'-cv --ploidy-file {PATH}/mpileup.ploidy -S {PATH}/mpileup.ped'); test_vcf_call($opts,in=>'mpileup.c.X',out=>'mpileup.c.X.2.out',args=>'-cv --ploidy-file {PATH}/mpileup.ploidy -S {PATH}/mpileup.2.samples'); +test_vcf_call($opts,in=>'call-G',out=>'call-G.1.out',args=>'-mv'); +test_vcf_call($opts,in=>'call-G',out=>'call-G.2.out',args=>'-mv -G AD:-'); +test_vcf_call($opts,in=>'call-G.2',out=>'call-G.2.1.out',args=>'-mv -F AN_POP,AC_POP'); +test_vcf_call($opts,in=>'call.af-fixation',out=>'call.af-fixation.1.out',args=>'-m'); +test_vcf_call($opts,in=>'call.af-fixation',out=>'call.af-fixation.2.out',args=>'-m -G {PATH}/call.af-fixation.txt'); +test_vcf_call($opts,in=>'call.af-fixation',out=>'call.af-fixation.3.out',args=>'-m -G {PATH}/call.af-fixation.txt -a GP,GQ'); test_vcf_filter($opts,in=>'view.filter',out=>'view.filter.6.out',args=>q[-S. -e'TXT0="text"'],reg=>''); test_vcf_filter($opts,in=>'view.filter',out=>'view.filter.7.out',args=>q[-S. -e'FMT/FRS[*:1]="BB"'],reg=>''); test_vcf_filter($opts,in=>'view.filter',out=>'view.filter.8.out',args=>q[-S. -e'FMT/FGS[*:0]="AAAAAA"'],reg=>''); diff --git a/test/trio-dnm.1.vcf b/test/trio-dnm.1.vcf deleted file mode 100644 index 7d15610dd..000000000 --- a/test/trio-dnm.1.vcf +++ /dev/null @@ -1,31 +0,0 @@ -##fileformat=VCFv4.2 -##FILTER= -##FORMAT= -##FORMAT= -##FORMAT= -##INFO= -##INFO= -##INFO= -##contig= -##contig= -##contig= -##reference=file:///lustre/scratch113/resources/ref/Homo_sapiens/1000Genomes_hs37d5/hs37d5.fa -#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT proband father mother -1 1 . G A,T . . TP GT:AD:PL 0/2:62,0,27:733,919,3060,0,2141,2060 0/0:35,0,0:0,99,1485,99,1485,1485 0/0:38,0,0:0,102,1530,102,1530,1530 -1 2 . G T,A . . TP GT:AD:PL 0/1:67,14,0:257,0,2105,458,2146,2605 0/0:36,0,0:0,99,1485,99,1485,1485 0/0:36,0,0:0,99,1485,99,1485,1485 -1 3 . G A . . TP GT:AD:PL 0/1:71,14:241,0,2314 0/0:45,0:0,101,1530 0/0:43,0:0,99,1485 -1 4 . C T . . TP GT:AD:PL 0/1:111,24:504,0,3776 0/0:35,0:0,99,1485 0/0:36,0:0,99,1485 -1 5 . C A . . TP GT:AD:PL 0/1:30,6:124,0,981 0/0:37,0:0,99,1485 0/0:32,0:0,90,1350 -1 8 . A G . . TP GT:AD:PL 0/1:434,52:859,0,18086 0/0:38,0:0,99,1485 0/0:38,0:0,99,1485 -2 1 . A G . . UN GT:AD:PL 0/1:2,5:179,0,55 0/0:7,0:0,0,180 0/0:4,0:0,12,148 -2 2 . A G . . UN GT:AD:PL 0/1:4,5:159,0,126 0/0:1,0:0,3,39 0/0:4,0:0,9,135 -2 3 . A G . . UN GT:AD:PL 0/1:4,4:137,0,107 0/0:6,0:0,18,213 0/0:8,0:0,0,232 -3 1 . A G . . FP GT:AD:PL 0/1:7,9:357,0,408 0/0:15,0:0,39,585 0/1:4,3:114,0,550 -3 2 . C A . . FP GT:AD:PL 0/1:13,15:453,0,442 0/1:29,30:913,0,1011 0/0:39,0:0,99,1485 -3 3 . A G . . FP GT:AD:PL 0/1:11,12:361,0,358 0/0:21,0:0,51,765 0/1:10,15:538,0,292 -3 4 . A G,C . . FP GT:PL:AD 0/0:0,255,255,255,255,255:306,11,0 0/0:0,255,255,255,255,255:328,1,1 0/0:0,255,255,255,255,255:318,0,0 -3 5 . A G . . FP GT:AD:PL 0/1:33,32:890,0,963 0/1:56,45:1328,0,1809 0/0:36,0:0,99,1485 -3 6 . A G . . FP GT:AD:PL 0/1:19,24:737,0,649 0/0:48,0:0,108,1620 0/1:25,22:644,0,836 -3 7 . A G . . FP GT:AD:PL 0/1:73,90:2864,0,2197 0/0:42,0:0,99,1485 0/1:69,74:2395,0,2064 -3 8 . A G . . FP GT:AD:PL 0/1:115,128:4130,0,3542 0/0:34,0:0,99,1360 0/1:137,89:2571,0,4411 -3 9 . A G . . FP GT:AD:PL 0/1:18,11:311,0,627 0/1:3,3:51,0,105 0/0:19,0:0,57,764 diff --git a/test/trio-dnm.2.vcf b/test/trio-dnm.2.vcf deleted file mode 100644 index a36f5d0ff..000000000 --- a/test/trio-dnm.2.vcf +++ /dev/null @@ -1,31 +0,0 @@ -##fileformat=VCFv4.2 -##FILTER= -##FORMAT= -##FORMAT= -##FORMAT= -##INFO= -##INFO= -##INFO= -##contig= -##contig= -##contig= -##reference=file:///lustre/scratch113/resources/ref/Homo_sapiens/1000Genomes_hs37d5/hs37d5.fa -#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT proband father mother -1 1 . G A,T . . TP GT:AD:PL 0/2:62,0,27,0:733,919,3060,0,2141,2060 0/0:35,0,0:0,99,1485,99,1485,1485 0/0:38,0,0:0,102,1530,102,1530,1530 -1 2 . G T,A . . TP GT:AD:PL 0/1:67,14,0,0:257,0,2105,458,2146,2605 0/0:36,0,0:0,99,1485,99,1485,1485 0/0:36,0,0:0,99,1485,99,1485,1485 -1 3 . G A . . TP GT:AD:PL 0/1:71,14,0:241,0,2314 0/0:45,0:0,101,1530 0/0:43,0:0,99,1485 -1 4 . C T . . TP GT:AD:PL 0/1:111,24,0:504,0,3776 0/0:35,0:0,99,1485 0/0:36,0:0,99,1485 -1 5 . C A . . TP GT:AD:PL 0/1:30,6,0:124,0,981 0/0:37,0:0,99,1485 0/0:32,0:0,90,1350 -1 8 . A G . . TP GT:AD:PL 0/1:434,52,0:859,0,18086 0/0:38,0:0,99,1485 0/0:38,0:0,99,1485 -2 1 . A G . . UN GT:AD:PL 0/1:2,5,0:179,0,55 0/0:7,0:0,0,180 0/0:4,0:0,12,148 -2 2 . A G . . UN GT:AD:PL 0/1:4,5,0:159,0,126 0/0:1,0:0,3,39 0/0:4,0:0,9,135 -2 3 . A G . . UN GT:AD:PL 0/1:4,4,0:137,0,107 0/0:6,0:0,18,213 0/0:8,0:0,0,232 -3 1 . A G . . FP GT:AD:PL 0/1:7,9,0,0:357,0,408 0/0:15,0:0,39,585 0/1:4,3:114,0,550 -3 2 . C A . . FP GT:AD:PL 0/1:13,15,0:453,0,442 0/1:29,30:913,0,1011 0/0:39,0:0,99,1485 -3 3 . A G . . FP GT:AD:PL 0/1:11,12,0:361,0,358 0/0:21,0:0,51,765 0/1:10,15:538,0,292 -3 4 . A G,C . . FP GT:PL:AD 0/0:0,255,255,255,255,255:306,11,0,0 0/0:0,255,255,255,255,255:328,1,1 0/0:0,255,255,255,255,255:318,0,0 -3 5 . A G . . FP GT:AD:PL 0/1:33,32,0:890,0,963 0/1:56,45:1328,0,1809 0/0:36,0:0,99,1485 -3 6 . A G . . FP GT:AD:PL 0/1:19,24,0:737,0,649 0/0:48,0:0,108,1620 0/1:25,22:644,0,836 -3 7 . A G . . FP GT:AD:PL 0/1:73,90,0:2864,0,2197 0/0:42,0:0,99,1485 0/1:69,74:2395,0,2064 -3 8 . A G . . FP GT:AD:PL 0/1:115,128,0:4130,0,3542 0/0:34,0:0,99,1360 0/1:137,89:2571,0,4411 -3 9 . A G . . FP GT:AD:PL 0/1:18,11,0:311,0,627 0/1:3,3:51,0,105 0/0:19,0:0,57,764 diff --git a/vcfcall.c b/vcfcall.c index 72387f105..0bbcadbec 100644 --- a/vcfcall.c +++ b/vcfcall.c @@ -740,6 +740,7 @@ static void destroy_data(args_t *args) free(args->samples_map); free(args->sample2sex); free(args->aux.ploidy); + free(args->aux.sample_groups); free(args->str.s); if ( args->gvcf ) gvcf_destroy(args->gvcf); bcf_hdr_destroy(args->aux.hdr); @@ -769,7 +770,20 @@ void parse_novel_rate(args_t *args, const char *str) else error("Could not parse --novel-rate %s\n", str); } -static int parse_format_flag(const char *str) +static void list_annotations(FILE *fp) +{ + fprintf(fp, + "\n" + "Optional INFO annotations available with -m (\"INFO/\" prefix is optional):\n" + " INFO/PV4 .. P-values for strand bias, baseQ bias, mapQ bias and tail distance bias (Number=4,Type=Float)\n" + "\n" + "Optional FORMAT annotations available with -m (\"FORMAT/\" prefix is optional):\n" + " FORMAT/GQ .. Phred-scaled genotype quality (Number=1,Type=Integer)\n" + " FORMAT/GP .. Phred-scaled genotype posterior probabilities (Number=G,Type=Float)\n" + "\n"); +} + +static int parse_output_tags(const char *str) { int flag = 0; const char *ss = str; @@ -777,8 +791,9 @@ static int parse_format_flag(const char *str) { const char *se = ss; while ( *se && *se!=',' ) se++; - if ( !strncasecmp(ss,"GQ",se-ss) ) flag |= CALL_FMT_GQ; - else if ( !strncasecmp(ss,"GP",se-ss) ) flag |= CALL_FMT_GP; + if ( !strncasecmp(ss,"GQ",se-ss) || !strncasecmp(ss,"FORMAT/GQ",se-ss) || !strncasecmp(ss,"FMT/GQ",se-ss) ) flag |= CALL_FMT_GQ; + else if ( !strncasecmp(ss,"GP",se-ss) || !strncasecmp(ss,"FORMAT/GP",se-ss) || !strncasecmp(ss,"FMT/GP",se-ss) ) flag |= CALL_FMT_GP; + else if ( !strncasecmp(ss,"PV4",se-ss) || !strncasecmp(ss,"INFO/PV4",se-ss) ) flag |= CALL_FMT_PV4; else { fprintf(stderr,"Could not parse \"%s\"\n", str); @@ -857,37 +872,41 @@ static void usage(args_t *args) fprintf(stderr, "Usage: bcftools call [options] \n"); fprintf(stderr, "\n"); fprintf(stderr, "File format options:\n"); - fprintf(stderr, " --no-version do not append version and command line to the header\n"); - fprintf(stderr, " -o, --output write output to a file [standard output]\n"); - fprintf(stderr, " -O, --output-type output type: 'b' compressed BCF; 'u' uncompressed BCF; 'z' compressed VCF; 'v' uncompressed VCF [v]\n"); - fprintf(stderr, " --ploidy [?] predefined ploidy, 'list' to print available settings, append '?' for details\n"); - fprintf(stderr, " --ploidy-file space/tab-delimited list of CHROM,FROM,TO,SEX,PLOIDY\n"); - fprintf(stderr, " -r, --regions restrict to comma-separated list of regions\n"); - fprintf(stderr, " -R, --regions-file restrict to regions listed in a file\n"); - fprintf(stderr, " -s, --samples list of samples to include [all samples]\n"); - fprintf(stderr, " -S, --samples-file PED file or a file with an optional column with sex (see man page for details) [all samples]\n"); - fprintf(stderr, " -t, --targets similar to -r but streams rather than index-jumps\n"); - fprintf(stderr, " -T, --targets-file similar to -R but streams rather than index-jumps\n"); - fprintf(stderr, " --threads use multithreading with worker threads [0]\n"); + fprintf(stderr, " --no-version Do not append version and command line to the header\n"); + fprintf(stderr, " -o, --output FILE Write output to a file [standard output]\n"); + fprintf(stderr, " -O, --output-type b|u|z|v Output type: 'b' compressed BCF; 'u' uncompressed BCF; 'z' compressed VCF; 'v' uncompressed VCF [v]\n"); + fprintf(stderr, " --ploidy ASSEMBLY[?] Predefined ploidy, 'list' to print available settings, append '?' for details\n"); + fprintf(stderr, " --ploidy-file FILE Space/tab-delimited list of CHROM,FROM,TO,SEX,PLOIDY\n"); + fprintf(stderr, " -r, --regions REGION Restrict to comma-separated list of regions\n"); + fprintf(stderr, " -R, --regions-file FILE Restrict to regions listed in a file\n"); + fprintf(stderr, " -s, --samples LIST List of samples to include [all samples]\n"); + fprintf(stderr, " -S, --samples-file FILE PED file or a file with an optional column with sex (see man page for details) [all samples]\n"); + fprintf(stderr, " -t, --targets REGION Similar to -r but streams rather than index-jumps\n"); + fprintf(stderr, " -T, --targets-file FILE Similar to -R but streams rather than index-jumps\n"); + fprintf(stderr, " --threads INT Use multithreading with INT worker threads [0]\n"); fprintf(stderr, "\n"); fprintf(stderr, "Input/output options:\n"); - fprintf(stderr, " -A, --keep-alts keep all possible alternate alleles at variant sites\n"); - fprintf(stderr, " -f, --format-fields output format fields: GQ,GP (lowercase allowed) []\n"); - fprintf(stderr, " -F, --prior-freqs use prior allele frequencies\n"); - fprintf(stderr, " -G, --group-samples group samples by population (file with \"sample\\tgroup\") or \"-\" for single-sample calling\n"); - fprintf(stderr, " -g, --gvcf ,[...] group non-variant sites into gVCF blocks by minimum per-sample DP\n"); - fprintf(stderr, " -i, --insert-missed output also sites missed by mpileup but present in -T\n"); - fprintf(stderr, " -M, --keep-masked-ref keep sites with masked reference allele (REF=N)\n"); - fprintf(stderr, " -V, --skip-variants skip indels/snps\n"); - fprintf(stderr, " -v, --variants-only output variant sites only\n"); + fprintf(stderr, " -A, --keep-alts Keep all possible alternate alleles at variant sites\n"); + fprintf(stderr, " -a, --annotate LIST Optional tags to output (lowercase allowed); '?' to list available tags\n"); +//todo? +// fprintf(stderr, " -a, --annots LIST Add annotations: GQ,GP,PV4 (lowercase allowed). Prefixed with ^ indicates a request for\n"); +// fprintf(stderr, " tag removal [^I16,^QS,^FMT/QS]\n"); + fprintf(stderr, " -F, --prior-freqs AN,AC Use prior allele frequencies, determined from these pre-filled tags\n"); + fprintf(stderr, " -G, --group-samples [TAG:]FILE|- Group samples by population (file with \"sample\\tgroup\") or \"-\" for single-sample calling.\n"); + fprintf(stderr, " This requires FORMAT/QS or other Number=R,Type=Integer tag such as FORMAT/AD\n"); + fprintf(stderr, " -g, --gvcf INT,[...] Group non-variant sites into gVCF blocks by minimum per-sample DP\n"); + fprintf(stderr, " -i, --insert-missed Output also sites missed by mpileup but present in -T\n"); + fprintf(stderr, " -M, --keep-masked-ref Keep sites with masked reference allele (REF=N)\n"); + fprintf(stderr, " -V, --skip-variants TYPE Skip indels/snps\n"); + fprintf(stderr, " -v, --variants-only Output variant sites only\n"); fprintf(stderr, "\n"); fprintf(stderr, "Consensus/variant calling options:\n"); - fprintf(stderr, " -c, --consensus-caller the original calling method (conflicts with -m)\n"); - fprintf(stderr, " -C, --constrain one of: alleles, trio (see manual)\n"); - fprintf(stderr, " -m, --multiallelic-caller alternative model for multiallelic and rare-variant calling (conflicts with -c)\n"); - fprintf(stderr, " -n, --novel-rate ,[...] likelihood of novel mutation for constrained trio calling, see man page for details [1e-8,1e-9,1e-9]\n"); - fprintf(stderr, " -p, --pval-threshold variant if P(ref|D) mutation rate (use bigger for greater sensitivity), use with -m [1.1e-3]\n"); + fprintf(stderr, " -c, --consensus-caller The original calling method (conflicts with -m)\n"); + fprintf(stderr, " -C, --constrain STR One of: alleles, trio (see manual)\n"); + fprintf(stderr, " -m, --multiallelic-caller Alternative model for multiallelic and rare-variant calling (conflicts with -c)\n"); + fprintf(stderr, " -n, --novel-rate FLOAT,[...] Likelihood of novel mutation for constrained trio calling, see man page for details [1e-8,1e-9,1e-9]\n"); + fprintf(stderr, " -p, --pval-threshold FLOAT Variant if P(ref|D)= 0) + while ((c = getopt_long(argc, argv, "h?o:O:r:R:s:S:t:T:ANMV:vcmp:C:n:P:f:a:ig:XYF:G:", loptions, NULL)) >= 0) { switch (c) { @@ -969,8 +989,12 @@ int main_vcfcall(int argc, char *argv[]) case 1 : ploidy = optarg; break; case 'X': ploidy = "X"; fprintf(stderr,"Warning: -X will be deprecated, please use --ploidy instead.\n"); break; case 'Y': ploidy = "Y"; fprintf(stderr,"Warning: -Y will be deprecated, please use --ploidy instead.\n"); break; - case 'G': args.aux.sample_groups = optarg; break; - case 'f': args.aux.output_tags |= parse_format_flag(optarg); break; + case 'G': args.aux.sample_groups = strdup(optarg); break; + case 'f': fprintf(stderr,"Warning: -f, --format-fields will be deprecated, please use -a, --annotate instead.\n"); + case 'a': + if (optarg[0]=='?') { list_annotations(stderr); return 1; } + args.aux.output_tags |= parse_output_tags(optarg); + break; case 'M': args.flag &= ~CF_ACGT_ONLY; break; // keep sites where REF is N case 'N': args.flag |= CF_ACGT_ONLY; break; // omit sites where first base in REF is N (the new default) case 'A': args.aux.flag |= CALL_KEEPALT; break;