In CNRGH/pintmf: Penalized Integrative Matrix Factorization for Multiomics Data

```r knitr::opts_chunk$set(echo = TRUE)

## Useful packages

```r
library(Matrix)
library(tidyverse)
library(PintMF)
library(dplyr)
library(ggplot2)
library(gridExtra)
library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)
library(matrixStats)

Run Bxd

library(PintMF)
data(bxd)

str(bxd)
remove_affy <- grep("Affy", colnames(bxd[[3]]))
bxd_filtered <- bxd
bxd_filtered[[3]] <- bxd_filtered[[3]][, -remove_affy]
true.clust.bxd <- rep(c("CD", "HFD"), times=c(33, 31))

library(future)
R_bxd <- lapply(2:5, function(p) SolveInt(Y=bxd_filtered, p=p,
                                          max.it=1,
                                          verbose=TRUE,
                                          init_flavor="snf",
                                          flavor_mod="glmnet"))

pves <- sapply(R_bxd, function (rr){
  rr$pve[length(rr$pve)]
})
df_pve <- data.frame(pve=pves, iteration=2:(length(R_bxd)+1))



df_bics <- mapply( function(res, Y) (computeLL(Y=bxd_filtered, res=res)), R_bxd, list(bxd_filtered)) %>% t %>% as.data.frame

coph <- sapply(R_bxd, function (rr) cor(hclust(dist(rr$W), method="ward.D2") %>% cophenetic, dist(rr$W)))
df_coph <- data.frame(cophenetic=coph, iteration=2:(length(coph)+1))

g_c <- ggplot(df_coph, aes(x=iteration, y=cophenetic))+geom_point()+theme_bw()+geom_line()+scale_x_continuous(breaks=1:10)+xlab("p")


g_bic <- df_bics %>% ggplot(aes(x=p, y=RSS)) +geom_point()+ theme(legend.position="bottom")+scale_x_continuous()+theme_bw()+theme_bw()+ theme(legend.position="bottom")+scale_x_continuous(
  breaks = 1:8)+xlab("p")+ylab("MSE")+geom_line()
g_bic


g_pve <- ggplot(df_pve, aes(x=iteration, y=pve))+geom_point()+theme_bw()+geom_line()+scale_x_continuous(breaks=1:8)+xlab("p")

g <- gridExtra::grid.arrange(g_bic, g_pve,g_c, ncol=3)
g
#ggsave("../../Figs/BXD_bics_pve.pdf", g, width =7, height=3.5)

my_meth <- SolveInt(Y=bxd_filtered, p=2, max.it=5,verbose=TRUE,init_flavor="snf",flavor_mod="glmnet")

library(RColorBrewer)
library(gplots)
cols <- brewer.pal(4, "Set2")
clust <- my_meth$W %>% dist %>% hclust(method="ward.D2") %>% cutree(2)
table(true.clust.bxd, clust)

pal <- rev(brewer.pal(5, "RdBu"))
heatmap.2(my_meth$W, scale="none", trace="none",RowSideColors = cols[as.factor(true.clust.bxd)],
          col = pal)

est_multiom <- lapply(1:length(bxd_filtered), function (ii) {
  rr <- my_meth$H[[ii]]
  rr%>% 
    apply(1,FUN=function(x) which(abs(x)>quantile(abs(x), 0.90)) %>% names)
})

ndat <- data %>% length
J <- sapply(data, ncol)

cluster <- apply(my_meth$W, 1, which.max)
H_transformed <- my_meth$H
H_coef <- lapply(H_transformed, function (hh) {
  s <- abs(hh[1,]-hh[2,])
  names(s) <- colnames(hh)
  s
})
H_coef <- lapply(H_transformed, function(hh) apply((hh), 2, sd))

H_sorted <- lapply(H_coef, sort,decreasing = TRUE) %>% lapply(names) %>% lapply(unique)
ht1 <- Heatmap(t(H_transformed[[2]]))

