```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
#If you want to knit the whole notebook, please specify below correct path to the Rmd-file
knitr::opts_knit$set(root.dir=getwd())
```


## Exploratory Data Analysis (EDA)

We will start with reading and imputing missing values using median imputation:

```{r,fig.width=10,fig.height=8}
#Uncomment this in case you can not change your working directory
#current_dir<-getwd()
#setwd(paste0(current_dir,"/session_ml/SupervisedOMICsIntegration"))

expr <- as.data.frame(t(read.delim("CLL_mRNA.txt", header = TRUE, sep="\t")))
for(i in 1:ncol(expr)){expr[,i][is.na(expr[,i])]<-median(expr[,i],na.rm=TRUE)}
expr[1:5,1:5]

mut <- as.data.frame(t(read.delim("CLL_Mutations.txt", header = TRUE, sep="\t")))
for(i in 1:ncol(mut)){mut[,i][is.na(mut[,i])]<-median(mut[,i],na.rm=TRUE)}
mut[1:5,1:5]

meth <- as.data.frame(t(read.delim("CLL_Methylation.txt", header = TRUE, sep="\t")))
for(i in 1:ncol(meth)){meth[,i][is.na(meth[,i])]<-median(meth[,i],na.rm=TRUE)}
meth[1:5,1:5]

drug <- as.data.frame(t(read.delim("CLL_Drugs.txt", header = TRUE, sep="\t")))
for(i in 1:ncol(drug)){drug[,i][is.na(drug[,i])]<-median(drug[,i],na.rm=TRUE)}
drug[1:5,1:5]
```

Now let us have a look at the phenotypic data, extract Gender and convert it to a variable Y that is going to be used later when running PLS-DA analysis:

```{r,fig.width=10,fig.height=8}
phen <- read.delim("CLL_Covariates.txt", header = TRUE, sep="\t")
head(phen)
Y<-factor(phen$Gender)
summary(Y)
```

Let us split the data set into train and test sub-sets. We select 60 samples (30%) for testing. Since we have more males than females (fraction of females is approximately 0.4) we select 36 males and 24 females.

```{r,fig.width=10,fig.height=8}
set.seed(1234)
female_sample<-rownames(expr)[as.character(phen$Gender)=="f"][sample(1:length(rownames(expr)[as.character(phen$Gender)=="f"]),24)]
female_sample
male_sample<-rownames(expr)[as.character(phen$Gender)=="m"][sample(1:length(rownames(expr)[as.character(phen$Gender)=="m"]),36)]
male_sample

expr_test<-expr[match(c(female_sample,male_sample),rownames(expr)),]
expr_test[1:5,1:5]
meth_test<-meth[match(c(female_sample,male_sample),rownames(meth)),]
meth_test[1:5,1:5]
mut_test<-mut[match(c(female_sample,male_sample),rownames(mut)),]
mut_test[1:5,1:5]
drug_test<-drug[match(c(female_sample,male_sample),rownames(drug)),]
drug_test[1:5,1:5]

Y.test<-Y[match(c(female_sample,male_sample),rownames(phen))]
Y.test
summary(Y.test)
length(Y.test)

expr<-expr[!rownames(expr)%in%rownames(expr_test),]
meth<-meth[!rownames(meth)%in%rownames(meth_test),]
mut<-mut[!rownames(mut)%in%rownames(mut_test),]
drug<-drug[!rownames(drug)%in%rownames(drug_test),]

Y<-Y[!rownames(phen)%in%c(female_sample,male_sample)]
Y
length(Y)
summary(Y)
```

Since mutations represent a binary data, there is always a lack of variation due to coding with 0 and 1. Therefore, we will prefilter the mutation matrix by excluding sites with variance across individuals close to zero:

```{r filter methylation matrix}
library("mixOmics")
library("matrixStats")
my_nearZeroVar<-nearZeroVar(mut)
head(my_nearZeroVar$Metrics)
dim(my_nearZeroVar$Metrics)
mut <- mut[,-which(colnames(mut)%in%rownames(my_nearZeroVar$Metrics))]
mut[1:5,1:5]
dim(mut)
mut_test<-subset(mut_test,select=colnames(mut))
```

Therefore we end up with just a few mutations which have enough variation for the PLS-DA model. Later perhaps it makes sense to include all of them and not select informative ones using the sparse PLS-DA algorithm.


### Feature Selection for OMICs Integration

Right now we can not integrate the OMICs data sets since at least expression and methylation data are high-dimensional so we need to perform a feature selection for those two OMICs. Here we use LASSO for doing feature selection for gene expression data:

