Abstract

Background

Tumour heterogeneity provides both challenges and opportunities in the development of cancer biomarkers. Tumours are mixed populations of multiple cell-types. Currently, the molecular profiles of interest – those of cancer cells or of distinct sub-populations of cancer cells – are blurred by the mixed signal from all cell types in a sample [1,2]. However, characterizing the heterogeneity of a patient’s tumour by identifying the sub-populations present, along with their proportions and molecular profiles, would provide a personalized cancer “fingerprint” that captures both cell-centred and whole-system information, opening up opportunities for targeted treatment [3,4]. The methods described in this article apply to mRNA expression data rather than to to DNA data (point mutations and copy number changes), which require other approaches [5-8].

As a first step, it is important to consider the two-population problem of normal and cancer cells. Even small fractions of contaminating normal cells can introduce noise in gene signatures [9,10], motivating the search for methods to deconvolve a mixed tumour profile by estimating the fraction of cancer cells and providing a personalized, purified mRNA abundance profile of the cancer cells.

Physical approaches for sample separation, such as laser capture micro-dissection [11] are costly, time-intensive, not always available and may degrade samples. Therefore, computational approaches to purification of tumour molecular profiles have become increasingly important. Table ​Table11 and Additional file 1 summarize the different methods currently available for deconvolving mRNA data.

Classical methods of profile deconvolution assume that a mixed profile is a linear combination of a predetermined number of pure constituent profiles. Written in matrix form, the measured, mixed profiles B are a product of A, a matrix of gene expression profiles of each constituent, weighted by the fractions X of each cell type in the mixture:

Equivalently, some algorithms start with the transpose of this equation. Different algorithms use different methods of deconvolution. Some assume that the fractions, X, are known [12-17]. Others use gene expression profiles [18], signatures [19-22] or markers [23-26] of the constituent profiles to recover X, and sometimes the expression profiles of other genes in the constituents [23]. (A gene marker is a set of genes assumed to be expressed solely in one cell subtype and in no other [4,27,28].) These approaches are limited in uncovering patient-specific variation in cancer, as they assume that all tumour profiles are mixtures of a small number of the same constituents. In addition, the expression data may be log-transformed [13,16], leading to a possible bias in the reconstruction of mixed tissue samples from constituent profiles [29].

Another approach to estimate both fraction of cancer content as well as patient-specific cancer profiles requires a matched normal profile for each patient [30,31]. In this case the normal profile is “electronically subtracted” from the bulk tumour profile. However, a matched normal profile may not always be available in existing datasets and may be difficult to obtain clinically. Furthermore, due to the biological variability of normal tissue, the provided normal profile may not match that of the tissue in the tumour sample.

Two algorithms, DeMix [32] and ISOpure [33] present statistical approaches for deconvolution of mixed tumour profiles given a set of unmatched normal samples. Cancer biomarkers generated from prostate and non-small cell lung cancer data purified using ISOpure were more effective at predicting survival relative to those generated using unpurified profiles [33].

Overview of ISOpure

In the following, we will use the term ‘ISOpure’ to refer to the algorithm in general, and ISOpureR to refer to the R package implementation. The next two sections will describe the statistical model and the algorithm in more detail, but we begin by providing a brief overview and example.

For the applications discussed here, the ISOpure algorithm is applied to microarray mRNA abundance data. The inputs to the model must be normalized (but not log-transformed) expression profiles. The following two sets of inputs are required: tumour mRNA abundance profiles, of the same cancer sub-type; and normal (i.e. healthy) mRNA abundance profiles, from the same tissue as the tumour.

The algorithm runs in two steps.

1. Cancer Profile Estimation (CPE). This step estimates and outputs an average cancer profile as well as a fraction of cancer in each tumour sample.

2. Patient Profile Estimation (PPE). This step estimates and outputs a cancer profile for each patient.These profiles are all similar to the average cancer profile, but contain patient-specific variations. The estimated cancer fraction for each tumour sample is fixed at the value calculated in the CPE step.

To run these two steps using ISOpureR, it suffices to apply the two functions ISOpure.step1.CPE and ISOpure.step2.PPE. The input expression data should be in matrix form, with samples along the columns and transcripts/features along the rows.

Open in a separate window

