BACKGROUND
practical application of genomic technologies, such as large-scale use of single nucleotide variant  arrays in animal and plant breeding, has become routine in many areas of the life sciences.
taking both polygenic additive  effects and genomic  effects into account, between  <dig> and 85% of the genetic variance observed in phenotypic traits of interest in cattle can be explained solely by snv effects
<cit>  and the number of genotyped animals in cattle populations worldwide is increasing steadily
<cit> .
fueled by decreasing costs, advances in next-generation sequencing  technologies enable identification of more complex forms of genetic variation , copy number variations , etc.).
these advances will inevitably foster our ability to partition the genetic variance underlying traits of interest.
while some applications of ngs require de novo sequencing of an individual organism , re-sequencing may also be possible if a reference genome for the species of interest is available.
the translation of raw ngs reads into tangible variants  via re-sequencing is a specific, delicate and computationally demanding task
<cit>  and comprises three steps.
first, short reads of dna are aligned to an existing reference genome .
second, sequence differences between the sample being sequenced and the reference genome are identified .
a myriad of alignment  and variant identification software programmes are available  or the haplotypecaller  of the genome analysis toolkit
<cit> ; platypus
<cit> ; samtools
<cit> ; etc.
), the majority of which can be obtained free of charge.
as a final step, variants are screened and filtered to remove potential false positives common in most ngs technologies .
as sequencing costs decline, reference genomes are becoming available for an increasing number of organisms, including agriculturally important species such as cattle .
the bos taurus reference genome umd <dig>  is similar in size to the human genome and contains ~ <dig>  billion base pairs, approximately 10% of which are not assigned to any chromosome
<cit> .
the n <dig> size can be used to compare the quality of genome assemblies of similar size: it represents contig size such that 50% of the genome is contained in contigs of length n or greater
<cit> .
because the n <dig> size for the umd <dig>  reference genome  is much shorter than that of the current human reference , the umd <dig>  reference genome will not likely allow the same accuracy in alignment, variant identification and further downstream analysis as the human reference allows.
nevertheless, algorithms and software developed for alignment and variant identification of human ngs data provide an excellent resource for translating ngs data of other non-human organisms, such as cattle, into genetic variants for application in genome wide association studies and genomic selection programmes.
several approaches to variant identification are possible.
the simplest variant detection methods identify variants on a per-sample basis, one position at a time.
once a variant locus is found, the most likely genotype for that locus is determined stochastically based on a consensus of aligned reads.
if multiple samples are analyzed simultaneously, an a priori likelihood of finding a variant locus given the observed data is derived, and the most likely genotype at a given position is determined.
either single or multi sample variant identification methods can be implemented in the ug
<cit> ; and sam
<cit> .
more advanced haplotype-based methods incorporate the correlation between adjacent variants within the variant detection procedure.
such methods use linkage disequilibrium between nearby variants to further enrich variant identification.
haplotype-based methods are implemented in pl
<cit>  and the gatk hc
<cit> .
the haplotype-based variant detection approach can also be conducted in either single or multi sample settings.
read chimerism, base pair tautomerisms and signal intensity issues can contribute to false positive variant detection by causing stochastic inaccuracies, general sequencing errors, and misalignments in ngs data
<cit> .
aside from the variant identification approach itself, a number of optional auxiliary steps have been recommended to improve the quality of ngs-derived variants.
these steps are conducted before  or after variant identification .
the first generally recommended preparatory step is the identification of falsely duplicated reads , which reduces bias in variant detection
<cit> .
secondly, local realignment around single or multiple bases that are either missing in the reference  or missing in the dna sequence being analyzed  is also commonly recommended
<cit> .
realignment cleans up spurious snvs that result from misalignment of reads around known alignment gaps and helps detect false negative snvs in the near vicinity of indels .
furthermore, the full alignment context is used to determine whether the reported divergence from the reference  actually exists
<cit> .
finally, flow cell lane, machine cycle , sequencing context  or other technical aspects may influence base quality scores, which help characterize the quality of the bases in the individual reads.
base quality score recalibration is recommended to lower the number of falsely identified snvs and to lower false confidence in identified bases
<cit> .
after variant identification, further filters can be implemented according to the individual dataset under consideration.
variant detection methods can be evaluated using genome simulations based on reference genomes
<cit> .
given that all mammalian reference genomes are incomplete and whole-genome alignment is imperfect, simulations may not provide realistic results, a difficulty that has prompted the use of real data
<cit> .
another common evaluation method in human studies is to compare variants from ngs data to those of array genotyping , to results from the  <dig> genomes project
<cit> , or to other existing resources of human sequence .
resources for cattle are not as comprehensive, however low, medium and high density snv arrays are available , thus allowing an estimation of relative accuracy between software .table  <dig>
variant genotypes compared for concordance between the array-based and sequence based methods to determine concordance, sensitivity and discrepancy between the two assays  and measures of concordance
a) ngs-based genotypes
array-based genotypes
aa
ab
bb
aa
a
b
c
ab
d
e
f
bb
g
h
i
–
k
l
m
b)
measures of concordance
array-based information from the illumina bovinehd beadchip®  was considered the “gold-standard” and compared to next-generation sequencing-based variants obtained using a illumina hiseq <dig> platform with various variant identification software, where genotypes are identified as:
a = homozygous reference in both ngs-based data and array-based data.
b = homozygous reference in ngs-based data, but as heterozygous in array-based data.
c = homozygous reference in ngs-based data, but as homozygous alternative in array-based data.
d = heterozygous in ngs-based data, but as homozygous reference in array-based data.
e = heterozygous in both ngs-based data and array-based data.
f = heterozygous in ngs-based data, but as homozygous alternative in array-based data.
g = homozygous reference in ngs-based data, but as homozygous reference in array-based data.
h = heterozygous in ngs-based data and array-based data, but as heterozygous in array-based data.
i = homozygous alternative in both ngs-based data and array-based data.
k = not found in ngs-based data, but as homozygous reference in array-based data.
l = not found in ngs-based data and array-based data, but as heterozygous in array-based data.
m = not found in ngs-based data, but as homozygous alternative in array-based data.
.
aside from measuring concordance between ngs data and array data, the ratio of transitions  to transversions  can be used as a convenient diagnostic to measure the quality of ngs data .
the genome-wide ti/tv ratio is reported between  <dig>  and  <dig>  in human whole-genome sequence data
<cit> , whereby this ratio is higher in exomes due to the increased presence of methelated cytosine in cpg dinucleotides in exonic regions
<cit> .
because this bias in favour of mutations between bases with similar biochemical properties  over those with dissimilar properties  is dependent on both cpg and gc content of the region, the ti/tv ratio is a useful diagnostic to measure quality across the genome
<cit> .
the objective of this study was to investigate which methods and software work best for detection of high quality genetic variants using ngs data in cattle, with a specific focus on single nucleotide variants.
using whole-genome sequence information from  <dig> individuals, we a) explore the implications of preparatory steps commonly recommended in human analysis, b) compare results of single and multi sample variant detection achieved using four publicly available variant detection software programmes, c) provide a comparison of computational processing time, and d) compare accuracy and completeness of snvs identified in ngs data by comparing them to genotypes from the same individuals generated with either high- or medium-density snv arrays, as well as to analyse genome-wide ti/tv ratios.
through benchmarking different variant detection methods in cattle, preliminary recommendations for variant identification in other organisms can be extrapolated.
our findings can serve as a reference for choosing variant identification software and can help assess the implication of preparatory steps in ngs pipelines for species with lower-quality or unfinished reference genomes.
RESULTS
alignment and coverage
approximately  <dig> billion paired-end reads were obtained for the  <dig> sequenced animals.
an average of  <dig> % of these reads  were mapped to  <dig> chromosomes  of the bovine reference genome assembly umd <dig>
<cit> .
approximately  <dig>  billion  duplicate reads  were marked and excluded from further analysis.
average coverage was  <dig>  reads per base; average coverage per-animal ranged from  <dig>  -  <dig> .
see additional file
1: table s <dig> for individual alignment and coverage information.
single sample variant detection
the results of single sample variant detection are in figure
<dig>  depending on software and preparatory steps, an average of  <dig> , <dig> -  <dig> , <dig> snvs ,  <dig>  –  <dig>  indels  and  <dig>  –  <dig>  multi-allelic sites  were identified per animal .
detailed information on the sequencing process  is given in additional file
2: table s <dig>  ug identified the most snvs, followed by pl and sam.
because pl identifies multi-nucleotide variants, pl originally had the lowest average number of snv per sample  because some variants were “hidden” in multi-snv replacements.
after splitting multi-nucleotide variants into their allelic primitives , a fair comparison was achieved and the average per-sample number of variants identified with pl_prim increased to  <dig> , <dig>  when single-sample results of all  <dig> animals were combined,  <dig> , <dig> -  <dig> , <dig> unique snv were identified .figure  <dig>
distributions of single nucleotide variant counts , insertion and deletion counts , and multi-allelic site counts  identified per animal.
for platypus results, multi-nucleotide variants were split into allelic primitives for fair comparison between software.
single nucleotide variant counts , insertion and deletion counts , and multi-allelic site counts  identified per animal  using single sample variant detection with platypus, samtools, and the unifiedgenotyper following three pre-calling approaches.
total number of single nucleotide variants , insertions and deletions , and transition/transversion ratios found using single and multi sample calling methods with haplotypecaller , platypus , platypus results after multi-nucleotide variants were split into allelic primitives , samtools , and the unifiedgenotyper  for  <dig> animals
single sample calling, combined
multi sample calling
single sample calling, combined
multi sample calling
single sample calling, combined
multi sample calling
the combined results of all  <dig> single samples represent single sample calling results.
indel realignment
indel realignment had a slight effect on the number of snvs and indels identified, with the largest effects observed in ug.
indel realignment reduced the number of snvs identified , increased the number of indels detected , and decreased the number of multi-allelic sites identified   in all samples.
indel realignment did not heavily reduce the number of snvs identified by sam compared to pl or ug, indicating that sam effectively removed false positives prior to ir with gatk .
indel realignment resulted in less snvs being identified, because spurious snvs caused by incorrect alignments in the close vicinity of real indels were no longer detected as variants.
during indel realignment, the realignertargetcreator of the gatk creates a list of regions in which indels are likely to occur depending on a set of known indels and snvs.
following this initial identification step, local realignment of reads spanning the indel occurs
<cit> .
however, if the indel is incorrect in the primary alignment , indel realignment may incorrectly realign the read segments surrounding the indel.
indel realignment therefore relies on a trustworthy set of known indels and snvs, which may not yet be available for all species.figure
2a shows the average ti/tv ratios of raw, ir and ir + bqsr for single sample results of pl, sam, ug and pl after multi-nucleotide variants were split into their allelic primitives .
the ti/tv ratios are in the same range as those in human data.
the benefit of indel realignment is most apparent in ug, where an increase in ti/tv ratio can be observed compared to raw.
for pl and sam the ti/tv ratio remained virtually unchanged after realignment.
splitting multi-nucleotide variants identified by pl into allelic primitives lowered the average ti/tv ratio, however a slight improvement after indel realignment was observed.
the use of indel realignment is therefore recommended for variant identification with ug and pl_prim, however sam and pl do not benefit markedly from additional realignment with gatk.figure  <dig>
average transition/transversion ratios over all animals using single sample variant identification  and transition/transversion ratios for variant identification with single and multi sample detection methods, as well as combined over all multi sample detection methods .
average transition/transversion ratios for variant identification with single sample detection methods using platypus, platypus primitives, samtools, unifiedgenotyper and haplotypecaller  are shown in .
transition / transversion ratios for variant identification with single and multi sample detection methods using platypus, platypus primitives, samtools, unifiedgenotyper and haplotypecaller  and a consensus data set  are shown in .
base quality score recalibration
base quality score recalibration reduced the average number of snvs identified in pl  and ug , but not in sam  .
the number of snvs identified after indel realignment and base quality score recalibration is expected to decrease, as quality scores are initially overestimated and more snvs fall below the cut-off after recalibration.
by lowering base quality scores through recalibration, confidence in weak variants should decrease and, as a result, the number of false positives is expected to drop
<cit> .
the observed increase in snvs identified with sam after base quality score recalibration likely corresponds with the effective removal of false positive snvs by sam prior to ir and bqsr and the less stringent default quality scores in sam, though the increase is minimal.
recalibration of base quality scores slightly increased the number of indels for pl and sam .the effect of ir + bqsr on the number of snvs identified was more obvious than the effect of ir alone for pl and sam.
depending on variant detection software used, the total average number of snvs identified per animal increased  or decreased  by around  <dig>  –  <dig> %.
in contrast, the number of indels identified after ir + bqsr did not change markedly compared to the number of indels identified after ir alone.
the ti/tv ratio decreased slightly for pl_prim, sam and ug  after bqsr, indicating that bqsr using default settings and the current available resources may actually decrease overall variant quality for these variant identification methods.
liu et al.
<cit>  analysed the effect of preparatory steps in whole exome sequencing data and found no clear effect of indel realignment or base quality score recalibration in five whole exome sequencing samples of breast cancer patients.
the authors state that the relative contribution of each preparatory step to the accuracy of variant identification is dependent on read depth; the lack of a sufficient number of reads in a low-coverage scenario limits the power of local multiple sequence alignment.
in contrast, local realignment can benefit from consistent alignment among reads when coverage is high, thus effectively reducing the number of false positive snvs.
li
<cit>  focused on deep illumina sequencing data from two human cell lines  and found “only” a  <dig> % difference in the number of variants before and after indel realignment and base quality score recalibration using both sam and gatk.
although the authors regarded this difference as negligible considering the increased computational costs, a  <dig> % difference in the number of variants may represent a good proportion of false positives being eliminated.
unfortunately, the resulting variants were not further analysed for quality  leaving the key question regarding the effect of indel realignment or base quality score recalibration unanswered.
both indel realignment and base quality score recalibration rely on a reference set of high-quality known indels and snvs.
existing resources for human sequence  provide qualitatively solid references, whereas the quality of bovine resources such as dbsnp
<cit>  is notably lower.
although there are currently close to  <dig> million bovine snvs included in dbsnp, only very few of them are validated .
in contrast, approximately  <dig> million human snv are validated.
high-density array information, such as that from the illumina  or affymetrix , provides higher-quality information, but for only a limited number of snvs .
as the number and quality of known indels and snvs in the bovine genome increase and reference information improves, we can expect better and more dependable effects of indel realignment and base quality score recalibration in variant identification pipelines.
multi sample variant detection
the differences between single and multi sample variant identification methods were generally slight in terms of number of snvs and ti/tv ratios, and depended on software.
results of multi sample variant detection  are shown in table
<dig>  variant detection with ug  resulted in the highest number of variants , followed by hc  and sam .
as observed in single sample variant identification, pl originally had the lowest number of variants , because some variants were “hidden” in multi-snv replacements as described previously.
after splitting multi-nucleotide variants into their allelic primitives, the number of variants identified with pl increased with pl_prim to  <dig> , <dig>  le roex et al.
<cit>  compared the number of snvs identified with sam and gatk in african buffalo and identified considerably more snvs with gatk than with sam using multi sample variant identification methods.
though not as pronounced, this agrees with both our single sample and multi sample results.
the total number of snvs identified by combining single sample results of all  <dig> animals was higher for pl and sam than when multi sample variant identification was carried out on all  <dig> animals simultaneously .
this was not the case for the ug, although the difference was very slight.
similarly, liu et al.
<cit>  found that the ug multi sample pipeline resulted in  <dig> % more raw snvs than single sample results, although they found no difference in the number of snvs found in single and multi sample sam pipelines.in terms of variant quality, ti/tv ratios for single sample calling with pl, pl_prim and ug were higher than those observed in multi sample calling , whereby the ti/tv ratio for single sample calling with sam was lower than that for multi sample calling .
this difference was most prominent in ug, in which the ti/tv ratio dropped by approximately 6% when multi sample calling was applied.
surprisingly, the total number of indels identified using single sample variant detection methods was also higher for pl  and sam  than when multi sample methods were used.
again, this was not the case for the ug .
in contrast, liu et al.
<cit>  analysed human data and found that multi sample analysis increased the number of indels identified considerably  compared to single sample variant detection.
it is possible that the number of indels is inherently lower in cattle populations because their effective population size is much smaller than it is in humans, however a more likely reason for this discrepancy could also be the quality of the reference genome.
the percentage of variants identified as indels in the human genome has been estimated at up to 18%
<cit>  whereas the number of indels in cattle has been estimated at only  <dig> %
<cit>  of the total variants identified, although this may also be only a difference in reference quality.
depending on software and methods used, we found  <dig>  –  <dig> % of the variation observed was due to indels.
consensus snvs and indels for multi sample results may provide a simple method to ensure higher variant accuracy, although this approach is computationally intensive.
analysis of variant sets produced in phase one of the human  <dig> genomes project
<cit>  showed that a consensus approach to identifying variants led to a higher quality data set; variants identified in all software applications were more accurate than those identified in any single variant set
<cit> .
li
<cit>  also recommended taking the intersection of raw variants from independent variant identification methods and applying a software-oblivious filter to derive a final variant set.
we used the vcf-intersect command of vcftools
<cit>  to count consensus snvs and indels .
in this study,  <dig>  million snv positions  were found by all four multi sample variant detection methods , and  <dig>  million indels were found by all four multi sample variant detection methods .
a consensus vcf file was created using the combinevariants walker of the gatk
<cit>  in which only snv found in all multi sample vcf files  were included.
the ti/tv ratio in the consensus vcf was higher than the ti/tv ratios in the individual single and multi sample variant sets , indicating higher variant quality in the consensus variant set.
this approach is, however, extremally computationally intensive, as variant identification must be conducted using multiple software applications and then combined.
an alternative could be to use only two variant sets , which also results in improved ti/tv ratios , but is computationally less intensive than including all possible variant sets.figure  <dig>
consensus single nucleotide variants  and insertions and deletions  identified using multi sample variant detection methods.
consensus single nucleotide variants  and insertions and deletions  identified from whole genome sequencing data using four multi sample variant detection methods .
by using default parameters, we did not fine-tune all possible options available in the individual software applications.
nevertheless, using default settings in both single and multi sample variant identification yielded good performance while maintaining output quality.
our goal was to provide an initial overview of methods using the default settings recommended; it should be noted that each dataset must be treated uniquely and alternative parameter settings may deliver more accurate results.
if possible, we recommend a consensus approach for variant identification using snvs identified by all software, which resulted in the highest snv quality and should be considered the “golden standard” for variant identification in organisms with lower-quality reference genomes.
if computational constraints do not allow a consensus approach to variant identification, the tradeoff between quality and quantity of snvs must be faced .
the ug identified the highest number of snv in both single and multi sample methods, however the number of false positive snvs was also highest.
snvs identified with pl had the highest quality of single sample methods, however the number of snv identified by pl may appear to be low if “hidden” snvs are not split into allelic primitives.
sam identified a good number of snvs, which were of comparable quality to those identified by pl.
both pl and sam are likely a good choice of software for organisms with lower-quality reference genomes, as the built-in indel realignment algorithm seems to efficiently remove false positives, making the use of lower quality resources  superfluous.
computation time
average computation time required for ir and ir + bqsr is depicted in figure
4; indel realignment was considerably faster than base quality score recalibration.
the time required for single sample variant detection varied considerably between software applications, whereby pl was fastest, followed by sam and ug, and hc was slowest.
to compare time required for single and multi sample variant detection, the time required for several single sample runs was summed and compared to multi sample runs for the same number of samples .
surprisingly, the difference in runtime between single and multi sample methods for the same number of samples was negligible with the exception of pl, in which a clear speed advantage of single over multi sample variant detection was observed.
as expected, increasing the length of the chromosomal region and the number of samples analysed also linearly increased computation time , except for ms hc.
this may have been caused by limited resources, which can cause unexpected behaviour.
for pl, sam and ug pipelines including indel realignment and base quality score recalibration, more than half of the computing time was for preparatory steps.
in pipelines involving hc, which was markedly slower than the other software, variant identification required more than double the amount of time necessary for preparatory steps.figure  <dig>
average per-sample wall clock computation time required for common preparatory steps indel realignment and base quality score recalibration .
wall clock computation time required for variant identification using platypus, haplotypecaller, samtools and unifiedgenotyper on a chromosomal region 5 mb in length with single  or multi  sample variant identification methods and varying numbers of samples .
average wall clock computation time required for multi sample variant identification with varying numbers of samples  and different lengths of chromosomal regions  using different software .
indel realignment and base quality score recalibration had only slight effects on the number of snvs, indels and multi-allelic sites identified.
the effect of indel realignment on ti/tv ratio was only positive for ug, and the effect of base quality score recalibration on ti/tv ratio was negligible  or even slightly negative .
given that computational costs in terms of time are very high, we recommend indel realignment only in combination with ug.
the use of bqsr for organisms with lower-quality resource information seems superfluous until better resources become available.
concordance
concordance, measured as non-reference sensitivity , non-reference discrepancy , snv concordance and genotype concordance, was calculated by comparing variants identified in ngs data to array information from the illumina bovinesnp <dig> v <dig> dna analysis beadchip®  and the illumina bovinehd beadchip® .
detailed results of ngs concordance with the illumina bovinesnp <dig> v <dig> dna analysis beadchip® are shown in additional file
3: table s3; results of ngs concordance with the illumina bovinehd beadchip® are given in additional file
4: table s <dig>  in this section we discuss concordance results with the high-density array .
an nrs of unity represents perfect concordance between the ngs variant set and the array.
jansen et al.
<cit>  compared nrs and nrd of array-derived genotypes with sequence-derived genotypes of  <dig> fleckvieh animals for bta <dig> using sam and found that low coverage  had a negative effect on both of these parameters.
in our study, the nrs was generally higher in multi sample methods than in single sample methods for all software , however this effect was most pronounced in pl and least pronounced in sam.
our nrs results for sam are similar to those of jansen et al.
<cit> .
liu et al.
<cit>  compared sensitivity of single and multi sample methods using whole exome sequence data of  <dig> individuals and observed only a slight improvement in sensitivity when multi sample methods were applied, with the exception of sam, in which a considerable drop  in sensitivity was observed.
for ug, liu et al.
<cit>  observed an increase in sensitivity of around 1%, whereas our results showed a slightly more pronounced improvement of nrs when multi sample methods were applied .
surprisingly, cheng et al.
<cit>  found slightly better sensitivity in single sample results of sam and ug compared to multi sample results in a population-based sample of  <dig> southeast asian malays with deep whole genome sequence information.figure  <dig>
non-reference sensitivity  and non-reference discrepancy  for single nucleotide variants identified using platypus primitives, samtools, unifiedgenotyper and haplotype caller  using variants identified with the illumina bovinehd beadchip® as a gold standard.
indel realignment and base quality score recalibration were conducted for both single and multi sample calling results.
nrd results for sam are slightly lower than those of jansen et al.
<cit> , which can be explained by the slightly higher coverage in our study.
the nrd is a measure of false positives; its importance depends on the purpose of the study .
generally, ratios should be as close to zero as possible.
figure
7b shows that single sample variant identification resulted in marginally lower nrd values than those observed using multi sample methods.
though slight, the trend is apparent in all software tested.
liu et al.
<cit>  found a higher number of false positive snvs in multi sample results , whereby this observation was most pronounced in sam.
cheng et al.
<cit>  observed that the number of false positives decreased with increasing read depth, while ug showed the lowest false positive rate of all tested software.the largest difference between single and multi sample methods was observed in snv concordance .
because homozygous reference loci show no difference to the reference genome, snv concordance alone is somewhat misleading, as such loci are inherently not identified .
this has a visible effect on the overall average concordance.
snv concordance by array genotype is therefore a better measure .
for homozygote reference loci, single sample methods provide no information whatsoever .
multi sample methods identify some of the homozygous reference loci as such, however concordance reached only  <dig> - <dig> %, which indicated that a sample size of n =  <dig> is too small to have all the loci covered .
once again, the “emit all” option of the gatk could be used, but the computational cost in a multi sample setting is even higher than that in single sample analysis.
there is also visible discrepancy between single and multi sample results with respect to heterozygous loci, though this discrepancy is considerably smaller than for homozygote reference loci.
sam showed the highest snv concordance for heterozygotes.
for homozygote alternative loci, snv concordance approached 100% for sam, ug and hc, with slightly lower concordance observed in pl.genotype concordance was very high in both single and multi sample results ; depending on software, genotype concordance for single sample results was between  <dig> - <dig> % and between  <dig> - <dig> % for multi sample methods.
once a polymorphic locus is determined as such, all software applications perform very well in deciphering the correct genotype.figure  <dig>
single nucleotide variant concordance  and single nucleotide variant concordance by array genotype  with variants identified using platypus primitives, samtools, unifiedgenotyper and haplotype caller  and variants identified with the illumina bovinehd beadchip® as a gold standard.
indel realignment and base quality score recalibration were conducted for both single and multi sample calling results.
genotype concordance between genotypes identified using platypus primitives, samtools, unifiedgenotyper and haplotypecaller  and genotypes identified with the illumina bovinehd beadchip® as a gold standard.
indel realignment and base quality score recalibration were conducted for both single and multi sample calling results.
the use of multi sample variant detection to identify snvs improved nrs but worsened nrd compared to single sample variant detection.
snv and genotype concordance improved when multi sample methods were applied.
this effect was most pronounced in snv concordance of homozygous reference genotypes and less pronounced in heterozygous genotypes, whereas both single and multi sample methods identified homozygote alternative genotypes equally well.
CONCLUSIONS
the objective of this study was to investigate which methods and software work best for detection of high quality genetic variants using ngs data in cattle.
we conclude that indel realignment and base quality score recalibration have only slight effects on the number and quality of variants identified with the currently available resources for cattle and are costly with respect to computation time.
the snv concordance between variants identified using ngs data and array-based data was higher for multi sample methods than for single sample calling methods, although this was due mainly to the lack of homozygous reference genotypes in single sample results.
the quality of snvs identified  using single sample methods was higher than that of multi sample calling for pl and ug and slightly lower for sam, whereby a consensus approach using results of different software generally provides the highest variant quality.
computation time for single and multi sample methods was similar when calculated on a per-sample basis.
these findings can serve as a reference for variant detection pipeline development in various organisms and help assess the value of preparatory steps in ngs pipelines for species with lower-quality reference genomes.

 practical application of genomic technologies, such as large-scale use of single nucleotide variant arrays in animal and plant breeding, has become routine in many areas of the life sciences. taking both polygenic additive effects and genomic effects into account, between and 85% of the genetic variance observed in phenotypic traits of interest in cattle can be explained solely by snv effects and the number of genotyped animals in cattle populations worldwide is increasing steadily. fueled by decreasing costs, advances in next-generation sequencing technologies enable identification of more complex forms of genetic variation, copy number variations, etc.). these advances will inevitably foster our ability to partition the genetic variance underlying traits of interest. while some applications of ngs require de novo sequencing of an individual organism, re-sequencing may also be possible if a reference genome for the species of interest is available. the translation of raw ngs reads into tangible variants via re-sequencing is a specific, delicate and computationally demanding task and comprises three steps. first, short reads of dna are aligned to an existing reference genome. second, sequence differences between the sample being sequenced and the reference genome are identified. a myriad of alignment and variant identification software programmes are available or the haplotypecaller of the genome analysis toolkit; platypus; samtools; etc.), the majority of which can be obtained free of charge. as a final step, variants are screened and filtered to remove potential false positives common in most ngs technologies. as sequencing costs decline, reference genomes are becoming available for an increasing number of organisms, including agriculturally important species such as cattle. the bos taurus reference genome umd is similar in size to the human genome and contains ~ billion base pairs, approximately 10% of which are not assigned to any chromosome. the n size can be used to compare the quality of genome assemblies of similar size: it represents contig size such that 50% of the genome is contained in contigs of length n or greater. because the n size for the umd reference genome is much shorter than that of the current human reference, the umd reference genome will not likely allow the same accuracy in alignment, variant identification and further downstream analysis as the human reference allows. nevertheless, algorithms and software developed for alignment and variant identification of human ngs data provide an excellent resource for translating ngs data of other non-human organisms, such as cattle, into genetic variants for application in genome wide association studies and genomic selection programmes. several approaches to variant identification are possible. the simplest variant detection methods identify variants on a per-sample basis, one position at a time. once a variant locus is found, the most likely genotype for that locus is determined stochastically based on a consensus of aligned reads. if multiple samples are analyzed simultaneously, an a priori likelihood of finding a variant locus given the observed data is derived, and the most likely genotype at a given position is determined. either single or multi sample variant identification methods can be implemented in the ug; and sam. more advanced haplotype-based methods incorporate the correlation between adjacent variants within the variant detection procedure. such methods use linkage disequilibrium between nearby variants to further enrich variant identification. haplotype-based methods are implemented in pl and the gatk hc. the haplotype-based variant detection approach can also be conducted in either single or multi sample settings. read chimerism, base pair tautomerisms and signal intensity issues can contribute to false positive variant detection by causing stochastic inaccuracies, general sequencing errors, and misalignments in ngs data. aside from the variant identification approach itself, a number of optional auxiliary steps have been recommended to improve the quality of ngs-derived variants. these steps are conducted before or after variant identification. the first generally recommended preparatory step is the identification of falsely duplicated reads, which reduces bias in variant detection. secondly, local realignment around single or multiple bases that are either missing in the reference or missing in the dna sequence being analyzed is also commonly recommended. realignment cleans up spurious snvs that result from misalignment of reads around known alignment gaps and helps detect false negative snvs in the near vicinity of indels. furthermore, the full alignment context is used to determine whether the reported divergence from the reference actually exists. finally, flow cell lane, machine cycle, sequencing context or other technical aspects may influence base quality scores, which help characterize the quality of the bases in the individual reads. base quality score recalibration is recommended to lower the number of falsely identified snvs and to lower false confidence in identified bases. after variant identification, further filters can be implemented according to the individual dataset under consideration. variant detection methods can be evaluated using genome simulations based on reference genomes. given that all mammalian reference genomes are incomplete and whole-genome alignment is imperfect, simulations may not provide realistic results, a difficulty that has prompted the use of real data. another common evaluation method in human studies is to compare variants from ngs data to those of array genotyping, to results from the genomes project, or to other existing resources of human sequence. resources for cattle are not as comprehensive, however low, medium and high density snv arrays are available, thus allowing an estimation of relative accuracy between software.table variant genotypes compared for concordance between the array-based and sequence based methods to determine concordance, sensitivity and discrepancy between the two assays and measures of concordance a) ngs-based genotypes array-based genotypes aa ab bb aa a b c ab d e f bb g h i – k l m b) measures of concordance array-based information from the illumina bovinehd beadchip® was considered the “gold-standard” and compared to next-generation sequencing-based variants obtained using a illumina hiseq platform with various variant identification software, where genotypes are identified as: a = homozygous reference in both ngs-based data and array-based data. b = homozygous reference in ngs-based data, but as heterozygous in array-based data. c = homozygous reference in ngs-based data, but as homozygous alternative in array-based data. d = heterozygous in ngs-based data, but as homozygous reference in array-based data. e = heterozygous in both ngs-based data and array-based data. f = heterozygous in ngs-based data, but as homozygous alternative in array-based data. g = homozygous reference in ngs-based data, but as homozygous reference in array-based data. h = heterozygous in ngs-based data and array-based data, but as heterozygous in array-based data. i = homozygous alternative in both ngs-based data and array-based data. k = not found in ngs-based data, but as homozygous reference in array-based data. l = not found in ngs-based data and array-based data, but as heterozygous in array-based data. m = not found in ngs-based data, but as homozygous alternative in array-based data.. aside from measuring concordance between ngs data and array data, the ratio of transitions to transversions can be used as a convenient diagnostic to measure the quality of ngs data. the genome-wide ti/tv ratio is reported between and in human whole-genome sequence data, whereby this ratio is higher in exomes due to the increased presence of methelated cytosine in cpg dinucleotides in exonic regions. because this bias in favour of mutations between bases with similar biochemical properties over those with dissimilar properties is dependent on both cpg and gc content of the region, the ti/tv ratio is a useful diagnostic to measure quality across the genome. the objective of this study was to investigate which methods and software work best for detection of high quality genetic variants using ngs data in cattle, with a specific focus on single nucleotide variants. using whole-genome sequence information from individuals, we a) explore the implications of preparatory steps commonly recommended in human analysis, b) compare results of single and multi sample variant detection achieved using four publicly available variant detection software programmes, c) provide a comparison of computational processing time, and d) compare accuracy and completeness of snvs identified in ngs data by comparing them to genotypes from the same individuals generated with either high- or medium-density snv arrays, as well as to analyse genome-wide ti/tv ratios. through benchmarking different variant detection methods in cattle, preliminary recommendations for variant identification in other organisms can be extrapolated. our findings can serve as a reference for choosing variant identification software and can help assess the implication of preparatory steps in ngs pipelines for species with lower-quality or unfinished reference genomes.  alignment and coverage approximately billion paired-end reads were obtained for the sequenced animals. an average of % of these reads were mapped to chromosomes of the bovine reference genome assembly umd. approximately billion duplicate reads were marked and excluded from further analysis. average coverage was reads per base; average coverage per-animal ranged from -. see additional file 1: table s for individual alignment and coverage information. single sample variant detection the results of single sample variant detection are in figure  depending on software and preparatory steps, an average of, -, snvs, – indels and – multi-allelic sites were identified per animal. detailed information on the sequencing process is given in additional file 2: table s ug identified the most snvs, followed by pl and sam. because pl identifies multi-nucleotide variants, pl originally had the lowest average number of snv per sample because some variants were “hidden” in multi-snv replacements. after splitting multi-nucleotide variants into their allelic primitives, a fair comparison was achieved and the average per-sample number of variants identified with pl_prim increased to, when single-sample results of all animals were combined,, -, unique snv were identified.figure distributions of single nucleotide variant counts, insertion and deletion counts, and multi-allelic site counts identified per animal. for platypus results, multi-nucleotide variants were split into allelic primitives for fair comparison between software. single nucleotide variant counts, insertion and deletion counts, and multi-allelic site counts identified per animal using single sample variant detection with platypus, samtools, and the unifiedgenotyper following three pre-calling approaches. total number of single nucleotide variants, insertions and deletions, and transition/transversion ratios found using single and multi sample calling methods with haplotypecaller, platypus, platypus results after multi-nucleotide variants were split into allelic primitives, samtools, and the unifiedgenotyper for animals single sample calling, combined multi sample calling single sample calling, combined multi sample calling single sample calling, combined multi sample calling the combined results of all single samples represent single sample calling results. indel realignment indel realignment had a slight effect on the number of snvs and indels identified, with the largest effects observed in ug. indel realignment reduced the number of snvs identified, increased the number of indels detected, and decreased the number of multi-allelic sites identified in all samples. indel realignment did not heavily reduce the number of snvs identified by sam compared to pl or ug, indicating that sam effectively removed false positives prior to ir with gatk. indel realignment resulted in less snvs being identified, because spurious snvs caused by incorrect alignments in the close vicinity of real indels were no longer detected as variants. during indel realignment, the realignertargetcreator of the gatk creates a list of regions in which indels are likely to occur depending on a set of known indels and snvs. following this initial identification step, local realignment of reads spanning the indel occurs. however, if the indel is incorrect in the primary alignment, indel realignment may incorrectly realign the read segments surrounding the indel. indel realignment therefore relies on a trustworthy set of known indels and snvs, which may not yet be available for all species.figure  2a shows the average ti/tv ratios of raw, ir and ir + bqsr for single sample results of pl, sam, ug and pl after multi-nucleotide variants were split into their allelic primitives. the ti/tv ratios are in the same range as those in human data. the benefit of indel realignment is most apparent in ug, where an increase in ti/tv ratio can be observed compared to raw. for pl and sam the ti/tv ratio remained virtually unchanged after realignment. splitting multi-nucleotide variants identified by pl into allelic primitives lowered the average ti/tv ratio, however a slight improvement after indel realignment was observed. the use of indel realignment is therefore recommended for variant identification with ug and pl_prim, however sam and pl do not benefit markedly from additional realignment with gatk.figure average transition/transversion ratios over all animals using single sample variant identification and transition/transversion ratios for variant identification with single and multi sample detection methods, as well as combined over all multi sample detection methods. average transition/transversion ratios for variant identification with single sample detection methods using platypus, platypus primitives, samtools, unifiedgenotyper and haplotypecaller are shown in. transition / transversion ratios for variant identification with single and multi sample detection methods using platypus, platypus primitives, samtools, unifiedgenotyper and haplotypecaller and a consensus data set are shown in. base quality score recalibration base quality score recalibration reduced the average number of snvs identified in pl and ug, but not in sam. the number of snvs identified after indel realignment and base quality score recalibration is expected to decrease, as quality scores are initially overestimated and more snvs fall below the cut-off after recalibration. by lowering base quality scores through recalibration, confidence in weak variants should decrease and, as a result, the number of false positives is expected to drop. the observed increase in snvs identified with sam after base quality score recalibration likely corresponds with the effective removal of false positive snvs by sam prior to ir and bqsr and the less stringent default quality scores in sam, though the increase is minimal. recalibration of base quality scores slightly increased the number of indels for pl and sam.the effect of ir + bqsr on the number of snvs identified was more obvious than the effect of ir alone for pl and sam. depending on variant detection software used, the total average number of snvs identified per animal increased or decreased by around – %. in contrast, the number of indels identified after ir + bqsr did not change markedly compared to the number of indels identified after ir alone. the ti/tv ratio decreased slightly for pl_prim, sam and ug after bqsr, indicating that bqsr using default settings and the current available resources may actually decrease overall variant quality for these variant identification methods. liu et al. analysed the effect of preparatory steps in whole exome sequencing data and found no clear effect of indel realignment or base quality score recalibration in five whole exome sequencing samples of breast cancer patients. the authors state that the relative contribution of each preparatory step to the accuracy of variant identification is dependent on read depth; the lack of a sufficient number of reads in a low-coverage scenario limits the power of local multiple sequence alignment. in contrast, local realignment can benefit from consistent alignment among reads when coverage is high, thus effectively reducing the number of false positive snvs. li focused on deep illumina sequencing data from two human cell lines and found “only” a % difference in the number of variants before and after indel realignment and base quality score recalibration using both sam and gatk. although the authors regarded this difference as negligible considering the increased computational costs, a % difference in the number of variants may represent a good proportion of false positives being eliminated. unfortunately, the resulting variants were not further analysed for quality leaving the key question regarding the effect of indel realignment or base quality score recalibration unanswered. both indel realignment and base quality score recalibration rely on a reference set of high-quality known indels and snvs. existing resources for human sequence provide qualitatively solid references, whereas the quality of bovine resources such as dbsnp is notably lower. although there are currently close to million bovine snvs included in dbsnp, only very few of them are validated. in contrast, approximately million human snv are validated. high-density array information, such as that from the illumina or affymetrix, provides higher-quality information, but for only a limited number of snvs. as the number and quality of known indels and snvs in the bovine genome increase and reference information improves, we can expect better and more dependable effects of indel realignment and base quality score recalibration in variant identification pipelines. multi sample variant detection the differences between single and multi sample variant identification methods were generally slight in terms of number of snvs and ti/tv ratios, and depended on software. results of multi sample variant detection are shown in table  variant detection with ug resulted in the highest number of variants, followed by hc and sam. as observed in single sample variant identification, pl originally had the lowest number of variants, because some variants were “hidden” in multi-snv replacements as described previously. after splitting multi-nucleotide variants into their allelic primitives, the number of variants identified with pl increased with pl_prim to, le roex et al. compared the number of snvs identified with sam and gatk in african buffalo and identified considerably more snvs with gatk than with sam using multi sample variant identification methods. though not as pronounced, this agrees with both our single sample and multi sample results. the total number of snvs identified by combining single sample results of all animals was higher for pl and sam than when multi sample variant identification was carried out on all animals simultaneously. this was not the case for the ug, although the difference was very slight. similarly, liu et al. found that the ug multi sample pipeline resulted in % more raw snvs than single sample results, although they found no difference in the number of snvs found in single and multi sample sam pipelines.in terms of variant quality, ti/tv ratios for single sample calling with pl, pl_prim and ug were higher than those observed in multi sample calling, whereby the ti/tv ratio for single sample calling with sam was lower than that for multi sample calling. this difference was most prominent in ug, in which the ti/tv ratio dropped by approximately 6% when multi sample calling was applied. surprisingly, the total number of indels identified using single sample variant detection methods was also higher for pl and sam than when multi sample methods were used. again, this was not the case for the ug. in contrast, liu et al. analysed human data and found that multi sample analysis increased the number of indels identified considerably compared to single sample variant detection. it is possible that the number of indels is inherently lower in cattle populations because their effective population size is much smaller than it is in humans, however a more likely reason for this discrepancy could also be the quality of the reference genome. the percentage of variants identified as indels in the human genome has been estimated at up to 18% whereas the number of indels in cattle has been estimated at only % of the total variants identified, although this may also be only a difference in reference quality. depending on software and methods used, we found – % of the variation observed was due to indels. consensus snvs and indels for multi sample results may provide a simple method to ensure higher variant accuracy, although this approach is computationally intensive. analysis of variant sets produced in phase one of the human genomes project showed that a consensus approach to identifying variants led to a higher quality data set; variants identified in all software applications were more accurate than those identified in any single variant set. li also recommended taking the intersection of raw variants from independent variant identification methods and applying a software-oblivious filter to derive a final variant set. we used the vcf-intersect command of vcftools to count consensus snvs and indels. in this study, million snv positions were found by all four multi sample variant detection methods, and million indels were found by all four multi sample variant detection methods. a consensus vcf file was created using the combinevariants walker of the gatk in which only snv found in all multi sample vcf files were included. the ti/tv ratio in the consensus vcf was higher than the ti/tv ratios in the individual single and multi sample variant sets, indicating higher variant quality in the consensus variant set. this approach is, however, extremally computationally intensive, as variant identification must be conducted using multiple software applications and then combined. an alternative could be to use only two variant sets, which also results in improved ti/tv ratios, but is computationally less intensive than including all possible variant sets.figure consensus single nucleotide variants and insertions and deletions identified using multi sample variant detection methods. consensus single nucleotide variants and insertions and deletions identified from whole genome sequencing data using four multi sample variant detection methods. by using default parameters, we did not fine-tune all possible options available in the individual software applications. nevertheless, using default settings in both single and multi sample variant identification yielded good performance while maintaining output quality. our goal was to provide an initial overview of methods using the default settings recommended; it should be noted that each dataset must be treated uniquely and alternative parameter settings may deliver more accurate results. if possible, we recommend a consensus approach for variant identification using snvs identified by all software, which resulted in the highest snv quality and should be considered the “golden standard” for variant identification in organisms with lower-quality reference genomes. if computational constraints do not allow a consensus approach to variant identification, the tradeoff between quality and quantity of snvs must be faced. the ug identified the highest number of snv in both single and multi sample methods, however the number of false positive snvs was also highest. snvs identified with pl had the highest quality of single sample methods, however the number of snv identified by pl may appear to be low if “hidden” snvs are not split into allelic primitives. sam identified a good number of snvs, which were of comparable quality to those identified by pl. both pl and sam are likely a good choice of software for organisms with lower-quality reference genomes, as the built-in indel realignment algorithm seems to efficiently remove false positives, making the use of lower quality resources superfluous. computation time average computation time required for ir and ir + bqsr is depicted in figure  4; indel realignment was considerably faster than base quality score recalibration. the time required for single sample variant detection varied considerably between software applications, whereby pl was fastest, followed by sam and ug, and hc was slowest. to compare time required for single and multi sample variant detection, the time required for several single sample runs was summed and compared to multi sample runs for the same number of samples. surprisingly, the difference in runtime between single and multi sample methods for the same number of samples was negligible with the exception of pl, in which a clear speed advantage of single over multi sample variant detection was observed. as expected, increasing the length of the chromosomal region and the number of samples analysed also linearly increased computation time, except for ms hc. this may have been caused by limited resources, which can cause unexpected behaviour. for pl, sam and ug pipelines including indel realignment and base quality score recalibration, more than half of the computing time was for preparatory steps. in pipelines involving hc, which was markedly slower than the other software, variant identification required more than double the amount of time necessary for preparatory steps.figure average per-sample wall clock computation time required for common preparatory steps indel realignment and base quality score recalibration. wall clock computation time required for variant identification using platypus, haplotypecaller, samtools and unifiedgenotyper on a chromosomal region 5 mb in length with single or multi sample variant identification methods and varying numbers of samples. average wall clock computation time required for multi sample variant identification with varying numbers of samples and different lengths of chromosomal regions using different software. indel realignment and base quality score recalibration had only slight effects on the number of snvs, indels and multi-allelic sites identified. the effect of indel realignment on ti/tv ratio was only positive for ug, and the effect of base quality score recalibration on ti/tv ratio was negligible or even slightly negative. given that computational costs in terms of time are very high, we recommend indel realignment only in combination with ug. the use of bqsr for organisms with lower-quality resource information seems superfluous until better resources become available. concordance concordance, measured as non-reference sensitivity, non-reference discrepancy, snv concordance and genotype concordance, was calculated by comparing variants identified in ngs data to array information from the illumina bovinesnp v dna analysis beadchip® and the illumina bovinehd beadchip®. detailed results of ngs concordance with the illumina bovinesnp v dna analysis beadchip® are shown in additional file 3: table s3; results of ngs concordance with the illumina bovinehd beadchip® are given in additional file 4: table s in this section we discuss concordance results with the high-density array. an nrs of unity represents perfect concordance between the ngs variant set and the array. jansen et al. compared nrs and nrd of array-derived genotypes with sequence-derived genotypes of fleckvieh animals for bta using sam and found that low coverage had a negative effect on both of these parameters. in our study, the nrs was generally higher in multi sample methods than in single sample methods for all software, however this effect was most pronounced in pl and least pronounced in sam. our nrs results for sam are similar to those of jansen et al.. liu et al. compared sensitivity of single and multi sample methods using whole exome sequence data of individuals and observed only a slight improvement in sensitivity when multi sample methods were applied, with the exception of sam, in which a considerable drop in sensitivity was observed. for ug, liu et al. observed an increase in sensitivity of around 1%, whereas our results showed a slightly more pronounced improvement of nrs when multi sample methods were applied. surprisingly, cheng et al. found slightly better sensitivity in single sample results of sam and ug compared to multi sample results in a population-based sample of southeast asian malays with deep whole genome sequence information.figure non-reference sensitivity and non-reference discrepancy for single nucleotide variants identified using platypus primitives, samtools, unifiedgenotyper and haplotype caller using variants identified with the illumina bovinehd beadchip® as a gold standard. indel realignment and base quality score recalibration were conducted for both single and multi sample calling results. nrd results for sam are slightly lower than those of jansen et al., which can be explained by the slightly higher coverage in our study. the nrd is a measure of false positives; its importance depends on the purpose of the study. generally, ratios should be as close to zero as possible. figure  7b shows that single sample variant identification resulted in marginally lower nrd values than those observed using multi sample methods. though slight, the trend is apparent in all software tested. liu et al. found a higher number of false positive snvs in multi sample results, whereby this observation was most pronounced in sam. cheng et al. observed that the number of false positives decreased with increasing read depth, while ug showed the lowest false positive rate of all tested software.the largest difference between single and multi sample methods was observed in snv concordance. because homozygous reference loci show no difference to the reference genome, snv concordance alone is somewhat misleading, as such loci are inherently not identified. this has a visible effect on the overall average concordance. snv concordance by array genotype is therefore a better measure. for homozygote reference loci, single sample methods provide no information whatsoever. multi sample methods identify some of the homozygous reference loci as such, however concordance reached only - %, which indicated that a sample size of n =  is too small to have all the loci covered. once again, the “emit all” option of the gatk could be used, but the computational cost in a multi sample setting is even higher than that in single sample analysis. there is also visible discrepancy between single and multi sample results with respect to heterozygous loci, though this discrepancy is considerably smaller than for homozygote reference loci. sam showed the highest snv concordance for heterozygotes. for homozygote alternative loci, snv concordance approached 100% for sam, ug and hc, with slightly lower concordance observed in pl.genotype concordance was very high in both single and multi sample results; depending on software, genotype concordance for single sample results was between - % and between - % for multi sample methods. once a polymorphic locus is determined as such, all software applications perform very well in deciphering the correct genotype.figure single nucleotide variant concordance and single nucleotide variant concordance by array genotype with variants identified using platypus primitives, samtools, unifiedgenotyper and haplotype caller and variants identified with the illumina bovinehd beadchip® as a gold standard. indel realignment and base quality score recalibration were conducted for both single and multi sample calling results. genotype concordance between genotypes identified using platypus primitives, samtools, unifiedgenotyper and haplotypecaller and genotypes identified with the illumina bovinehd beadchip® as a gold standard. indel realignment and base quality score recalibration were conducted for both single and multi sample calling results. the use of multi sample variant detection to identify snvs improved nrs but worsened nrd compared to single sample variant detection. snv and genotype concordance improved when multi sample methods were applied. this effect was most pronounced in snv concordance of homozygous reference genotypes and less pronounced in heterozygous genotypes, whereas both single and multi sample methods identified homozygote alternative genotypes equally well.  the objective of this study was to investigate which methods and software work best for detection of high quality genetic variants using ngs data in cattle. we conclude that indel realignment and base quality score recalibration have only slight effects on the number and quality of variants identified with the currently available resources for cattle and are costly with respect to computation time. the snv concordance between variants identified using ngs data and array-based data was higher for multi sample methods than for single sample calling methods, although this was due mainly to the lack of homozygous reference genotypes in single sample results. the quality of snvs identified using single sample methods was higher than that of multi sample calling for pl and ug and slightly lower for sam, whereby a consensus approach using results of different software generally provides the highest variant quality. computation time for single and multi sample methods was similar when calculated on a per-sample basis. these findings can serve as a reference for variant detection pipeline development in various organisms and help assess the value of preparatory steps in ngs pipelines for species with lower-quality reference genomes. 
