Detecting copy number status and uncovering subclonal markers in heterogeneous tumor biopsies

De-mixing of samples

In a homogeneous sample s, the exact genotype of each cell in the sample can be represented by the pair of vectors, A and B, containing respectively the number of copies of the A- and B-allele at each locus. For the sake of simplicity, we assume that each locus is at most bi-allelic, namely that only two variants, A and B, can occur. In the case of tumor biopsies, or studies in which cells from several individuals are pooled together, the measured and vectors correspond to the average copy number of the A- and B-alleles at each locus across all cells in the sample. When studying heterogeneous samples, it is very unlikely to find any cell in the sample whose genotype corresponds exactly to the measured profile.

The "de-mixing problem" is the task of recovering the fraction x _i , called the mixing coefficient, and the genotype g _i = {A _i , B _i } of each i-th homogeneous group of cells in the sample (e.g. the different subclones). This problem corresponds to solving the linear system

(1)

,

where x _i is an element of the vector X, j is the locus index, a _ji and b _ji are elements of the matrices A and B respectively, and N ₀ is the set of all non-negative integers.

Deconvoluting genotyping data signal from heterogeneous cell population

The de-mixing problem is trivially solved when only two independent components (e.g. from two different individuals) are mixed together, and their copy number, a _ji + b _ji , is measured using independent experimental approaches in at least two loci. However, in general, deconvolution of CNA signals from a heterogeneous sample to its homogeneous components is an ill-posed mathematical problem. Linear algebra solution approaches will encounter non-unique solutions corresponding to a non-unique combination of possible discrete vectors representing different subclonal aberration profiles.

Regularization methods (such as L1 minimization for recovery of sparse signals) are often employed to reduce the number of feasible solutions in underdetermined problems. The underlying assumption of applying regularization in this context is that a sample will not generally contain all theoretical subclones, defined by all possible combinations of aberrations. Instead, a few subclones arising during cancer progression will be selected. This assumption is also consistent with the possibility that subclones' viability may be hindered by the lethal synergy of some of their aberrations. Donoho and Tanner reported key results on regularization methods from counting faces of randomly projected polytopes when the projection radically lowers dimension [9]. It is possible to show that even in the simplified problem with a _ji ∈ {0,1} b _ji ∈ {0,1}, the required condition for sparse reconstruction,

(2)

is not fulfilled for any s > 2, where 2^s-1 is the size of the search space, with s equal to the number of loci measured, k is the number of components (e.g. observable subclones) and s+1 is the maximum number of equations, corresponding to the number s of unique measurements (e.g. SNPs), including the additional constraint that the sum of mixing coefficients is 1. Thus, a mixture of 3 or more components, which is typical in tumor biopsies, cannot be uniquely deconvolved solving the linear system, or using sparse reconstruction methods. In the presence of heterogeneous samples no inference can be made of the exact copy number from the aggregate signal of the subclones.

To illustrate these points, we provide a detailed example of the conditions and solution of copy number inference for a mixture of a tumor subclone component with its matched stromal component. In this scenario, de-mixing is elicited by additional constraints on the copy number of the stromal and of the tumor components. The stromal component is assumed to have no CNA, thus its number of copies is constant and equal to 2. The linear system in Eq.1 of the measured aggregate genotypes at the j-th locus is reduced to:

(3)

,

where x is the mixing coefficient of the stromal component, a _j1 and b _j1 are the A- and B-allele copy number of the tumor, a _j2 is the A-allele copy number of the stromal component and N ₀ is the set of all non-negative integers. In addition, we include the equation a _j1 +b _j1 = m _j in the system in order to restrict the copy number of the tumor component at the j-th locus to an integer m _j . We consider only the loci harboring deletion events with an aggregate copy number smaller than two, where m_j - the copy number of the tumor component - can take one of two values: one if the deletion in the tumor cells is hemizygous, or zero if it is homozygous. We then solve the system at each locus independently for the two possible values of m _j and determine the two corresponding possible values of x. Since there is only one tumor component, there must exist one value of x to solve Eq.3 for all loci simultaneously. The mode of all values of x across all loci corresponds to this value. The corresponding value of m _j at the j-th locus is then the exact copy number of the tumor component.

A simple framework for CNA analysis

CNA analysis encompasses the tasks of detecting and classifying copy number aberrations. Current HMM based applications for CNA analysis rely on the assumption that increasing the number of hidden states, thereby increasing the complexity of the underlying model, will eventually result in higher CNA classification accuracy [10–12]. As detailed above, determination of the exact number of copies is an ill-posed problem in the presence of heterogeneous samples. Therefore, we propose to detect substantial deviation from the normal diploid state by classifying the aberrations as gains or losses inferred from the dominant component in the mixture. The goal of CNA analysis then becomes to infer the most recurring aberration that would be observed in a hypothetical, high-resolution karyotyping of the same heterogeneous sample, if it were available.

