SHARP.Clustering.LIHC.Rmd

---
title: "SHARP: TCGA Liver Cancer"
author: "ks"
date: "`r Sys.Date()`"
output: html_document
---

# {.tabset}

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
library(corrplot)
library(mogsa)
library(survival)
```

## TCGA liver cancer

This lab is going to compare several joint dimension reduction approaches using multiomic data. This time we're going to analyze liver cancer data from the Cancer Genome Atlas (TCGA). I grabbed data from this website http://acgt.cs.tau.ac.il/multi_omic_benchmark/download.html referenced  as the benchmark data used in https://pubmed.ncbi.nlm.nih.gov/33402734/ and https://pubmed.ncbi.nlm.nih.gov/33402734/. 

The molecular data are gene expression, DNA methylation, and micro RNA expression. I subset the samples to the 410 samples that were analyzed on each platform.  

The expression data were a count matrix from RNA sequencing and there was considerable variation in total counts (library size) across the samples. I normalized the counts by sample (library) using counts per million and removed genes that were not expressed in at least 10% of the samples.  Then I log2-transformed the data to stabilize the variances of the genes. On the log2 scale, the variance of gene expression is approximately independent of the mean expression level.  This is not true on the count scale.

I filtered the provided DNA methylation data (5000 features) to remove target CpGs on the X,Y chromosome (n=158) and those that do not pass the MASK_quality filter from https://zwdzwd.github.io/InfiniumAnnotation (n=259). We are left with 4583 features for the analysis.

The microRNA data had similar library sizes so I applied the same filter I applied to the gene expression data, and log2 transformed those remaining (470 features).

Now we're ready to load the subset of processed LIHC data for multi-omic analysis.  We also have survival time for the patients with multi-omic analysis. Our goal will be to see if we can estimate latent variables that are associated with patient's survival time.
```{r load-data}
load("mbliver.rda")
sapply(mbliver,dim)
```

The  above shows the different number of features in each omic platform (row 1)  and the number of patients (row 2).

```{r survar}
str(survar)
```
There are some missing survival times, but not many.

Try the following analyses:

1. latent variables to predict survival (Liver data)
a. concatenated PCA   
b. consensus PCA (CPCA) 
c. regularized generalized CCA  (RGCCA)
d. multiple co-inertia analysis (MCIA)  

2. PCA of individual data blocks  
a. expression  
b. DNAm  
c. microRNA  

At the end we'll compare our discoveries across the 7 analyses.

## concatenated PCA 

I don't need the same number of variables in each block to run this, but let's arbitrarily filter on the top 400 varying features for each data type and explore changing the filter later. As it is, this file currently takes ~3 min. to run start-to-finish on my laptop.

```{r topvarfeatures}
# pick 400 most variable features
fmad  <- matrixStats::rowMads(mbliver$lcpm)
rfmad <- rank(-fmad)
fidE <- which(rfmad <= 400)

fmad  <- matrixStats::rowMads(mbliver$dnam)
rfmad <- rank(-fmad)
fidM <- which(rfmad <= 400)

fmad  <- matrixStats::rowMads(mbliver$lmirna)
rfmad <- rank(-fmad)
fidmir <- which(rfmad <= 400)

mbmat <- t(rbind(mbliver$lcpm[fidE,],
               mbliver$dnam[fidM,],
               mbliver$lmirna[fidmir,]))
dim(mbmat)
```

Now run PCA on the concatenated data set.

```{r pcall}
sdztall <- scale(mbmat)
concatpca <- prcomp(sdztall,retx=TRUE)
pcall <- as.data.frame(concatpca$x)

res.cox <- coxph(Surv(Survival, Death) ~ pcall$PC1, data = survar)
res.cox
```

Do any of the top 10 PCs predict death?
```{r coxphp}
pccvec <- NULL
for (i in 1:10)
pccvec[i] <- 1-pchisq(coxph(Surv(Survival, Death) ~ pcall[,i], 
      data = survar)$score,1)
