BACKGROUND
the decreasing cost of producing sequence data has made the sequencing of genomes routine.
researchers use genome resequencing to identify how genomic changes are related to phenotype in their organism of interest.
in the case of humans and certain other genomes , resequencing projects aim to associate genetic changes to disease risk, medical treatment efficacy or other traits of interest.
in some applications it is desirable to detect de novo somatic mutations, which may affect disease progression, prognosis and therapy.
in other applications like genomic medicine for cancer, genomic variants in normal tissue can be compared to genomic variants of the tumor to identify relevant somatic mutations.
however, the accurate detection of single nucleotide variants  and small insertions or deletions  is not trivial.
there is no standard protocol for detecting snp predictions with the highest sensitivity and specificity.
each algorithm used in snp detection creates a different balance of sensitivity and specificity, to either increase the number of true positives at the cost of additional false positives or decrease the number of false positives at the cost of reducing the number of true positives.
additionally, many variant calling algorithms do not provide quantitative values for filtering the vcf output file, or if they do provide users with numerical values for quality score filtering, it often remains unclear to the naïve user what is an appropriate filter.
variant calling programs like gatk and atlas provide only qualitative values for filtering, with language like “pass” or “lowqual”.
in addition, some algorithms, e.g.
gatk, recommend the user include many samples in order to recalibrate quality scores or classify snps with distinctions between pass and lowqual, and thereby increase variant call accuracy.
maximal sensitivity is desirable to minimize false negative calls and therefore avoid missing true mutations.
the consequences of not detecting real variation include: 1) failure to diagnosis or detect real disease and correctly predict elevated or reduced risk for medical problems or potential drug effects, and 2) excess mortality or suffering because of nonintervention or non-optimal treatment.
maximal specificity is also essential to minimize false positive calls and thereby avoid erroneous over-diagnosis and the time, patient distress and expense of confirmatory testing and potential morbidity from unneeded overtreatment.
unfortunately, any classifier performing a nontrivial detection operation on real-world data achieves improved sensitivity only by accepting some elevated rate of false positives, and thus reduced specificity.
this detection error tradeoff  is an essential feature of detection task performance  <cit> .
because any single classifier has an inherent sensitivity versus specificity tradeoff, we hypothesized that more sophisticated methods for systematically integrating the output from multiple independent classifiers  – some with higher inherent sensitivity, some with higher intrinsic specificity - would result in overall improvement in the receiver operating characteristics of the baysic integrated call set compared to the input call sets.
managing sensitivity and specificity of variant calls is critical in projects using genomic data for clinical care  <cit> .
variant call accuracy may be affected by multiple factors, including systematic sequencing error, sequence read depth, allele variant fraction and position-specific error rate, among others.
while there have been recent descriptions of other methods to improve variant call accuracy, including means of combining read mapping and variant call algorithms, these methods typically require training on a gold standard dataset considered to be the truth  <cit> .
by contrast, baysic is a completely unsupervised machine learning method.
baysic does not depend upon training and discordant call arbitration with validated data, yet still achieves gains in sensitivity and specificity over input call sets.
moreover, clinical genome sequencing often involves small sample numbers and/or variant calls in genomic regions with low sequencing coverage.
for example, many clinical applications involve only trios of exomes or genomes, comparing snps differential between two healthy parents and a sick child for diagnosis and treatment selection.
in other clinical cases, real snps could be missed in low read depth regions where the number of reads containing a snp do not meet a strict a priori evidence threshold for inclusion in lists of putative clinically relevant variants  <cit> .
here we describe baysic , a novel algorithm that uses a bayesian statistical method based on latent class analysis to combine variant sets produced by different bioinformatic packages  into a high-confidence set of genome variants.
the strengths of baysic are several.
first, baysic integrates data produced from multiple snp callers, each with differing biases, and produces a call set with a posterior probability that is intuitive and can be used for quantitative filtering.
equally important, baysic is a completely unsupervised method of clustering or classification and requires no training on a “gold standard” or validated data sets.
third, baysic performance improves along with the sensitivity or specificity gains of the input call sets.
if new calling methods yield improved specificity and sensitivity, then baysic will reap those rewards too.
for example, in applications in which sensitivity is a priority, the baysic posterior probability cutoff can be set low to minimize false negatives, and for applications in which specificity is a priority it can be set high to minimize false positives.
baysic run with a posterior probability threshold of  <dig>  produces more sensitive and specific snp prediction than any individual caller used as input.
implementation
baysic algorithm
the user provides variant calls from one or more variant calling programs of their choice in vcf format and, optionally, a posterior probability cutoff .
while not required, the user may also provide a vcf file containing the contents of third party database  as an additional source of variant information for baysic.
false positive and false negative error rates for each evidence source  are estimated using a latent class analysis  approach similar to the approach previously used to combine sets of gene prediction  <cit>  and to infer orthologous genes from different genomes  <cit> .
briefly, this approach assumes a multinomial probability model that uses the number of observed counts for each possible combination of evidence sources that detect a given snp to calculate the underlying parameters for each evidence source: the background frequency of true cases , and the independent and identically distributed  false positive and false negative error rates of each evidence source.
this lca model is implemented using a fully bayesian markov chain monte carlo  simulation using the r2jags r package .
for each of the three parameters to be estimated , a random value is selected from a beta distribution with shape parameters a of  <dig> and b of  <dig> for  <dig>  iterations to yield an estimated value for each of these three parameters.
the posterior probability for each possible combination of agreement amongst the evidence sources  are then calculated as:
θ∏i=1rβi1−xi1−βixiθ∏i=1rβi1−xi1−βixi+1−θ∏i=1rαixi1−αi1−xi
where r is the number of evidence sources used, αi is the false positive rate for the ith program, βi is the false negative rate for the ith evidence source, and θ is the estimate of rate of overall snp occurrence, xi is  <dig> or  <dig> depending on whether the ith evidence source called a snp at the given location.
for each variant, a posterior probability is determined based on which evidence source detected the variant, and the posterior probability cutoff is applied to yield a set of integrated variant calls.
baysic is implemented as a perl script that receives input parameters from the user .
the perl script invokes a separate r script, which computes the α, β and θ parameters and the posterior probabilities for each possible combination of programs.
the perl script then determines the posterior probability for each snp variant based on which callers detected the variant, and writes out to a vcf those variants whose posterior probability is greater than the posterior probability cutoff specified by the user or a default value of  <dig>  if no cutoff was specified.
methods
detection of genome variants using samples from the  <dig> genome project
to detect genome variants, gatk version  <dig> - <dig>  <cit> , atlas version v <dig> .
<dig>  <cit> , samtools version  <dig> .
<dig>  and freebayes version  <dig> .
<dig>  <cit>  were used.
bam files for the following ten samples were downloaded and used as input for the four variant calling programs above: na <dig>  na <dig>  na <dig>  na <dig>  na <dig>  na <dig>  na <dig>  na <dig>  na <dig> and na <dig>  vcf files output by these programs as well as a vcf for dbsnp build version  <dig> were used as input for baysic.
measurement of sensitivity and specificity using data from the  <dig> genome project
sensitivity of each variant detection program was measured as the percent of snps detected by the given program that were confirmed by orthogonal technology  detected by each program.
specificity for each program was measured as the ratio of transitions to transversion  for the set of snp variants produced by each program using vcftools  <cit> .
detection of clinically associated genome variants in a previously verified sample
peripheral blood was taken from a male patient diagnosed with vanishing white matter leukodystrophy, as well as from the unaffected father, mother and sister.
genomic dna was extracted from each sample using standard protocols, and exome capture was carried out using illumina’s truseq technology according to the manufacturer’s protocols.
enriched exome libraries were then subjected to next generation sequencing using standard truseq sample preparation protocols from the manufacturer , and paired end sequencing was carried out on an illumina hiseq.
image analysis and base calling was carried out using casava  <dig> .
bwa was used to align sequence reads to reference genome hg <dig> with subsequent processing by samtools  and picard  to ensure proper file formatting.
alignments were then recalibrated and realigned using gatk  <cit> .
detection of somatic mutations and measurement of sensitivity and specificity in tumor versus normal pair data
using sequencing data from tumor and normal pair from a single patient available in cosmic , we produced somatic mutation calls using mutect  <cit> , varscan <dig>  <cit> , shimmer  <cit>  and strelka  <cit> .
these four sets of somatic mutation calls were combined using baysic with a posterior probability cutoff of  <dig> .
sensitivity was approximated as the overall number of somatic mutations detected by the program, and specificity was measured as percent of somatic mutation calls produced by the program that were present in cosmic version v <dig>  <cit> .
RESULTS
overview of baysic algorithm
several programs exist for the detection of genome variants such as snps and insertions and deletions   <cit> .
previous studies have demonstrated that the agreement between sets of genome variants produced by these methods is poor  <cit> .
the impact of this disagreement among callers on the analytical validity and clinical utility of genomic sequencing is obvious.
baysic allows users to combine two or more sets of genome variants.
the user supplies one or more vcf files containing the sets to be combined and a posterior probability cutoff based on the user’s tolerance for false positive and false negative errors .
optionally, the user may also supply a set of known variants from third party databases in order to increase accuracy, such as dbsnp or cosmic.
the rate of false positive and false negative errors for each set of variant calls are estimated based on the input data using a mcmc simulation, and the posterior probability for each possible combination of agreement between the sets of calls is determined .
the posterior probability cutoff is then applied, and each variant that passes the cutoff is written out to a new vcf file containing the integrated set of variant calls.
sensitivity and specificity of baysic algorithm
to evaluate baysic, we first detected genome variants using ten samples from the  <dig> genomes project  <cit>  using gatk version  <dig>  freebayes, atlas and samtools.
as seen previously  <cit> , there was alarming disagreement among the four sets of genome variant calls.
many snps were present only in one set  .
further, only  <dig> %  of calls were present in all four sets , and only  <dig> %  of snps were present in two or more sets.
we next combined these four sets of variant calls using baysic.
we used as input to baysic the vcf files generated from gatk, freebayes, atlas and samtools as well as a vcf containing variants from dbsnp version  <dig>  the number of positions and posterior probabilities for each possible combination of variant callers and dbsnp are shown in figure  <dig>  for this particular dataset, snps detected by any two prediction methods  would have passed the  <dig>  posterior probability threshold with the exception of a prediction by atlas and dbsnp.
to evaluate the performance of baysic in comparison to existing variant calling programs, we measured the sensitivity and specificity of each method.
sensitivity was measured as the percent of snps detected using an orthogonal technology – array based genotyping   <cit> .
specificity was measured as the ratio of transitions and transversions , previously demonstrated to be approximately  <dig> in coding regions and approximately  <dig> in non-coding regions for true positive snps  <cit> , but  <dig>  for false positive snps  <cit> .
contamination of snp call sets with many false positives results in a ti/tv closer to  <dig> , while fewer false positives will result in a value close to the normal value of ti/tv:  <dig> or  <dig> for coding regions and non-coding regions, respectively.
ti/tv may therefore be used as a measure of specificity since it is proportional to the rate of false positive snp detection.
as expected, as a more stringent posterior probability cutoff was applied, the specificity of the variant set produced by baysic improved at the expense of a slight decrease in sensitivity.
the sensitivity of baysic  in detecting snps in coding regions was identical to the union of the set of all snps detected by these same four programs, which represents maximal sensitivity .
similarly, the sensitivity of baysic calls in noncoding regions was nearly identical to the union of all variant calls.
for snps in coding regions, baysic, using a default posterior probability cutoff of  <dig> , detected  <dig>   of the omni chip-confirmed snps in coding regions, and the union of the call sets, combining calls from each of the callers detected  <dig>   of the omnichip-confirmed snps in coding regions.
for snps in non-coding region, baysic detected  <dig>   of omnichip-confirmed snps, while the union of all snps detected by these four programs yielded  <dig>   of these omni chip-confirmed snps.
as the posterior probability cutoff applied to the baysic set was increased from  <dig>  to  <dig> , the specificity increased while the sensitivity decreased .
at a posterior probability cutoff of  <dig> , the specificity approached the specificity of the intersection of the set of snps called by all four programs with default parameters, which represents maximal specificity .
baysic improved the sensitivity and specificity of the snp detection programs used as input to baysic.
in detecting snps in non-coding regions, baysic with the default posterior probability cutoff of  <dig>  was more sensitive than freebayes, samtools and atlas <dig> and gatk with no filter applied, and more specific than freebayes, gatk and atlas <dig>  and comparable in specificity to samtools .
in detecting snps in coding regions, baysic with the default posterior probability cutoff of  <dig>  was more sensitive than freebayes, samtools and atlas <dig> and gatk with the low quality filter applied, and higher in specificity than gatk and samtools in non-coding regions.
freebayes, atlas <dig> and gatk with low quality filter applied, however, were higher in specificity in coding regions than baysic with the default posterior probability cutoff of  <dig> .
when the baysic posterior probability threshold was increased to  <dig> , the specificity of baysic in coding regions exceeded samtools, freebayes and gatk with low quality filter, and the specificity of baysic in non-coding regions exceeded all  <dig> input call sets.
samtools sensitivity was slightly higher than baysic with a posterior probability cutoff of  <dig> , and atlas <dig> coding region specificity is slightly higher than baysic with posterior probability set to  <dig> .
since other variant calling programs offer filtering options to increase the specificity of snp detection at the expense of sensitivity similar to the posterior probability cutoff available in baysic, we compared the performance of these filtering options with those of baysic.
baysic performed favorably compared with gatk snp call sets filtered using the tranche and vqslod options, and also with freebayes snp call sets filtered using the qual score.
in snps occurring in non-coding regions, baysic  with increasing posterior probability cutoffs described a curve that was above and to the right of curves for gatk with increasingly stringent tranche and vqslod filtering, and freebayes with increasingly stringent qual score filtering .
in snps occurring in coding regions, baysic  with increasingly stringent filtering described a curve that was above and to the right of freebayes using qual filtering, and more sensitive and specific than gatk using tranche filtering when baysic was run with a posterior probability p >  <dig>  .
at p >  <dig>  and p =  <dig> , baysic was slightly more sensitive but less specific than gatk tranche  <dig> and tranche  <dig> , and less sensitive and specific than tranche  <dig>  compared with gatk using with vqslod filtering, baysic  was generally more sensitive, but less specific.
to assess the importance of each individual variant calling program used as input to baysic, we next investigated the effect of leaving out one of the four variant calling programs  on the sensitivity and specificity of baysic.
overall, the sensitivity and specificity of baysic using as input any three of variant caller programs were comparable to the specificity using all four variant calling programs , although the sensitivity of two sets  dropped markedly when applying the highest posterior probability cutoff of  <dig> .
also, applying the same posterior probability cutoff to sets of variant produced by different programs resulted in different sensitivity and specificity .
it is also noteworthy that the inclusion of the gatk set as input into baysic had little effect on the sensitivity and specificity of the resulting integrated snp set produced by baysic.
the sensitivity and specificity of baysic using atlas, freebayes and samtools as input was comparable to that of baysic using atlas, freebayes, samtools and gatk as input .
baysic analysis of exome data from a subject with a previously detected known clinically relevant mutation
to evaluate the performance of baysic in a clinical setting, we tested the ability of these programs to detect a compound heterozygous mutation  in the eif2b <dig> gene, mutations in which have been shown to be causative for vanishing white matter leukodystrophy  <cit> .
the genotypes at these two snps were characterized in a family consisting of one child affected with vanishing white matter leukodystrophy, two healthy parents and one healthy child using standard clinical genetic testing by sanger sequencing.
.
whole exome sequence data was obtained for these same four individuals, and sequence variants were detected using gatk, atlas, samtools, fasd  and snp calls from these four programs were then integrated using baysic.
baysic with a default posterior probability cutoff  correctly identified the genotypes of both snps in all family members.
both snps were identified by three of the four callers, and the posterior probability calculated by baysic for each snp was greater than  <dig> , far above the default cutoff of  <dig> .
gatk and samtools were able to detect the risk allele in all four family members, but gatk did so only when operating with the most relaxed filtering option .
atlas <dig> only detected one  of the two causal snps .
fasd did not predict genotypes for healthy family members; failing to detect the variants as candidates.
therefore, results from samtools or gatk with low quality filter would have lead researchers to correctly identify the compound heterozygous snps.
however, if researchers were to use any of the more stringent filtering options commonly used in gatk, e.g., tranche  <dig>   <dig>   <dig>  or a pass filter, these snps would not have been detected.
because one of the two snps is rare, our integrated approach provides additional confidence for researchers interested in detecting rare snps.
snp variants were detected in samples from four individuals from a single family with one individual affected with white matter leukodytrophy  <dig> healthy parents and a healthy sibling.
for each of two causative snps in the eif2b <dig> gene associated with white matter leukodystrophy , the ability of atlas, gatk, samtools, fasd or baysic  to detect each variant is shown.
for gatk, the most stringent filter that could be applied and still detect each variant is shown in parentheses.
“detected” indicates that the snp was predicted for one of the four subjects.
using baysic to combine sets of somatic mutation calls produced with tumor/normal pair data
a common application of genome sequencing is to sequence samples taken from normal and tumorous tissue and detect somatic mutations that may be involved in cancer  <cit> .
many programs exist to detect somatic mutations, but as with programs for detecting snp variants, the agreement of these programs is poor  <cit> .
the problem of combining these sets of somatic mutations is analogous to the problem of combining disparate sets of snps produced by different snp detection programs.
we applied baysic to this related problem of combining disparate sets of somatic mutation calls.
using sequencing data from tumor and normal pair from a single patient available in a catalog of previously observed somatic mutations , we produced somatic mutation calls using mutect  <cit> , varscan <dig>  <cit> , shimmer  <cit>  and strelka  <cit> , and then combined these four sets of somatic mutation calls using baysic with a default posterior probability cutoff of  <dig> .
as a measure of specificity, we determined the overall number of somatic mutations detected by each program that were present in cosmic .
mutect, varscan <dig>  shimmer, strelka and baysic, using as input the sets of somatic mutation produced by all four callers, detected  <dig>   <dig>   <dig>   <dig> and  <dig> somatic mutations that were present in cosmic, respectively ; this translates to  <dig> %,  <dig> %,  <dig> %,  <dig> % and  <dig> % of the total snps.
using this as a measure, baysic therefore improves the specificity of all four callers used as input.
baysic predicted a lower number of somatic mutations  compared to mutect, varscan and shimmer , but more than strelka .
if the overall number of somatic mutations is taken as a measure of sensitivity, the sensitivity of baysic is lower than mutect and shimmer, much lower than varscan and higher than strelka.
in general, each of the four individual callers is highly concordant with the genotype predicted by the genotype chip, one measure of accuracy.
for all programs, all somatic mutations that occurred at positions present on the genotype chip predicted genotypes in agreement with the genotype chip .
somatic mutations from exome data from a single patient were predicted using mutect, strelka, varscan <dig> and shimmer and these four sets of somatic mutation calls were combined using baysic with a posterior probability cutoff of  <dig> .
the total number of somatic mutations detected by each program, as well as the agreement of these somatic mutations with those cosmic and with sets determined by chip are shown.
CONCLUSIONS
clinical applications of genomics demand reliable detection of real variants and discrimination and rejection of false alarms due to sequencing error, low sequence coverage or low allelic variant fraction.
accordingly, the utility of genomic medicine will be improved by better methods for accurately identifying snps and other genomic variants.
our analyses support our initial hypothesis: baysic variant calls demonstrated improved variant detection accuracy and superior receiver operating characteristics compared to the variant call methods used as input for baysic.
importantly, baysic will accept as input any number of alternative variant detection algorithms, allowing the user to combine methods that emphasize sensitivity with methods that enhance specificity and achieve overall gains in detection accuracy.
as the sensitivity or specificity of input call sets improve, the sensitivity and specificity of baysic variant calls also increases.
likewise, baysic may be used to focus on specific types of variant detection problems such as somatic mutations in cancer, and achieves similar gains in receiver operating characteristics compared to the individual somatic variant calling algorithms used as input.
another program was recently described to combine somatic mutation calls  <cit> .
future work will determine the relative performance of baysic compared with this program, and assess how inclusion of improved somatic mutation call sets, as input to baysic, affects baysic’s overall performance in somatic mutation detection.
it is possible that the degree of improvement offered by baysic in combining sets of germline snp variant calls compared to somatic mutation calls is explainable by the different error rates in these two different experimental settings.
that is, germline snp discovery has very low false positive and low false negative rates relative to somatic mutation calls, with generally good sensitivity and specificity  <cit> .
therefore, producing a consensus set of germline snp variants with baysic provides marginal but noticeable improvements to both sensitivity and specificity .
in contrast, somatic mutation discovery has very high false positive  rates, with poor specificity .
producing a consensus snp set using baysic therefore makes dramatic improvements to specificity without losing sensitivity .
it is possible that correlations between the errors in the sets of variant calls used as input to baysic could result in false positive errors in the integrated variant set produced by baysic.
to address this, future versions of baysic will measure the bivariate residuals after latent class analysis is performed, and will penalize the significance of input sets that are highly correlated  <cit> .
baysic currently only integrates sets of snp variant calls.
future work will expand this to include other sorts of variants such as insertions/deletions , and additional modifications to facilitate improved performance in somatic mutation detection.
availability and requirements
project name: baysic
project home page:http://genformatic.com/baysic
operating systems: linux, os x, windows
programming languages: perl, r
other requirements: jags, json file::temp getopt::long list::util file::next test::warn file::slurp perlio::gzip file::which local::lib
license: free for academic use, license needed for commercial use
competing interests
bc, dw and jr own an interest in genformatic, llc, a company which performs sequencing and sequence data analysis as a service.
this does not alter the authors’ adherence to the journal’s policies on sharing data and materials.
author’s contributions
jr, am, bc, dw, nm and jz conceived, designed and executed experiments.
nm and jz provided access to clinical data.
am and jr designed the algorithm and wrote software.
all authors contributed to preparation of the manuscript.

 the decreasing cost of producing sequence data has made the sequencing of genomes routine. researchers use genome resequencing to identify how genomic changes are related to phenotype in their organism of interest. in the case of humans and certain other genomes, resequencing projects aim to associate genetic changes to disease risk, medical treatment efficacy or other traits of interest. in some applications it is desirable to detect de novo somatic mutations, which may affect disease progression, prognosis and therapy. in other applications like genomic medicine for cancer, genomic variants in normal tissue can be compared to genomic variants of the tumor to identify relevant somatic mutations. however, the accurate detection of single nucleotide variants and small insertions or deletions is not trivial. there is no standard protocol for detecting snp predictions with the highest sensitivity and specificity. each algorithm used in snp detection creates a different balance of sensitivity and specificity, to either increase the number of true positives at the cost of additional false positives or decrease the number of false positives at the cost of reducing the number of true positives. additionally, many variant calling algorithms do not provide quantitative values for filtering the vcf output file, or if they do provide users with numerical values for quality score filtering, it often remains unclear to the naïve user what is an appropriate filter. variant calling programs like gatk and atlas provide only qualitative values for filtering, with language like “pass” or “lowqual”. in addition, some algorithms, e.g. gatk, recommend the user include many samples in order to recalibrate quality scores or classify snps with distinctions between pass and lowqual, and thereby increase variant call accuracy. maximal sensitivity is desirable to minimize false negative calls and therefore avoid missing true mutations. the consequences of not detecting real variation include: 1) failure to diagnosis or detect real disease and correctly predict elevated or reduced risk for medical problems or potential drug effects, and 2) excess mortality or suffering because of nonintervention or non-optimal treatment. maximal specificity is also essential to minimize false positive calls and thereby avoid erroneous over-diagnosis and the time, patient distress and expense of confirmatory testing and potential morbidity from unneeded overtreatment. unfortunately, any classifier performing a nontrivial detection operation on real-world data achieves improved sensitivity only by accepting some elevated rate of false positives, and thus reduced specificity. this detection error tradeoff is an essential feature of detection task performance. because any single classifier has an inherent sensitivity versus specificity tradeoff, we hypothesized that more sophisticated methods for systematically integrating the output from multiple independent classifiers – some with higher inherent sensitivity, some with higher intrinsic specificity - would result in overall improvement in the receiver operating characteristics of the baysic integrated call set compared to the input call sets. managing sensitivity and specificity of variant calls is critical in projects using genomic data for clinical care. variant call accuracy may be affected by multiple factors, including systematic sequencing error, sequence read depth, allele variant fraction and position-specific error rate, among others. while there have been recent descriptions of other methods to improve variant call accuracy, including means of combining read mapping and variant call algorithms, these methods typically require training on a gold standard dataset considered to be the truth. by contrast, baysic is a completely unsupervised machine learning method. baysic does not depend upon training and discordant call arbitration with validated data, yet still achieves gains in sensitivity and specificity over input call sets. moreover, clinical genome sequencing often involves small sample numbers and/or variant calls in genomic regions with low sequencing coverage. for example, many clinical applications involve only trios of exomes or genomes, comparing snps differential between two healthy parents and a sick child for diagnosis and treatment selection. in other clinical cases, real snps could be missed in low read depth regions where the number of reads containing a snp do not meet a strict a priori evidence threshold for inclusion in lists of putative clinically relevant variants. here we describe baysic, a novel algorithm that uses a bayesian statistical method based on latent class analysis to combine variant sets produced by different bioinformatic packages into a high-confidence set of genome variants. the strengths of baysic are several. first, baysic integrates data produced from multiple snp callers, each with differing biases, and produces a call set with a posterior probability that is intuitive and can be used for quantitative filtering. equally important, baysic is a completely unsupervised method of clustering or classification and requires no training on a “gold standard” or validated data sets. third, baysic performance improves along with the sensitivity or specificity gains of the input call sets. if new calling methods yield improved specificity and sensitivity, then baysic will reap those rewards too. for example, in applications in which sensitivity is a priority, the baysic posterior probability cutoff can be set low to minimize false negatives, and for applications in which specificity is a priority it can be set high to minimize false positives. baysic run with a posterior probability threshold of produces more sensitive and specific snp prediction than any individual caller used as input. implementation baysic algorithm the user provides variant calls from one or more variant calling programs of their choice in vcf format and, optionally, a posterior probability cutoff. while not required, the user may also provide a vcf file containing the contents of third party database as an additional source of variant information for baysic. false positive and false negative error rates for each evidence source are estimated using a latent class analysis approach similar to the approach previously used to combine sets of gene prediction and to infer orthologous genes from different genomes. briefly, this approach assumes a multinomial probability model that uses the number of observed counts for each possible combination of evidence sources that detect a given snp to calculate the underlying parameters for each evidence source: the background frequency of true cases, and the independent and identically distributed false positive and false negative error rates of each evidence source. this lca model is implemented using a fully bayesian markov chain monte carlo simulation using the r2jags r package. for each of the three parameters to be estimated, a random value is selected from a beta distribution with shape parameters a of and b of for iterations to yield an estimated value for each of these three parameters. the posterior probability for each possible combination of agreement amongst the evidence sources are then calculated as: θ∏i=1rβi1−xi1−βixiθ∏i=1rβi1−xi1−βixi+1−θ∏i=1rαixi1−αi1−xi where r is the number of evidence sources used, αi is the false positive rate for the ith program, βi is the false negative rate for the ith evidence source, and θ is the estimate of rate of overall snp occurrence, xi is or depending on whether the ith evidence source called a snp at the given location. for each variant, a posterior probability is determined based on which evidence source detected the variant, and the posterior probability cutoff is applied to yield a set of integrated variant calls. baysic is implemented as a perl script that receives input parameters from the user. the perl script invokes a separate r script, which computes the α, β and θ parameters and the posterior probabilities for each possible combination of programs. the perl script then determines the posterior probability for each snp variant based on which callers detected the variant, and writes out to a vcf those variants whose posterior probability is greater than the posterior probability cutoff specified by the user or a default value of if no cutoff was specified. methods detection of genome variants using samples from the genome project to detect genome variants, gatk version -, atlas version v., samtools version. and freebayes version. were used. bam files for the following ten samples were downloaded and used as input for the four variant calling programs above: na na na na na na na na na and na vcf files output by these programs as well as a vcf for dbsnp build version were used as input for baysic. measurement of sensitivity and specificity using data from the genome project sensitivity of each variant detection program was measured as the percent of snps detected by the given program that were confirmed by orthogonal technology detected by each program. specificity for each program was measured as the ratio of transitions to transversion for the set of snp variants produced by each program using vcftools. detection of clinically associated genome variants in a previously verified sample peripheral blood was taken from a male patient diagnosed with vanishing white matter leukodystrophy, as well as from the unaffected father, mother and sister. genomic dna was extracted from each sample using standard protocols, and exome capture was carried out using illumina’s truseq technology according to the manufacturer’s protocols. enriched exome libraries were then subjected to next generation sequencing using standard truseq sample preparation protocols from the manufacturer, and paired end sequencing was carried out on an illumina hiseq. image analysis and base calling was carried out using casava. bwa was used to align sequence reads to reference genome hg with subsequent processing by samtools and picard to ensure proper file formatting. alignments were then recalibrated and realigned using gatk. detection of somatic mutations and measurement of sensitivity and specificity in tumor versus normal pair data using sequencing data from tumor and normal pair from a single patient available in cosmic, we produced somatic mutation calls using mutect, varscan, shimmer and strelka. these four sets of somatic mutation calls were combined using baysic with a posterior probability cutoff of. sensitivity was approximated as the overall number of somatic mutations detected by the program, and specificity was measured as percent of somatic mutation calls produced by the program that were present in cosmic version v.  overview of baysic algorithm several programs exist for the detection of genome variants such as snps and insertions and deletions. previous studies have demonstrated that the agreement between sets of genome variants produced by these methods is poor. the impact of this disagreement among callers on the analytical validity and clinical utility of genomic sequencing is obvious. baysic allows users to combine two or more sets of genome variants. the user supplies one or more vcf files containing the sets to be combined and a posterior probability cutoff based on the user’s tolerance for false positive and false negative errors. optionally, the user may also supply a set of known variants from third party databases in order to increase accuracy, such as dbsnp or cosmic. the rate of false positive and false negative errors for each set of variant calls are estimated based on the input data using a mcmc simulation, and the posterior probability for each possible combination of agreement between the sets of calls is determined. the posterior probability cutoff is then applied, and each variant that passes the cutoff is written out to a new vcf file containing the integrated set of variant calls. sensitivity and specificity of baysic algorithm to evaluate baysic, we first detected genome variants using ten samples from the genomes project using gatk version freebayes, atlas and samtools. as seen previously, there was alarming disagreement among the four sets of genome variant calls. many snps were present only in one set. further, only % of calls were present in all four sets, and only % of snps were present in two or more sets. we next combined these four sets of variant calls using baysic. we used as input to baysic the vcf files generated from gatk, freebayes, atlas and samtools as well as a vcf containing variants from dbsnp version the number of positions and posterior probabilities for each possible combination of variant callers and dbsnp are shown in figure  for this particular dataset, snps detected by any two prediction methods would have passed the posterior probability threshold with the exception of a prediction by atlas and dbsnp. to evaluate the performance of baysic in comparison to existing variant calling programs, we measured the sensitivity and specificity of each method. sensitivity was measured as the percent of snps detected using an orthogonal technology – array based genotyping. specificity was measured as the ratio of transitions and transversions, previously demonstrated to be approximately in coding regions and approximately in non-coding regions for true positive snps, but for false positive snps. contamination of snp call sets with many false positives results in a ti/tv closer to, while fewer false positives will result in a value close to the normal value of ti/tv: or for coding regions and non-coding regions, respectively. ti/tv may therefore be used as a measure of specificity since it is proportional to the rate of false positive snp detection. as expected, as a more stringent posterior probability cutoff was applied, the specificity of the variant set produced by baysic improved at the expense of a slight decrease in sensitivity. the sensitivity of baysic in detecting snps in coding regions was identical to the union of the set of all snps detected by these same four programs, which represents maximal sensitivity. similarly, the sensitivity of baysic calls in noncoding regions was nearly identical to the union of all variant calls. for snps in coding regions, baysic, using a default posterior probability cutoff of, detected of the omni chip-confirmed snps in coding regions, and the union of the call sets, combining calls from each of the callers detected of the omnichip-confirmed snps in coding regions. for snps in non-coding region, baysic detected of omnichip-confirmed snps, while the union of all snps detected by these four programs yielded of these omni chip-confirmed snps. as the posterior probability cutoff applied to the baysic set was increased from to, the specificity increased while the sensitivity decreased. at a posterior probability cutoff of, the specificity approached the specificity of the intersection of the set of snps called by all four programs with default parameters, which represents maximal specificity. baysic improved the sensitivity and specificity of the snp detection programs used as input to baysic. in detecting snps in non-coding regions, baysic with the default posterior probability cutoff of was more sensitive than freebayes, samtools and atlas and gatk with no filter applied, and more specific than freebayes, gatk and atlas and comparable in specificity to samtools. in detecting snps in coding regions, baysic with the default posterior probability cutoff of was more sensitive than freebayes, samtools and atlas and gatk with the low quality filter applied, and higher in specificity than gatk and samtools in non-coding regions. freebayes, atlas and gatk with low quality filter applied, however, were higher in specificity in coding regions than baysic with the default posterior probability cutoff of. when the baysic posterior probability threshold was increased to, the specificity of baysic in coding regions exceeded samtools, freebayes and gatk with low quality filter, and the specificity of baysic in non-coding regions exceeded all input call sets. samtools sensitivity was slightly higher than baysic with a posterior probability cutoff of, and atlas coding region specificity is slightly higher than baysic with posterior probability set to. since other variant calling programs offer filtering options to increase the specificity of snp detection at the expense of sensitivity similar to the posterior probability cutoff available in baysic, we compared the performance of these filtering options with those of baysic. baysic performed favorably compared with gatk snp call sets filtered using the tranche and vqslod options, and also with freebayes snp call sets filtered using the qual score. in snps occurring in non-coding regions, baysic with increasing posterior probability cutoffs described a curve that was above and to the right of curves for gatk with increasingly stringent tranche and vqslod filtering, and freebayes with increasingly stringent qual score filtering. in snps occurring in coding regions, baysic with increasingly stringent filtering described a curve that was above and to the right of freebayes using qual filtering, and more sensitive and specific than gatk using tranche filtering when baysic was run with a posterior probability p > . at p >  and p = , baysic was slightly more sensitive but less specific than gatk tranche and tranche, and less sensitive and specific than tranche compared with gatk using with vqslod filtering, baysic was generally more sensitive, but less specific. to assess the importance of each individual variant calling program used as input to baysic, we next investigated the effect of leaving out one of the four variant calling programs on the sensitivity and specificity of baysic. overall, the sensitivity and specificity of baysic using as input any three of variant caller programs were comparable to the specificity using all four variant calling programs, although the sensitivity of two sets dropped markedly when applying the highest posterior probability cutoff of. also, applying the same posterior probability cutoff to sets of variant produced by different programs resulted in different sensitivity and specificity. it is also noteworthy that the inclusion of the gatk set as input into baysic had little effect on the sensitivity and specificity of the resulting integrated snp set produced by baysic. the sensitivity and specificity of baysic using atlas, freebayes and samtools as input was comparable to that of baysic using atlas, freebayes, samtools and gatk as input. baysic analysis of exome data from a subject with a previously detected known clinically relevant mutation to evaluate the performance of baysic in a clinical setting, we tested the ability of these programs to detect a compound heterozygous mutation in the eif2b gene, mutations in which have been shown to be causative for vanishing white matter leukodystrophy. the genotypes at these two snps were characterized in a family consisting of one child affected with vanishing white matter leukodystrophy, two healthy parents and one healthy child using standard clinical genetic testing by sanger sequencing.. whole exome sequence data was obtained for these same four individuals, and sequence variants were detected using gatk, atlas, samtools, fasd and snp calls from these four programs were then integrated using baysic. baysic with a default posterior probability cutoff correctly identified the genotypes of both snps in all family members. both snps were identified by three of the four callers, and the posterior probability calculated by baysic for each snp was greater than, far above the default cutoff of. gatk and samtools were able to detect the risk allele in all four family members, but gatk did so only when operating with the most relaxed filtering option. atlas only detected one of the two causal snps. fasd did not predict genotypes for healthy family members; failing to detect the variants as candidates. therefore, results from samtools or gatk with low quality filter would have lead researchers to correctly identify the compound heterozygous snps. however, if researchers were to use any of the more stringent filtering options commonly used in gatk, e.g., tranche or a pass filter, these snps would not have been detected. because one of the two snps is rare, our integrated approach provides additional confidence for researchers interested in detecting rare snps. snp variants were detected in samples from four individuals from a single family with one individual affected with white matter leukodytrophy healthy parents and a healthy sibling. for each of two causative snps in the eif2b gene associated with white matter leukodystrophy, the ability of atlas, gatk, samtools, fasd or baysic to detect each variant is shown. for gatk, the most stringent filter that could be applied and still detect each variant is shown in parentheses. “detected” indicates that the snp was predicted for one of the four subjects. using baysic to combine sets of somatic mutation calls produced with tumor/normal pair data a common application of genome sequencing is to sequence samples taken from normal and tumorous tissue and detect somatic mutations that may be involved in cancer. many programs exist to detect somatic mutations, but as with programs for detecting snp variants, the agreement of these programs is poor. the problem of combining these sets of somatic mutations is analogous to the problem of combining disparate sets of snps produced by different snp detection programs. we applied baysic to this related problem of combining disparate sets of somatic mutation calls. using sequencing data from tumor and normal pair from a single patient available in a catalog of previously observed somatic mutations, we produced somatic mutation calls using mutect, varscan, shimmer and strelka, and then combined these four sets of somatic mutation calls using baysic with a default posterior probability cutoff of. as a measure of specificity, we determined the overall number of somatic mutations detected by each program that were present in cosmic. mutect, varscan shimmer, strelka and baysic, using as input the sets of somatic mutation produced by all four callers, detected and somatic mutations that were present in cosmic, respectively; this translates to %, %, %, % and % of the total snps. using this as a measure, baysic therefore improves the specificity of all four callers used as input. baysic predicted a lower number of somatic mutations compared to mutect, varscan and shimmer, but more than strelka. if the overall number of somatic mutations is taken as a measure of sensitivity, the sensitivity of baysic is lower than mutect and shimmer, much lower than varscan and higher than strelka. in general, each of the four individual callers is highly concordant with the genotype predicted by the genotype chip, one measure of accuracy. for all programs, all somatic mutations that occurred at positions present on the genotype chip predicted genotypes in agreement with the genotype chip. somatic mutations from exome data from a single patient were predicted using mutect, strelka, varscan and shimmer and these four sets of somatic mutation calls were combined using baysic with a posterior probability cutoff of. the total number of somatic mutations detected by each program, as well as the agreement of these somatic mutations with those cosmic and with sets determined by chip are shown.  clinical applications of genomics demand reliable detection of real variants and discrimination and rejection of false alarms due to sequencing error, low sequence coverage or low allelic variant fraction. accordingly, the utility of genomic medicine will be improved by better methods for accurately identifying snps and other genomic variants. our analyses support our initial hypothesis: baysic variant calls demonstrated improved variant detection accuracy and superior receiver operating characteristics compared to the variant call methods used as input for baysic. importantly, baysic will accept as input any number of alternative variant detection algorithms, allowing the user to combine methods that emphasize sensitivity with methods that enhance specificity and achieve overall gains in detection accuracy. as the sensitivity or specificity of input call sets improve, the sensitivity and specificity of baysic variant calls also increases. likewise, baysic may be used to focus on specific types of variant detection problems such as somatic mutations in cancer, and achieves similar gains in receiver operating characteristics compared to the individual somatic variant calling algorithms used as input. another program was recently described to combine somatic mutation calls. future work will determine the relative performance of baysic compared with this program, and assess how inclusion of improved somatic mutation call sets, as input to baysic, affects baysic’s overall performance in somatic mutation detection. it is possible that the degree of improvement offered by baysic in combining sets of germline snp variant calls compared to somatic mutation calls is explainable by the different error rates in these two different experimental settings. that is, germline snp discovery has very low false positive and low false negative rates relative to somatic mutation calls, with generally good sensitivity and specificity. therefore, producing a consensus set of germline snp variants with baysic provides marginal but noticeable improvements to both sensitivity and specificity. in contrast, somatic mutation discovery has very high false positive rates, with poor specificity. producing a consensus snp set using baysic therefore makes dramatic improvements to specificity without losing sensitivity. it is possible that correlations between the errors in the sets of variant calls used as input to baysic could result in false positive errors in the integrated variant set produced by baysic. to address this, future versions of baysic will measure the bivariate residuals after latent class analysis is performed, and will penalize the significance of input sets that are highly correlated. baysic currently only integrates sets of snp variant calls. future work will expand this to include other sorts of variants such as insertions/deletions, and additional modifications to facilitate improved performance in somatic mutation detection. availability and requirements project name: baysic project home page:http://genformatic.com/baysic operating systems: linux, os x, windows programming languages: perl, r other requirements: jags, json file::temp getopt::long list::util file::next test::warn file::slurp perlio::gzip file::which local::lib license: free for academic use, license needed for commercial use competing interests bc, dw and jr own an interest in genformatic, llc, a company which performs sequencing and sequence data analysis as a service. this does not alter the authors’ adherence to the journal’s policies on sharing data and materials. author’s contributions jr, am, bc, dw, nm and jz conceived, designed and executed experiments. nm and jz provided access to clinical data. am and jr designed the algorithm and wrote software. all authors contributed to preparation of the manuscript. 