Peiffer et al. [13] employed a convenient transformation of A- and B- allele copy numbers to B- allele frequency β and the total-DNA enrichment ρ. For clarity of discussion, we use a simpler version of this transformation for each SNP defined as:

(4)

,

where N is the total number of cells in the sample, are the allele counts estimated from the SNP-arrays. Peiffer and et al. suggest estimating this quantity from the signal of a normal population (e.g. HapMap samples [5]). In general, occurrence of CNA is reflected in the total-DNA enrichment. In the presence of heterogeneous samples, these changes can be difficult to detect due to the high level of noise. As expected, in heterogeneous samples, aberrations also lead to allelic imbalance at positions of heterozygous SNPs (Figure 1). We thus propose to base CNA detection in heterogeneous samples on measures of allelic imbalance.

One such measure of allelic imbalance is

(5)

,

where W corresponds to a window of appropriate size. This quantity is an aberration score for each SNP, based on the measure of allelic imbalance in its neighborhood. We denote this measure as the M-measure and use it for the remainder of the analyses with W = 20. W = 20 is an arbitrary choice of a window large enough to have a reasonably robust estimate of the M-Measure. The M-measure is insensitive to hemizygous deletion events if they occur in the majority of cells in the sample. We address this issue by testing whether the total-DNA enrichment ρ, also known as R-ratio, is significantly below its expected value of 1. We use the M-measure to classify SNP CNA profiles into three states: gain, loss and neither. The M-measure is a measure of deviation from the expected state of normal non-tumor sample. In order to have a measure that is robust to noise, we chose trigonometric functions such that their product is close to zero in a wide region around the value of the expected state of a normal, non-tumor sample. These trigonometric functions have a non-linear rapid change of the deviation score once the threshold of having a true gain or loss is crossed. As previously discussed, in this context there are numerous solutions to the problem of inferring the exact number of copies consistent with the measured aggregate cell population signal (Eq.1). As none of the possible solutions can be preferred to the others, characterization of CNAs in terms of gain and loss is a practical alternative.

The three-state M-measure classification can also be used in a three-state Viterbi algorithm whose hidden states are normal, deletion, and duplication segments. We name this classification approach 3SMM. We use 3SMM only in our simulated scenarios. In order to use the best parameters for our simulations, we determined transition and emission probabilities from the underlying true copy number status of the simulated data for each run.

Simulated data

We randomly generated two haplotypes of 10,000 loci to simulate a mixture of a homogeneous tumor subclone and its matched stromal cells. A haplotype is a vector of zeros and ones, representing the two alleles respectively. Our random generation process ensures that at each locus, the germline had a 90.5% probability of being homozygous and a 9.5% probability of being heterozygous. The stromal (non aberrant) component of the mixture was obtained by combining the two haplotypes as described in Eq.1. The tumor (aberrant) component was obtained by adding aberrations to the two haplotypes. To have a unique profile of aberration status, we excluded simultaneous occurrence of different kinds of aberrations (e.g. deletion and duplication) at the same locus. We mixed the normal and aberrant components with proportions of α and (1-α) respectively. Gaussian noise with a signal-to-noise ratio of 30 was used to simulate experimental noise.

We then simulated samples comprising two tumoral subclones and their matched stromal component. As in the case with one subclone, we excluded simultaneous occurrence of different kind of aberrations (e.g. deletion and duplication) at the same locus in the same subclone (e.g. neutral LOH); however, the two tumor components were generated such that their aggregate signal contained all combinations of aberrations (e.g. deletion and duplication). The three components were then mixed by a linear combination with proportions of α, 2(1-α)/3 and (1-α)/3 corresponding to the normal component and the remaining two tumoral subclones respectively. As in the case with two components, Gaussian noise with a signal-to-noise ratio of 30 was used to simulate experimental noise.

We tested four different types of aberrations: homozygous deletion, hemizygous deletion, a gain of one copy, and a gain of two copies. We did not test neutral LOH since it is not an event that would affect the copy number. Each aberration covers 1,000 loci (~20 Mbps). Aberrations were separated by 1,000 aberration-free loci (~20 Mbps). Our aim was to test the classification performance of different approaches to heterogenous-like data rather than testing the sensitivity to detect short focal events. This motivated our choice of aberrations of 1000 loci each.