col.clust <- brewer.pal(7, "Set2")[c(4,5)]
clust_col = structure(names = c("1", "2"),col.clust)
col.true <- brewer.pal(6, "Set3")[c(5,6)]
true_col = structure(names = c("CD", "HFD"),col.true)
mat1 <- bxd_filtered[[1]][, H_sorted[[1]][1:10]] %>% t
f21 = colorRamp2(seq(min(mat1), max(mat1), length = 3), c("blue", "#EEEEEE", "red"),
                 space = "RGB")
mat2 <- bxd_filtered[[2]][, H_sorted[[2]][1:10]] %>% t
f22 = colorRamp2(seq(min(mat2), max(mat2), length = 3), c("blue", "#EEEEEE", "red"),
                 space = "RGB")
mat3 <- bxd_filtered[[3]][, H_sorted[[3]][1:10]] %>% t
f23 = colorRamp2(seq(min(mat3), max(mat3), length = 3), c("blue", "#EEEEEE", "red"),
                 space = "RGB")

ha = HeatmapAnnotation(PintMF = clust,
                       Truth = true.clust.bxd,
                       col = list( PintMF=clust_col, Truth=true_col),
                       show_annotation_name = FALSE
)


ht1 <- Heatmap(t(mat1),col = f21,column_title = "Metabolites",
               name="Metabolites",
               show_row_dend = FALSE,
               show_column_dend = FALSE,
               cluster_rows=FALSE,
               column_names_gp = gpar(fontsize =10),
               row_names_gp = gpar(col=rep(col.true,c(33,31)), fontsize =7),
               show_row_names = FALSE,heatmap_legend_param = list(direction = "horizontal"))
rownames(mat2) <- rownames(mat2) %>% stringr::str_remove(pattern=".*;")
ht2 <- Heatmap(t(mat2),col = f22,name="Proteins",,column_title = "Proteins",
               show_row_dend = FALSE,
               show_column_dend = FALSE,
               cluster_columns=FALSE,
               column_names_gp = gpar(fontsize =8),
               row_names_gp = gpar(col=rep(col.true,c(33,31)), fontsize =7),
               show_row_names = ,heatmap_legend_param = list(direction = "horizontal")
)
ht3 <- Heatmap(t(mat3), col = f23,name="RNA",column_title = "RNA",
               show_row_dend = FALSE,
               show_column_dend = FALSE,
               cluster_rows=FALSE,
               column_names_gp = gpar(fontsize =8),
               row_names_gp = gpar(col=rep(col.true,c(33,31)), fontsize =7),
               show_row_names = FALSE,heatmap_legend_param = list(direction = "horizontal"))
hc = Heatmap(clust,name="PIntMF", col = clust_col,heatmap_legend_param = list(direction = "horizontal"))
hc2 = Heatmap(true.clust.bxd,name="Truth", col = true_col,heatmap_legend_param = list(direction = "horizontal"))
ht_list <- ht1+ht2+ht3+hc+hc2

draw(ht_list, merge_legend = TRUE, heatmap_legend_side = "bottom", 
     annotation_legend_side = "bottom")

pdf(sprintf("../../Figs/PintMF_heatmap_all.pdf"), width = 6, height = 3)
draw(ht_list, merge_legend = TRUE, heatmap_legend_side = "top", 
     annotation_legend_side = "top")
dev.off()


pdf(sprintf("../../Figs/PintMF_heatmap_%s.pdf", names(bxd_filtered)[1]), width = 6, height = 6)

draw(ht1, merge_legend = TRUE, heatmap_legend_side = "bottom", 
     annotation_legend_side = "bottom")
dev.off()
pdf(sprintf("../../Figs/PintMF_heatmap_%s.pdf", names(bxd_filtered)[2]), width = 6, height = 6)

draw(ht2, merge_legend = TRUE, heatmap_legend_side = "bottom", 
     annotation_legend_side = "bottom")
dev.off()
pdf(sprintf("../../Figs/PintMF_heatmap_%s.pdf", names(bxd_filtered)[3]), width = 6, height = 6)

draw(ht3, merge_legend = TRUE, heatmap_legend_side = "bottom", 
     annotation_legend_side = "bottom")
dev.off()