We evaluate and demonstrate GeTallele's functionality using matched whole exome and RNA sequencing datasets from paired normal and tumor tissue obtained from 72 female patients with breast invasive carcinoma (BRCA) from TCGA. Each dataset contains four matched sequencing sets: normal exome (Nex), normal transcriptome (Ntr), tumor exome (Tex), and tumor transcriptome (Ttr) (see Supplementary Table 1). The raw sequencing data were processed as previously described (Movassagh et al., 2016) to generate the inputs for GeTallele.

In short, all datasets were generated through paired-end sequencing on an Illumina HiSeq platform. The human genome reference (hg38)-aligned sequencing reads (Binary Alignment Maps, bams) were downloaded from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/) and processed downstream through an in-house pipeline. After variant calling (Li, 2011), the RNA-seq and whole exome sequencing (WES) alignments, together with their respective variant calls, were processed through the read count module of the package RNA2DNAlign (Movassagh et al., 2016), to produce variant and reference sequencing read counts for all the variant positions in all four sequencing signals (normal exome, normal transcriptome, tumor exome and tumor transcriptome). Selected read count assessments were visually examined using the Integrative Genomics Viewer (Thorvaldsdóttir et al., 2013).

For each sample, to select SNV positions for analysis, we start with heterozygous SNV calls in the normal exome (Li et al., 2009). In each of these positions, we estimate the counts of the variant and reference reads (n_VAR and n_REF, respectively) across the 4 matching datasets, and retain positions covered by a minimum total (variant + reference) read depth for further analyses. This threshold is flexible and is required to ensure that only sufficiently covered positions will be analyzed; it is set to 3 in the herein presented results. For further analysis (without loss of generality), we transform all the original VAF values to VAF = |VAF−0.5|+0.5. We introduce this transformation due to the symmetric nature of the VAF distributions.

In addition, we required each tumor sample to have at least three of the following five purity estimates—Estimate, Absolute, LUMP, IHC, and the consensus purity estimate (CPE) (Katkovnik et al., 2002; Pagès et al., 2010; Carter et al., 2012; Yoshihara et al., 2013; Zheng et al., 2014; Aran et al., 2015). On the same datasets, we applied THetA (Oesper et al., 2013, 2014)—a popular tool for assessing CNA and admixture from sequencing data—was also applied to the datasets.

To test statistical significance, GeTallele uses parametric and non-parametric methods and statistical tests (Hollander et al., 2013; Corder and Foreman, 2014). Namely, to compare distributions of the variant allele frequencies (VAF) we use the Kolmogorov–Smirnov test (examples of VAF distributions are depicted in Figures 2, 3). To study concurrence of windows, we use permutation/bootstrap tests. To test relations between v_PR and copy number alterations (CNA), we use Pearson's correlation coefficient.

To account for multiple comparisons, we set the probability for rejecting the null hypothesis at p < 1e−5, which corresponds to Bonferroni (Dunn, 1961) family-wise error rate (FWER) correction against 100,000 comparisons. We use a fixed value, rather than other approaches, to ensure better consistency and reproducibility of the results. Alternatively, we apply Benjamini and Hochberg (Benjamini and Hochberg, 1995) false discovery rate (FDR) correction with a probability of accepting false positive results p_FDR <0.05. We specify the method used in the text when reporting the results.

To analyse variant allele frequencies (VAF) at genome-wide level, GeTallele first divides the VAF sequence into a set of non-overlapping segments along the chromosomes. To partition the data into segments, GeTallele uses a parametric global method, which detects the breakpoints in a signal using its mean, as implemented in the Matlab function findchangepts (Lavielle, 2005; Killick et al., 2012) In each segment, the VAFs of the chosen signal must have a different mean than that of the adjacent segment. In the Matlab implementation, sensitivity of breakpoint detection can be controlled using parameter MinThreshold; with a default setting of 0.2. Segments containing fewer than 10 data points were merged with the preceding segment. For analysis of matched signals, segmentation is based on one signal, and then applied to the others. In the presented analysis, segmentation is based mainly on Tex dataset, for comparison we also use Ttr dataset (dataset used for segmentation is specified in the description of the results presented in Section Results).