Measure of performance

To assess algorithmic performance we used the Area Under the Receiver-Operator Characteristic Curve (AUCROC) measure. In the present study, which is different from the usual formulation, in which the classification task is to distinguish between a set of positive instances and a set of negative instances, we had three class labels: gain, loss and normal. In a typical machine learning scenario in which there is no parameter to vary, the AUCROC in the operative point (which is also termed balanced accuracy) is computed as

(6)

,

where TP is the number of true positives, P is the number of positives, TN is the number of true negatives and N is the number of negatives. Similarly, in our three-classes scenario, the AUCROC can be computed as

(7)

,

where TP _gain is the number of SNPs that are amplified and correctly inferred as such by the algorithm, P _gain is the total number of SNPs that are amplified in the sample, TN _gain is the number of SNPs that are not amplified and that are correctly inferred as such by the algorithm, N _gain is the total number of SNPs that are not amplified in the sample; TP _loss is the number of SNPs that are deleted and correctly inferred as such by the algorithm, P _loss is the total number of SNPs that are deleted in the sample, TN _loss is the number of SNPs that are not deleted and that are correctly inferred as such by the algorithm, N _loss is the total number of SNPs that are not deleted in the sample; TP _normal is the number of SNPs that are not aberrated and that are correctly inferred as such by the algorithm, P _normal is the total number of SNPs that are not aberrated in the sample, TN _normal is the number of SNPs that are aberrated and correctly inferred as such by the algorithm, and N _normal is the total number of SNPs that are aberrated in the sample.

SNP profiling using microarrays

DNA from 30 melanoma cell lines were hybridized to Illumina's Human1M BeadChip (Illumina Inc. San Diego, CA). We generated B-allele frequencies and Log-R ratios using standard procedures included in the Illumina BeadStudio package. We normalized with respect to the population of western European ancestry (CEU) from the HapMap project that was analyzed on the Illumina Human1M BeadChip.

Design, probe annotation, and data processing of the arrays for detection of genome-wide gene expression

We used NimbleGen genome-wide human expression arrays (2005-04-20_Human_60mer_1in2) with a total of < 400,000 probes for < 30,000 transcripts and < 20,000 known genes, as of NimbleGen annotations. NimbleGen provides design and probe annotation. Standard methods for one-channel and two-channel microarrays from the R statistical software were used as previously described [14].

Transcriptome profiling using next-generation sequencing

We re-analyzed the RNA-seq sample MeWo from a recent study [15, 16]. Namely, the reads were aligned to the reference genome using Bowtie [17] with standard parameters. Nucleotide variations were determined after pileup using Samtools [18], and the frequency of the variant, β, was calculated as in Eq.4.

Fluorescent in situ hybridization (FISH)

Fluorescence in situ hybridization (FISH) was performed using probes from the bacterial artificial chromosome (BAC) clones (RPCI-11 human BAC library) containing the selected genes E2F8 (248D22 and 80B10 at 11p15.1), ETV4 (100E5 and 147C10 at 17q21.31), EZH2 (140E16 and 24N19 at 7q36.1) and FAM84B (455K11 and 90G11 at 8q24.21). All BAC clones were cultured in 100 ml LB media supplemented with chloramphenicol at 37°C shaker incubator overnight, and cell pellets collected by centrifuge were used for DNA extraction using the large-construct kit (Qiagen, Valencia, CA). Two BAC clones for the 5'-end or the 3'-end of each gene were labeled differently by SpectrumGreen-dUTP or SpectrumOrange-dUTP using the nick translation kit (Abbott Molecular, Des Plaines, IL). Probe hybridization on slides of interphase cells was performed following the laboratory's standardized protocol. Hybridization signals were visualized and captured using an Olympus BX60 fluorescence microscope with CytoVision software version 4.5.2 (Genetix, San Jose, CA). In each sample, 200 nuclei were inspected and the signal patterns were documented.

The measure of allelic imbalance (M-measure) is robust to heterogeneity

We performed simulations to study the behavior of the allelic imbalance M-measure in presence of aberrations (Figure 1C and Figure 1D). A simple threshold procedure can be used to identify high confidence CNA loci. An arbitrary choice of 0.1 for the cutoff of the M-measure and window size W = 20 is sufficient to achieve satisfactory accuracy (Figure 1). Remarkably, the M-measure is robust in detecting aberrations even when the aberrated component is present at low concentrations (Figure 1).

