Comparative Identification and Assessment of Single Nucleotide Variants Through Shifts in Base Call Frequencies

Starting with a set of aligned reads we observe the sequencing depth \(N^{S}_{i}\) and the number of non-consensus base calls \(K^{S}_{i}\) for a sample \(S\) at position \(i\). Non-consensus is here defined as base calls that differ from the consensus sequence, which can be either the reference sequence or, in a comparative setting, the most abundant base call in the control sample. The non-consensus rate \(p^{S}_{i}\), which we assume to be binomially distributed, can then be estimated as

$$\hat{p}^{S}_{i} = \frac{K^{S}_{i}}{N^{S}_{i}}$$

for sites with \(N^{S}_{i} > 0\). Since \(K^{S}_{i} \leq N^{S}_{i}\), the rates are bound to the range \([0,1]\).

The true non-consensus rate

$$p^{S}_{i} = v^{S}_{i} + e^{S}_{i}$$

comprises the presence of a putative variant with a frequency \(v^{S}_{i}\) and a technical error rate \(e^{S}_{i}\). In order to detect and describe the change in the variant frequency, we focus on the shift \(d_{i}\) in non-consensus rates as the difference of the rates between the test and control samples, which we estimate as

$$\hat{d}_{i} = \hat{p}^{T}_{i} - \hat{p}^{C}_{i}.$$

If we assume that the true site-specific technical error rates are identical between the two matched samples [10], the difference of the rates yields an unbiased estimate for the change in the variant frequency. Thus, positions not haboring biological alterations will result in \(\hat{d}_{i} \approx 0\).

Distinguishing biological variants from noise requires knowledge about the variance of the point estimate \(\hat{d_{i}}\). By constructing a confidence interval (CI) for \(d_{i}\) with confidence level \(\beta\) [14], we assess the certainty of the estimated shift in non-consensus frequencies. The probability of the true value being outside the confidence interval is less than \(\alpha = 1 - \beta\). This is in concordance with the type I or \(\alpha\) error definition in statistical testing.

Under the assumption that the non-consensus counts \(K^{S}_{i}\) in our samples follow binomial distributions with parameters \(p^{S}_{i}\) and \(N^{S}_{i}\), several methods have been established for estimating confidence intervals for the difference of two rate parameters [12] [7]. The performance of an approach is generally described in terms of its coverage probability indicating the probability of a confidence interval to cover the true value (see 4.2). Coverage probabilities greater and less than the confidence level \(\beta\) describe conservative and liberal behaviors, respectively. Due to the conservative coverage probabilities and high computational effort of exact confidence interval estimates, approximate methods are generally preferred [1] [7].

The Agresti-Caffo (AC) confidence interval [2]

$$\tilde{p}^{T} - \tilde{p}^{C} \pm z \sqrt{ \frac{\tilde{p}^{T} (1 - \tilde{p}^{T})} {\tilde{N}^{T}} + \frac{\tilde{p}^{C}(1 - \tilde{p}^{C})} {\tilde{N}^{C}} }$$

with

$$\tilde{p}^{X} = \frac{K^{X}+\zeta}{N^{X}+2\zeta},$$

$$\tilde{N}^{X} = N^{X} + 2\zeta,$$

$$\zeta = \frac{1}{4} z^2,$$

and \(z = z_{(1-\beta)/2}\) as the upper \((1-\beta)/2\) percentile of the standard normal distribution), is an approximation of the score test-based confidence interval. Several publications emphasize the usefulness and advantages of the AC method over related approaches [7] [3] [5].

Focusing on the comparative shift of non-consensus frequencies allows us to detect and distinguish different types of events. Since Rariant does not make explicit assumptions about the abundance of a potential variant in the control sample, we are further able to find clonal shifts, for example between different tumor samples, or losses of heterozygocity. Generally, gains and losses of variant alleles are characterized by positive and negative values of \(d\), respectively. For a differentiated interpretation of the results, we classify a variant into one of four classes:

somatic

A somatic variant that does not occur in the control sample

hetero/LOH

A shift away from heterozygous SNP in the control sample

undecided

Both of the 'somatic' or 'hetero' are possible

powerless

A distinction between the two classes cannot be made due to a lack of power

The classification is based on two binomial tests for each position:

1. Somatic variants where the variant allele is not present in the control sample, rejecting a binomial test with the alternative hypothesis \(H_{1}: p^{C} > 0\).
2. Sites with a loss of heterozygosity with a shift away from a heterozygous variant in the control sample, rejecting a binomial test with the alternative hypothesis \(H_{1}: p^{C} \neq \frac{1}{2}\).

The method that we have described before is suited for detecting variant positions efficiently in large sequencing datasets, including whole-exome and whole-genome sequencing. For this purpose, we test for a shift in non-consensus frequencies between two samples at each genomic position individually:

1. Form the base counts table for four bases A, C, G, T from the aligned reads. In order to reduce the number of false counts, we can optionally exclude reads with low base calling quality and clip the head of each read.
2. Determine the consensus sequence: In a comparative setting, we will use the most abundant base call.
3. Calculate the sequencing depth \(N^{S}_{i}\), mismatch counts \(K^{S}_{i}\), and derived statistics for both samples, based on the consensus sequence (see 2.1).
4. Find potential variant sites with a Fisher's Exact Test, comparing the number of mismatching and total bases between the samples: \({K^{T}_{i}, N^{T}_{i}, K^{C}_{i}, N^{C}_{i}}\). The p-values are corrected for multiple testing according to the Benjamini-Hochberg procedure. Only positions rejecting the null hypothesis at a significance level \(\alpha\) are further on considered as potential variants.
5. Calculate Agresti-Caffo confidence intervals with confidence level \(\beta\), in order to evaluate presence or absence of the variant (see 2.2).
6. Classify variant sites into the groups: somatic, LOH, undecided, and powerless (see 2.3).

We want to further demostrate the usage and abilities of the Rariant package on a real-life data set. Due to legal and privacy issues, most human cancer sequencing data is not publicly accessible and therefore cannot serve as an example data set here. Alternatively, we conduct an analysis to mimic the characteristics of current cancer sequencing studies.

For the purpose of the analysis, we compare three samples from the 1000 Genomes project [6], serving as a control/normal (control) and two related test/tumor samples (test and test2). Further, we simulate a clonal mixture (mix) of the two test samples by combining their reads.

library(Rariant)

library(GenomicRanges)
library(ggbio)

tp53_region = GRanges("chr17", IRanges(7565097, 7590856))

control_bam = system.file("extdata", "platinum", "control.bam", package = "Rariant", mustWork = TRUE)
test1_bam = system.file("extdata", "platinum", "test.bam", package = "Rariant", mustWork = TRUE)
test2_bam = system.file("extdata", "platinum", "test2.bam", package = "Rariant", mustWork = TRUE)
mix_bam = system.file("extdata", "platinum", "mix.bam", package = "Rariant", mustWork = TRUE)

v_test1 = rariant(test1_bam, control_bam, tp53_region, select = FALSE)

   Warning: 'rbind_all' is deprecated.
   Use 'bind_rows()' instead.
   See help("Deprecated")

v_test2 = rariant(test2_bam, control_bam, tp53_region, select = FALSE)

   Warning: 'rbind_all' is deprecated.
   Use 'bind_rows()' instead.
   See help("Deprecated")

v_mix = rariant(mix_bam, control_bam, tp53_region, select = FALSE)

   Warning: 'rbind_all' is deprecated.
   Use 'bind_rows()' instead.
   See help("Deprecated")

v_all = GenomicRangesList(T1 = v_test1, T2 = v_test2, M = v_mix)

   Warning: 'GenomicRangesList' is deprecated.
   Use 'GRangesList(..., compress=FALSE)' instead.
   See help("Deprecated")

v_all = endoapply(v_all, updateCalls)

To better understand the evidence for the presence or absence of particular variants across samples, we plot the confidence intervals, colored according to the predicted event type, and abundance shifts for all sites of interest, colored according to the sign of the shift.

t_ci = tracks(lapply(v_all, plotConfidenceIntervals, color = "verdict")) + verdictColorScale()

print(t_ci)

   Coordinate system already present. Adding new coordinate system, which will replace the existing one.
   Coordinate system already present. Adding new coordinate system, which will replace the existing one.
   Coordinate system already present. Adding new coordinate system, which will replace the existing one.
   Coordinate system already present. Adding new coordinate system, which will replace the existing one.

   Warning: Removed 6 rows containing missing values (geom_pointrange).

Confidence intervals for simulation study

t_ci = tracks(lapply(v_all, plotConfidenceIntervals, color = "eventType"))

print(t_ci)

   Coordinate system already present. Adding new coordinate system, which will replace the existing one.
   Coordinate system already present. Adding new coordinate system, which will replace the existing one.
   Coordinate system already present. Adding new coordinate system, which will replace the existing one.
   Coordinate system already present. Adding new coordinate system, which will replace the existing one.

   Warning: Removed 6 rows containing missing values (geom_pointrange).

Confidence intervals for simulation study

t_rates = tracks(lapply(v_all, plotAbundanceShift))

print(t_rates)

   Coordinate system already present. Adding new coordinate system, which will replace the existing one.
   Coordinate system already present. Adding new coordinate system, which will replace the existing one.
   Coordinate system already present. Adding new coordinate system, which will replace the existing one.
   Coordinate system already present. Adding new coordinate system, which will replace the existing one.

   Warning: Removed 6 rows containing missing values (geom_linerange).

   Warning: Removed 6 rows containing missing values (geom_point).

Abundance shifts for simulation study

In the following, we look at positions which showed a significant effect in at least one sample. This gives us 12 positions to consider in the following.