Variant probability is a biologically interpretable quantitative descriptor of the VAF distribution. It is the common probability of observing a variant allele at any site in a given chromosomal segment. The v_PR is a measure describing the genomic event that, through the sequencing process, was transformed into an observed distribution of VAFs. For example, in VAF_DNA from a diploid genome, we assume variant probability v_PR = 0.5 (meaning that both alleles are equally probable) corresponds to a true allelic ratio of 1:1 for heterozygous sites. The value might differ from 0.5 due to reference mapping biases (Degner et al., 2009). For heterozygous sites in the DNA from a diploid monoclonal samples, the corresponding tumor VAF_DNA is expected to have the following interpretations: v_PR = 1 or v_PR = 0 corresponding to a monoallelic status resulting from a deletion, and v_PR = 0.8 (or 0.2), 0.75 (or 0.25), 0.67 (or 0.33) corresponding to allele-specific tetra-, tri-, and duplication of the variant-bearing allele, respectively.

The v_PR of the VAF_RNA is interpreted as follows. In positions corresponding to heterozygote sites in DNA, alleles not preferentially targeted by regulatory traits are expected to have expression rates with variant probability v_PR = 0.5, which (by default) scale with the DNA allele distribution. Differences between VAF_DNA and VAF_RNA values are observed in special cases of transcriptional regulation where one of the alleles is preferentially transcribed over the other. In the absence of allele-preferential transcription, VAF_DNA, and VAF_RNA are anticipated to have similar v_PR across both diploid (normal) and copy number altered genomic regions. Consequently, VAF_DNA, and VAF_RNA are expected to synchronously switch between allelic patterns along the chromosomes, with the switches indicating breakpoints of DNA deletions or amplifications.

Since we observed that DNA and RNA signals have different distributions of total reads and also that the distributions of total reads vary between participants, the synthetic VAF distributions are generated individually for each sequencing signal and each participant.

To estimate v_PR in the signals, GeTallele first generates synthetic VAF distributions and then uses the earth mover's distance (EMD) (Kantorovich and Rubinstein, 1958; Levina and Bickel, 2001) to fit them to the data. To generate a synthetic VAF distribution with a given variant probability, v_PR, GeTallele, bootstraps 10,000 values of the total reads (sum of the variant and reference reads; n_VAR + n_REF) from the analyzed signal in the dataset. It then uses binomial pseudorandom number generator to get number of successes for given number of total reads and a given value of v_PR (implemented in the Matlab function binornd). The v_PR is the common value of the probability of success and generated number of successes is interpreted as an n_VAR. Since the v_PR of the synthetic sample can take any value, it can correspond to a single genomic event as well as any combination of genomic events in any mixture of normal and tumor populations (See section v_PR Values in Mixtures of Normal and Tumour Populations).

The analysis presented in the paper uses 51 synthetic VAF distributions with v_PR values that vary from 0.5 to 1 with step (increment of) 0.01. The synthetic VAF distributions are parametrized using only v_PR ≥ 0.5, however, to generate them we use v_PR and its symmetric counterpart 1-v_PR. The process of generating synthetic VAF distributions along with examples of synthetic and real VAF distributions with different values of v_PR are illustrated in Figure 3.

To estimate v_PR, we compute the Earth mover's distance between the distribution of VAF values in the considered window and the 51 synthetic VAF distributions (i.e., observed vs. synthetic VAF). The estimate is given by the v_PR of the synthetic VAF distribution that is closest to the VAF distribution in the segment.