nsign  <- sum(pccvec<0.005)
idx <- which(pccvec<0.005)
```

`r nsign` of the top 10 PCs predict death (bonferroni adj p < 0.05). It was the `r idx`th PC.

Do the samples cluster based on these PCs?
Est. # of clusters using the gap statistic.  I should increase nstart and B (number of replicates), but I want this to  run a little faster in lab.

```{r gapstatpcall, echo=FALSE}
set.seed(500)
gap_stat <- cluster::clusGap(pcall[,1:10], FUN = kmeans, nstart = 30, K.max = 12, B = 50)
plot(gap_stat,main=c("Gap Statistic"), xlab="Number of Clusters")
```

There must be a function that provides the number of clusters from this output but here are my calculations. I'm looking for the first time 'the gap' does not exceed the previous gap + 1 SE.
```{r}
ng <- function(gs) {
  nc <- nrow(gs$Tab)
  firstSEmax <- cbind.data.frame(nclust = 1:nc,
                    gap = gs$Tab[,"gap"],
                    SEmax = c(gs$Tab[1,"gap"],
                              rowSums(gs$Tab[-12,3:4]))
  )
  nclust <- min(which(firstSEmax$gap<firstSEmax$SEmax))
 return(nclust)
}
nClust <- ng(gap_stat)
```

This picks `r nClust` clusters.
```{r kmeansMultStartkm8}
set.seed(99)
km10 <- stats::kmeans(pcall[,1:10],centers  = nClust,nstart = 200)
coxph(Surv(Survival, Death) ~ factor(km10$cluster), data = survar)
```

```{r kmplot}
kmfit = survfit(Surv(Survival, Death) ~ km10$cluster, data = survar)
plot(kmfit, col = c(1:10), xlab = "Survival Time In Days", ylab = "Survival Probabilities")
```

Group 7 (yellow line) has the worst survival.
```{r}
table(km10$cluster)
```

We have `r table(km10$cluster)[7]` people in group 7.


```{r pcplot}
ggplot(data.frame(pcall),  aes(x=PC1, y=PC2,
              color = factor(km10$cluster))) + 
          geom_point() + labs(color="Cluster")
```

They are not identified by the first 2 PCs.

EXTRA CREDIT:  Which PC is most associated with group 7?
```{r}

```


## CPCA

Now let's combine the blocks of data using consensus PCA. The mbpca() function expects a list object with each data type as either a matrix or data.frame object, with genes on the rows and sample in the columns. I don't need the same number of variables in each block to run this, but we'll start with the top 400 varying features for each data type again and explore alternatives later.

```{r mblist}
mblist <- list(    E = mbliver$lcpm[fidE,],
                dnam = mbliver$dnam[fidM,],
               mirna = mbliver$lmirna[fidmir,])
sapply(mblist,dim)
```

Option = "lambda1"  will scale each block by it's 1st eigenvalue, to give equal weight to each data block. Type ?mbpca to see other "option" choices.

```{r mbpca}
moa <- mbpca(mblist, ncomp = 10, k = "all", method = "globalScore", option = "lambda1",
center=TRUE, scale=TRUE, moa = TRUE, svd.solver = "fast", maxiter = 1000)

```

```{r}
plot(moa, value="eig", type=2)
```

This next step takes a while to run and can be skipped until later. It estimates the coherence of different datasets, and estimates the number of components to include in the analysis.
```{r bootmbpcaE}
#set.seed(90)
#r <- bootMbpca(moa, mc.cores = 1, B=20, replace = FALSE, resample = "sample")
```

I'm going to use the top 10 PCs for all the analyses that follow. The results suggest that using this criterion may be more strict than we want for the purpose of discovering latent variables associated with survival.

The next function will fit our  same latent variable model with a sparsity constraint. k=0.3 says we'll estimate 30\% of our coefficients and set the others to 0.

```{r moas}
moas <- mbpca(mblist, ncomp = 2, k = 0.3, method = "globalScore", option = "lambda1",
center=TRUE, scale=TRUE, moa = TRUE, svd.solver = "fast", maxiter = 1000)
```

Let's get the resulting latent variables from the output and see how correlated the originals are to those estimated using the sparsity constraint.
```{r scorecor}
scr1 <- moaScore(moa)
scrs <- moaScore(moas)
diag(cor(scr1[, 1:2], scrs))
```

If we make the vectors too sparse, they don't estimate the original latent variables very well. We can play with above 2 code chunks, and try different k, but we seem to lose correlation if we fit too sparse of a model.  

Let's do a Cox regression analysis and see if any of our 10 latent variables predict survival time.  Do any of the top 10 common PCs predict death?
```{r coxphpcpca}
pcpca <- NULL
for (i in 1:10)
pcpca[i] <- 1-pchisq(coxph(Surv(Survival, Death) ~ scr1[,i], 
      data = survar)$score,1)