While most of the variants are somatic, i.e. they do not appear in the control sample, the last variant position shows a loss of a heterozygous SNP. Looking for example in more detail into the group of 5 variant sites around 7.85 Mbp: We can identify them as consistent with a heterozygous somatic variant in the first sample, since their 95% CIs overlap the value of 0.5. In contrast, we can show the absence of the same variants in the second sample. The third sample again shows the presence of the variants, as seen in the first case, but with lower abundance. Such a result could be expected in a mixture of subclones, in which some clones carry a somatic variant and others not. Further, we can also see the case of the next variant which consistently exists in all three samples with the same abundance.

z = filterCalls(v_all, verdict %in% c("present", "inbetween", "dontknow"))

elementNROWS(z)

   T1 T2  M 
   49 49 49

evidenceHeatmap(z, fill = "d", color = "verdict") + verdictColorScale()

By comparing the confidence intervals between the two strands, we can further detect and characterize effects such as variations in sequencing depth and strand biases. We illustrate this with a set of hypothetical cases for confidence intervals for two strands. The upper row (cases 4-7) corresponds to sites with overlapping CIs, whereas the lower row (cases 1-3) shows cases of disagreements between the CIs indicative of strand biases. When analyzing the probability for the overlap of confidence intervals, an adjustment of the confidence level has to be taken into account [8].

Illustrative cases of confidence intervals for somatic variant frequency estimates for two strands

Motivated by the analysis of different Illumina genome and exome sequencing, we consider strand-biases, in which the non-consensus base call rates differ significantly between strands at sites with sufficient sequencinq depth, a neglectable problem with current data sets and analysis pipelines (see also Best practices for short-read processing). In the presence of strand biases, pooling the counts of both plus and minus strand may be not desirable. A possible solution may be to perform a strand-specific analysis, and later combine the resulting statistics. Gerstung and colleagues discuss different approaches for combining p-values [9], in particular taking the minimum, maximum, average, or Fisher combination. These can be also applied for confidence intervals, with Fisher's method being equivalent to taking the sum of both strands.

As outlined before, an important property for assessing confidence intervals is given by their coverage probabilities. Ideally, we would expect a method to have coverage probabilities close to the nominal confidence level β over a wide range in the parameter space. Previous publications analyzing the performance focus on parameter settings that deviate from those of sequencing data sets [7]. Therefore, we perform a simulation that demonstrates the behavior of the Agresti-Caffo methods for a whole-genome sequencing study. For a fixed sequencing depth of 30 in both test and control sample, the coverage probability of 95% AC confidence intervals is computed for all possible combinations of mismatch counts \(K^{T}\) and \(K^{C}\).

## WGS
n1 = 30
n2 = 30
k1 = 0:(n1-1)
k2 = 0:(n2-1)
cl = 0.95
n_sample = 1e4

pars = expand.grid(k1 = k1, k2 = k2, n1 = n1, n2 = n2, conf_level = cl)

cp_ac = coverageProbability(pars, fun = acCi, n_sample = n_sample)

p_ac = ggplot(cp_ac) + geom_tile(aes(x = k1, y = k2, fill = cp)) +
    scale_fill_gradient2(midpoint = 0.95, limits = c(0.9, 1)) +
    theme_bw() + xlab("kT") + ylab("kC")

print(p_ac)

Coverage probabilities for whole-genome setting

For mismatch rates close to 0 or 1 in both samples, the Agresti-Caffo method shows a conservative perfomance.

For an analysis of two matched human tumor samples, we performed a benchmark to assess the computational time and memory usage on a standard laptop (Thinkpad X220 built in 2011). Both samples contain about 95M reads mapped to the 1000genomes reference sequence reads that are considered in the analysis. For the analysis of chromosome 22, the analysis with default parameters required ~873s and 600MB of RAM. For an analysis of all linear toplevel chromosomes (autosomes and allosomes), this would require ~15h of time. Please consider that the current version of Rariant is under active development and computational efficiency will increase with newer versions.

We welcome emails with questions or suggestions about our software, and want to ensure that we eliminate issues if and when they appear. We have a few requests to optimize the process:

• All emails and follow-up questions should take place over the Bioconductor mailing list, which serves as a repository of information.
• The subject line should contain Rariant and a few words describing the problem. First search the Bioconductor mailing list, for past threads which might have answered your question.
• If you have a question about the behavior of a function, read the sections of the manual page for this function by typing a question mark and the function name, e.g. ?rariant. Additionally, read through the vignette to understand the interplay between different functions of the package. We spend a lot of time documenting individual functions and the exact steps that the software is performing.
• Include all of your R code related to the question you are asking.
• Include complete warning or error messages, and conclude your message with the full output of sessionInfo().

Before you want to install the Rariant package, please ensure that you have the latest version of R and Bioconductor installed. For details on this, please have a look at the help packages for R and Bioconductor. Then you can open R and run the following commands which will install the latest release version of Rariant:

source("http://bioconductor.org/biocLite.R")
biocLite("Rariant")