############################################################################################
#### MULTIOMICS WORKSHOP 2026                                                           ####
#### LESSON 3: MULTI-OMICS INTEGRATION WITH DIABLO                                      ####
#### =============================================                                      ####
#### This worksheet contains an incomplete version of the code presented in the lesson. ####
#### Fill in all ... with code.                                                         ####
############################################################################################
############################################################################################

### Setup

library(mixOmics)
library(Mus.musculus)

custom_colors <- c("salmon", "darkolivegreen3", "turquoise3")

# Import pre-transformed data.
transcriptomics_log2_z <- read.csv('data/transcriptomics_log2.csv', row.names = 1)
proteomics_log2_z <- read.csv('data/proteomics_log2.csv', row.names = 1)
metabolomics_log2_z <- read.csv('data/metabolomics_filtered_log2.csv', row.names = 1)
load('data/diablo_tune_multi.rda')

# If you are running the code on your own computer, uncomment these lines to read the files from the GitHub Pages URLs.
# transcriptomics_log2_z <- read.csv('https://usda-ree-ars.github.io/SEAStatsData/omics2026/transcriptomics_log2.csv', row.names = 1)
# proteomics_log2_z <- read.csv('https://usda-ree-ars.github.io/SEAStatsData/omics2026/proteomics_log2.csv', row.names = 1)
# metabolomics_log2_z <- read.csv('https://usda-ree-ars.github.io/SEAStatsData/omics2026/metabolomics_filtered_log2.csv', row.names = 1)
# load(file = url('https://usda-ree-ars.github.io/SEAStatsData/omics2026/diablo_tune_multi.rda'))

### Set up data object

data <- list(transcripts = ...,
             metabolites = t(na.omit(metabolomics_log2_z)),
             proteins = t(na.omit(proteomics_log2_z)))

lapply(data, ...)

lapply(data, ...)

Y <- factor(rep(c('0h', '4h', '24h'), each = 6), levels = c('0h', '4h', '24h'))

### Get gene annotations

genes <- AnnotationDbi::select(
  Mus.musculus, 
  keys = colnames(data$transcripts), 
  columns = c("SYMBOL", "TXCHROM", "ENTREZID"),
  keytype = "ENSEMBL"
)

genes <- genes[!duplicated(genes$ENSEMBL), ]
colnames(data$transcripts)[!duplicated(genes$SYMBOL) & !is.na(genes$SYMBOL)] <- genes$SYMBOL[!duplicated(genes$SYMBOL) & !is.na(genes$SYMBOL)]

### Pairwise PLS to get correlations

pls1 <- ...(data..., data..., ncomp = ...)
pls2 <- pls(data[["transcripts"]], data[["proteins"]], ncomp = 1)
pls3 <- pls(data[["metabolites"]], data[["proteins"]], ncomp = 1)

# calculate correlation of transcripts and metabolites
print(...(pls1$..., pls1$...))
# calculate correlation of transcripts and proteins
print(cor(pls2$variates$X, pls2$variates$Y))
# calculate correlation of metabolites and proteins
print(cor(pls3$variates$X, pls3$variates$Y))

### Generate design matrix

design <- matrix(..., ... = length(data), ... = length(data), 
                dimnames = list(names(data), names(data)))
...(design) <- ... 

design

### Fit initial DIABLO model and tune for number of components

basic.diablo.model <- ...(X = ..., Y = ..., ncomp = ..., design = ...)

perf.diablo <- ...(basic.diablo.model, validation = ..., 
                    folds = ..., nrepeat = ..., seed = 1234)

plot(perf.diablo) 

perf.diablo$choice.ncomp$...

### Tune for number of features per component
### NOTE: Do not run the code between the two dashed lines!

# ---------------------------------------------------------------------------- #

test.keepX <- list(
  ... = c(1, 2, 3, 4, 5, 10, 15, 20, 30, 40, 50),
  ... = c(1, 2, 3, 4, 5, 10, 15, 20, 30, 40, 50),
  ... = c(1, 2, 3, 4, 5, 10, 15, 20, 30, 40, 50)
)

tune.MULTI <- ...(X = data, Y = Y, ncomp = ..., 
                               test.keepX = test.keepX, design = design,
                               validation = 'Mfold', folds = 3, nrepeat = ...,
                               dist = "...", seed = 123)

# ---------------------------------------------------------------------------- #

list.keepX <- tune.MULTI$...
list.keepX

### Fit final model

final.diablo.model <- ...(X = data, Y = Y, ncomp = 2, 
                          keepX = list.keepX, design = design)

### Show selected features

...(final.diablo.model, block = '...', comp = ...)$transcripts$... 
selectVar(final.diablo.model, block = 'metabolites', comp = 1)$metabolites$name 
selectVar(final.diablo.model, block = 'proteins', comp = 1)$proteins$name

selectVar(final.diablo.model, block = 'transcripts', comp = 2)$transcripts$name 
selectVar(final.diablo.model, block = 'metabolites', comp = 2)$metabolites$name 
selectVar(final.diablo.model, block = 'proteins', comp = 2)$proteins$name 

### Diagnostic correlation plot

...(final.diablo.model, ncomp = ..., col.per.group = custom_colors)

### Scatter plots

...(..., ind.names = FALSE, legend = TRUE, 
          title = 'DIABLO Sample Plots', col = custom_colors)

### Arrow plot

...(..., ind.names = FALSE, legend = TRUE, 
          title = 'DIABLO', col = custom_colors)

### Correlation circle plot

...(final.diablo.model, var.names = FALSE, 
        style = 'lattice', legend = TRUE,
        pch = c(16, 17, 15), cex = c(0.8, 0.8, 0.8), 
        col = custom_colors)

### Circos plots

...(final.diablo.model, cutoff = ..., line = TRUE,
           color.blocks = custom_colors,
           color.cor = c("mediumpurple1","grey20"), size.labels = 1,
           size.variables = 0.4, linkWidth = 0.3, showIntraLinks = FALSE, var.adj = -0.33)

### Loadings plots

...(final.diablo.model, 
             comp = ..., block = '...',
             contrib = '...', method = '...', 
             size.name = 0.6, legend.color = custom_colors)

plotLoadings(final.diablo.model, 
             comp = 1, block = '...',
             contrib = 'max', method = 'median', 
             size.name = 0.6, legend.color = custom_colors)

### Heat map

color_palette <- colorRampPalette(c("turquoise4", "turquoise3", "paleturquoise1", "lemonchiffon", "peachpuff", "salmon", "darkred"))(100)
...(..., margins = c(2, 18), size.legend = 1, color.blocks = custom_colors, color.Y = c("gray60", "gray90", "gray70"), color = color_palette)

### Cross-validation to evaluate final model performance

perf.diablo <- ...(final.diablo.model, validation = '...',
                    folds = ..., nrepeat = ..., 
                    dist = 'centroids.dist', seed = 321) 

perf.diablo$WeightedVote.error.rate