BACKGROUND
next-generation sequencing  technology is a powerful and cost-effective approach for large-scale dna sequencing  <cit> .
it has significantly propelled the sequence-based genetics and genomics research and its downstream applications which include, but are not limited to, de novo sequencing  <cit> , quantifying expression level s <cit> , providing a genome-scale look at transcription-factor binding  <cit> , creating a foundation for understanding human disease  <cit>  and systematically investigating of human variation  <cit> .
a number of projects based on ngs technology are underway.
for example,  <dig> genomes project http://www.1000genomes.org/ aims to provide a comprehensive resource of human genetic variation as a foundation for understanding the relationship between genotype and phenotype  <cit> .
the nhlbi go exome sequencing project  http://evs.gs.washington.edu/evs/ focuses on protein coding regions to discover novel genes and mechanisms contributing to heart, lung and blood disorders.
tcga  http://cancergenome.nih.gov/ has been sequencing a large number of tumor/normal pairs to provide insights into the landscape of somatic mutations and the great genetic heterogeneity that defines the unique signature of individual tumor  <cit> .
the ability to discover a comprehensive list of human genetic variation and to search for causing variation or mutation underlying diseases depends crucially on the accurate calling of snps and genotypes  <cit> .
translating the raw sequencing data into the final snp and genotype calls requires two essential steps: read mapping and snp/genotype inference.
first, reads are aligned onto an available reference genome, then variable sites are identified and genotypes at those sites are determined.
snp and genotype calling suffers from high error rates that are due to the following factors.
poor quality or low-quality tails prevent reads from being properly mapped.
each read is aligned independently, causing many reads that span indels are misaligned  <cit> .
the raw base-calling quality scores often co-vary with features like sequence technology, machine cycle and sequence context and, thus, cannot reflect the true base-calling error rates  <cit> .
these alignment and base-calling errors propagate into snp and genotype inference and lead to false variant detection.
moreover, low-coverage sequencing always introduces considerable uncertainty into the results and makes accurate snp and genotype calling difficult.
to obtain high quality snp and genotype data, most contemporary algorithms use a probabilistic framework to quantify the uncertainty and to model errors introduced in alignment and base calling  <cit> .
in addition, a number of optional steps are recommended.
some are prior to variant calling, including raw reads preprocessing, duplicate marking, local realignment, and base quality score recalibration <cit> .
others are posterior to variant calling, including linkage-based genotype refining  <cit>  and snp filtering  <cit>  or variant quality score recalibration  <cit> .
here we focused on those optional steps preceding variant calling.
we assessed their relative contributions and evaluated the effect of their orders on the accuracy of snp and genotype calling with data generated on illumina sequencing platform, which is currently the most widely used sequencing technology.
besides, we also compared the performance of three popular multi-sample snp callers, samtools  <cit> , gatk  <cit> , and glfmultiples  <cit> , in terms of dbsnp rate, transition to transversion ratio , and concordance rate with snp arrays .
our findings can serve as a general guide for choosing appropriate steps for snp and genotype calling from illumina sequencing data with different coverage.
methods
sequencing data and snp calling
five samples were selected for whole exome sequencing.
all samples were taken from women with very early-onset  breast cancer or early-onset  plus a first-degree family history of breast cancer  <cit> .
genomic dna from buffy coat was extracted using qiamp dna kit  following the manufacture's protocol.
exonic regions were captured using illumina truseq exome enrichment kit.
it targeted  <dig>  regions , covering  <dig> % of consensus coding sequence database .
an illumina hiseq  <dig> was used to generate 100-bp paired-end reads .
reads were mapped to the ncbi build  <dig> reference genome with bwa  <cit> , sorted and indexed with samtools  <cit> .
those reads were classified into three categories by their mapped locations on the genome, inside target regions, outside target regions with ≤  <dig> bp distance and outside target regions >  <dig> bp distance.
for these five samples, there was an average of  <dig> % bases  mapped to target regions,  <dig> %  mapped to outside ≤  <dig> bp regions, and  <dig> %  mapped to outside >  <dig> bp regions.
as expected, the depth of coverage was the highest for inside target regions  and lowest for outside >  <dig> bp regions  .
<dig> % target regions,  <dig> % of outside ≤  <dig> bp regions and  <dig> % of outside >  <dig> bp regions are accessed by sequencing data .
high
medium
low
poor-quality tails of reads were dynamically trimmed off by the bwa parameter .
duplicated reads were marked by picard.
base quality recalibration and local realignment were carried out using genome analysis toolkit   <cit> .
snps were called simultaneously on five samples by gatk unified genotyper, samtools mpileup and glfmultiples using bases with base quality≥ <dig> and reads with mapping quality ≥ <dig>
definition of performance metrics
dbsnp rate
the percentage of variants found in dbsnp database  <cit>  is used to measure an approximate false-positive rate of snp calling.
here dbsnp  <dig> was used, which contains approximately  <dig> million snp entries  <cit> .
it excludes the impact of the  <dig> genomes project and is useful for evaluation.
multi-sample snp calling is able to find more rare variants than single sample calling, thus the aggregate dbsnp rate is lower.
of ~ <dig> k variants discovered from these five samples, about 77% were already catalogued in dbsnp  <dig> .
it should be noted that dbsnp rate is not an absolute measurement of which variant calls are better, but the same number of variants with higher dbsnp rate may reasonably suggest lower false-positive rates.
call set
raw: without any preprocessing steps; filtery: removing those reads that fail the illumina chastity filter; trim: trimming off low-quality tails from reads with the bwa parameter ; filtery&trim: removing those reads that fail the illumina chastity filter and trimming off low quality tails.
snps were called for five samples together by gatk using bases with base quality≥ <dig> and reads with mapping quality ≥ <dig>  only sites with qual > =  <dig> were considered as potentially variable sites.
transition/transversion ratio
the variants are observed either as transitions  or transversions .
the ratio of the number of transitions to the number of transversions is particularly helpful for assessing the quality of snp calls  <cit> .
ti/tv ratios are often calculated for known and novel snps separately.
the expected ti/tv ratios in whole-genome sequencing are  <dig>  and  <dig>  for known and novel variants, respectively, and in the exome target regions are  <dig>  and  <dig> , respectively  <cit> .
the higher ti/tv ratio generally indicates higher accuracy.
when detected variants demonstrate a ratio closer to the expected ratio for random substitutions , low-quality variant calling or data is implied.
genotype concordance
all five samples have been genotyped using the affymetrix snp  <dig>  array in a previous genome-wide association study  <cit> .
detailed genotyping methods and stringent quality control criteria were described in zheng et al.,  <cit> .
the original scan included three quality control samples in each 96-well plate, and the snp calls showed a very high concordance rate  for the quality control samples.
genotypes obtained from the sequencing data were compared with those from the snp array.
the non-reference discrepancy  rate was used to measure the accuracy of genotype calls, which reported the percent of discordant genotype calls at commonly called on-reference sites on the snp array and exome-sequencing.
the mathematical definition of nrd can be found in depristo et al.,  <cit> .
the lower nrd generally indicates higher accuracy of genotype calls.
RESULTS
effects of data preprocessing
using high-quality reads is expected to identify true variants.
generally, there are two ways to extract high-quality reads from illumina sequencing data: removing reads that fail the illumina chastity filter  and trimming off low-quality ends from reads .
the trim step obtained the largest number of mapped reads, while the filtery produced the fewest number of mapped reads resulting from lots of low-quality reads being discarded .
although the trim step helped align more reads and identify slightly more variants , it obtained a lower dbsnp rate  and a lower novel transition/transversion ratio   compared with those using raw sequencing data  .
trimming low-quality tails added  <dig>  novel variants, representing about 8% of all novel calls, with a ti/tv ratio of  <dig> , while it eliminated  <dig>  novel variants with a ti/tv ratio of  <dig>  from the raw call set .
the novel variants unique to the trim call set had a much lower ti/tv ratio  compared with the ti/tv ratio  of those unique to the raw call set, which suggested that more false positive variants were introduced by the trim step.
results from applying both filtery and trim steps  compared with those from performing filtery step alone also revealed that trim step would increase the number of false positives .
the filtery step identified fewer variants ; however, those variants showed the similar dbsnp rate  and ti/tv ratio  compared with the raw call set.
removing poor-quality reads from raw data  added  <dig> known variants with a ti/tv ratio of  <dig> , while it eliminated  <dig> known variants with a ti/tv ratio of  <dig>  from the raw call set .
that is, filtery step dropped more than  <dig>  known variants, representing about 2% of all known calls.
these results suggested that throwing out those poor quality reads which failed the chastity filter might not be necessary for further snp calling.
comparison results from applying both filtery and trim steps  with those from performing trim step alone also revealed the useless of filtery step on improving snp calling performance .
a comprehensive comparison using variable quality thresholds for high-coverage data , medium-coverage data  and low-coverage data  came to the same conclusion, that these two preprocessing step, filtery and trim, could not improve the performance of snp calling, a conclusion contrary to the usual expectation.
application of the trim step might even introduce false positives, especially for high-coverage data.
compared with low coverage data, the problem of introducing false positives caused by the trim step is more serious for high coverage data .
effects of duplicate marking, realignment and recalibration
among the three optional steps, local realignment, marking duplication and base quality recalibration, local realignment obtained the highest dbsnp rate  and novel ti/tv ratio  for high-coverage data  .
local realignment eliminated  <dig> novel variants from the initial call set, representing more than 7% of all novel calls, with a ti/tv ratio of  <dig> , which indicated that about 90% of these novel calls were false-positives .
in contrast, base quality recalibration eliminated only  <dig> novel variants with a ti/tv ratio of  <dig>  but added  <dig> novel variants with a ti/tv ratio of  <dig>  from the initial call set .
marking duplication removed  <dig> novel variants with a ti/tv ratio of  <dig>  but it added  <dig> novel variants with a ti/tv ratio of  <dig>  from the initial call set .
these results suggested that local realignment was efficient in reducing the false-positive rate, while the effect of recalibration and marking duplications was limited for deep-sequencing data.
snps were called for  <dig> samples together by gatk using bases with base quality≥ <dig> and reads with mapping quality ≥ <dig>  only sites with qual >  <dig> for deep-coverage or qual >  <dig> for shallow coverage were considered as potentially variable sites.
for low-coverage sequencing , however, the ability of these three steps to eliminate false-positive variants changed.
marking duplication obtained the highest performance with  <dig> % dbsnp rate and a novel ti/tv ratio of  <dig>  .
marking duplication removed  <dig> novel variants from the initial call set, representing more than 10% of all novel calls, with a ti/tv ratio of  <dig>  .
in contrast, local realignment only eliminated  <dig> novel variants with a ti/tv ratio of  <dig>   and recalibration only removed  <dig> novel variants with a ti/tv ratio of  <dig>  .
these results suggested that marking duplication was more efficient in reducing false-positive rates than other two optional steps for low-coverage sequencing data.
a comprehensive comparison using variable quality thresholds also suggested that realignment was more efficient in removing false positives than base call recalibration and marking duplication for high-coverage data, whereas marking duplication was more efficient than the other two for low-coverage data .
the effect of orders of the optional steps on snp calling was also evaluated.
we obtained the same accuracy of snp and genotype calling using different order arrangements, suggesting that the order of steps had no effect on the calling performance .
comparing the performance of gatk, samtools and glfmultiples
samtools and gatk obtained higher known and novel ti/tv ratios than glfmultiples for deep-sequencing data , while they produced a lower dbsnp rate and known and novel ti/tv ratios than glfmultiples for low-sequencing data  when the same number of snps were identified .
for those data with medium-coverage, these three multi-sample calling tools produced similar dbsnp rate, known and novel ti/tv ratios .
all of these three tools produced a similar genotype concordance with snp chip data for all regions .
these results suggested that samtools and gatk had better performance than glfmultiples for high-coverage data, while glfmultiples were superior to samtools and gatk for low-coverage data.
discussion
intriguingly, we found that the read preprocessing steps before mapping were not necessary.
trimming off low-quality tails from reads even worsen the power of variant calling, although it helps align more reads with high error rate in the tail.
a possible explanation is that although the quality of tails is not good enough, they are still helpful for reads mapping.
thus trimming off low-quality tails would lead to more alignment artifacts than using raw reads and, in turn, cause false-positive variants discovery.
it should be noted that trimming reads is somehow a question of trial and error and a balance between the number of mapped reads and mapping accuracy.
if the decrease of the quality of the 3' end is acceptable and the loss of coverage is affordable, trimming is not necessary.
in contrast, if there is a dramatic quality decrease at the tail and poor quality was observed at very earlier sequencing cycle, trimming might be helpful by increasing the number of mapped reads greatly but without reducing the mapping accuracy much.
for the steps after read mapping, including marking duplication, realignment and recalibration, the relative contribution of each step to the accuracy of variant calling depends on the sequencing depth.
when the sequencing depth is high, read mapping can benefit from finding consistent alignment among all reads and thus reduce the number of false-positives effectively.
when the sequencing depth is low, however, the lack of sufficient reads mapping to the locus limits the power of local multiple sequence alignment and thus it cannot improve the quality of variant calls much.
in such circumstances, marking duplication plays a more important role in reducing false positives than realignment and recalibration.
moreover, the performances of three popular multi-sample calling tools, samtools, gatk and glfmultiples, also depend on the sequencing depth.
they use the same genotype likelihood model, but glfmultiples not only takes into account the maximized likelihood but also an overall prior for each type of polymorphism.
for example, they favor sites with transition polymorphisms over those with transversion  <cit> .
thus, incorporating such additional information helps reduce the uncertainty associated with shallow-sequencing data.
however, the additional information will disturb the identification of variants when enough evidence is already involved with deep-sequencing data.
the steps posterior to variant calling, including linkage-based genotype refining and snp filtering or variant quality score recalibration, also contribute a lot to the accurate snp and genotype calling.
the use of ld  patterns can substantially improve genotype calling when multiple samples have been sequenced  <cit> .
because not all information regarding errors can be fully incorporated into the statistical framework, the proper snp filtering strategies are recommended to reduce the error rates  <cit> .
besides, the consensus of multiple call sets from different methods provide higher quality than any of individual call sets  <cit> .
even with the best pipelines, however, we are still far from obtaining a complete and accurate picture of snps and genotypes in the human genome.
the most challenging task is to distinguish rare variants from sequencing errors.
snp and genotype calling for rare variants, which would not be represented in any reference panel, may not improve much by the use of ld information.
to identify rare variants, a direct and more powerful approach is to sequence a large number of individuals  <cit> .
in addition to using the proper sequencing strategies, developing more accurate snp detection methods is needed.
more research is also needed in other areas, including longer read depths, improved protocols for generate paired ends, advances in sequencing technology with lower base calling error rates, and more powerful alignment methods.
CONCLUSIONS
here, we evaluated the effect of a number of computational steps on the accuracy of snp and genotype calling from illumina sequencing data with different coverage.
to our knowledge, no other study has made a systematic assessment of whether each step is valuable and how it affects the quality of variant detection.
our findings can serves as the general guideline for choosing snp calling strategies.
competing interests
the authors declare that they have no competing interests.
authors' contributions
ys led the project and oversaw the analysis.
ql and yg designed and performed the research.
jl participated in the data analysis.
jrl guided the experiments and provided the ngs data.
bz guided the analysis and revised the manuscript.
ql wrote the manuscript.
all authors have read and approved of the final manuscript.
supplementary material
additional file 1
comparison of effect of different preprocessing steps.
a detailed comparison of calling results with different preprocessing steps in terms of dbsnp rate, ti/tv ratio, novel ti/tv ratio and nrd for all regions, inside target regions, outside ≤  <dig> bp regions, and outside >  <dig> bp regions from illumina whole-exome sequencing data.
raw , filtery , trim  and filtery&trim .
click here for file
additional file 2
comparison of effect of marking duplication, realignment and recalibration.
a detailed comparison of results using different steps, marking duplication, realignment and recalibration, in terms of dbsnp rate, ti/tv ratio, novel ti/tv ratio and nrd for all regions, inside target regions, outside ≤  <dig> bp regions, and outside >  <dig> bp regions from illumina whole-exome sequencing data.
initial alignment , marking duplication , realignment , recalibration , marking duplication followed by realignment , marking duplication followed by realignment and recalibration .
click here for file
additional file 3
comparison of effect of different arrangements of marking duplication, realignment and recalibration.
a detailed comparison of results by arranging three steps, marking duplication, realignment and recalibration, in different orders in terms of dbsnp rate, ti/tv ratio, novel ti/tv ratio and nrd for all regions, inside target regions, outside ≤  <dig> bp regions, and outside >  <dig> bp regions from illumina whole-exome sequencing data.
marking duplication followed by realignment and recalibration , marking duplication followed by recalibration and realignment , realignment followed by recalibration and marking duplication .
click here for file
acknowledgements
the authors wish to thank peggy schuyler for editorial work on this manuscript and wei zheng for his support.
this work was supported by national cancer institute grants u <dig> ca <dig>  p <dig> ca <dig>  p <dig> ca <dig>  p <dig> ca <dig> and p <dig> ca <dig>  and the national institutes of health grants r01gm <dig> .
subject recruitment and exome sequencing is supported by ca <dig>  and ca <dig> .
ql's work was partially supported by the national natural science foundation of china  <dig> .
this article has been published as part of bmc genomics volume  <dig> supplement  <dig>  2012: proceedings of the international conference on intelligent biology and medicine : genomics.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/13/s <dig>

 next-generation sequencing technology is a powerful and cost-effective approach for large-scale dna sequencing. it has significantly propelled the sequence-based genetics and genomics research and its downstream applications which include, but are not limited to, de novo sequencing, quantifying expression level s, providing a genome-scale look at transcription-factor binding, creating a foundation for understanding human disease and systematically investigating of human variation. a number of projects based on ngs technology are underway. for example, genomes project http://www.1000genomes.org/ aims to provide a comprehensive resource of human genetic variation as a foundation for understanding the relationship between genotype and phenotype. the nhlbi go exome sequencing project http://evs.gs.washington.edu/evs/ focuses on protein coding regions to discover novel genes and mechanisms contributing to heart, lung and blood disorders. tcga http://cancergenome.nih.gov/ has been sequencing a large number of tumor/normal pairs to provide insights into the landscape of somatic mutations and the great genetic heterogeneity that defines the unique signature of individual tumor. the ability to discover a comprehensive list of human genetic variation and to search for causing variation or mutation underlying diseases depends crucially on the accurate calling of snps and genotypes. translating the raw sequencing data into the final snp and genotype calls requires two essential steps: read mapping and snp/genotype inference. first, reads are aligned onto an available reference genome, then variable sites are identified and genotypes at those sites are determined. snp and genotype calling suffers from high error rates that are due to the following factors. poor quality or low-quality tails prevent reads from being properly mapped. each read is aligned independently, causing many reads that span indels are misaligned. the raw base-calling quality scores often co-vary with features like sequence technology, machine cycle and sequence context and, thus, cannot reflect the true base-calling error rates. these alignment and base-calling errors propagate into snp and genotype inference and lead to false variant detection. moreover, low-coverage sequencing always introduces considerable uncertainty into the results and makes accurate snp and genotype calling difficult. to obtain high quality snp and genotype data, most contemporary algorithms use a probabilistic framework to quantify the uncertainty and to model errors introduced in alignment and base calling. in addition, a number of optional steps are recommended. some are prior to variant calling, including raw reads preprocessing, duplicate marking, local realignment, and base quality score recalibration. others are posterior to variant calling, including linkage-based genotype refining and snp filtering or variant quality score recalibration. here we focused on those optional steps preceding variant calling. we assessed their relative contributions and evaluated the effect of their orders on the accuracy of snp and genotype calling with data generated on illumina sequencing platform, which is currently the most widely used sequencing technology. besides, we also compared the performance of three popular multi-sample snp callers, samtools, gatk, and glfmultiples, in terms of dbsnp rate, transition to transversion ratio, and concordance rate with snp arrays. our findings can serve as a general guide for choosing appropriate steps for snp and genotype calling from illumina sequencing data with different coverage. methods sequencing data and snp calling five samples were selected for whole exome sequencing. all samples were taken from women with very early-onset breast cancer or early-onset plus a first-degree family history of breast cancer. genomic dna from buffy coat was extracted using qiamp dna kit following the manufacture's protocol. exonic regions were captured using illumina truseq exome enrichment kit. it targeted regions, covering % of consensus coding sequence database. an illumina hiseq was used to generate 100-bp paired-end reads. reads were mapped to the ncbi build reference genome with bwa, sorted and indexed with samtools. those reads were classified into three categories by their mapped locations on the genome, inside target regions, outside target regions with ≤ bp distance and outside target regions > bp distance. for these five samples, there was an average of % bases mapped to target regions, % mapped to outside ≤ bp regions, and % mapped to outside > bp regions. as expected, the depth of coverage was the highest for inside target regions and lowest for outside > bp regions. % target regions, % of outside ≤ bp regions and % of outside > bp regions are accessed by sequencing data. high medium low poor-quality tails of reads were dynamically trimmed off by the bwa parameter. duplicated reads were marked by picard. base quality recalibration and local realignment were carried out using genome analysis toolkit. snps were called simultaneously on five samples by gatk unified genotyper, samtools mpileup and glfmultiples using bases with base quality≥ and reads with mapping quality ≥ definition of performance metrics dbsnp rate the percentage of variants found in dbsnp database is used to measure an approximate false-positive rate of snp calling. here dbsnp was used, which contains approximately million snp entries. it excludes the impact of the genomes project and is useful for evaluation. multi-sample snp calling is able to find more rare variants than single sample calling, thus the aggregate dbsnp rate is lower. of ~ k variants discovered from these five samples, about 77% were already catalogued in dbsnp. it should be noted that dbsnp rate is not an absolute measurement of which variant calls are better, but the same number of variants with higher dbsnp rate may reasonably suggest lower false-positive rates. call set raw: without any preprocessing steps; filtery: removing those reads that fail the illumina chastity filter; trim: trimming off low-quality tails from reads with the bwa parameter; filtery&trim: removing those reads that fail the illumina chastity filter and trimming off low quality tails. snps were called for five samples together by gatk using bases with base quality≥ and reads with mapping quality ≥ only sites with qual > = were considered as potentially variable sites. transition/transversion ratio the variants are observed either as transitions or transversions. the ratio of the number of transitions to the number of transversions is particularly helpful for assessing the quality of snp calls. ti/tv ratios are often calculated for known and novel snps separately. the expected ti/tv ratios in whole-genome sequencing are and for known and novel variants, respectively, and in the exome target regions are and, respectively. the higher ti/tv ratio generally indicates higher accuracy. when detected variants demonstrate a ratio closer to the expected ratio for random substitutions, low-quality variant calling or data is implied. genotype concordance all five samples have been genotyped using the affymetrix snp array in a previous genome-wide association study. detailed genotyping methods and stringent quality control criteria were described in zheng et al.,. the original scan included three quality control samples in each 96-well plate, and the snp calls showed a very high concordance rate for the quality control samples. genotypes obtained from the sequencing data were compared with those from the snp array. the non-reference discrepancy rate was used to measure the accuracy of genotype calls, which reported the percent of discordant genotype calls at commonly called on-reference sites on the snp array and exome-sequencing. the mathematical definition of nrd can be found in depristo et al.,. the lower nrd generally indicates higher accuracy of genotype calls.  effects of data preprocessing using high-quality reads is expected to identify true variants. generally, there are two ways to extract high-quality reads from illumina sequencing data: removing reads that fail the illumina chastity filter and trimming off low-quality ends from reads. the trim step obtained the largest number of mapped reads, while the filtery produced the fewest number of mapped reads resulting from lots of low-quality reads being discarded. although the trim step helped align more reads and identify slightly more variants, it obtained a lower dbsnp rate and a lower novel transition/transversion ratio compared with those using raw sequencing data. trimming low-quality tails added novel variants, representing about 8% of all novel calls, with a ti/tv ratio of, while it eliminated novel variants with a ti/tv ratio of from the raw call set. the novel variants unique to the trim call set had a much lower ti/tv ratio compared with the ti/tv ratio of those unique to the raw call set, which suggested that more false positive variants were introduced by the trim step. results from applying both filtery and trim steps compared with those from performing filtery step alone also revealed that trim step would increase the number of false positives. the filtery step identified fewer variants; however, those variants showed the similar dbsnp rate and ti/tv ratio compared with the raw call set. removing poor-quality reads from raw data added known variants with a ti/tv ratio of, while it eliminated known variants with a ti/tv ratio of from the raw call set. that is, filtery step dropped more than known variants, representing about 2% of all known calls. these results suggested that throwing out those poor quality reads which failed the chastity filter might not be necessary for further snp calling. comparison results from applying both filtery and trim steps with those from performing trim step alone also revealed the useless of filtery step on improving snp calling performance. a comprehensive comparison using variable quality thresholds for high-coverage data, medium-coverage data and low-coverage data came to the same conclusion, that these two preprocessing step, filtery and trim, could not improve the performance of snp calling, a conclusion contrary to the usual expectation. application of the trim step might even introduce false positives, especially for high-coverage data. compared with low coverage data, the problem of introducing false positives caused by the trim step is more serious for high coverage data. effects of duplicate marking, realignment and recalibration among the three optional steps, local realignment, marking duplication and base quality recalibration, local realignment obtained the highest dbsnp rate and novel ti/tv ratio for high-coverage data. local realignment eliminated novel variants from the initial call set, representing more than 7% of all novel calls, with a ti/tv ratio of, which indicated that about 90% of these novel calls were false-positives. in contrast, base quality recalibration eliminated only novel variants with a ti/tv ratio of but added novel variants with a ti/tv ratio of from the initial call set. marking duplication removed novel variants with a ti/tv ratio of but it added novel variants with a ti/tv ratio of from the initial call set. these results suggested that local realignment was efficient in reducing the false-positive rate, while the effect of recalibration and marking duplications was limited for deep-sequencing data. snps were called for samples together by gatk using bases with base quality≥ and reads with mapping quality ≥ only sites with qual > for deep-coverage or qual > for shallow coverage were considered as potentially variable sites. for low-coverage sequencing, however, the ability of these three steps to eliminate false-positive variants changed. marking duplication obtained the highest performance with % dbsnp rate and a novel ti/tv ratio of. marking duplication removed novel variants from the initial call set, representing more than 10% of all novel calls, with a ti/tv ratio of. in contrast, local realignment only eliminated novel variants with a ti/tv ratio of and recalibration only removed novel variants with a ti/tv ratio of. these results suggested that marking duplication was more efficient in reducing false-positive rates than other two optional steps for low-coverage sequencing data. a comprehensive comparison using variable quality thresholds also suggested that realignment was more efficient in removing false positives than base call recalibration and marking duplication for high-coverage data, whereas marking duplication was more efficient than the other two for low-coverage data. the effect of orders of the optional steps on snp calling was also evaluated. we obtained the same accuracy of snp and genotype calling using different order arrangements, suggesting that the order of steps had no effect on the calling performance. comparing the performance of gatk, samtools and glfmultiples samtools and gatk obtained higher known and novel ti/tv ratios than glfmultiples for deep-sequencing data, while they produced a lower dbsnp rate and known and novel ti/tv ratios than glfmultiples for low-sequencing data when the same number of snps were identified. for those data with medium-coverage, these three multi-sample calling tools produced similar dbsnp rate, known and novel ti/tv ratios. all of these three tools produced a similar genotype concordance with snp chip data for all regions. these results suggested that samtools and gatk had better performance than glfmultiples for high-coverage data, while glfmultiples were superior to samtools and gatk for low-coverage data. discussion intriguingly, we found that the read preprocessing steps before mapping were not necessary. trimming off low-quality tails from reads even worsen the power of variant calling, although it helps align more reads with high error rate in the tail. a possible explanation is that although the quality of tails is not good enough, they are still helpful for reads mapping. thus trimming off low-quality tails would lead to more alignment artifacts than using raw reads and, in turn, cause false-positive variants discovery. it should be noted that trimming reads is somehow a question of trial and error and a balance between the number of mapped reads and mapping accuracy. if the decrease of the quality of the 3' end is acceptable and the loss of coverage is affordable, trimming is not necessary. in contrast, if there is a dramatic quality decrease at the tail and poor quality was observed at very earlier sequencing cycle, trimming might be helpful by increasing the number of mapped reads greatly but without reducing the mapping accuracy much. for the steps after read mapping, including marking duplication, realignment and recalibration, the relative contribution of each step to the accuracy of variant calling depends on the sequencing depth. when the sequencing depth is high, read mapping can benefit from finding consistent alignment among all reads and thus reduce the number of false-positives effectively. when the sequencing depth is low, however, the lack of sufficient reads mapping to the locus limits the power of local multiple sequence alignment and thus it cannot improve the quality of variant calls much. in such circumstances, marking duplication plays a more important role in reducing false positives than realignment and recalibration. moreover, the performances of three popular multi-sample calling tools, samtools, gatk and glfmultiples, also depend on the sequencing depth. they use the same genotype likelihood model, but glfmultiples not only takes into account the maximized likelihood but also an overall prior for each type of polymorphism. for example, they favor sites with transition polymorphisms over those with transversion. thus, incorporating such additional information helps reduce the uncertainty associated with shallow-sequencing data. however, the additional information will disturb the identification of variants when enough evidence is already involved with deep-sequencing data. the steps posterior to variant calling, including linkage-based genotype refining and snp filtering or variant quality score recalibration, also contribute a lot to the accurate snp and genotype calling. the use of ld patterns can substantially improve genotype calling when multiple samples have been sequenced. because not all information regarding errors can be fully incorporated into the statistical framework, the proper snp filtering strategies are recommended to reduce the error rates. besides, the consensus of multiple call sets from different methods provide higher quality than any of individual call sets. even with the best pipelines, however, we are still far from obtaining a complete and accurate picture of snps and genotypes in the human genome. the most challenging task is to distinguish rare variants from sequencing errors. snp and genotype calling for rare variants, which would not be represented in any reference panel, may not improve much by the use of ld information. to identify rare variants, a direct and more powerful approach is to sequence a large number of individuals. in addition to using the proper sequencing strategies, developing more accurate snp detection methods is needed. more research is also needed in other areas, including longer read depths, improved protocols for generate paired ends, advances in sequencing technology with lower base calling error rates, and more powerful alignment methods.  here, we evaluated the effect of a number of computational steps on the accuracy of snp and genotype calling from illumina sequencing data with different coverage. to our knowledge, no other study has made a systematic assessment of whether each step is valuable and how it affects the quality of variant detection. our findings can serves as the general guideline for choosing snp calling strategies. competing interests the authors declare that they have no competing interests. authors' contributions ys led the project and oversaw the analysis. ql and yg designed and performed the research. jl participated in the data analysis. jrl guided the experiments and provided the ngs data. bz guided the analysis and revised the manuscript. ql wrote the manuscript. all authors have read and approved of the final manuscript. supplementary material additional file 1 comparison of effect of different preprocessing steps. a detailed comparison of calling results with different preprocessing steps in terms of dbsnp rate, ti/tv ratio, novel ti/tv ratio and nrd for all regions, inside target regions, outside ≤ bp regions, and outside > bp regions from illumina whole-exome sequencing data. raw, filtery, trim and filtery&trim. click here for file additional file 2 comparison of effect of marking duplication, realignment and recalibration. a detailed comparison of results using different steps, marking duplication, realignment and recalibration, in terms of dbsnp rate, ti/tv ratio, novel ti/tv ratio and nrd for all regions, inside target regions, outside ≤ bp regions, and outside > bp regions from illumina whole-exome sequencing data. initial alignment, marking duplication, realignment, recalibration, marking duplication followed by realignment, marking duplication followed by realignment and recalibration. click here for file additional file 3 comparison of effect of different arrangements of marking duplication, realignment and recalibration. a detailed comparison of results by arranging three steps, marking duplication, realignment and recalibration, in different orders in terms of dbsnp rate, ti/tv ratio, novel ti/tv ratio and nrd for all regions, inside target regions, outside ≤ bp regions, and outside > bp regions from illumina whole-exome sequencing data. marking duplication followed by realignment and recalibration, marking duplication followed by recalibration and realignment, realignment followed by recalibration and marking duplication. click here for file acknowledgements the authors wish to thank peggy schuyler for editorial work on this manuscript and wei zheng for his support. this work was supported by national cancer institute grants u ca p ca p ca p ca and p ca and the national institutes of health grants r01gm. subject recruitment and exome sequencing is supported by ca and ca. ql's work was partially supported by the national natural science foundation of china. this article has been published as part of bmc genomics volume supplement 2012: proceedings of the international conference on intelligent biology and medicine: genomics. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/13/s 
