inst/doc/ConsensusOPLS.R

In ConsensusOPLS: Consensus OPLS for Multi-Block Data Fusion

## ----setup, class.source='fold-hide', warning=FALSE---------------------------
#install.packages("knitr")
library(knitr)
opts_chunk$set(echo = TRUE)

# To ensure reproducibility
set.seed(12)

## ----packages_installation, eval=FALSE, class.source='fold-hide'--------------
# pkgs <- c("ggplot2", "ggrepel", "DT", "psych", "ConsensusOPLS")
# sapply(pkgs, function(x) {
#   if (!requireNamespace(x, quietly = TRUE)) {
#     install.packages(x)
#   }
# })
# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# if (!requireNamespace("ComplexHeatmap", quietly = TRUE)) {
#   BiocManager::install("ComplexHeatmap")
# }

## ----packages_load, warning=FALSE, message=FALSE, class.source='fold-hide'----
library(ggplot2)        # to make beautiful graphs
library(ggrepel)        # to annotate ggplot2 graph
library(DT)             # to make interactive data tables
library(psych)          # to make specific quantitative summaries
library(ComplexHeatmap) # to make heatmap with density plot
library(ConsensusOPLS)  # to load ConsensusOPLS

## ----theme_ggplot2------------------------------------------------------------
require(ggplot2)
theme_graphs <- theme_bw() + theme(strip.text = element_text(size=14),
                                   axis.title = element_text(size=16),
                                   axis.text = element_text(size=14),
                                   plot.title = element_text(size=16),
                                   legend.title = element_text(size=14))

## ----import_demo_3_Omics------------------------------------------------------
data(demo_3_Omics, package = "ConsensusOPLS")

## ----check_dims, class.source='fold-hide'-------------------------------------
# Check dimension
BlockNames <- c("MetaboData", "MicroData", "ProteoData")
nbrBlocs <- length(BlockNames)
dims <- lapply(X=demo_3_Omics[BlockNames], FUN=dim)
names(dims) <- BlockNames
dims

# Remove unuseful object for the next steps
rm(dims)

## ----check_orders_and_names, class.source='fold-hide'-------------------------
# Check rows names in any order
row_names <- lapply(X=demo_3_Omics[BlockNames], FUN=rownames)
rns <- do.call(cbind, row_names)
rns.unique <- apply(rns, 1, function(x) length(unique(x)))
if (max(rns.unique) > 1) {
  stop("Rows names are not identical between blocks.")
}

# Check order of samples
check_row_names <- all(sapply(X=row_names, FUN=identical, y = row_names[[1]]))
if (!check_row_names && max(rns.unique) == 1) {
  print("Rows names are not in the same order for all blocks.")
}

# Remove unuseful object for the next steps
rm(row_names, rns, rns.unique, check_row_names)

## ----describe_data_by_Y_function, class.source='fold-hide'--------------------
require(psych)
require(DT)
describe_data_by_Y <- function(data, group) {
  bloc_by_Y <- describeBy(x = data, group = group,
                          mat = TRUE)[, c("group1", "n", "mean", "sd",
                                          "median", "min", "max", "range",
                                          "se")]
  bloc_by_Y[3:ncol(bloc_by_Y)] <- round(bloc_by_Y[3:ncol(bloc_by_Y)], 
                                        digits = 2)
  return (datatable(bloc_by_Y))
}

## ----describe_data_by_Y_bloc1-------------------------------------------------
describe_data_by_Y(data = demo_3_Omics[[BlockNames[1]]],
                   group = demo_3_Omics$ObsNames[,1])

## ----describe_data_by_Y_bloc2-------------------------------------------------
describe_data_by_Y(data = demo_3_Omics[[BlockNames[2]]],
                   group = demo_3_Omics$ObsNames[,1])

## ----describe_data_by_Y_bloc3-------------------------------------------------
describe_data_by_Y(data = demo_3_Omics[[BlockNames[3]]],
                   group = demo_3_Omics$ObsNames[,1])

## ----scale_data, class.source='fold-hide'-------------------------------------
# Save not scaled data
demo_3_Omics_not_scaled <- demo_3_Omics

