-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathIntLIM-BreastCancerAmbsVignette.Rmd
198 lines (139 loc) · 11 KB
/
IntLIM-BreastCancerAmbsVignette.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
---
title: "Running the linear model on breast cancer study"
author: "Jalal K. Siddiqui"
date: "3/29/2018"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
##Introduction
A previous study had conducted both gene expression and metabolomics profiling of tissue samples from breast cancer patients (1). This vigenette will highlight the analysis we conduct on the breast cancer data. More details can be found in <a href="https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2085-6" target="_blank"> our publication "IntLIM: integration using linear models of metabolomics and gene expression data"</a> (2).
##Loading in IntLIM and Files
IntLIM, available from Github, can be installed as in the documentation. Once IntLIM is installed, what is necessary is loading in the library. First clear the workspace.
```{r}
rm(list = ls())
library(IntLIM)
```
For breast cancer study, both gene expression (available on Gene Expression Omnibus Accession Number: GSE37751) and metabolomics(1) (http://www.jci.org/articles/view/71180/sd/2) data are available online. Much of this data had been processed as previously described (1). Probes from the Affymetrix data not mapping to a gene symbol were removed. Additionally, as with NCI-60 data, only the probe corresponding to the highest mean expression was used for analysis when multiple probes corresponded to a single gene. This resulted in a total of 20,254 genes for 108 patient samples. The Metabolon data did not need to be filtered by coefficient of variation, as there were no technical replicates. The resulting data consisted of 536 metabolites with 132 patient samples.
This data has been formatted for IntLIM. We load in the data as follows. The bc.ambs.csv meta file contains a list of phenotypic data file, metabolite data file, gene expression data file, metabolite meta file, and gene expression meta file. The function OutputStats will give a summary of the NCI-60 data
Note: The data is found in the BC.data.vignette.zip file. This file has to be de-zipped prior to loading.
Load in the breast cancer data using the ReadData() function. ShowStats allow us to
```{r}
inputData <- IntLIM::ReadData('input.csv',metabid='id',geneid='id')
IntLIM::ShowStats(inputData)
```
From the OutputStats, we find that we have gene expression data involving 20,254 genes with 108 patient samples and metabolite abundance data involving 536 metabolites with 132 patient samples.
##Filtering Gene Expression and Metabolite Data
The __FilterData()__ function is used to filter the data. We remove genes with mean belows below the 10th percentile. Furthermore, we remove metabolites with more than 80% missing values. This results in gene expression data involving 18,228 genes and 108 patient samples and metabolite abundance data involving 379 metabolites and 132 patient samples.
```{r}
inputDatafilt <- IntLIM::FilterData(inputData,geneperc=0.10, metabmiss = 0.80)
IntLIM::ShowStats(inputDatafilt)
```
We can obtain boxplot distributions of the data as follows. This is used to make figures.
```{r}
IntLIM::PlotDistributions(inputDatafilt, palette = c("black", "black"))
```
##Principal Component Analysis
The principal component analysis is performed on filtered metabolite and gene expression data to obtain visual representations showing how different sub-sections of the data could be grouped into different clusters. Common samples patient samples (either tumor or adjacent non-tumor samples). Blue samples indicate tumor samples and red samples indicate non-tumor samples. Note the clear delineation of samples.
```{r}
PlotPCA(inputDatafilt, stype = "DIAG", common = F)
```
##Running Linear Model
The linear model is for integrating transcriptomic and metabolomics data is: E(m|g,t) = β1 + β2 g + β3 p + β4 (g:p) + ε where ‘m’ and ‘g’ are log-transformed metabolite abundances and gene levels respectively, ‘p’ is phenotype (cancer type, patient diagnosis, treatment group, etc), ‘(g:p)’ is the association between gene expression and phenotype, and ‘ε’ is the error term that is normally distributed. A statistically significant p-value of the ‘(g:p)’ association term indicates that the slope relating gene expression and metabolite abundance is different from one phenotype compared to another. We run a linear model on tumor (n = 61) and non-tumor samples (n = 47) that included 18,228 genes and 379 metabolites (total of 6,908,412 possible associations and hence models). For genes and metabolites that had standard deviations of 0 in one of the groups, we assign a p-value of NA. The model is run as below by calling __RunIntLim()__. __DistPvalues()__ allows us to obtain a distribution of p-values for the (g:p) term. The __pvalCorrVolcano()__ function allows us to observe how p-values vary with correlation differences.
```{r}
myres <- IntLIM::RunIntLim(inputDatafilt, stype="DIAG")
IntLIM::DistPvalues(myres)
IntLIM::pvalCorrVolcano(myres, inputDatafilt, diffcorr = 0.5, pvalcutoff = 0.05)
```
The next step is to process the results of this model by filtering the results of the linear model by FDR-adjusted p-value cutoff (0.10 selected here) for the (g:p) association coefficient and calculate the correlations of the gene-metabolite pairs in each group (tumor and non-tumor) for the filtered results. We further may only interested in results that have an absolute correlation value difference (0.50 selected here). This is done with the __ProcessResults()__ function. In addition we also develop a heatmap of the gene-metabolite association correlations for the selected groups.
```{r}
myres10 <- IntLIM::ProcessResults(myres, inputDatafilt, diffcorr = 0.50, pvalcutoff = 0.10)
IntLIM::CorrHeatmap(myres10)
OutputResults(myres10, filename = "bcresults10.csv")
dim([email protected])
```
We find that we obtain 10,861 correlations. We try to lessen the number by setting a more stringent p-value cutoff of 0.05 as done below.
```{r}
myres <- IntLIM::ProcessResults(myres, inputDatafilt, diffcorr = 0.50, pvalcutoff = 0.05)
#colnames([email protected])[3] <- "NORMAL"
#colnames([email protected])[4] <- "TUMOR"
IntLIM::CorrHeatmap(myres, top_pairs = 3000)
OutputResults(myres, filename = "bcresults5.csv")
dim([email protected])
```
From this model we find 2842 gene-metabolite correlations that have an association FDR-adjusted p-value of 0.05 and an absolute value correlation difference of 0.5 or greater. The top pairs are shown below.
```{r}
corr.table <- [email protected]
abs.corrdiff <- abs([email protected]$NORMAL_cor - [email protected]$TUMOR_cor)
sort.table <- corr.table[order(-abs.corrdiff),]
sort.table[1:20,]
```
We can show some example plots of some of these pairs. The first example is the PDLIM4 vs. adrenate (22:4n6). There appears to be an error with the scatterplot.
```{r}
IntLIM::PlotGMPair(inputDatafilt, stype = "DIAG", geneName = "PDLIM4", metabName = "adrenate (22:4n6)")
```
We do find that there is one pair with 2-hydroxyglutarate at FDR adjusted p-value of 0.05 and correlation difference of over 0.5 with GPT2.
```{r}
IntLIM::PlotGMPair(inputDatafilt, stype = "DIAG", geneName = "GPT2", metabName = "2-hydroxyglutarate")
```
GPT2 and MYC are not linked as was in the Ambs paper.
```{r}
IntLIM::PlotGMPair(inputDatafilt, stype = "DIAG", geneName = "MYC", metabName = "2-hydroxyglutarate")
```
##Analyzing Clusters
We will cut the heatmap as below. We find 1038 gene-metabolite pairs in the Tumor-correlated cluster and 1804 pairs in the Tumor anti-correlated cluster. The summary of the cluster.1 shows that the gene-metabolite correlations are positive for the tumor samples making this into the tumor-correlated cluster. For cluster.2, the gene-metabolite correlations are negative for tumor samples making it the tumor anti-correlated cluster.
```{r}
hc.rows<- hclust(dist([email protected][,c(3,4)]))
ct<- cutree(hc.rows, k=2)
cluster.1 <- [email protected][which(ct == 1), ]
cluster.2 <- [email protected][which(ct == 2), ]
dim(cluster.1)
dim(cluster.2)
summary(cluster.1$TUMOR_cor)
summary(cluster.2$TUMOR_cor)
```
We have 288 unique genes in the Tumor-correlated cluster and 479 unique genes in the Tumor-anti-correlated cluster.
```{r}
tumor.corr.uniqgene <- unique(cluster.1$gene)
tumor.anti.corr.uniqgene <- unique(cluster.2$gene)
write.csv(tumor.corr.uniqgene, "tumor.corr.uniqgene.bc.csv")
write.csv(tumor.anti.corr.uniqgene, "tumor.anti.corr.uniqgene.bc.csv")
length(tumor.corr.uniqgene)
length(tumor.anti.corr.uniqgene)
```
We also find out how many unique metabolites we have. We have 155 in the tumor-correlated cluster and 188 in the tumor-anti-correlated cluster.
```{r}
tumor.corr.uniqmetab <- unique(cluster.1$metab)
tumor.anti.corr.uniqmetab <- unique(cluster.2$metab)
length(tumor.corr.uniqmetab)
length(tumor.anti.corr.uniqmetab)
```
We need to fit the metabolite data to Human Metabolome Database (HMDB) IDs
The Metabolon fData file is imported as follows. This is used to convert the list of unique metabolites into HMDB IDs for the IPA analysis.
```{r}
fData.metab <- read.csv("fData.metab.csv")
hmdb.match <- fData.metab[,c('id', 'HMDB_ID')]
rownames(hmdb.match) <- hmdb.match$id
write.csv(hmdb.match, "hmdb.match.csv")
tumor.corr.uniqmetab.hmdbid <- hmdb.match[as.character(tumor.corr.uniqmetab),]
tumor.corr.hmdb.list <- tumor.corr.uniqmetab.hmdbid[is.na(tumor.corr.uniqmetab.hmdbid$HMDB_ID) == FALSE,'HMDB_ID']
tumor.anti.corr.uniqmetab.hmdbid <- hmdb.match[as.character(tumor.anti.corr.uniqmetab),]
tumor.anti.corr.hmdb.list <- tumor.anti.corr.uniqmetab.hmdbid[is.na(tumor.anti.corr.uniqmetab.hmdbid$HMDB_ID) == FALSE,'HMDB_ID']
write.csv(tumor.corr.hmdb.list, "tumor.corr.uniqmetab.hmdbid.bc.csv")
write.csv(tumor.anti.corr.hmdb.list, "tumor.anti.corr.uniqmetab.hmdbid.bc.csv")
```
The following 5 genes are found to be in both clusters.
```{r}
intersect(tumor.corr.uniqgene, tumor.anti.corr.uniqgene)
length(intersect(tumor.corr.uniqgene, tumor.anti.corr.uniqgene))
```
These metabolites are found in both clusters
```{r}
intersect(tumor.corr.uniqmetab, tumor.anti.corr.uniqmetab)
length(intersect(tumor.corr.uniqmetab, tumor.anti.corr.uniqmetab))
```
The list of unique genes and metabolites for each cluster are used to conduct a pathway enrichment analysis using Ingenuity Pathway Analysis (IPA) (https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/).
###Reference
1. Terunuma A, Putluri N, Mishra1 P, Mathé EA, Dorsey TH, Yi M, Wallace TA, Issaq HJ, Zhou M, Killian JK, Stevenson HS, Karoly ED, Chan K, Samanta S, Prieto D, Hsu TY.T., Kurley SJ, Putluri V, Sonavane R, Edelman DC, Wulff J, Starks AM, Yang Y, Kittles RA, Yfantis HG, Lee DH, Ioffe OB, Schiff R, Stephens RM, Meltzer PS, Veenstra TD, Westbrook TF, Sreekumar A, and Stefan Ambs S. MYC-driven 2-hydroxyglutarate associates with poor prognosis in breast cancer. J Clin Invest. 2014 Jan 2;124(1):398-412.
2. Siddiqui JK, Baskin E, Liu M, Cantemir-Stone CZ, Zhang B, Bonneville R, McElroy JP, Coombes KR, Mathé EA. IntLIM: integration using linear models of metabolomics and gene expression data. BMC bioinformatics. 2018 Dec;19(1):81.