although trimming off low-quality tails helped align more reads, it introduced lots of false positives. our findings clarify the necessity and effectiveness of computational steps for improving the accuracy of snp and genotype calls from illumina sequencing data and can serve as a general guideline for choosing snp calling strategies for data with different coverage. in this study, we made a systematic assessment of the relative contribution of each step to the accuracy of variant calling from illumina dna sequencing data. however, whether each step does contribute to the performance of variant calling and how it affects the accuracy still remain unclear, making it difficult to select and arrange appropriate steps to derive high quality variants from different sequencing data. rearranging these steps did not affect the results. 22- april the international conference on intelligent biology and medicine nashville, tn, usa a number of computational steps are required or recommended when translating the raw sequencing data into the final calls. the relative performance of three popular multi-sample snp callers, samtools, gatk, and glfmultiples, also varied with the sequencing depth. the ability of markup duplication, local realignment and recalibration, to help eliminate false positive variants depended on the sequencing depth.
 accurate calling of snps and genotypes from next-generation sequencing data is an essential prerequisite for most human genetics studies. a number of computational steps are required or recommended when translating the raw sequencing data into the final calls. however, whether each step does contribute to the performance of variant calling and how it affects the accuracy still remain unclear, making it difficult to select and arrange appropriate steps to derive high quality variants from different sequencing data. in this study, we made a systematic assessment of the relative contribution of each step to the accuracy of variant calling from illumina dna sequencing data.  we found that the read preprocessing step did not improve the accuracy of variant calling, contrary to the general expectation. although trimming off low-quality tails helped align more reads, it introduced lots of false positives. the ability of markup duplication, local realignment and recalibration, to help eliminate false positive variants depended on the sequencing depth. rearranging these steps did not affect the results. the relative performance of three popular multi-sample snp callers, samtools, gatk, and glfmultiples, also varied with the sequencing depth.  our findings clarify the necessity and effectiveness of computational steps for improving the accuracy of snp and genotype calls from illumina sequencing data and can serve as a general guideline for choosing snp calling strategies for data with different coverage. 22- april the international conference on intelligent biology and medicine nashville, tn, usa 