Earth mover's distance (EMD) is a metric for quantifying differences between probability distributions (Kantorovich and Rubinstein, 1958; Levina and Bickel, 2001) and in the case of univariate distributions it can be computed as: EMD(PDF1,PDF2)=∫Z|CDF1(z)CDF2(z)|dz. Here, PDF₁ and PDF₂ are two probability density functions, and CDF₁ and CDF₂ are their respective cumulative distribution functions. Z is the support of the PDFs (i.e., set of all the possible values of the random variables described by them). Because VAFs are defined as simple fractions with values between 0 and 1, their support is given by a Farey sequence (Hardy and Wright, 2008) of order n; n is the highest denominator in the sequence. For example, Farey sequence of order 2 is 0, 1/2, 1, and Farey sequence of order 3 is 0, 1/3, 1/2, 2/3, 1. We use a Farey sequence of order 1,000 as the support Z for estimating the v_PR.

Examples of VAF distributions with fitted synthetic VAF distributions are shown in Figures 3A,D. The dependence of the confidence intervals of the estimation on the number of VAF values in a segment is illustrated in Figure 4, which clearly demonstrates that the accuracy of the estimate is positively correlated with the number of VAFs in the chosen segment.

Since the v_PR can take any value between 0.5 and 1 it can correspond to a single genomic event as well as any combination of genomic events in any mixture of normal and tumor populations. A mixture v_PR value that corresponds to a combination of genomic events can be computed using the following expression:

Where e_VAR and e_REF are the multiplicities of variant and reference alleles and p_PL is a proportion of one of the populations. For heterozygote sites e_VAR = 1 and e_REF = 1, for deletions e_VAR = 0 or e_REF = 0, for du-, tri- and tetraplications e_VAR or e_REF can be equal to 2, 3 or 4, respectively. The sum of proportions p_PL over the populations is equal 1. For example, for a mixture of 1 normal (N, p_N = 0.44) and 2 tumor populations (T1, p_T1 = 0.39 and T2, p_T2 = 0.17), T1 with deletion and T2 with deletion the mixture v_PR value can be computed as follows:

By comparing the v_PR values estimated from data with possible mixture v_PR values we propose to estimate sample purity and its clonal composition. To this end, we first generate a full set of proportions of all the population in the mixture with step (increment of) 0.01 and compute all the possible v_PR values that each of the mixtures could produce. For step 0.01: two populations (1 tumor) give 99 proportions, three populations (2 tumors) give 4,851 proportions, four populations (3 tumors) give 156,849 proportions. The matrices with mixture v_PR values for each proportion, vary from 2 × 2, for two populations with deletions, to 35 × 35 for four populations with all events up to tetra-plications. Then, we run an exhaustive approximate search over all the matrices with mixture v_PR values over all the proportions. The search is approximate because the estimated v_PR values have limited accuracy and because we consider only discrete values of proportions. In the analysis we define a match between estimated and mixture v_PR values if they differ by <0.009 (we chose a value that is smaller than the smallest difference between possible v_PR estimates). The search returns a large number of admissible mixtures that could produce the estimated v_PR values. This process is illustrated in Figure 5.

To visualize the admissible mixtures, we use ternary plots, which allow us to illustrate composition of three components in two dimensions. The composition, represented by ratios of the three components, which sum to a constant, is depicted as point inside or on the edge of an equilateral triangle. If the point is on the edges, the composition has only two components. To help interpretation of the ternary plots, we also plot the grid lines that are parallel to the sides of the triangle. These gridlines indicate the directions of constant ratios of the components. Along such direction the ratio of one of the components is fixed and only the other two ratios vary. Examples of visualization of admissible mixtures on ternary plots are shown in Figures 5, 6.

To facilitate analysis of the admissible mixtures returned by the search procedure we introduce mixture complexity. Mixture complexity is a measure that increases with number of populations as well as with variety of genetic events. From the simplest mixture of 1 normal and 1 tumor population in which only deletions are possible to a model with 1 normal and multiple tumor populations where each can have deletions, and any level of multiplications. In practice, we set the limit at 3 tumor populations and tetra-plications. Mixture complexity helps to group and visualize admissible mixtures. Mixtures with higher complexity allow more possible v_PR values, meaning that it is easier to find the match with the estimated v_PR values but that the number of admissible mixtures increases (see Figure 6). We, further, observe that proportion of normal population, p_N, increases with a number of clonal tumor populations included in the model mixture and that, generally, p_N stays constant with increasing variety of genetic events, for a fixed number of clonal tumor populations. We note that this is just one of many possible ways of deciding which solution should be chosen.