```{r,fig.width=10,fig.height=8}
library("glmnet")
lasso_fit <- cv.glmnet(as.matrix(expr), Y, family = "binomial", alpha = 1)
plot(lasso_fit)
coef <- predict(lasso_fit, s = "lambda.min", type = "nonzero")
colnames(expr)[unlist(coef)]
result_expr <- data.frame(GENE = names(as.matrix(coef(lasso_fit, s = "lambda.min"))
                                [as.matrix(coef(lasso_fit, s = "lambda.min"))[,1]!=0, 1])[-1], 
                   SCORE = as.numeric(as.matrix(coef(lasso_fit, s = "lambda.min"))
                                      [as.matrix(coef(lasso_fit, 
                                                      s = "lambda.min"))[,1]!=0, 1])[-1])
result_expr <- result_expr[order(-abs(result_expr$SCORE)),]
print(head(result_expr,10))
print(as.character(result_expr$GENE))
expr <- subset(expr, select = as.character(result_expr$GENE))
expr_test<-subset(expr_test,select=colnames(expr))

ens2gs<-read.delim("ENSEMBLE_TO_GENE_SYMBOL_COORD.txt",header=TRUE,sep="\t")
ens2gs<-ens2gs[match(colnames(expr),as.character(ens2gs$ensembl_gene_id)),]
colnames(expr)<-ens2gs$external_gene_name
colnames(expr_test)<-ens2gs$external_gene_name
```

We have also converted the Ensembl gene IDs into gene names on the fly for clarity. Next we proceed with LASSO feature selection for methylation data:

```{r,fig.width=10,fig.height=8}
library("glmnet")
lasso_fit <- cv.glmnet(as.matrix(meth), Y, family = "binomial", alpha = 1)
plot(lasso_fit)
coef <- predict(lasso_fit, s = "lambda.min", type = "nonzero")
colnames(meth)[unlist(coef)]
result_meth <- data.frame(CPG = names(as.matrix(coef(lasso_fit, s = "lambda.min"))
                                      [as.matrix(coef(lasso_fit, s = "lambda.min"))[,1]!=0, 1])[-1], 
                          SCORE = as.numeric(as.matrix(coef(lasso_fit, s = "lambda.min"))
                                             [as.matrix(coef(lasso_fit, s = "lambda.min"))[,1]!=0, 1])[-1])
result_meth <- result_meth[order(-abs(result_meth$SCORE)),]
print(head(result_meth,10))
print(as.character(result_meth$CPG))
meth <- subset(meth, select = as.character(result_meth$CPG))
meth_test <- subset(meth_test, select = as.character(result_meth$CPG))
```

We can see that we dramatically decreased the number of dimensions / features in the gene expression and methylation data sets. Now we can proceed with OMICs integration via PLS-DA algorithm.

### Integrative OMICs Analysis

Now we will start integrating the four OMICs: 1) gene expression, 2) methylation and 3) mutations and 4) drug response. For this purpose we will concatenate gene expression, methylation, mutation and drug response matrices into X matrix and use the Gender as Y variable, so it is a typical Machine Learning setup: y=f(x), where x is the input, y is the class labels of individuals and the f-function is learnt from the data. Note that the f-function is a-priori linear for PLS-DA and non-linear for e.g. artificial neural networks.

```{r,fig.width=10,fig.height=8}
data<-list(expr=expr,mut=mut,meth=meth,drug=drug)
names(data)
lapply(data, dim) 
```

First, we fit a DIABLO model without variable selection to assess the global performance and choose the number of components for the final DIABLO model. The function "perf" will perform M-fold cross validation (number of folds is specified by "folds" parameter) repeated "nrepeat" times. The design matrix sets expected correlation between the OMICs. The values may range between 0 (no correlation) to 1 (strong correlation), the design can be chosen based on prior knowledge inferred by e.g. MOFA (see Unsupervised OMICs integration session). Here due to the lack of prior knowledge we assume a strong correlation 1 between the OMICs.

```{r perf,fig.width=10,fig.height=8}
library("mixOmics")
design=matrix(1,ncol=length(data),nrow=length(data),dimnames=list(names(data),names(data)))
diag(design)=0
design

splsda.res = block.splsda(X = data, Y = Y, ncomp = 8, design = design)
perf.diablo = perf(splsda.res, validation = 'Mfold', folds = 2, nrepeat = 5, progressBar=FALSE, cpus=4)
perf.diablo
plot(perf.diablo,overlay='dist',sd=TRUE)
```