nsign  <- sum(pcpca<0.005)
idx <- which(pcpca<0.005)
```

`r nsign` of the top 10 PCs predict death (bonferroni adj p < 0.05). It was the `r idx`th PC.



Compare PC1 from concatenated data and consensus PCA
```{r pc-compare}
plot(scr1[,c("PC1")],pcall$PC1,pch=16)
```

They are similar.

Let's try clustering the top 10 latent variables from consensus PCA.

```{r gapstatcpca, echo=FALSE}
set.seed(500)
gap_stat <- cluster::clusGap(scr1[,1:10], FUN = kmeans, nstart = 30, K.max = 12, B = 50)
plot(gap_stat,main=c("Gap Statistic"), xlab="Number of Clusters")
```

```{r}
nClust <- ng(gap_stat)
```

This solution supports `r nClust` clusters.
```{r kmcpca}
set.seed(99)
km10 <- stats::kmeans(scr1[,1:10],centers = nClust,nstart = 200)
coxph(Surv(Survival, Death) ~ factor(km10$cluster), data = survar)
```

This says there is differential survival across our 11 groups (p=0.004).

What do the clusters look like in our latent variable space?
```{r  spcpca}
ggplot(data.frame(scr1),  aes(x=PC1, y=PC2,
                              color = factor(km10$cluster))) + 
                geom_point() + labs(color="Cluster")
```


## RGCAA

Now let's try using this same function to run generalized CCA  (>2 blocks). For this we change the method to "blockScore".
```{r mbrgcca}
moa <- mbpca(mblist, ncomp = 10, k = "all", method = "blockScore", option = "lambda1",
center=TRUE, scale=TRUE, moa = TRUE, svd.solver = "fast", maxiter = 1000)
```

```{r}
plot(moa, value="eig", type=2)
```

Now it looks like the DNA methylation data aren't contributing after the first latent variable.  Again, for the sake of time I'm going to skip this bootstrap procedure to pick the number of latent variables.

```{r}
#set.seed(90)
#r <- bootMbpca(moa, mc.cores = 1, B=20, replace = FALSE, resample = #"sample")
```

Let's estimate 3 latent variables with sparse loadings.

```{r moasrgcca}
moas <- mbpca(mblist, ncomp = 3, k = 0.2, method = "blockScore", option = "lambda1",
center=TRUE, scale=TRUE, moa = TRUE, svd.solver = "fast", maxiter = 1000)
```

```{r cor-rgcca-slvs}
scr2 <- moaScore(moa)
scrs <- moaScore(moas)
diag(cor(scr2[, 1:3], scrs))
```

Running the sparsity model saving 20\% of the features gives a result that is relatively correlated still with the original PCs.

Est. # of clusters using the gap statistic. 
```{r gapstatrgcca, echo=FALSE}
set.seed(500)
gap_stat <- cluster::clusGap(scr2[,1:10], FUN = kmeans, nstart = 30, K.max = 12, B = 50)
plot(gap_stat,main=c("Gap Statistic"), xlab="Number of Clusters")
```


```{r}
nClust <- ng(gap_stat)
```

This solution supports `r nClust` clusters.

```{r kmrgcca}
set.seed(99)
km10rgcaa <- stats::kmeans(scr2[,1:10],centers = nClust,nstart = 200)
coxph(Surv(Survival, Death) ~ factor(km10rgcaa$cluster), data = survar)
```

These groups also explain differences in survival times.

Let's look at the Survival curves.
```{r kmplotrgcca}
kmfit = survfit(Surv(Survival, Death) ~ km10rgcaa$cluster, data = survar)
plot(kmfit, col = c(1:10), xlab = "Survival Time In Days", ylab = "Survival Probabilities")
```


```{r}
ggplot(data.frame(scr2),  aes(x=PC1, y=PC2,
                              color = factor(km10rgcaa$cluster))) + 
                geom_point() + labs(color="Cluster")