# Scaling data
demo_3_Omics[BlockNames] <- lapply(X = demo_3_Omics[BlockNames], 
                                   FUN = function(x){
                                     scale(x, center = TRUE, scale = TRUE)
                                   }
)

## ----heatmap_function, message = FALSE, class.source='fold-hide'--------------
heatmap_data <- function(data, bloc_name, factor = NULL){
  if(!is.null(factor)){
    ht <- Heatmap(
      matrix = data, name = "Values",
      row_dend_width = unit(3, "cm"),
      column_dend_height = unit(3, "cm"),
      column_title = paste0("Heatmap of ", bloc_name),
      row_split = factor,
      row_title = "Y = %s",
      row_title_rot = 0
    )
  } else{
    ht <- Heatmap(
      matrix = data, name = "Values",
      row_dend_width = unit(3, "cm"),
      column_dend_height = unit(3, "cm"),
      column_title = paste0("Heatmap of ", bloc_name)
    )
  }
  return(ht)
}

## ----heatmap_no_scale, message = FALSE, class.source='fold-hide'--------------
# Heat map for each data block
lapply(X = 1:nbrBlocs,
       FUN = function(X){
         bloc <- BlockNames[X]
         heatmap_data(data = demo_3_Omics_not_scaled[[bloc]],
                      bloc_name = bloc,
                      factor = demo_3_Omics_not_scaled$Y[,1])})

## ----heatmap_scale, message = FALSE, class.source='fold-hide'-----------------
# Heat map for each data block
lapply(X = 1:nbrBlocs,
       FUN = function(X){
         bloc <- BlockNames[X]
         heatmap_data(data = demo_3_Omics[[bloc]],
                      bloc_name = bloc,
                      factor = demo_3_Omics$Y[,1])})

## ----heatmap_density, class.source='fold-hide'--------------------------------
# Heatmap with density for each data bloc
lapply(X = 1:nbrBlocs,
       FUN = function(X){
         bloc <- BlockNames[X]
         factor <- demo_3_Omics$Y[, 1]
         densityHeatmap(t(demo_3_Omics[[bloc]]),
                        ylab = bloc,
                        column_split  = factor,
                        column_title = "Y = %s")})

## ----rm_unscale_data, class.source='fold-hide'--------------------------------
# Remove unscaled data
rm(demo_3_Omics_not_scaled)

## ----define_cv_parameters-----------------------------------------------------
# Number of predictive component(s)
LVsPred <- 1

# Maximum number of orthogonal components
LVsOrtho <- 3

# Number of cross-validation folds
CVfolds <- nrow(demo_3_Omics[[BlockNames[[1]]]])
CVfolds

## ----run_consensusOPLSmodel---------------------------------------------------
copls.da <- ConsensusOPLS(data = demo_3_Omics[BlockNames],
                          Y = demo_3_Omics$Y,
                          maxPcomp = LVsPred,
                          maxOcomp  = LVsOrtho,
                          modelType = "da",
                          nperm = 1000,
                          cvType = "nfold",
                          nfold = 14,
                          nMC = 100,
                          cvFrac = 4/5,
                          kernelParams = list(type = "p", 
                                              params = c(order = 1)),
                          mc.cores = 1)

## ----outputs_model_summary, class.source='fold-hide'--------------------------
copls.da

## ----outputs_model_attributes, class.source='fold-hide'-----------------------
summary(attributes(copls.da))

## ----print_main_results_R2, class.source='fold-hide'--------------------------
position <- copls.da@nOcomp

paste0('R2: ', round(copls.da@R2Y[paste0("po", position)], 4))

## ----print_main_results_Q2, class.source='fold-hide'--------------------------
paste0('Q2: ', round(copls.da@Q2[paste0("po", position)], 4))

## ----print_main_results_DQ2, class.source='fold-hide'-------------------------
paste0('DQ2: ', round(copls.da@DQ2[paste0("po", position)], 4))

## ----extract_VIP, class.source='fold-hide'------------------------------------
# Compute the VIP
VIP <- copls.da@VIP