Here all BER (balanced error rate, which we should use in our case since we have more males than females) distances seem to reach their minimum / plateau at ncomp=2, and do not seem to change further. Therefore, we will use ncomp=2 as an optimal number of PLS components to keep in the further downstream analysis. To keep more PLS components would lead to a risk of overfitting the model, therefore we should keep as few PLS components as possible. In the BER figure, different colors correspond to different distances between the reference train samples and the test samples projected on the reference during the cross-validation procedure. In every cross-validation iteration we decide whether a test sample belongs to the Male or Female clusters, and the distance measuring the proximity to a cluster can be max.dist, mahalanobis.dist and centroids.dist. After a test sample has been assigned to one of the reference clusters (Males or Females), and got a predicted label (Male or Female), it is validated against its true label, and the accuracy (or error rate in this case) demonstrate how generalizable the model is.

After we have selected the optimal number of PLS components, let us now perform tuning of the model which implies selecting most informative variables in all layers of the data by LASSO algorithm. For this purpose, we will again need to provide the design matrix which shows a-priori correlation between the phenotype of interest and the OMICs. We assume strong correlation, but one can in principle play with this parameter and find out how it influences the final result. Further, we will concentrate on the first two PLS components based on the output of "perf"-function above.

```{r tune block splsda,fig.width=10,fig.height=8}
library("mixOmics")
test.keepX=list("expr"=c(1:5),"mut"=c(1:5),"meth"=c(1:5),"drug"=c(1:5))
ptm<-proc.time()
BPPARAM <- BiocParallel::MulticoreParam(workers = parallel::detectCores()-1)
tune.omics=tune.block.splsda(X=data,Y=Y,ncomp=2,test.keepX=test.keepX,design=design, BPPARAM=BPPARAM,progressBar=FALSE,validation="Mfold",folds=2,nrepeat=5,near.zero.var=FALSE,dist = "mahalanobis.dist")
my_time<-proc.time()-ptm
print(paste0("Thus it takes ",as.numeric(my_time["elapsed"])/3600,"h to run this piece of code"))
```

LASSO selected the following numbers of most informative variables from the loading vectors for each OMICs data set for each of the two principal components:

```{r,fig.width=10,fig.height=8}
tune.omics$choice.keepX
```

Now let us perform the final sPLS-DA modelling and display PCA plots and loadings. Here, however, for simplicity we will use all available mutations due to its special binary structure. Genetic variation data is notoriously difficult to handle by traditional statistical approaches that assume normality, therefore PLS (as well as many statistical other tools) has unfortunately problems working with this type of data.

```{r,fig.width=10,fig.height=8}
#list.keepX=list("expr"=tune.omics$choice.keepX$expr,"mut"=tune.omics$choice.keepX$mut,"meth"=tune.omics$choice.keepX$meth,"drug"=tune.omics$choice.keepX$drug)
list.keepX=list("expr"=tune.omics$choice.keepX$expr,"mut"=c(dim(mut)[2],dim(mut)[2]),"meth"=tune.omics$choice.keepX$meth,"drug"=tune.omics$choice.keepX$drug)
res=block.splsda(X=data,Y=Y,ncomp=2,keepX=list.keepX,design=design,near.zero.var=FALSE)
selectVar(res,block='expr',comp=1)$expr$name
plotIndiv(res,legend=TRUE,title="CLL Omics",ellipse=FALSE,ind.names=FALSE,cex=2)
```

Here we can visualize the samples for the 4 OMICs using their features selected via the integration procedure. Male and female samples seem to be clearly separated for gene expression and perhaps methylation, while no obvious separation is present for the mutation and drug OMICs. Let us look at the loadings in order to see what features drive the scattering of the samples in the plots above:

```{r,fig.width=10,fig.height=8}
plotLoadings(res,comp=1,contrib='max',method='mean')
plotLoadings(res,comp=2,contrib='max',method='mean')
```

Now, we will diplay the data points using PCA principal components from individual OMICs. This can be an informative way to to visualise correlation between components from different data sets. In our case mutation and drug OMICs seem to have the strongest correlation.

```{r,fig.width=10,fig.height=8}
plotDiablo(res,ncomp=1)
plotDiablo(res,ncomp=2)
```

Now let us display so-called “arrow plot” which demonstrates the samples (individuals) for all the OMICS superimposed. Here, each sample will be indicated using arrows. The start of the arrow indicates the centroid / consensus between all OMICS for a given individual, and the tips of the arrows the location of that individual in each individual OMICs. Short arrows indicate a strong agreement between the matching OMICs, long arrows a disagreement between the matching OMICs. 

```{r,fig.width=10,fig.height=8}
plotArrow(res,ind.names=FALSE,legend=TRUE,title="CLL Omics Integration")
```