```

Do any of the top 10 latent variables predict death?
```{r coxphprgcca}
prgcca <- NULL
for (i in 1:10)
prgcca[i] <- 1-pchisq(coxph(Surv(Survival, Death) ~ scr2[,i], 
      data = survar)$score,1)
nsign  <- sum(prgcca<0.005)
idx <- which(prgcca<0.005)
nsign
idx
```

Interesting. The patients show survival differences by group, but none of the individual latent variables are associated with survival.

## MCIA

Now I'll use this same function to run MCIA (>2 blocks). For this we change the method to "blockLoading".

```{r mcia}
moa <- mbpca(mblist, ncomp = 10, k = "all", method = "blockLoading", option = "lambda1",
center=TRUE, scale=TRUE, moa = TRUE, svd.solver = "fast", maxiter = 1000)
```

```{r}
plot(moa, value="eig", type=2)
```


```{r sparsemcia}
moas <- mbpca(mblist, ncomp = 3, k = 0.2, method = "blockLoading", option = "lambda1",
center=TRUE, scale=TRUE, moa = TRUE, svd.solver = "fast", maxiter = 1000)
```


```{r}
scr3 <- moaScore(moa)
scrs <- moaScore(moas)
diag(cor(scr3[, 1:3], scrs))
```

The sparse solutions don't explain these very well at all!


Est. # of clusters using the gap statistic.
```{r gapstatmcia, echo=FALSE}
set.seed(500)
gap_stat <- cluster::clusGap(scr3[,1:10], FUN = kmeans, nstart = 30, K.max = 12, B = 50)
plot(gap_stat,main=c("Gap Statistic"), xlab="Number of Clusters")
```


```{r}
nClust <- ng(gap_stat)
```

This solution supports `r nClust` clusters.

```{r kmmcia}
set.seed(99)
km10mcia <- stats::kmeans(scr3[,1:10],centers = nClust,nstart = 200)
coxph(Surv(Survival, Death) ~ factor(km10mcia$cluster), data = survar)
```

These groups also explain differences in survival times (p=0.02).

Let's look at the Survival curves.
```{r kmplotmcia}
kmfit = survfit(Surv(Survival, Death) ~ km10mcia$cluster, data = survar)
plot(kmfit, col = c(1:10), xlab = "Survival Time In Days", ylab = "Survival Probabilities")
```


```{r}
ggplot(data.frame(scr3),  aes(x=PC1, y=PC2,
                              color = factor(km10mcia$cluster))) + 
              geom_point() + labs(color="Cluster")
```

Do any of the top 10 latent variables predict death?
```{r coxphmcia}
pmcia <- NULL
for (i in 1:10)
pmcia[i] <- 1-pchisq(coxph(Surv(Survival, Death) ~ scr3[,i], 
      data = survar)$score,1)
nsign  <- sum(pmcia<0.005)
idx <- which(pmcia<0.005)
```

`r nsign` of the top 10 PCs predict death (bonferroni adj p < 0.05). It was PC:`r idx`.

## Compare Latent Vars

**concatenated PCA vs consensus PCA** 
Compare the first 3 latent variables from each method.

```{r cormat}
scrmat <- cbind(pcall[,1:3],scr1[,1:3])
corrplot::corrplot(cor(scrmat),is.corr=FALSE)
```

**CPCA vs. RGCCA**
```{r cormat2}
scrmat <- cbind(scr1[,1:3],scr2[,1:3])
corrplot::corrplot(cor(scrmat),is.corr=FALSE)
```

**CPCA vs. MCIA**
```{r cormat3}
scrmat <- cbind(scr1[,1:3],scr3[,1:3])
corrplot::corrplot(cor(scrmat),is.corr=FALSE)
```

## Gene Expression only: PCA

Select the 400 most variable features.

```{r pca}
# pick 400 most variable features
# subset and transpose the expression matrix
xE <- t(mbliver$lcpm[fidE,])
sdztE <- scale(xE)
#pca
my.pcaE <- prcomp(sdztE,retx=TRUE)
epcs <- as.data.frame(my.pcaE$x)