# Multiply VIP * sign(loadings for predictive component)
VIP_plot <- lapply(X = 1:nbrBlocs,
                   FUN = function(X){
                     sign_loadings <- sign(copls.da@loadings[[X]][, "p_1"])
                     result <- VIP[[X]][, "p"]*sign_loadings
                     return(sort(result, decreasing = TRUE))})
names(VIP_plot) <- BlockNames

## ----plot_VIP, class.source='fold-hide'---------------------------------------
# Metabo data
ggplot(data = data.frame(
  "variables" = factor(names(VIP_plot[[1]]),
                       levels=names(VIP_plot[[1]])[order(abs(VIP_plot[[1]]), 
                                                         decreasing=T)]), 
  "valeur" = VIP_plot[[1]]), 
  aes(x = variables, y = valeur)) +
  geom_bar(stat = "identity") +
  labs(title = paste0("Barplot of ", names(VIP_plot)[1])) +
  xlab("Predictive variables") +
  ylab("VIP x loading sign") +
  theme_graphs +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) 

# Microarray data
ggplot(data = data.frame(
  "variables" = factor(names(VIP_plot[[2]]),
                       levels=names(VIP_plot[[2]])[order(abs(VIP_plot[[2]]), 
                                                         decreasing=T)]), 
  "valeur" = VIP_plot[[2]]), 
  aes(x = variables, y = valeur)) +
  geom_bar(stat = "identity") +
  labs(title = paste0("Barplot of ", names(VIP_plot)[2])) +
  xlab("Predictive variables") +
  ylab("VIP x loading sign") +
  theme_graphs +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) 

# Proteo data
ggplot(data = data.frame(
  "variables" = factor(names(VIP_plot[[3]]),
                       levels=names(VIP_plot[[3]])[order(abs(VIP_plot[[3]]), 
                                                         decreasing=T)]), 
  "valeur" = VIP_plot[[3]]), 
  aes(x = variables, y = valeur)) +
  geom_bar(stat = "identity") +
  labs(title = paste0("Barplot of ", names(VIP_plot)[3])) +
  xlab("Predictive variables") +
  ylab("VIP x loading sign") +
  theme_graphs +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) 

## ----ggplot_score_data, class.source='fold-hide', warning=FALSE---------------
ggplot(data = data.frame("p_1" = copls.da@scores[, "p_1"],
                         "o_1" = copls.da@scores[, "o_1"],
                         "Labs" = as.matrix(unlist(demo_3_Omics$ObsNames[, 1]))),
       aes(x = p_1, y = o_1, label = Labs, 
           shape = Labs, colour = Labs)) +
  xlab("Predictive component") +
  ylab("Orthogonal component") +
  ggtitle("ConsensusOPLS Score plot")+
  geom_point(size = 2.5) + 
  geom_text_repel(size = 4, show.legend = FALSE) + 
  theme_graphs+
  scale_color_manual(values = c("#7F3C8D", "#11A579"))

## ----ggplot_data_pred_compo, class.source='fold-hide'-------------------------
ggplot(
  data = data.frame("Values" = copls.da@blockContribution[, "p_1"],
                    "Blocks" = as.factor(labels(demo_3_Omics[1:nbrBlocs]))),
  aes(x = Blocks, y = Values,
      fill = Blocks, labels = Values)) +
  geom_bar(stat = 'identity') + 
  geom_text(aes(label = round(Values, 2), y = Values), 
            vjust = 1.5, color = "black", fontface = "bold") +
  ggtitle("Block contributions to the predictive component")+
  xlab("Data blocks") +
  ylab("Weight") +
  theme_graphs +
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3"))

## ----ggplot_data_1st_ortho_compo, class.source='fold-hide'--------------------
ggplot(
  data = 
    data.frame("Values" = copls.da@blockContribution[, "o_1"],
               "Blocks" = as.factor(labels(demo_3_Omics[1:nbrBlocs]))),
  aes(x = Blocks, y = Values,
      fill = Blocks, labels = Values)) +
  geom_bar(stat = 'identity') + 
  geom_text(aes(label = round(Values, 2), y = Values), 
            vjust = 1.5, color = "black", fontface = "bold") +
  ggtitle("Block contributions to the first orthogonal component") +
  xlab("Data blocks") +
  ylab("Weight") +
  theme_graphs +
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3"))