Segmentation of the data, based on the tumor exome signal, resulted in 2,697 chromosomal segments across the 72 datasets. We excluded from further analysis 294 chromosomal segments where either tumor exome or transcriptome had v_PR > = 0.58 but their VAF distribution could not be differentiated from the model VAF distributions with v_PR = 0.5 (p > 1e−5, Kolmogorov Smirnov test, equivalent to Bonferroni FWER correction for 100,000 comparisons). The 294 excluded chromosomal segments, corresponding to 4% of the total length of the data in base pairs and 4% of all the available data points. This implies these short segments containing few VAF values. In the remaining 2,403 chromosomal segments, we systematically examined the similarity between corresponding VAF_TEX (tumor exome), VAF_TTR (tumor transcriptome), and CNA. We obtained several distinct patterns of coordinated RNA-DNA allelic behavior as well as correlations with CNA data.

In 60% of all analyzed chromosomal segments the distributions of VAF_TEX and VAF_TTR were statistically concordant (P > 1e−5, Kolmogorov Smirnov test), and in 40% they were statistically discordant (P < 1e−5, Kolmogorov Smirnov test). In two chromosomal segments, VAF_TEX and VAF_TTR, had the same v_PR, while having statistically different VAF distributions (P < 1e−5, Kolmogorov Smirnov test). We consider such chromosomal segments as concordant. The v_PR robustly characterizes VAF sample while the Kolmogorov-Smirnov test is very sensitive for differences between distributions that might be caused by to technical variance. In the vast majority of the discordant chromosomal segments v_PR of the VAF_TTR, v_PR,TTR, was higher than v_PR of the VAF_TEX, v_PR,TEX, (only in 21 out of 959 discordant chromosomal segments v_PR,TTR was lower than v_PR,TEX).

We next analyzed the concurrence between chromosomal segments resulting from independent segmentations of the tumor exome (VAF_TEX) and transcriptome (VAF_TTR) datasets (2,697 and 3,605 chromosomal segments, respectively, across all the samples). We first assessed chromosome-wise alignment of the start and end points of the chromosomal segments. In 45% of the chromosomes both VAF_TEX and VAF_TTR signals produce a single segment that contains the whole chromosome. In 33% of chromosomes both signals produced multiple chromosomal segments. These chromosomal segments are well aligned, with 90% of the breakpoints differing <7% of data points in the chromosome, e.g., they are <70 points apart if the chromosome contains 1,000 data points; Q50 = 0.02%, Q75 = 2% of data points in the chromosome. The probability of observing such an alignment by chance is smaller than p = 1e−5 (100,000 bootstrap samples with breakpoints assigned randomly in all the individual chromosomes where both signals produced multiple chromosomal segments). In 22% of the chromosomes, segments based on VAF_TEX and VAF_TTR signals were positionally discordant—one signal produced a single segment containing whole chromosome while the other produced multiple chromosomal segments.