The Arrow Plot can be viewed as a Consensus Plot between all the OMICs on the sample level, the X and Y coordiantes correspond to common latent variables between the OMICs, i.e. we projected all the OMICs into some common latent space where they loose information about their technological origin and hence can be superimposed onto each other.

Now we will display the so-called Correlation Circle Plot, where the top loadings variables from each of the OMICs are superimposed. Clustering of variables around the poles of the circle implies strong correlation between the variables from different OMICs. Variables on the opposite poles of the correlation circle plot imply strong anti-correlation. While the Arrow Plot serves as a main integrative OMICs visualization on sample level, the Correlation Circle Plots can be considered as a main visualization of OMICs integration on the feature level.

```{r,fig.width=10,fig.height=8}
plotVar(res,var.names=TRUE,style='ggplot2',legend=TRUE,cex=c(3,3,3,3),col=c('blue','red2','darkgreen','darkorange'))
```

For further visualization of the results of OMICs integration, we will calculate the so-called Circos Plot that diaplays correlations between features from different OMICs dat sets. Please note that here for the Circos Plot, as well as for the Correlation Circle Plot above, the features were selected simultaneously from all the OMICs when performing integration, i.e. they are not equavivalent to those obtained from each individual OMIC separately.

```{r,fig.width=10,fig.height=8}
circosPlot(res,cutoff=0.7,line=FALSE,size.variables=0.5)
```

Correlation network is another handy way to demostrate correlations between top loadings of the OMICs data sets in a pairwise fashion. Here we can choose a pair of OMICs and make a network comprising most informative features (selected via the integration) from both OMICs. The color of the edges corresponds to the strength of the correlation between the OMICs. One can use a "cutoff" paremeter to display one the the edges above a specified threshold.

```{r,fig.width=10,fig.height=8}
network(res,blocks=c(1,2),cex.node.name=0.6,color.node=c('blue','red2'),breaks=NULL)
network(res,blocks=c(1,3),cex.node.name=0.6,color.node=c('blue','darkgreen'),breaks=NULL)
network(res,blocks=c(1,4),cex.node.name=0.6,color.node=c('blue','darkorange'),breaks=NULL)
network(res,blocks=c(2,3),cex.node.name=0.6,color.node=c('red2','darkgreen'),breaks=NULL)
network(res,blocks=c(2,4),cex.node.name=0.6,color.node=c('red2','darkorange'),breaks=NULL)
network(res,blocks=c(3,4),cex.node.name=0.6,color.node=c('darkgreen','darkorange'),breaks=NULL)
```

Finally, the correlation heatmap displays strongly correlated blocks of gene expression, methylation, mutations and drug markers. This is an unsupervised visualization equivalent to hierarchical clustering. So if the DIABLO supervised integration worked well, one should ideally see females and males clusters on the y-axis and the corresponding "mixtures" of features from different OMICs that provide the clustering of the samples.

```{r,fig.width=10,fig.height=8}
cimDiablo(res,margins=c(11,18))
```
Here the color of a cell on a heatmap demonstrates abundance of the feature in the sample that are "coordinates" of the cell. In our case we can see somewhat noticable separation between males and females and can clearly see that it is not a single OMIC that drives this separation but the features from different OMICs are well mixed on the x-axis.

### Predict Gender from OMICs Integration

Now it is time for prediction. Once we have trained the PLS-DA model, we can use it and utilize the 60 test samples for making prediction of their gender and accessing the accuracy of the prediction:

```{r,fig.width=10,fig.height=8}
data.test<-list(expr=expr_test,mut=mut_test,meth=meth_test,drug=drug_test)
lapply(data.test, dim) 

predict.diablo=predict(res,newdata=data.test,dist='mahalanobis.dist')
#auroc.diablo=auroc(res,newdata=data.test,outcome.test=Y.test,plot=TRUE,roc.comp=c(1),roc.block=c(1,2,3,4))
data.frame(predict.diablo$class,Truth=Y.test)

table(predict.diablo$MajorityVote$mahalanobis.dist[,1],Y.test)
round((sum(diag(table(predict.diablo$MajorityVote$mahalanobis.dist[,1],Y.test)))/sum(table(predict.diablo$MajorityVote$mahalanobis.dist[,1],Y.test)))*100)
```

Therefore the the success rate of the first predictive component is `r round((sum(diag(table(predict.diablo$MajorityVote$mahalanobis.dist[,1],Y.test)))/sum(table(predict.diablo$MajorityVote$mahalanobis.dist[,1],Y.test)))*100)`%, it is quite high and hopefully provides new candidate bio-markers (see the variable plot) from different OMICs and connections between them for understanding of mechanisms of CLL pathogenesis.

### RSession
```{r,fig.width=10,fig.height=8}
sessionInfo()
```