## ----ggplot_data_pred_ortho, message = FALSE, class.source='fold-hide'--------
data_two_plots <- data.frame("Values" = copls.da@blockContribution[, "p_1"],
                             "Type" = "Pred",
                             "Blocks" = labels(demo_3_Omics[1:nbrBlocs]))
data_two_plots <- data.frame("Values" = c(data_two_plots$Values,
                                          copls.da@blockContribution[, "o_1"]),
                             "Type" = c(data_two_plots$Type,
                                        rep("Ortho", times = length(copls.da@blockContribution[, "o_1"]))),
                             "Blocks" = c(data_two_plots$Blocks,
                                          labels(demo_3_Omics[1:nbrBlocs])))

ggplot(data = data_two_plots,
       aes(x = factor(Type), 
           y = Values, 
           fill = factor(Type))) +
  geom_bar(stat = 'identity') + 
  ggtitle("Block contributions to each component")+
  geom_text(aes(label = round(Values, 2), y = Values), 
            vjust = 1.5, color = "black", fontface = "bold") +
  xlab("Data blocks") +
  ylab("Weight") +
  facet_wrap(. ~ Blocks)+
  theme_graphs+
  scale_fill_discrete(name = "Component")+
  scale_fill_manual(values = c("#7F3C8D", "#11A579"))

## ----plot_bloc_PredVSOrtho_bis, message = FALSE, class.source='fold-hide'-----
ggplot(data = data_two_plots,
       aes(x = Blocks, 
           y = Values, 
           fill = Blocks)) +
  geom_bar(stat = 'identity') +
  geom_text(aes(label = round(Values, 2), y = Values), 
            vjust = 1.5, color = "black", fontface = "bold") +
  ggtitle("Block contributions to each component") +
  xlab("Components") +
  ylab("Weight") +
  facet_wrap(. ~ factor(Type, levels = c("Pred", "Ortho"))) +
  theme_graphs +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        plot.title = element_text(hjust = 0.5, 
                                  margin = margin(t = 5, r = 0, b = 0, l = 100))) +
  scale_fill_manual(values = c("#1B9E77", "#D95F02", "#7570B3"))

## ----ggplot_data_pred_vs_ortho, message = FALSE, warning = FALSE, class.source='fold-hide'----
ggplot(data = data.frame("Pred" = copls.da@blockContribution[, "p_1"],
                         "Ortho" = copls.da@blockContribution[, "o_1"],
                         "Labels" = labels(demo_3_Omics[1:nbrBlocs])),
       aes(x = Pred, y = Ortho, label = Labels, 
           shape = Labels, colour = Labels)) +
  xlab("Predictive component") +
  ylab("Orthogonal component") +
  ggtitle("Block contributions predictive vs. orthogonal") +
  geom_point(size = 2.5) + 
  geom_text_repel(size = 4, show.legend = FALSE) + 
  theme_graphs +
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3"))

## ----create_data_loadings-----------------------------------------------------
loadings <- copls.da@loadings
data_loads <- sapply(X = 1:nbrBlocs,
                     FUN = function(X){
                       data.frame("Pred" = 
                                    loadings[[X]][, grep(pattern = "p_",
                                                         x = colnames(loadings[[X]]),
                                                         fixed = TRUE)],
                                  "Ortho" = 
                                    loadings[[X]][, grep(pattern = "o_",
                                                         x = colnames(loadings[[X]]),
                                                         fixed = TRUE)],
                                  "Labels" = labels(demo_3_Omics[1:nbrBlocs])[[X]])
                     })
data_loads <- as.data.frame(data_loads)

## ----ggplot_data_loadings, class.source='fold-hide'---------------------------
ggplot() +
  geom_point(data = as.data.frame(data_loads$V1),
             aes(x = Pred, y = Ortho, colour = Labels), 
             size = 2.5, alpha = 0.5) + 
  geom_point(data = as.data.frame(data_loads$V2),
             aes(x = Pred, y = Ortho, colour = Labels),
             size = 2.5, alpha = 0.5) +
  geom_point(data = as.data.frame(data_loads$V3),
             aes(x = Pred, y = Ortho, colour = Labels),
             size = 2.5, alpha = 0.5) +
  xlab("Predictive component") +
  ylab("Orthogonal component") +
  ggtitle("Loadings plot on first orthogonal and predictive component")+
  theme_graphs+
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3"))