Comparison of the performance of the M-measure to that of state-of-the-art HMM-based CNA methods requires data in which the subclonal composition and copy number of each subclone are known. Currently, there is no practical solution to catalog all aberrations in all clones in a given sample, and we therefore used simulated data to test the accuracy of CNA classification of complex mixtures. Following the binary mixture procedure described in the Methods section, we generated 200 independent datasets for selected values of the mixing coefficient α to simulate a scenario of contaminating a homogeneous tumor sample (composed of only one subclone) with stromal cells. This binary mixture scenario should not be confused with the one used by Peiffer et al [13]. Peiffer et al. mixed tumor samples of one individual with normal cells from another individual, generating a mixed sample that was not reflective of clinical settings. In our simulations, we mixed components from the same individual, reflecting stromal contamination, which is often present in tumor samples. Deconvolution strategies to analyze mixtures comprising one normal and one homogenous tumor component have been proposed by others [10, 13].

We tested PennCNV [11], genoCNA [10], the M-measure and the three-state Viterbi algorithm based on the M-measure (3SMM, see Methods) using the mixture scenario of a stromal contamination with a "pure" tumor sample. We used default parameters for both PennCNV and genoCNA. As a measure of performance, we calculated the average of the areas under the receiver-operator characteristic curve (AUCROC) at the operative points (Eq. 7), reducing the classification task to the correct classification of gain, loss and normal states (Figure 2B). As expected, PennCNV performed the worst. PennCNV is not designed for tumor copy number inference. Its poor performance indicates that algorithms that do not take into account the possibility that the data could be a mixture of more than one component will perform worse than algorithms that do not ignore this possibility. Given the low performance, PennCNV was not included in any further analysis. Overall, the M-measure performed as well as, or better than, the state-of-the-art HMM-based methods. Remarkably, it exhibited higher robustness to the mixing coefficient. Interestingly, 3SMM shows no evident improvement compared to the M-measure (Figure 2A, blue and green AUCROCs). We also compared the actual value of the mixing coefficient to the value inferred by solving the linear system, as described in the Methods, and to the value returned by genoCNA, as described in the genoCNA documentation. As shown in Figure 2B, our solution obtained using linear algebra is close to the actual solution (green and blue curves), even when the normal component was overwhelmingly abundant. The main effect of increasing stromal contamination seems to be reflected in a slightly larger variance of the estimated mixing coefficient, which is particularly visible in the inference based upon the 3SMM algorithm. This increase in the variance reflects the overall decrease in performance in terms of AUCROC (Figure 2A green and blue curves). The surprisingly weak performance of genoCNA in inferring the correct mixing coefficient is probably linked to the non-uniqueness of solutions in the problem of inferring exact copy numbers. As the use of Viterbi decoding did not significantly improve algorithmic performance, we did not use 3SMM for further analyses.

To test classification performance in more heterogeneous conditions, we generated 200 independent datasets of a tumor sample with two subclones contaminated by a stromal component for selected values of the normal component mixing coefficient α, as described above. In this scenario, the M-measure clearly shows its superior robustness in the inference of CNA status from heterogeneous samples with an AUCROC above 0.9 for mixtures in which the tumor components together constitute at least half of the sample (Figure 3). The classification performance of the M-measure decreases as the total fraction of stromal cells increased, suggesting that the decreasing signal-to-noise ratio of allelic imbalance, as already shown in Figure 1, is the main cause for misclassification. Further comparison using real SNP-array tumor data revealed higher robustness of the M-Measure in identifying de novo aberrations during tumor progression (Additional file 1, Additional file 1, Figure.S1, Additional file 1, Figure.S2 and Additional file 1, Table S1).

In general, the linear relationship between probe signal and allelic content shown in Equation 1 recapitulates that pooling and counting the number of occurrences of each allele from a biological sample follows basic algebraic relationships. This is the case in karyotyping, FISH experiments, exome sequencing, DNA-seq and many other experimental procedures. Equation 1 does not refer to a specific technology; instead, it reflects basic statistics notions and is, in its theoretical formulation, exact. It has been shown that SNP-array data may require preprocessing steps involving non-linear transformations [19]. We have shown that copy numbers cannot be inferred at the linear level and that a fortiori cannot be inferred after applying the non-linear transformations required to faithfully analyze SNP-array data.

Experimental validation of robustness of the allelic imbalance measure to heterogeneity