unlike other classification methods, the underlying baysic model does not require training using a “gold standard” of true positives. we demonstrated baysic’s superior variant-calling accuracy, both for somatic mutation detection and germline variant detection. here we describe baysic, which combines snp variant calls produced by different methods into a more accurate set of variant calls. baysic differs from majority voting, consensus or other ad hoc intersection-based schemes for combining sets of genome variant calls. baysic provides a method for combining sets of snp variant calls produced by different variant calling programs. the integrated set of snp variant calls produced by baysic improves the sensitivity and specificity of the variant calls used as input. the user specifies a posterior probability threshold according to the user’s tolerance for false positive and false negative errors; lowering the posterior probability threshold allows the user to trade specificity for sensitivity while raising the threshold increases specificity in exchange for sensitivity. in addition to combining sets of germline variants, baysic can also be used to combine sets of somatic mutations detected in the context of tumor/normal sequencing experiments.
 accurate genomic variant detection is an essential step in gleaning medically useful information from genome data. however, low concordance among variant-calling methods reduces confidence in the clinical validity of whole genome and exome sequence data, and confounds downstream analysis for applications in genome medicine. here we describe baysic, which combines snp variant calls produced by different methods into a more accurate set of variant calls. baysic differs from majority voting, consensus or other ad hoc intersection-based schemes for combining sets of genome variant calls. unlike other classification methods, the underlying baysic model does not require training using a “gold standard” of true positives. rather, with each new dataset, baysic performs an unsupervised, fully bayesian latent class analysis to estimate false positive and false negative error rates for each input method. the user specifies a posterior probability threshold according to the user’s tolerance for false positive and false negative errors; lowering the posterior probability threshold allows the user to trade specificity for sensitivity while raising the threshold increases specificity in exchange for sensitivity.  we assessed the performance of baysic in comparison to other variant detection methods using ten low coverage samples from the genomes project, a tumor/normal exome pair, and exome sequences from positive control samples previously identified to contain clinically relevant snps. we demonstrated baysic’s superior variant-calling accuracy, both for somatic mutation detection and germline variant detection.  baysic provides a method for combining sets of snp variant calls produced by different variant calling programs. the integrated set of snp variant calls produced by baysic improves the sensitivity and specificity of the variant calls used as input. in addition to combining sets of germline variants, baysic can also be used to combine sets of somatic mutations detected in the context of tumor/normal sequencing experiments. snpgenome variantsbayesianlatent class analysiscancersomatic mutation 