Multiple Omics Data Integration to Identify Long Noncoding RNA Responsible for Breast Cancer–Related Mortality

Abstract

Long non-coding RNAs (lncRNAs) are a large and diverse class of transcribed RNAs, which have been shown to play a significant role in developing cancer. In this study, we apply integrative modeling framework to integrate the DNA copy number variation (CNV), lncRNA expression, and downstream target protein expression to predict patient survival in breast cancer. We develop a 3-stage model combining a mechanical model (lncRNA regressed on CNV and target proteins regressed on lncRNA) and a clinical model (survival regressed on estimated effects from the mechanical models). Using lncRNAs (such as HOTAIR and MALAT1) along with their CNV, target protein expressions, and survival outcomes from The Cancer Genome Atlas (TCGA) database, we show that predicted mean square error and integrated Brier score (IBS) are both lower for the proposed 3-step integrated model than that of 2-step model. Therefore, the integrative model has better predictive ability than the 2-step model not considering target protein information.

Introduction

Several evidences highlight the emerging impact of long noncoding RNAs (lncRNAs) in cancer progression.^1-4 The aim of this study is to identify the predictive capability of some oncogenic lncRNAs in tumor progression and prognosis of breast cancer.

Breast cancer is the most common malignancy and the leading cause of cancer death in women. By focusing on a single type of genetic alteration such as copy number variation (CNV), scientists have identified significant genes that may contribute to cancer progression.^5-8 Due to its complexity, the study of cancer should focus on incorporating data from multiple platforms ranging from genes, transcripts, and proteins found in cancer cells, ⁹ to whole biological systems, represented by molecular pathways and cell populations. ¹⁰ The integration, where multiple levels of omics data (ie, CNV, methylation, and gene expression) are gathered from the same subjects and analyzed, is known as vertical integration.^10-12

In this study, we introduce an easy and simplified way to integrate multiple omics data to show that the survival prediction due to the presence of lncRNAs increases significantly in breast cancer. We consider the genomic platform such as CNV, mRNA expression, proteomic platform such as protein expression, and the phenotype such as the survival of the patients. This study focuses only on the lncRNA expressions from The Cancer Genome Atlas (TCGA) breast cancer data. We consider the target protein expressions as proteomics data.

An Integrative Model

We consider a 3-stage model here. Suppose that n is the number of patients, p is the number of lncRNAs, and L is the number of CNV expressions.

The mechanistic model for each lncRNA can be expressed as

$$lnc\mathit{R}\mathit{N}{\mathit{A}}_k=\sum _{l=1}^L{c}_{kl}{\alpha }_l+{{\rm O}}_k,k=1,\dots ,p$$

(1)

where $lnc\mathit{R}\mathit{N}{\mathit{A}}_k$ is the level of gene expression for gene $k,k=1,\dots ,p$ , and is of dimension $n\times 1$ ; $c_{kl}$ is part of the $lncRN{A}_k$ expression that is attributed to the lth CNV; $O_k$ is the other (remaining) part of the gene expression which is not regulated by CNV and is of dimension $n\times 1$ ; and $\mathrm{\alpha }$ is the regression coefficient vector.

Next, the downstream target protein of each specific lncRNA was identified from PubMed articles, TCGA RNA-Seq database, and other extensive analyses such as differential expression analysis. The mechanistic model for each protein (for every lncRNA) can be expressed as

$$Protein=C{\gamma }_1+O{\gamma }_2+O^{*}$$

(2)

where $C={(c_{kl})}_{p\times L}$ and $O^{*}$ represents the “other” part of the protein expression that is not regulated by lncRNA. ${\gamma }_1$ and ${\gamma }_2$ are the regression coefficients corresponding to the CNV expressions and the error part from equation (1), respectively.

The clinical component part models the effect of the mechanistic parts of the genes on a clinical outcome of interest and can be written as

$$\log t=lnc\mathit{R}\mathit{N}{\mathit{A}}_1{\mathrm{\beta }}_1+{\mathit{O}}_2^{*}{\mathrm{\beta }}_2+ϵ$$

(3)

where $\mathrm{t}$ is the survival outcome, $ϵ$ is the error term, and ${\mathrm{\beta }}_1$ and ${\mathrm{\beta }}_2$ are the usual regression coefficients corresponding to lncRNA and the estimated error part from equation (2), respectively.

The variable $\mathit{P}\mathit{r}\mathit{o}\mathit{t}\mathit{e}\mathit{i}\mathit{n}$ represents the vectorized downstream gene effects attributed to protein expressions and is estimated from the second-stage mechanistic model. Therefore, the clinical component additively models the effects of all the gene expressions and their components—derived from different sources (gene expression, CNV) in a unified manner.

