BACKGROUND
in the last decade, sequencing technologies, the so-called next generation sequencing , have delivered a step change in the ability to sequence genome, leading to a state of permanent evolution.
ngs platforms allow to detect mutations significantly reducing time and costs  <cit> .
ion torrent personal genome machine  started to be distributed in  <dig>  <cit>  and thus to be used for the identification of genetic variants associated to human diseases  <cit> .
indel detection, in particular, in homopolymer region have a high positive rate  <cit>  which have to be lowered for clinical applications  <cit> .
thus the major challenge in ngs regards the correct manipulation of output data  <cit>  assembling appropriate pipeline, including aligner and variant caller.
diverse algorithms for alignment have been compared in many studies  <cit> .
caboche et al  <cit>  compared different mapping algorithms with ion torrent data, in terms of computational requirement, mapper robustness, ability to map reads in repeated regions and behavior with mutated reference genome.
they were able to optimize a benchmark procedure from whole genome sequencing of small genomes, highlighting the importance to evaluate mappers or to optimize parameters of a chosen mapper for a specific application.
moreover, variant calling pipelines have been compared in relation to different applications and platforms.
yeo et al.
<cit>  optimized an indel detection workflow for brca1/ <dig> pgm sequencing panel.
they compared the proprietary software torrent suite and two open source variant callers, gatk and samtools.
their results showed that snv detection was less problematic than indel identification using torrent suite.
moreover, they demonstrated how the combination of bwa or tmap mappers and samtools is able to improve indel detection.
as demonstrated by these studies, the bioinformatic challenge on ngs data and, in particular, ion torrent data from targeted sequencing requires a lot of efforts in order to correctly identify the best analysis pipeline.
gatk is a well-known toolbox for ngs data analysis.
variant quality score recalibration  step generate an adaptive model based on metrics, such as strand bias, from true variants.
thus it could be possible to calculate if a variant is true or false.
however, this step could be used only for whole genome data or for dataset including more than  <dig> exomes.
for targeted gene panels, gatk’s best practices suggest to set up hard filters specific for the study.
in the present study, we compared variant calling results of gatk pipeline including the use of hard filtering, suggested by gatk’s best practices, and the proprietary torrent suite variant caller regarding a custom panel including  <dig> genes.
then, we focused on two simulated datasets , with high and low coverage, and then we processed the raw variants called by gatk to set parameters of quality in order to increase the number of true variants and decrease the number of false variants.
methods
real dataset
a dataset of  <dig> metastastic melanoma formalin-fixed paraffin embedded  samples was studied.
exonic regions of a panel of  <dig> genes were sequenced with ion torrent pgm.
sequencing data were analyzed using  <dig> standard pipelines including either the torrent suite  with the parameters “somatic” and “high stringency” switched on) or the gatk suite for few variants .
simulated datasets
the above reported panel of  <dig> genes was devised to identify sequence variants in the coding regions.
when the present study was conceived we prepared a catalog of known variants and their allele frequency using the information from the real  <dig> sequenced individuals.
the simulation-based study was then set up using a reference dataset of  <dig> individuals with randomly assigned variants according to the catalog of known variants.
the simulated genotypic profiles of each individual were recorded.
for every sample it was prepared by simulation a file  containing the amplicon sequences that were modified to introduce the assigned variants, repeated  <dig> times  in both the forward and in the reverse sequence.
the dataset of  <dig> fasta files was then processed with the art simulator  <cit>  to generate files  similar to those produced by a sequencer of new generation, mimicking its features and biases.
for every sample the art simulator was launched twice under the hypothesis of a sequencing depth of 20x  and 100x , respectively.
the  <dig> files were then processed using a standard pipeline for variant calling that includes the aligner bwa and the gatk suite of programs .
all the still unfiltered variants were tagged as true  or false  variants.
thus one hundred independent datasets of  <dig> individuals were simulated for a total of  <dig> simulated samples.
by this approach we tried to simulate several times the most likely scenario  resembling to the real one.
gatk pipeline
we followed the toolkit for genome analysis   <cit>  recommendations of dnaseq best practices for calling variants.
hence the following software were used: bwa-mem  for sequence alignment  <cit>  and gatk  <dig>  software for the later steps.
in more detail, sequences underwent the following steps: 1) alignment to the human genome reference version hg <dig>  2) realignment around indels, 3) base recalibration and 4) variant discovery  without been marked for duplicates.
the discovered variants were hard filtered after having selected the rules to setup the filters from the classification trees as described in the “filters for the variant calling” and “classification trees” sections.
filters for the variant calling
hard filtering evaluated  <dig> standard gatk filters  and  <dig> filters  that were not present in the standard gatk vcf output files.
fs is the p value from the contingency table of the number of reads calling the alleles at the variant site on either the dna strands .
adt and adtl evaluate imbalances in calling the reference and the alternative allele , also depending on the amount of reads that map on the variant locus .
adt and adtl are defined as follows: adt=ad1−ad2ad1+ad2adtl=log10ad1+ad2*adt where ad <dig> and ad <dig> are the number of unfiltered reads calling the reference and alternate allele, respectively.
descriptive statistics were performed using r  <dig> .
<dig> and the kruskal.test function, and rocr library version  <dig> - <dig>
classification trees
every filter was included into classification trees to target the filter rules that better discriminate the true variants  from false variants .
these rules have been then used for the hard filtering.
eight classification trees were generated to investigate separately snv and indel, homozygotes and heterozygotes  under a simulated coverage of 20x or 100x.
for every tree we extracted the selected filters and their threshold values in order to use them in the phase of hard-filtering.
we extracted the filters that were listed starting from the root node of the classification tree to the next  <dig> consecutive daughter nodes .
the filters selection for each classification tree was made as follows.
comparing every daughter node with the root node we targeted the nodes that 1) contained at least 10% of the true calls and 2) where the ratio of the number of true calls and false calls was greater than  <dig>  to minimize the number and complexity of the filter rules we then considered for each targeted node l  the upper connected node u  closer to the root node.
we then selected the node u instead of the node l when 1) the node l contained less than 80% of the true calls or 2) when containing more than 70% of false calls of node u , respectively.
once the relevant nodes were selected we extracted the filter rules starting from the root node toward each of the finally selected daughter nodes.
classification trees were produced using the library rpart  of the r package
RESULTS
results in the real dataset: gatk vs tvc
in the first part of this study we looked at the results obtained using  <dig> different standard pipelines on the same set of sequences.
tvc is the software use for targeted sequencing in bundle with the ion torrent sequencer.
gatk is considered the “gold standard” in managing large ngs data  and can be used for targeted sequencing.
for this reason, we compared tvc calls with those produced by gatk  <dig> .
tvc called  <dig> variants in the entire dataset,  <dig> of which were shared with gatk that detected  <dig> snvs.
then we performed the vcf files, which are the output files of both tvc and gatk, focusing on some parameters of sequencing quality.
in particular, dp  and af  tags were shared by vcf outputs.
tvc calls showed mean, dp and af values of  <dig>  and  <dig> , respectively.
the  <dig> snvs called by both tvc and gatk showed higher mean.
this observation could suggest that shared snvs are true positive calls, having higher coverage and quality by depth values.
it is noteworthy that mean poq values of the variants called only by gatk were similar to tvc calls .
snvs identified only by tvc showed lower mean dp and af values than shared variations .
such results suggest that snvs called by tvc may be enriched of many false positives.
indel calling is a highly debated problem when we refer to ion torrent data.
in a very preliminary way, we evaluated mean dp and af values of  <dig> indels identified by gatk, which displayed similar values to those of shared snvs .
no indel was detected by tvc.
high-coverage and low-coverage simulated datasets: descriptive statistics
to better explore gatk variant calling and to try to tune the hard filtering parameters , we performed a simulation-based study, as described in the “methods” section.
in the high-coverage  dataset, gatk identified  <dig> unfiltered snvs and  <dig> unfiltered indels.
it is noteworthy that  <dig> % and only  <dig> % of snvs and indels were true variants, respectively.
in the low-coverage  dataset,  <dig> and  <dig> unfiltered snvs and unfiltered indels were respectively called.
as expected, we found that the percentage of true variants was lower in the lc dataset  .table  <dig> overall gatk unfiltered alterations identified in hc and lc dataset
indels
snvs
indels
snvs
tv true variants, fv false variants
we observed that  <dig> % of snvs and  <dig> % of indels were homozygous in the hc dataset whereas  <dig> % of snvs and  <dig> % of indels were homozygous in the lc dataset.
we investigated the distribution of the values of the individual gatk filters, namely baseqranksum, readposranksum, clippingranksum, dp, mq, mqranksum, and gq  between true and false variants.
table  <dig> reports the descriptive statistics for the hc dataset and shows that baseqranksum, readposranksum and dp display a statistically significant difference both in snvs and indels.
mqrs showed a statistically significant difference in indels only.table  <dig> descriptive statistics of gatk filters in the hc dataset, stratifying calls by type .
data are displayed as mean ± sd
tv
mean ± sd
fv
mean ± sd
p-value
tv
mean ± sd
fv
mean ± sd
p-value
the mean value is the mean value of the median value from each of the  <dig> replicates
bqrs baseqranksum, rprs readposranksum, crs clippingranksum, dp depth of coverage, mq mappingquality, mqrs mappingqualityranksum, gq genotype quality, tv true variants, fv false variants
table  <dig> reports the descriptive statistics for the lc dataset and shows a statistically significant difference for all the gatk filters with the exception of mq both in snvs and indels subsets.
of note that gq always  presented a value equal to  <dig> for the true variants.
it is intriguing that in the case of snv subset, dp presents the lowest values in the tvs compared to fvs even the difference is relatively small.table  <dig> descriptive statistics of gatk filters in the lc dataset, stratifying calls by type .
data are displayed as mean ± sd
tv
mean ± sd
fv
mean ± sd
p-value
tv
mean ± sd
fv
mean ± sd
p-value
-
-
-
-
the mean value is the mean value of the median value from each of the  <dig> replicates
bqrs baseqranksum, rprs readposranksum, crs clippingranksum, dp depth of coverage, mq mappingquality, mqrs mappingqualityranksum, gq genotype quality, tv true
the performance of all individual filters to discriminate between true and false variants was summarized by estimating the area under the roc curve .
table  <dig> reports that adt and gq showed the best performance.
additionally every filter showed a better performance for snvs rather than indels .table  <dig> perfomance of the individual filters evaluated to discriminate between true and false variants by the auc values from roc curve, grouped by type of variants  and status of the genotype call  according to the depth of sequencing
homo
het
homo
het
it could be noticed that correctly called alterations showed a higher coverage than false variants, highlighting the importance of this parameter.
classification trees
we performed the analyses of either snvs or indels subsets stratified by genotype, in hc and lc datasets.
classification trees  allowed to set a series of filter rules for each of the  <dig> type of sequence alteration.
table  <dig> shows the parameters and threshold values to be used for hard-filtering been extracted from the classification trees.table  <dig> parameters and their thresholds selected by regression trees
: no reliable filters were selected by classification trees
notably, the classification tree did not select any reliable filter for homozygous indels in the lc dataset .table  <dig> results by the application of selection parameters and their thresholds on simulated datasets
overall
selected
overall
selected
overall
selected
overall
selected
overall
selected
overall
selected
overall
selected
overall
selected
% have to be intended as the percentage of unfiltered variants for “overall”calls and as the percentage of alterations which were not filtered out in the hard filtering process for “selected”calls; % of selection indicates the amount of variants selected from the total callset.
tv true variants, fv false variants
we then explored the sequence of the flanking regions of each type of alterations, in particular we observed that short homopolymeric strings are recurrent and therefore partly responsible for false positive calls .table  <dig> homopolymeric sequences flanking false positive variants
variant locus is on the 6th nucletide  of the  <dig> nucleotide string
discussion
it is well known that there are many technical challenges involved in getting an accurate variant calling procedure of ngs data including the bioinformatic analysis.
a number of tools based on complex statistical models has been developed but many concerns related to their performance remain still open.
since the number of the called variants varies from software to software, typically more than one computer program is then used.
if the variant is actually called by all the programs then its support increases.
however, the problem occurs when the variant is called only by some programs, raising the suspicion that it is not true.
ngs is now applied in many fields.
we were interested in studying the case of targeted sequencing of small set of genes when using a common ngs platform such as ion torrent.
when analyzing a few variants  the gatk guides suggest to use filters that must be set by the user  rather than the adaptive filtering that needs a high number of variants to work properly.
under these conditions  it becomes important to tailor the filters to call the true variants on the specific design.
it is likely that different panels of genes and even different designs for the same panel of genes require a different setup of filters.
we therefore tried to explore through a simulation-based study the outcomes that the pipeline for the variant calling may encounter.
we were interested in the study of a specific scenario made up a group of individuals with a specific sequencing design.
we measured the performance inthe calling true variants for each of the filters that can be set when working with hard-filtering.
hence we used classification trees on a large data simulated dataset of true and false variants.
in the present paper, we studied several standard and non standard gatk filters to be used for hard-filtering in the context of a targeted gene panel sequencing.
firstly, we analyzed a real dataset coming from the sequencing of an ion torrent targeted gene panel observing a high discrepancy between tvc and gatk, particularly for indels, suggesting that such type variants are even difficult to be detected by the present bionformatic tools.
in fact the importance to define a “gold standard” dataset to test variant calling methods is a very hot topic.
recently, “synthetic” matched tumor/normal samples was created for comparing performances of popular variant callers in detection of “somatic” snvs.
however, even if they had the advantage to refer to nist-giab  <cit>  as gold standard, authors could not discriminate “somatic” snvs from germline background, an important issue when studying tumors, and moreover the batch that they purchased was not the same used for nist-giab  <cit> .
in the present study, we decided to simulate two datasets, each with a different coverage and carrying alterations found in real data.
notwithstanding this investigation did not simulate tumor heterogeneity, gatk variant calling was tested both in a relative high coverage and low coverage conditions.
recently, vanni et al  <cit>  highlighted the discrepancy between tvc and gatk, which was also observed in our study, excluding indel calls from comparison.
they considered phred score ranging 5– <dig> to mark low-quality variants.
our results show that such approach could not be enough to have a high quality gatk call set.
in detail, we evidenced that different parameters could be tuned depending on type of mutations and genotypes suggested.
in a previous study, authors focused on the detection of parameters that could allow to improve indel detection  <cit> .
they focused on two parameters regarding the frequency of reference and alternate alleles and the variance of the width of inserted/deleted sequences.
the first parameters is similar to our adt and adtl filters, which were involved in the step of selection for the reduction of false positives.
they observed that the numbers of false positive regarded in particular indel in homopolymeric regions.
in a similar way we observed that flanking regions, were homopolymeric in a high number of false positive calls.
it is important to reduce errors in these regions because they occur in genomic regions where the occurrence of true alterations is also higher  <cit> .
variant calling of tvc is improving but indels are still a problem and thus parallel pipeline with opportune set up and filters could be helpful to solve this question, with particular attention on type of platform used for sequencing and on type of design .
carson et al  <cit>  recently demonstrated how dp and gq filters could be able to enhance sensitivity and specificity in whole exome sequencing data.
they tested different thresholds and showed that over a certain threshold accuracy reached a plateau and notably they demonstrated that vqsr is not enough to improve variant calling.
indeed, they concluded that, also when vqsr could be applied, opportune hard filtering strategies need to be set up.
our intent was not to target the precise hard filter parameter values and our results have to be intended as suggestion in handling data coming from targeted gene panel sequencing in order to optimize gatk variant calling outputs.
we observed that hard filtering was able to reduce the number of unfiltered false positives, with a different efficiency between snvs  and indels .
true indels were hard to be filtered and the performance of filtering was generally lower than in the case of the snvs .
we also observed that the very majority to the unfiltered indels called at the homozygous status is represented by false variants.
hence there are some regions of the dna reference sequence that are prone to be recognized as carrying indels by the bioinformatic pipeline.
our preliminary investigations show that these regions often contain low complexity sequences .
a good strategy would be to train in advance the program that operates in the indel parameter recalibration phase to recognize these regions but this hypothesis needs to be investigated in more detail.
the reader should note that even if some ngs technologies are known to read with difficulty the dna regions having low complexity, we are here asserting that the call errors of the variants in the low complexity regions are due to the bioinformatic analysis and not to sequencing errors as we worked with data produced by the simulator and not by a ngs sequencer.
so, in the real world, regions with low complexity sequences are doubly condemned to possible abundant errors due to both sequencing and the following bioinformatics analysis.
some results were quite different from the expected.
in particular, we observed that they were detected  more snvs in the case of low coverage  than in high coverage.
of note that in such cases the rate of false variant is very small even for the unfiltered variants.
however, the rate of false variant is about  <dig> folds greater in the lc dataset.
therefore we hypothesize that variant calling process is more sensitive because of less specific when analyzing dataset with a low coverage.
in general terms, as already known, the indels are more difficult to be analyzed than the snvs and a deeper sequencing helps to improve the performance of filtering the true variants.
however, such a performance could vary since complexity of the sequence changes along the sequence itself.
in the case of targeted sequencing, our suggestion is to study in advance the region that we will be sequenced in order to evaluate the performance of the variant calling procedure over such regions in order to figure out the most problematic areas to be treated with caution.
we also suggest the use of simulations based on the specific target region which can help to calibrate the filters for the specific problem.
we argue that it could be useful to set specific filters for different regions and for different known variants.
CONCLUSIONS
the results of our study showed that filters could be correctly tuned according to coverage and type of alterations.
moreover, it could be useful to test by appropriate simulations the design of amplicon gene panels to gain a priori knowledge of the possible issues in variant calling by gatk.
additional files
additional file 1: technical definitions of parameters of sequencing quality.
additional file 2: performance of all individual filters to discriminate between true and false variants estimated by the area under the roc curve .
additional file 3: classification trees to set filter rules.
abbreviations
snvsingle nucleotide variants
acknowledgment
fondazione cariverona
funding
publication of this article was funded by italian ministry of health .
availability of data and materials
the real datasets are available upon request from the author.
the simulated datasets can be accessed at http://medgen.univr.it/downloads/bmcbioinf_desumma/index.html
authors’ contributions
sds, gm: conceived and designed the study; rp: experiments for real dataset data; am, vm: bioinformatic analyses; st: revision of the manuscript.
all authors have read and approved the final manuscript.
competing interests
the authors declare that they have no competing interests.
ethics approval and consent to participate
the study was approved by the local ethics committee of the istituto tumori “giovanni paolo ii” of bari  and was conducted in accordance with the international standards of good clinical practice.
all patients signed an informed consent.
consent for publication
not applicable.
about this supplement
this article has been published as part of bmc bioinformatics volume  <dig> supplement  <dig>  2017: selected works from the joint 15th network tools and applications in biology international workshop and 11th integrative bioinformatics international symposium .
the full contents of the supplement are available online at https://bmcbioinformatics.biomedcentral.com/articles/supplements/volume-18-supplement- <dig>

 in the last decade, sequencing technologies, the so-called next generation sequencing, have delivered a step change in the ability to sequence genome, leading to a state of permanent evolution. ngs platforms allow to detect mutations significantly reducing time and costs. ion torrent personal genome machine started to be distributed in and thus to be used for the identification of genetic variants associated to human diseases. indel detection, in particular, in homopolymer region have a high positive rate which have to be lowered for clinical applications. thus the major challenge in ngs regards the correct manipulation of output data assembling appropriate pipeline, including aligner and variant caller. diverse algorithms for alignment have been compared in many studies. caboche et al compared different mapping algorithms with ion torrent data, in terms of computational requirement, mapper robustness, ability to map reads in repeated regions and behavior with mutated reference genome. they were able to optimize a benchmark procedure from whole genome sequencing of small genomes, highlighting the importance to evaluate mappers or to optimize parameters of a chosen mapper for a specific application. moreover, variant calling pipelines have been compared in relation to different applications and platforms. yeo et al. optimized an indel detection workflow for brca1/ pgm sequencing panel. they compared the proprietary software torrent suite and two open source variant callers, gatk and samtools. their results showed that snv detection was less problematic than indel identification using torrent suite. moreover, they demonstrated how the combination of bwa or tmap mappers and samtools is able to improve indel detection. as demonstrated by these studies, the bioinformatic challenge on ngs data and, in particular, ion torrent data from targeted sequencing requires a lot of efforts in order to correctly identify the best analysis pipeline. gatk is a well-known toolbox for ngs data analysis. variant quality score recalibration step generate an adaptive model based on metrics, such as strand bias, from true variants. thus it could be possible to calculate if a variant is true or false. however, this step could be used only for whole genome data or for dataset including more than exomes. for targeted gene panels, gatk’s best practices suggest to set up hard filters specific for the study. in the present study, we compared variant calling results of gatk pipeline including the use of hard filtering, suggested by gatk’s best practices, and the proprietary torrent suite variant caller regarding a custom panel including genes. then, we focused on two simulated datasets, with high and low coverage, and then we processed the raw variants called by gatk to set parameters of quality in order to increase the number of true variants and decrease the number of false variants. methods real dataset a dataset of metastastic melanoma formalin-fixed paraffin embedded samples was studied. exonic regions of a panel of genes were sequenced with ion torrent pgm. sequencing data were analyzed using standard pipelines including either the torrent suite with the parameters “somatic” and “high stringency” switched on) or the gatk suite for few variants. simulated datasets the above reported panel of genes was devised to identify sequence variants in the coding regions. when the present study was conceived we prepared a catalog of known variants and their allele frequency using the information from the real sequenced individuals. the simulation-based study was then set up using a reference dataset of individuals with randomly assigned variants according to the catalog of known variants. the simulated genotypic profiles of each individual were recorded. for every sample it was prepared by simulation a file containing the amplicon sequences that were modified to introduce the assigned variants, repeated times in both the forward and in the reverse sequence. the dataset of fasta files was then processed with the art simulator to generate files similar to those produced by a sequencer of new generation, mimicking its features and biases. for every sample the art simulator was launched twice under the hypothesis of a sequencing depth of 20x and 100x, respectively. the files were then processed using a standard pipeline for variant calling that includes the aligner bwa and the gatk suite of programs. all the still unfiltered variants were tagged as true or false variants. thus one hundred independent datasets of individuals were simulated for a total of simulated samples. by this approach we tried to simulate several times the most likely scenario resembling to the real one. gatk pipeline we followed the toolkit for genome analysis recommendations of dnaseq best practices for calling variants. hence the following software were used: bwa-mem for sequence alignment and gatk software for the later steps. in more detail, sequences underwent the following steps: 1) alignment to the human genome reference version hg 2) realignment around indels, 3) base recalibration and 4) variant discovery without been marked for duplicates. the discovered variants were hard filtered after having selected the rules to setup the filters from the classification trees as described in the “filters for the variant calling” and “classification trees” sections. filters for the variant calling hard filtering evaluated standard gatk filters and filters that were not present in the standard gatk vcf output files. fs is the p value from the contingency table of the number of reads calling the alleles at the variant site on either the dna strands. adt and adtl evaluate imbalances in calling the reference and the alternative allele, also depending on the amount of reads that map on the variant locus. adt and adtl are defined as follows: adt=ad1−ad2ad1+ad2adtl=log10ad1+ad2*adt where ad and ad are the number of unfiltered reads calling the reference and alternate allele, respectively. descriptive statistics were performed using r. and the kruskal.test function, and rocr library version - classification trees every filter was included into classification trees to target the filter rules that better discriminate the true variants from false variants. these rules have been then used for the hard filtering. eight classification trees were generated to investigate separately snv and indel, homozygotes and heterozygotes under a simulated coverage of 20x or 100x. for every tree we extracted the selected filters and their threshold values in order to use them in the phase of hard-filtering. we extracted the filters that were listed starting from the root node of the classification tree to the next consecutive daughter nodes. the filters selection for each classification tree was made as follows. comparing every daughter node with the root node we targeted the nodes that 1) contained at least 10% of the true calls and 2) where the ratio of the number of true calls and false calls was greater than to minimize the number and complexity of the filter rules we then considered for each targeted node l the upper connected node u closer to the root node. we then selected the node u instead of the node l when 1) the node l contained less than 80% of the true calls or 2) when containing more than 70% of false calls of node u, respectively. once the relevant nodes were selected we extracted the filter rules starting from the root node toward each of the finally selected daughter nodes. classification trees were produced using the library rpart of the r package  results in the real dataset: gatk vs tvc in the first part of this study we looked at the results obtained using different standard pipelines on the same set of sequences. tvc is the software use for targeted sequencing in bundle with the ion torrent sequencer. gatk is considered the “gold standard” in managing large ngs data and can be used for targeted sequencing. for this reason, we compared tvc calls with those produced by gatk. tvc called variants in the entire dataset, of which were shared with gatk that detected snvs. then we performed the vcf files, which are the output files of both tvc and gatk, focusing on some parameters of sequencing quality. in particular, dp and af tags were shared by vcf outputs. tvc calls showed mean, dp and af values of and, respectively. the snvs called by both tvc and gatk showed higher mean. this observation could suggest that shared snvs are true positive calls, having higher coverage and quality by depth values. it is noteworthy that mean poq values of the variants called only by gatk were similar to tvc calls. snvs identified only by tvc showed lower mean dp and af values than shared variations. such results suggest that snvs called by tvc may be enriched of many false positives. indel calling is a highly debated problem when we refer to ion torrent data. in a very preliminary way, we evaluated mean dp and af values of indels identified by gatk, which displayed similar values to those of shared snvs. no indel was detected by tvc. high-coverage and low-coverage simulated datasets: descriptive statistics to better explore gatk variant calling and to try to tune the hard filtering parameters, we performed a simulation-based study, as described in the “methods” section. in the high-coverage dataset, gatk identified unfiltered snvs and unfiltered indels. it is noteworthy that % and only % of snvs and indels were true variants, respectively. in the low-coverage dataset, and unfiltered snvs and unfiltered indels were respectively called. as expected, we found that the percentage of true variants was lower in the lc dataset.table overall gatk unfiltered alterations identified in hc and lc dataset   indels   snvs   indels   snvs tv true variants, fv false variants we observed that % of snvs and % of indels were homozygous in the hc dataset whereas % of snvs and % of indels were homozygous in the lc dataset. we investigated the distribution of the values of the individual gatk filters, namely baseqranksum, readposranksum, clippingranksum, dp, mq, mqranksum, and gq between true and false variants. table  reports the descriptive statistics for the hc dataset and shows that baseqranksum, readposranksum and dp display a statistically significant difference both in snvs and indels. mqrs showed a statistically significant difference in indels only.table descriptive statistics of gatk filters in the hc dataset, stratifying calls by type. data are displayed as mean ± sd tv mean ± sd fv mean ± sd p-value tv mean ± sd fv mean ± sd p-value the mean value is the mean value of the median value from each of the replicates bqrs baseqranksum, rprs readposranksum, crs clippingranksum, dp depth of coverage, mq mappingquality, mqrs mappingqualityranksum, gq genotype quality, tv true variants, fv false variants table  reports the descriptive statistics for the lc dataset and shows a statistically significant difference for all the gatk filters with the exception of mq both in snvs and indels subsets. of note that gq always presented a value equal to for the true variants. it is intriguing that in the case of snv subset, dp presents the lowest values in the tvs compared to fvs even the difference is relatively small.table descriptive statistics of gatk filters in the lc dataset, stratifying calls by type. data are displayed as mean ± sd tv mean ± sd fv mean ± sd p-value tv mean ± sd fv mean ± sd p-value - - - - the mean value is the mean value of the median value from each of the replicates bqrs baseqranksum, rprs readposranksum, crs clippingranksum, dp depth of coverage, mq mappingquality, mqrs mappingqualityranksum, gq genotype quality, tv true the performance of all individual filters to discriminate between true and false variants was summarized by estimating the area under the roc curve. table  reports that adt and gq showed the best performance. additionally every filter showed a better performance for snvs rather than indels.table perfomance of the individual filters evaluated to discriminate between true and false variants by the auc values from roc curve, grouped by type of variants and status of the genotype call according to the depth of sequencing homo het homo het it could be noticed that correctly called alterations showed a higher coverage than false variants, highlighting the importance of this parameter. classification trees we performed the analyses of either snvs or indels subsets stratified by genotype, in hc and lc datasets. classification trees allowed to set a series of filter rules for each of the type of sequence alteration. table  shows the parameters and threshold values to be used for hard-filtering been extracted from the classification trees.table parameters and their thresholds selected by regression trees: no reliable filters were selected by classification trees notably, the classification tree did not select any reliable filter for homozygous indels in the lc dataset.table results by the application of selection parameters and their thresholds on simulated datasets    overall    selected    overall    selected    overall    selected    overall    selected    overall    selected    overall    selected    overall    selected    overall    selected % have to be intended as the percentage of unfiltered variants for “overall”calls and as the percentage of alterations which were not filtered out in the hard filtering process for “selected”calls; % of selection indicates the amount of variants selected from the total callset. tv true variants, fv false variants we then explored the sequence of the flanking regions of each type of alterations, in particular we observed that short homopolymeric strings are recurrent and therefore partly responsible for false positive calls.table homopolymeric sequences flanking false positive variants variant locus is on the 6th nucletide of the nucleotide string discussion it is well known that there are many technical challenges involved in getting an accurate variant calling procedure of ngs data including the bioinformatic analysis. a number of tools based on complex statistical models has been developed but many concerns related to their performance remain still open. since the number of the called variants varies from software to software, typically more than one computer program is then used. if the variant is actually called by all the programs then its support increases. however, the problem occurs when the variant is called only by some programs, raising the suspicion that it is not true. ngs is now applied in many fields. we were interested in studying the case of targeted sequencing of small set of genes when using a common ngs platform such as ion torrent. when analyzing a few variants the gatk guides suggest to use filters that must be set by the user rather than the adaptive filtering that needs a high number of variants to work properly. under these conditions it becomes important to tailor the filters to call the true variants on the specific design. it is likely that different panels of genes and even different designs for the same panel of genes require a different setup of filters. we therefore tried to explore through a simulation-based study the outcomes that the pipeline for the variant calling may encounter. we were interested in the study of a specific scenario made up a group of individuals with a specific sequencing design. we measured the performance inthe calling true variants for each of the filters that can be set when working with hard-filtering. hence we used classification trees on a large data simulated dataset of true and false variants. in the present paper, we studied several standard and non standard gatk filters to be used for hard-filtering in the context of a targeted gene panel sequencing. firstly, we analyzed a real dataset coming from the sequencing of an ion torrent targeted gene panel observing a high discrepancy between tvc and gatk, particularly for indels, suggesting that such type variants are even difficult to be detected by the present bionformatic tools. in fact the importance to define a “gold standard” dataset to test variant calling methods is a very hot topic. recently, “synthetic” matched tumor/normal samples was created for comparing performances of popular variant callers in detection of “somatic” snvs. however, even if they had the advantage to refer to nist-giab as gold standard, authors could not discriminate “somatic” snvs from germline background, an important issue when studying tumors, and moreover the batch that they purchased was not the same used for nist-giab. in the present study, we decided to simulate two datasets, each with a different coverage and carrying alterations found in real data. notwithstanding this investigation did not simulate tumor heterogeneity, gatk variant calling was tested both in a relative high coverage and low coverage conditions. recently, vanni et al highlighted the discrepancy between tvc and gatk, which was also observed in our study, excluding indel calls from comparison. they considered phred score ranging 5– to mark low-quality variants. our results show that such approach could not be enough to have a high quality gatk call set. in detail, we evidenced that different parameters could be tuned depending on type of mutations and genotypes suggested. in a previous study, authors focused on the detection of parameters that could allow to improve indel detection. they focused on two parameters regarding the frequency of reference and alternate alleles and the variance of the width of inserted/deleted sequences. the first parameters is similar to our adt and adtl filters, which were involved in the step of selection for the reduction of false positives. they observed that the numbers of false positive regarded in particular indel in homopolymeric regions. in a similar way we observed that flanking regions, were homopolymeric in a high number of false positive calls. it is important to reduce errors in these regions because they occur in genomic regions where the occurrence of true alterations is also higher. variant calling of tvc is improving but indels are still a problem and thus parallel pipeline with opportune set up and filters could be helpful to solve this question, with particular attention on type of platform used for sequencing and on type of design. carson et al recently demonstrated how dp and gq filters could be able to enhance sensitivity and specificity in whole exome sequencing data. they tested different thresholds and showed that over a certain threshold accuracy reached a plateau and notably they demonstrated that vqsr is not enough to improve variant calling. indeed, they concluded that, also when vqsr could be applied, opportune hard filtering strategies need to be set up. our intent was not to target the precise hard filter parameter values and our results have to be intended as suggestion in handling data coming from targeted gene panel sequencing in order to optimize gatk variant calling outputs. we observed that hard filtering was able to reduce the number of unfiltered false positives, with a different efficiency between snvs and indels. true indels were hard to be filtered and the performance of filtering was generally lower than in the case of the snvs. we also observed that the very majority to the unfiltered indels called at the homozygous status is represented by false variants. hence there are some regions of the dna reference sequence that are prone to be recognized as carrying indels by the bioinformatic pipeline. our preliminary investigations show that these regions often contain low complexity sequences. a good strategy would be to train in advance the program that operates in the indel parameter recalibration phase to recognize these regions but this hypothesis needs to be investigated in more detail. the reader should note that even if some ngs technologies are known to read with difficulty the dna regions having low complexity, we are here asserting that the call errors of the variants in the low complexity regions are due to the bioinformatic analysis and not to sequencing errors as we worked with data produced by the simulator and not by a ngs sequencer. so, in the real world, regions with low complexity sequences are doubly condemned to possible abundant errors due to both sequencing and the following bioinformatics analysis. some results were quite different from the expected. in particular, we observed that they were detected more snvs in the case of low coverage than in high coverage. of note that in such cases the rate of false variant is very small even for the unfiltered variants. however, the rate of false variant is about folds greater in the lc dataset. therefore we hypothesize that variant calling process is more sensitive because of less specific when analyzing dataset with a low coverage. in general terms, as already known, the indels are more difficult to be analyzed than the snvs and a deeper sequencing helps to improve the performance of filtering the true variants. however, such a performance could vary since complexity of the sequence changes along the sequence itself. in the case of targeted sequencing, our suggestion is to study in advance the region that we will be sequenced in order to evaluate the performance of the variant calling procedure over such regions in order to figure out the most problematic areas to be treated with caution. we also suggest the use of simulations based on the specific target region which can help to calibrate the filters for the specific problem. we argue that it could be useful to set specific filters for different regions and for different known variants.  the results of our study showed that filters could be correctly tuned according to coverage and type of alterations. moreover, it could be useful to test by appropriate simulations the design of amplicon gene panels to gain a priori knowledge of the possible issues in variant calling by gatk. additional files additional file 1: technical definitions of parameters of sequencing quality. additional file 2: performance of all individual filters to discriminate between true and false variants estimated by the area under the roc curve. additional file 3: classification trees to set filter rules. abbreviations snvsingle nucleotide variants acknowledgment fondazione cariverona funding publication of this article was funded by italian ministry of health. availability of data and materials the real datasets are available upon request from the author. the simulated datasets can be accessed at http://medgen.univr.it/downloads/bmcbioinf_desumma/index.html authors’ contributions sds, gm: conceived and designed the study; rp: experiments for real dataset data; am, vm: bioinformatic analyses; st: revision of the manuscript. all authors have read and approved the final manuscript. competing interests the authors declare that they have no competing interests. ethics approval and consent to participate the study was approved by the local ethics committee of the istituto tumori “giovanni paolo ii” of bari and was conducted in accordance with the international standards of good clinical practice. all patients signed an informed consent. consent for publication not applicable. about this supplement this article has been published as part of bmc bioinformatics volume supplement 2017: selected works from the joint 15th network tools and applications in biology international workshop and 11th integrative bioinformatics international symposium. the full contents of the supplement are available online at https://bmcbioinformatics.biomedcentral.com/articles/supplements/volume-18-supplement- 