The vector ISOpureS1model$alphapurities produced by the first step contains the proportion of cancer for each patient. The matrix ISOpureS2 model $cc_cancerprofiles produced by the second step contains the patient-specific profiles, rescaled to be of the same scale as the tumour expression data. That is, these profiles are estimated as probabilities within the algorithm, but are scaled to represent microarray signal intensity data. A detailed example is given in Section 3 of the ISOpureR package vignette, included as Additional file 2.

In the statistical and algorithm descriptions of ISOpure, the notation t _n describes the matrix of tumour profiles, where the index n denotes the n-th patient. The fractions α _n represent the cancer fractions. The b _r and c _n represent normal and purified cancer profiles, but they will be interpreted as probabilities, as described in more detail below. An overview of the algorithm workflow is given in Figure ​Figure11.

Open in a separate window

Figure 1

ISOpure workflow. An overview of the ISOpure algorithm, illustrating the inputs and the most important outputs of the Cancer Profile Estimation (CPE) and the Patient Profile Estimation (PPE) steps. The CPE step estimates an average cancer profile over all patients and the proportion of cancer in each tumour, as well as the patient healthy profiles (as weights of input profiles). The healthy profile weights are re-estimated in the PPE step. This second step estimates the purified cancer profile for each patient. All estimated parameters from the CPE and PPE steps are output by the ISOpureR functions ISOpure.step1.CPE and ISOpure.step2.PPE.

While we focus on microarray mRNA abundance data, ISOpure is generic when it comes to different species of RNA. For example, it is able to deconvolve both mRNA (non-coding RNAs inclusive) and microRNAs. This is true for both microarray as well as next-generation sequencing (RNA-Seq) data. Input data matrices can represent genes, isoforms, exons or microRNAs of samples. The data matrices should represent approximate counts of molecules (normalised but not log converted) for a given RNA profile, e.g. gene level mRNA. However, the input should not contain a mixture of mRNA and microRNA data. ISOpure has been applied to the RNA-Seq The Cancer Genome Atlas (TCGA) PRAD dataset as part of that project (manuscript in preparation), helping to demonstrate its wider applicability.

Implementation

Results and discussion

To verify the numerical equivalence of the MATLAB and R implementations of ISOpure, their performances on four datasets were compared. Two sets were of lung adenocarcinoma from Bhattacharjee [39] and Beer [40] and two were of prostate cancer from Wallace [41] and Wang [42] (Table ​(Table2).2). These four datasets are among the tumour datasets purified using ISOpure in [33], chosen because the cancer types are not yet known to have established subtypes. The array data processing was detailed in [33].

Each dataset was purified using both the MATLAB and R implementations of ISOpure and we compared the resulting parameters. The algorithms were run 50 times for each dataset, with different initial conditions. To minimize differences due to random number generation implementations, the initial values of parameters were loaded from a file, and the extra optimizations of parameters ν, ω, and k [33] in the CPE step, which included some random initializations, were omitted. The results of these models are very similar to results generated by the full, randomized version of ISOpureR; the motivation for minimizing randomness was simply to reduce differences in the performance comparisons between MATLAB and R.