Assumptions such as $\mathrm{O}\sim N(0,{\mathrm{\sigma }}_1^2{I}_n)$ and ${\mathrm{O}}^{*}\sim N(0,{\mathrm{\sigma }}_2^2{I}_n)$ give rise to the usual linear model, whereas we obtain the log-normal accelerated failure time (AFT) model when we assume $ϵ\sim N(0,{\mathrm{\sigma }}_3^2{I}_n)$ .

In the presence of right censoring, we observe the tuple $(t_i^{*},{\delta }_i),i=1,\dots ,n$ , where ${\mathrm{\delta }}_i=1$ if the event is observed (death in this case), and 0 otherwise; $t_i^{*}=\min (t_i,c_i)$ with $c_i$ being the censoring time. A standard statistical software can be used to fit a log-normal AFT model and the other linear regression models.

To quantify the prediction accuracy, we consider a standard comparative predictive approach Brier score (BS) ¹³ which uses the predicted survival times

$$\begin{array}{c}\mathrm{BS}(t)=n^{-1}\sum _{i=1}^n\frac{\widehat{S}{(t|x_i)}^{2}I(t_i\le t\wedge {\delta }_i=1)}{\widehat{K}(t_i)}\\ +\frac{{(1-\widehat{S}(1|x_i))}^{2}I(t_i>t)}{\widehat{K}(t)}\end{array}$$

where $\widehat{K}(.)$ denotes the Kaplan-Meier estimate of the censoring distribution which is based on the observations $(t_i,1-{\delta }_i)$ , and $\widehat{S}(.)$ stands for the estimated survival function. As the mathematical form suggests, BS provides a numerical comparison between the observed and estimated survival functions. Brier score is defined for each time point $t$ and hence can be added for the entire time range to obtain IBS, $\mathrm{IBS}=\max {(t_i)}^{-1}{\int }_0^{\max (t_i)}\mathrm{BS}(t)dt$ . We can see that models with smaller scores are preferred. We compute integrated Brier score (IBS) using ipred package. ¹⁴

Nevertheless, we also compute the prediction square error by comparing the observed data and their posterior predicted values.

From TCGA database, we consider the information of 222 breast tumor samples with their survival data. We observe that at least 82% data are right censored.

Along with the clinical observations, we also collected measurements of 12 lncRNA expressions (Table 1). Among those, we found the CNV information available for 9 genes (or lncRNAs). We also consider 64 target protein expressions for these genes.

Table 1.

The lncRNA considered for our experiment.

Gene	Function
BCAR4 a	Oncogenic, promotes invasion and metastasis 15
BCYRN1	Oncogenic, promotes tumor progression 16
GAS5 a	Tumor suppressor 17
H19 a	Oncogenic, promotes proliferation and metastasis 18
HOTAIR a	Oncogenic, promotes EMT, proliferation, and metastasis 19
MALAT1 a	Oncogenic, promotes proliferation, invasion, and migration 20
MEG3 a	Tumor suppressor, induces accumulation of p53 21
PVT1 a	Oncogenic, promotes tumor progression 22
SOX2OT	Oncogenic, promotes tumor growth and metastasis 23
SRA1 a	Oncogenic 24
UCA1 a	Oncogenic, promotes cell growth, suppresses the tumor suppressor p27 25
XIST	Tumor suppressor 26

Abbreviation: EMT, epithelial-mesenchymal transition.

a

The copy number variation available (among those lncRNAs, SRA1 transcribes both long noncoding and protein-coding RNAs which are produced by alternative splicing).

We apply the integrative modeling in these data and obtain the results shown in Table 2. We notice that the mean squared prediction error and IBS are both lower for the proposed model than for the 2-stage model after omitting the protein expressions from the analysis.

Table 2.

MSPE and IBS for fitted models in TCGA breast cancer data.

Models	MSPE	IBS
2-stage	1.903	0.488
3-stage	1.196	0.395

Abbreviations: IBS, integrated Brier score; MSPE, mean squared prediction error; TCGA, The Cancer Genome Atlas.

In this article, we have shown that when the contribution of lncRNA’s target protein expression measurement is not ignored, then the survival prediction has improved dramatically. Toward this, we have developed a simple yet integrative modeling strategy which borrows strengths from all 3 platforms such as DNA CNV, mRNA expressions for the long noncoding genes, and their target protein expressions to predict the survival of the subjects. We have shown that this integrated model outperforms its closest competitor.