To compare the v_PR values in the 55% of chromosomes where at least one signal produced more than one chromosomal segment, we computed chromosome-wise mean absolute error (MAE) between the v_PR in two sets of chromosomal segments. To account for different start and end points of the segments we interpolated the v_PR values (nearest neighbor interpolation) at each data point in the chromosome. We separately compared the v_PR,TEX and v_PR,TTR values. Assessment of alignment using MAE showed strong concordance: v_PR,TEX agreed perfectly in 11% of the chromosomes and had the percentiles of MAE equal to Q50 = 0.012, Q75 = 0.022 and Q97.5 = 0.047, while v_PR,TTR agreed perfectly in 8% but had slightly higher percentiles of MAE Q50 = 0.019, Q75 = 0.034 and Q97.5 = 0.07. v_PR,TEX and v_PR,TTR values had MAE = 0 simultaneously in 4% of the chromosomes. Probability of observing such values of MAE by chance is smaller than p = 1e−3 (1,000 random assignments of v_PR,TEX and v_PR,TTR values to windows in the 873 chromosomes where at least one signal had more than one chromosomal segment). It is noteworthy that MAE Q97.5 < 0.07 is comparable with the confidence interval of single v_PR estimate based on 50 VAF values. In other words, both signals in a sample (Tex and Ttr) give very similar results in terms of segmentation and estimated values of the v_PR. Albeit, segmentation of VAF_TTR generates a higher number of chromosomal segments. The higher number of VAF_TTR chromosomal segments indicates that transcriptional regulation occurs at a smaller scale than alterations in DNA. Figure 7 shows examples of concurrence between chromosomal segments based on VAF_TEX and VAF_TTR signals in a positionally concordant chromosome (both signals produced multiple segments).

Finally, we assess the correlations between v_PR and CNA in the individual samples. We separately computed correlations for deletions and amplifications. In order to separate deletions and amplifications, for each data set we found CNA_MIN, value of the CNA in the range −0.3 to 0.3 that had the smallest corresponding v_PR,TEX. To account for observed variability of the CNA values near the CNA_MIN, we set the threshold for amplifications to CNA_{AMPLIFICATION} = CNA_MIN-0.05, and for deletions we set it to CNA_DELETION = CNA_MIN + 0.05 (each data set had a different threshold).