The numerical differences in the parameter estimates produced by MATLAB and R are small enough to have no biological significance. We calculated the means of the parameters (${\overline{\mathbf{x}}}_{\mathit{\text{MATLAB}}}$, ${\overline{\mathbf{x}}}_R$, for each of the estimated parameters such as ν,α, etc.) from the MATLAB and R models over the 50 iterations. A comparison of the entries of these mean parameter values produced by R and MATLAB shows that the two implementations are numerically analogous (Figure ​(Figure22 and Additional file 6). In particular, β ₁ and Spearman’s ρ are 1 for both the estimated fraction of cancer cells in the tumour, α, and for the log of the individual cancer profiles c _n.

Open in a separate window

Figure 2

A comparison of parameters estimated by the MATLAB and R implementations of ISOpure for the Beer dataset. Each plot shows the entries of a parameter estimated using ISOpureR plotted against the corresponding entries estimated using the MATLAB code. The parameter is an average over 49 models run with different initial conditions (one MATLAB model for the Beer dataset resulted in a zero θ value, and was dropped). The line y=x is indicated in black, and the linear regression line, or robust regression line for θ, is dashed orange. (A) Parameters from the Cancer Profile Estimation step of ISOpure: (i) ν, the hyper-parameter for the Dirichlet distribution over θ, (ii) θ, the proportion of a patient sample from a known healthy-tissue profile, (iii) m, the average mRNA abundance cancer profile, (iv) α, the fraction of cancer cells for every patient sample, (v) ω a hyper-parameter for the Dirichlet distribution over m. (B) Parameters from the Patient Profile Estimation step of ISOpure: (i) ν, the hyper-parameter for the Dirichlet distribution over θ, (ii) θ, the proportion of a patient sample from a known healthy-tissue profile, (iii) c _n, the purified mRNA abundance cancer profile for each patient.

For the all datasets, the mean and median of the fractional difference between the entries of the parameter means are very close to zero for six of the nine parameters (Additional file 7). For instance, for both the Beer and the Bhattacharjee datasets the means and medians for the fractional difference of α and logc _n are between 10⁻⁵ and 10⁻⁴. The three parameters having larger differences, ω in the CPE step and θ in both CPE and PPE steps, contain very small entries which are not biologically significant. The entries of θ represent the weights of known normal-tissue profiles adding up to a patient’s particular normal profile; a weight of 10⁻⁵ essentially means that particular known normal profile is not present for that patient. When entries larger than a certain threshold, such as 10⁻⁵, are compared, the fractional differencedecreases.

Furthermore, the differences between the MATLAB and R algorithms are similar to the differences within each implementation, under different initial conditions (Additional file 8). In particular, smaller numbers are not as precisely predicted. Computationally, operations with smaller numbers are susceptible to floating point rounding errors; however differences in small numbers do not alter the model, as they are biologically unimportant.

Finally, the mean and median of the vector components of $\left({\overline{\mathbf{x}}}_{\mathit{\text{MATLAB}}}-{\overline{\mathbf{x}}}_R\right)$ are very close to 0, and are sometimes positive and sometimes negative. The R algorithm does seem to under-estimate parameters compared to MATLAB, but this bias disappears when we compare parameter values larger than 10⁻⁵.

Conclusions

The first stage in the development of ISOpureR focused on establishing the numerical equivalence of the results produced by the R and MATLAB code. Future steps include testing backward compatibility with ISOLATE [43], the precursor to ISOpure.

The translation of MATLAB code into R code was surprisingly challenging. Debugging took five times as long as the translation and initial testing, as differences in function performance appeared only for certain input values during the execution of the full algorithm rather than in the testing of individual functions. Perhaps some of the differences in basic MATLAB and R operators mentioned can be of help for others tackling a similar project. A list of key issues encountered is in Additional file 9.

One of the recommendations for the comparison of implementations would be to eliminate all sources of randomness in the algorithms and postpone the improvement of running time only once numerical results are consistent, and consistently reproducible given different initial random seeds.

Most importantly, the contribution of the R implementation of the ISOpure algorithm is that it can now be readily integrated into existing analysis pipelines. ISOpureR can easily be included in benchmark comparisons of different deconvolution algorithms. A promising upcoming project is the ICGC-TCGA DREAM Somatic Mutation Calling - Tumour Heterogeneity Challenge [44], a benchmarking competition of computational methods for determining the best sub-clonal reconstruction algorithms. The collaborative-competitive framework of the DREAM projects encourages transparency in algorithm development and the availability of open-source code. The R implementation of ISOpure provides increased flexibility and ease of parallelization, so that the algorithm may be easily modified, extended and tested by thecommunity.

Availability and Requirements

The ISOpureR package is submitted to the Comprehensive R Archive Network (CRAN) which maintains an active package homepage for ISOpureR. The version of the code at the time of publication is included in Additional file 10. The package is written an implemented in the R programming language (version ≥ 3.1.1), with some C++ code incorporated using Rcpp [45] and RcppEigen [38]. It is platform-independent, but has been primarily tested on Linux, and is available under the GLP-2 license.

Acknowledgements

Footnotes

Contributor Information

Catalina V Anghel, Email: [email protected].

Gerald Quon, Email: ude.tim@nouqg.

Syed Haider, Email: [email protected].

Francis Nguyen, Email: moc.liamtoh@neyugN-F.

Amit G Deshwar, Email: [email protected].

Quaid D Morris, Email: [email protected].

Paul C Boutros, Email: [email protected].

References