the average quality of snvs identified, measured as the ratio of transitions to transversions, was higher using single sample methods than multi sample methods. we explored the consequences of commonly recommended preparatory steps and the effects of single and multi sample variant identification methods using four publicly available software applications on whole genome sequence data of key ancestors of swiss dairy cattle populations. preparatory realignment of insertions and deletions and subsequent base quality score recalibration had only minor effects on the number and quality of snvs identified by different software, but increased computing time considerably. a consensus approach using results of different software generally provided the highest variant quality in terms of transition/transversion ratio. the total number of snvs identified varied by software and method, with single sample results ranging from to million variants. average concordance for single sample results with high-density chip data was % and average genotype concordance in correctly identified snvs was % across software. the challenges of variant identification in organisms with lower quality reference genomes are less well documented. our findings serve as a reference for variant identification pipeline development in non-human organisms and help assess the implication of preparatory steps in next-generation sequencing pipelines for organisms with incomplete reference genomes. next-generation sequencing analysissingle nucleotide variant identificationpipelineissue-copyright-statement© the author 2014
 advances in human genomics have allowed unprecedented productivity in terms of algorithms, software, and literature available for translating raw next-generation sequence data into high-quality information. the challenges of variant identification in organisms with lower quality reference genomes are less well documented. we explored the consequences of commonly recommended preparatory steps and the effects of single and multi sample variant identification methods using four publicly available software applications on whole genome sequence data of key ancestors of swiss dairy cattle populations. accuracy of calling next-generation sequence variants was assessed by comparison to the same loci from medium and high-density single nucleotide variant arrays.  the total number of snvs identified varied by software and method, with single sample results ranging from to million variants. computing time varied considerably between software. preparatory realignment of insertions and deletions and subsequent base quality score recalibration had only minor effects on the number and quality of snvs identified by different software, but increased computing time considerably. average concordance for single sample results with high-density chip data was % and average genotype concordance in correctly identified snvs was % across software. the average quality of snvs identified, measured as the ratio of transitions to transversions, was higher using single sample methods than multi sample methods. a consensus approach using results of different software generally provided the highest variant quality in terms of transition/transversion ratio.  our findings serve as a reference for variant identification pipeline development in non-human organisms and help assess the implication of preparatory steps in next-generation sequencing pipelines for organisms with incomplete reference genomes. benchmarking this information should prove particularly useful in processing next-generation sequencing data for use in genome-wide association studies and genomic selection. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords next-generation sequencing analysissingle nucleotide variant identificationpipelineissue-copyright-statement© the author 2014 