For VAF_TEX, we observed significant correlations with negative trend between v_PR,TEX and CNA ≤ CNA_DELETION in 57 datasets and with a positive trend between v_PR,TEX and CNA ≥ CNA_{AMPLIFICATION} in 39 datasets (p_FDR < 0.05, Pearson's correlation with Benjamini Hochberg multiple comparison correction for 72 samples). For VAF_TTR, we observed significant correlations with a negative trend between v_PR,TTR and CNA ≤ CNA_DELETION in 62 datasets and with positive trend between v_PR,TTR and CNA ≥ CNA_{AMPLIFICATION} in 33 datasets (p_FDR < 0.05, Pearson correlation with Benjamini Hochberg correction). These correlations indicate that the segmentation and the estimated v_PR values are concordant with CNA calls. However, the v_PR values (estimated at the level of chromosomal segments) do not differentiate between positive and negative values of the CNA, meaning it is not possible to use v_PR alone to call amplifications and deletions.

Figure 8 shows four typical patterns of correlation between the CNA and v_PR values observed in the data. In Figure 8A, all the values of CNA are close to CNA_MIN. In Figure 8B, the relationship between CNA and v_PR is noisy, only correlations between v_PR,TTR and CNA ≤ CNA_DELETION are statistically significant (r_TEX,CNA,DEL = −0.29, p_FDR = 0.063; r_TEX,CNA,DEL = −0.38, p_FDR = 0.012; r_TEX,CNA,AMPL = 0.14, p_FDR = 0.58; r_TEX,CNA,AMPL = 0.19, p_FDR = 0.47; Pearson's correlation with Benjamini Hochberg multiple comparison correction for 72 samples). In Figure 8C all the correlations are statistically significant, v_PR,TTR values (circles) follow closely the v_PR,TEX (squares) indicating that in most of the windows distributions of the VAF_TEX and VAF_TTR are concordant (r_TEX,CNA,D = −0.91, p_FDR < 1e−10; r_TEX,CNA,DEL = −0.96, p_FDR < 1e−10; r_TEX,CNA,AMPL = 0.92, p_FDR < 1e−10; r_TEX,CNA,AMPL = 0.95, p_FDR < 1e−10). In Figure 8D correlations between v_PR,TEX, v_PR,TTR and CNA ≤ CNA_D are statistically significant, but there is a large difference (with median of 0.18) between v_PR,TEX and v_PR,TTR values, indicating that in most of the windows the distributions of the VAF_TEX and VAF_TTR in this dataset are discordant (r_TEX,CNA,DEL = −0.44, p_FDR = 0.047; r_TEX,CNA,DEL = −0.64, p_FDR = 0.0017; r_TEX,CNA,AMPL = 0.44 p_FDR = 0.16; r_TEX,CNA,AMPL = 0.28, p_FDR = 0.41). In many of the datasets we observe that the v_PR,TTR values are higher than the corresponding v_PR,TEX values (median v_PR,TTR-v_PR,TEX = 0.03), likely indicative of preferential transcription of some alleles in the chromosomal segment. Correlations between v_PR and CNA in all datasets are shown in the Supplementary Figure 1.

To demonstrate a practical application of the v_PR values we use them to estimate tumor purity of the samples. To this end we compared the v_PR based purity (VBP) estimates with ESTIMATE, ABSOLUTE, LUMP, IHC, and the Consensus Purity Estimation (CPE) (Katkovnik et al., 2002; Pagès et al., 2010; Carter et al., 2012; Yoshihara et al., 2013; Zheng et al., 2014; Aran et al., 2015).

To obtain the VBP estimate we used v_PR,TEX values. We, first, selected the v_PR,TEX values that: 1. are estimated with high confidence, i.e., are based on at least 50 VAF values; 2. are most likely heterozygous in normal exome, i.e., have a corresponding v_PR value in normal exome v_PR,NEX <0.58; 3. most likely have v_PR,TEX > 0.5, i.e., their p-value for comparison with v_PR,TEX = 0.5 is very small p < 1e−5 (Kolmogorov-Smirnov test).

Next, we used the selected v_PR,TEX values to find all admissible mixtures (with 1–3 tumor populations and allowing for all events, from deletions to tetraplications). To estimate the VBP, out of all the admissible mixtures we chose these with lowest mixture complexity and among these mixtures we take one with the highest p_N (proportion of the normal population). The VBP, percentage of tumor populations in the sample, is then given as 1-p_N. Such approach provides rather conservative estimates of VBP (the smallest 1-p_N). However, GetAllele can be extended to offer alternative methods of employing the admissible mixtures to estimate VBP. Development, analysis and comparison of alternative VBP estimation methods is beyond scope of the current paper.

Figure 9A shows violin plots of all considered 1-p_N values and (x) indicates the smallest value taken as a VBP estimate. In two of the datasets we could not estimate the purity due to lack of suitable v_PR,TEX values. The VBP estimates shows the best agreement with ABSOLUTE method (y = 0.86 x + 0.02, r = 0.76, p < 3.4e−14, Pearson's correlation, Figure 9B2). We suppose that this is because the ABSOLUTE method is based on copy number distributions, and our analysis (Section Correlation Between v_PR and CNA) revealed high correlations between the CNAs and v_PR values. Similar, to the ABSOLUTE method, VBP estimates are generally lower than the other purity estimates (ESTIMATE, LUMP, IHC, CPE); see Figures 9B1–5.

The approach presented in this section differs from other methods for inferring genomic mixture composition in that it is based on chromosomal segments with at least 50 VAF values which can extend over millions of base pairs. In contrast, PyClone (Roth et al., 2014) is based on sets of carefully selected individual deeply sequenced VAF values, while SciClone (Miller et al., 2014) and TPES (Locallo et al., 2019) are based on analysis of selected VAF values aggregated from the genome-wide sequences (multiple chromosomes). By using chromosomal segments, v_PR allows for a more granular description of the VAF distributions than aggregating genome-wide VAF values. At the same time, basing purity estimation on v_PR values allows for the use of SNVs with a low sequencing depth (3 in the presented analysis). Rigorous comparison of the performance of the different methods is beyond the scope of this demonstration of potential practical applications of v_PR.