res.cox <- coxph(Surv(Survival, Death) ~ epcs$PC1, data = survar)
res.cox
```

EXTRA CREDIT: How many PCs correlate with survival time?  
```{r}

```

EXTRA CREDIT:  Repeat this without the filter of the top 400 varying genes. Now is PC1 a stronger or weaker predictor of survival time?
```{r}

```


## DNA methylation

Now let's see what happens when we do a similar analysis with DNA methylation data.

```{r pca-dnam}
# pick 400 most variable features
# transpose the dnam matrix
tM <- t(mbliver$dnam[fidM,])
sdztM <- scale(tM)
my.mpca <- prcomp(sdztM,retx=TRUE)
mpca <- as.data.frame(my.mpca$x)

res.cox <- coxph(Surv(Survival, Death) ~ mpca$PC1, data = survar)
res.cox
```

This one is not significant. 

EXTRA CREDIT: Do any of the top 10 PCs predict survival?
```{r}

```



## micro RNA

Let's try micro-RNA.

```{r pca-mirna}
# pick 400 most variable features
# transpose the dnam matrix
tM <- t(mbliver$lmirna[fidmir,])
sdztM <- scale(tM)
my.mirnapca <- prcomp(sdztM,retx=TRUE)
mirnapca <- as.data.frame(my.mirnapca$x)

res.cox <- coxph(Surv(Survival, Death) ~ mirnapca$PC1, data = survar)
res.cox
```

That one is highly statistically significant.

## Compare jDR vs DR

Let's create the plot of -log10 pvalue of the 10 PCs for each method, as we saw was done in the 2021 Nature Communications paper.

```{r  pmat}
results <- cbind(ccPCA = pccvec,
                  CPCA = pcpca,
                 RGCCA = prgcca,
                  MCIA = pmcia,
                   gex = rep(NA,10),
                  dnam = rep(NA,10),
                  mirna = rep(NA,10)
                 )

for (i in 1:10) {
  results[i,"gex"] <- 1-pchisq(coxph(Surv(Survival, Death) ~ epcs[,i], 
      data = survar)$score,1)
  results[i,"dnam"] <- 1-pchisq(coxph(Surv(Survival, Death) ~ mpca[,i], 
      data = survar)$score,1)
  results[i,"mirna"] <- 1-pchisq(coxph(Surv(Survival, Death) ~ mirnapca[,i], 
      data = survar)$score,1)
}
```

Apply a Bonferroni correction for performing 10 tests.

```{r}
adjp <- apply(results,2,p.adjust,"bonferroni")
```



```{r plot, echo=FALSE}
xm <- cbind.data.frame(x = rep(1:7,each=10),
            y = as.vector(unlist(-log10(adjp))))
sign <- xm[xm$y > -log10(0.05),]
plot(xm$x,xm$y,pch=16,col=1,
        axes=F,xlab="method",ylab="bonf adj p (-log10)")
axis(1,at=1:7,colnames(adjp))
axis(2)
box()
abline(-log10(0.05),0,lty=3)
points(sign,col=2,pch=16)
```

CONCLUSIONS:

1. Most methods found a signal that correlated with survival, either through the latent variables, a cluster analysis of the latent variables, or both.
2. The most significant signal came from PCA on the concatenated data blocks.
3. These results are not consistent with Suppl Fig 3. (2021) NatComm. Is it because I used a different R function? Those authors used 3 different packages: JIVE for CPCA, RGCCA for RGCCA, and omicade4 for MCIA. Or, because I chose to process the data differently?  I do think there are differences related to the data processing because even my PCA results on only the gene expression data (1 block) turned up a significant PC and theirs did not. 


## SessionInfo

```{r SessionInfo}
sessionInfo()
```