## ----create_data_loadings_VIP-------------------------------------------------
loadings <- do.call(rbind.data.frame, copls.da@loadings)
loadings$block <- do.call(c, lapply(names(copls.da@loadings), function(x) 
  rep(x, nrow(copls.da@loadings[[x]]))))
loadings$variable <- gsub(paste(paste0(names(copls.da@loadings), '.'), 
                                collapse='|'), '', 
                          rownames(loadings))

VIP <- do.call(rbind.data.frame, copls.da@VIP)
VIP$block <- do.call(c, lapply(names(copls.da@VIP), function(x) 
  rep(x, nrow(copls.da@VIP[[x]]))))
VIP$variable <- gsub(paste(paste0(names(copls.da@VIP), '.'), 
                           collapse='|'), '', 
                     rownames(VIP))

loadings_VIP <- merge(x = loadings[, c("p_1", "variable")], 
                      y = VIP[, c("p", "variable")], 
                      by = "variable", all = TRUE)
colnames(loadings_VIP) <- c("variable", "loadings", "VIP")
loadings_VIP <- merge(x = loadings_VIP, 
                      y = loadings[, c("block", "variable")], 
                      by = "variable", all = TRUE)
loadings_VIP$label <- ifelse(loadings_VIP$VIP > 1, loadings_VIP$variable, NA)

## ----ggplot_data_loadings_VIP, class.source='fold-hide'-----------------------
ggplot(data = loadings_VIP,
       aes(x=loadings, y=VIP, col=block, label = label)) +
  geom_point(size = 2.5, alpha = 0.5) + 
  xlab("Predictive component") +
  ylab("Variable Importance in Projection") +
  ggtitle("VIP versus loadings on predictive components")+
  theme_graphs+
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3"))

## ----run_permutations, warning=FALSE------------------------------------------
PermRes <- copls.da@permStats

## ----plot_R2_perm, warning=FALSE, class.source='fold-hide'--------------------
ggplot(data = data.frame("R2Yperm" = PermRes$R2Y),
       aes(x = R2Yperm)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 color="grey", fill="grey") +
  geom_density(color = "blue", linewidth = 0.5) +
  geom_vline(aes(xintercept=PermRes$R2Y[1]), 
             color="blue", linetype="dashed", size=1) +
  xlab("R2 values") +
  ylab("Frequency") +
  ggtitle("R2 Permutation test")+
  theme_graphs

## ----plot_Q2_perm, warning=FALSE, class.source='fold-hide'--------------------
ggplot(data = data.frame("Q2Yperm" = PermRes$Q2Y),
       aes(x = Q2Yperm)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 color="grey", fill="grey") +
  geom_density(color = "blue", size = 0.5) +
  geom_vline(aes(xintercept=PermRes$Q2Y[1]), 
             color="blue", linetype="dashed", size=1) +
  xlab("Q2 values") +
  ylab("Frequency") +
  ggtitle("Q2 Permutation test")+
  theme_graphs

## ----plot_DQ2_perm, warning=FALSE, class.source='fold-hide'-------------------
ggplot(data = data.frame("DQ2Yperm" = PermRes$DQ2Y),
       aes(x = DQ2Yperm)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 color="grey", fill="grey") +
  geom_density(color = "blue", size = 0.5) +
  geom_vline(aes(xintercept=PermRes$DQ2Y[1]), 
             color="blue", linetype="dashed", size=1) +
  xlab("DQ2 values") +
  ylab("Frequency") +
  ggtitle("DQ2 Permutation test")+
  theme_graphs

## ----prediction---------------------------------------------------------------
reprediction <- predict(copls.da, newdata = demo_3_Omics[BlockNames])

## ----prediction_output--------------------------------------------------------
reprediction$Y
reprediction$class

## ----reproducibility----------------------------------------------------------
sessionInfo()