simulated sequences were aligned with standard protocols and then regression trees were built up to identify the most reliable parameters and cutoff values to discriminate true and false variant calls. moreover, we analyzed flanking sequences of region presenting a high rate of false positive calls observing that such sequences present a low complexity make up. this was pursued by using classification trees on true and false variants from simulated sequences of a real dataset data. our results showed that gatk hard filtering parameter values can be tailored through a simulation study based-on the dna region of interest to ameliorate the accuracy of the variant calling. we simulated two datasets, with different coverages, including all the sequence alterations identified in a real dataset according to their observed frequencies. ngsvariant callingvariant filteringtargeted gene panelsnvindeljoint 14th network tools and applications in biology workshop and 11th integrative bioinformatics symposium nettab/ib 2015bari, italy 14- october issue-copyright-statement© the author 2017 the aim of the present paper was to set up a procedure to assess the best parameters to be used in the hard filtering of gatk.
 ngs technology represents a powerful alternative to the standard sanger sequencing in the context of clinical setting. the proprietary software that are generally used for variant calling often depend on preset parameters that may not fit in a satisfactory manner for different genes. gatk, which is widely used in the academic world, is rich in parameters for variant calling. however the self-adjusting parameter calibration of gatk requires data from a large number of exomes. when these are not available, which is the standard condition of a diagnostic laboratory, the parameters must be set by the operator. the aim of the present paper was to set up a procedure to assess the best parameters to be used in the hard filtering of gatk. this was pursued by using classification trees on true and false variants from simulated sequences of a real dataset data.  we simulated two datasets, with different coverages, including all the sequence alterations identified in a real dataset according to their observed frequencies. simulated sequences were aligned with standard protocols and then regression trees were built up to identify the most reliable parameters and cutoff values to discriminate true and false variant calls. moreover, we analyzed flanking sequences of region presenting a high rate of false positive calls observing that such sequences present a low complexity make up.  our results showed that gatk hard filtering parameter values can be tailored through a simulation study based-on the dna region of interest to ameliorate the accuracy of the variant calling. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords ngsvariant callingvariant filteringtargeted gene panelsnvindeljoint 14th network tools and applications in biology workshop and 11th integrative bioinformatics symposium nettab/ib 2015bari, italy 14- october issue-copyright-statement© the author 2017 