We sought to demonstrate that the allelic imbalance measure can detect CNA in highly heterogeneous tumor samples. For this purpose, we selected loci from melanoma tumor samples (Table 1 for details on the samples) that exhibited deviations relative to normal melanocytes in terms of CNA and mRNA expression. We identified 299 genes that were over-expressed in more than 24 melanoma samples (> 80% of the 30 samples in the cohort) relative to normal melanocytes. We selected the 20 most frequently over-expressed genes for further analyses by intersecting these genes with all aberrant genes whose M-measure and log-R were indicative of gain in copy number. Fluorescent in situ hybridization (FISH) analysis of four genes, E2F8, ETV4, EZH2 and FAM84B, which show the most recurring gains, revealed remarkably complex mixtures with varying gains and number of subclones (Figure 4). The prevalence of copy number gains of ETV4 and FAM84B suggests that aberrations involving these genes may be functionally involved in maintaining the disease, or that they occurred early in tumor genesis and inherited by most subclones during tumor expansion. This finding is consistent with melanoma cytogenetic reports of recurring early aberrations of their hosting chromosomes [3].

Sample_ID	Normal/Nevus/Melanoma	Stage	BRAF status	NRAS status
HFSC	Normal	Normal	NA	NA
Nbmel	Normal	Normal	NA	NA
YULOVY	Melanoma	I, primary	WT	Q61L
YUPLA	Melanoma	II	WT	WT
YUGOE	Melanoma	III	WT	WT
YUKIM	Melanoma	III	WT	Q61R
YUROL	Melanoma	III	WT	WT
YUPAO	Melanoma	III, acral	WT	WT
YUCAS	Melanoma	IV	WT	WT
YUCHER	Melanoma	IV	WT	Q61R
YUMAG	Melanoma	IV	WT	Q61R / WT
YUROB	Melanoma	IV	WT	WT
YUSIV	Melanoma	IV	WT	WT
YUTUR	Melanoma	IV	WT	WT
YUZOR	Melanoma	IV	WT	WT
YUWERA	Melanoma	IV, acral	WT	WT
YUHOIN	Melanoma	IV, primary	WT	WT
YUDOSO	Melanoma	llb, primary	WT	Q61K / WT
YUHEIK	Melanoma	primary	WT	WT
YUFULO	Melanoma	primary	WT	Q61L / WT
YUSTE	Melanoma	III	V600E	WT
YUCAL	Melanoma	IV	V600E	WT
YUSAC	Melanoma	IV	V600E	WT
YUGEN8	Melanoma	IV	V600E	WT
YUCLIR	Giant nevus	Giant nevus	V600E / WT	WT
YUSIK	Melanoma	III+	V600E / WT	WT
YUNIBO	Melanoma	IIb, primary	V600K	WT
YUKSI	Melanoma	IV	V600K	WT
YULAC	Melanoma	IV	V600K	WT
YUMAC	Melanoma	IV	V600K	WT
YURIF	Melanoma	IV	V600K	WT
YUSIT	Melanoma	IV	V600K / WT	WT

Samples have been characterized by disease status, stage and BRAF or NRAS mutation.

High-resolution profiling of heterogeneity using next-generation sequencing

Allelic imbalance can also be used to detect heterogeneity using next-generation sequencing, in particular genomic DNA sequencing, mRNA sequencing (RNA-seq) and whole-exome sequencing. Similar to the CNA analysis, allelic imbalance at the single nucleotide level can give insights into tumor heterogeneity, provided that there is sufficient sequencing coverage to reliably identify variations distinctive of subclonal populations. A recent study produced an in-depth sequencing profile of melanoma samples [15, 16]. We re-analyzed sequencing data from the MeWo melanoma cell line, which had the highest sequencing depth, to identify patterns associated with subclonal heterogeneity. Inferring heterogeneity from RNA-Seq is very difficult due to the co-occurrence of numerous effects, such as allele-specific expression, varying expression levels across subclones, and alternative splicing, all of which can alter the allelic balance. To overcome the presence of confounding effects, we focused on the variation of single nucleotide allelic imbalance along exons. We argue that in the absence of subclonal heterogeneity and of CNA at the level of DNA, allelic imbalance is expected to be constant along the entire exon. On the other hand, in the presence of subclonal heterogeneity, acquired variations in some subclones will introduce allelic imbalance at a population level with significant fluctuations between neighboring loci (Figure 5A). This imbalance is not associated with differential expression between paternal and maternal alleles, but rather with the co-occurrence of subclonal variants whose observed expression level may reflect their relative abundance. In addition, large fluctuations between contiguous nucleotides can exclude the presence of focal CNA, which would instead result in highly correlated fluctuations at each nucleotide along the aberration. Given the lack of independent CNA information for the MeWo sample, we studied heterogeneity in RNA-seq using single-nucleotide mutations. However, the proposed approach is general enough to be extended to the study of CNA heterogeneity in next-generation sequencing.