############################################################################################
#### MULTIOMICS WORKSHOP 2026                                                           ####
#### LESSON 1: WORKING WITH TRANSCRIPTOME DATA                                          ####
#### =========================================                                          ####
#### This worksheet contains an incomplete version of the code presented in the lesson. ####
#### Fill in all ... with code.                                                         ####
############################################################################################
############################################################################################

### Setup

# CRAN packages
library(ggplot2)
library(RColorBrewer)
library(gplots)

# Bioconductor packages
library(limma)
library(edgeR)
library(Mus.musculus)
library(ComplexHeatmap)

theme_set(
  theme_bw(base_size = 14) + 
    theme(panel.grid = element_blank(),
          strip.background = element_blank(),
          axis.text = element_text(color = 'black'))
)

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

# Import data. If running code locally, uncomment the GitHub Pages URL line and import from there.
load(file = 'data/transcriptome.rda')
# load(file = url('https://usda-ree-ars.github.io/SEAStatsData/omics2026/transcriptome.rda'))

### Inspect and process data

dim(transcriptome$counts)

groups <- as.factor(c("tp_0h", "tp_0h", "tp_4h", "tp_24h", "tp_0h", "tp_4h", "tp_24h", "tp_0h", "tp_4h", "tp_24h", "tp_4h", "tp_24h", "tp_0h", "tp_4h", "tp_24h", "tp_0h", "tp_4h", "tp_24h"))
transcriptome$samples$group <- groups
transcriptome$samples

head(transcriptome$counts)

geneid <- rownames(transcriptome)
genes <- AnnotationDbi::select(
  Mus.musculus, 
  keys = geneid, 
  columns = c("SYMBOL", "TXCHROM", "ENTREZID"),
  keytype = "ENSEMBL"
)
dim(genes)

genes <- genes[!duplicated(genes$ENSEMBL), ]
transcriptome$genes <- genes

### Transform and filter data

cpm <- cpm(transcriptome)

lcpm <- cpm(transcriptome, log = TRUE)

summary(lcpm)

table(rowSums(transcriptome$counts) == 0)

# Makes a copy of the unfiltered data, so that we can make the figure comparing filtered vs. unfiltered later.
transcriptome_unfiltered <- transcriptome 

# Decide which genes to keep and then filter the transcriptome data by that criterion
keep.exprs <- filterByExpr(transcriptome, group = transcriptome$samples$group)
transcriptome <- transcriptome[keep.exprs, , keep.lib.sizes = FALSE]
dim(transcriptome)

### Code to draw figure illustrating the filtering

L <- mean(transcriptome_unfiltered$samples$lib.size) * 1e-6
M <- median(transcriptome_unfiltered$samples$lib.size) * 1e-6

lcpm.cutoff <- log2(10/M + 2/L)

lcpm_unfiltered <- cpm(transcriptome_unfiltered, log = TRUE)
lcpm_filtered <- cpm(transcriptome, log = TRUE)

nsamples <- ncol(transcriptome_unfiltered)
col <- brewer.pal(nsamples, "Set3")

par(mfrow = c(1, 2))
plot(density(lcpm_unfiltered[, 1]), col = col[1], lwd = 2, ylim = c(0,0.7), las = 2, main="", xlab="")
title(main = "A. Raw data", xlab = "Log-cpm")
abline(v = lcpm.cutoff, lty = 3)
for (i in 2:nsamples) {
  den <- density(lcpm_unfiltered[, i])
  lines(den$x, den$y, col = col[i], lwd = 2)
}

plot(density(lcpm_filtered[, 1]), col = col[1], lwd=2, ylim = c(0,0.26), las = 2, main="", xlab="")
title(main = "B. Filtered data", xlab = "Log-cpm")
abline(v = lcpm.cutoff, lty = 3)
for (i in 2:nsamples) {
  den <- density(lcpm_filtered[, i])
  lines(den$x, den$y, col = col[i], lwd = 2)
}

### Normalize the data

transcriptome <- calcNormFactors(transcriptome, method = "TMM")
transcriptome$samples$norm.factors

### Multidimensional scaling plot

par(mfrow = c(1, 1))
plotMDS(
  lcpm_filtered, 
  labels = groups, 
  col = as.character(factor(groups, labels = custom_colors[1:3]))
)

### Create design and contrast matrix for differential expression analysis

design <- model.matrix(~ 0 + groups)
colnames(design) <- gsub("groups", "", colnames(design))
design

contr.matrix <- makeContrasts(
  Comp_0h_v_4h = tp_0h - tp_4h,
  Comp_0h_v_24h = tp_0h - tp_24h,
  Comp_4h_v_24h = tp_4h - tp_24h,
  levels = colnames(design))
contr.matrix

### Run voom() to correct mean-variance relationship

v <- voom(transcriptome, design, plot=TRUE)

head(v$weights)

### Fit linear models and use empirical Bayes to adjust for multiple testing

vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts = contr.matrix)
efit <- eBayes(vfit)
plotSA(efit)

summary(decideTests(efit))

### Use log fold change cutoff instead of empirical Bayes

tfit <- treat(vfit, lfc = 1)
dt <- decideTests(tfit)
summary(dt)

### Look at overlapping genes and construct Venn diagrams

de.common <- which(dt[, 1] != 0 & dt[, 2] != 0)

length(de.common)
tfit$genes$SYMBOL[de.common]

vennDiagram(dt[, 1:2], circle.col = custom_colors[c(3, 1)])

de.common <- which(dt[, 1] != 0 & dt[, 3] != 0)
tfit$genes$SYMBOL[de.common]

vennDiagram(dt[, c(1, 3)], circle.col = custom_colors[c(3, 1)])

de.common <- which(dt[, 1] != 0 & dt[, 2] != 0 & dt[, 3] != 0)
length(de.common)

vennDiagram(dt, circle.col = custom_colors[c(3, 1, 2)])

### UpSet plot (alternative to Venn diagram)

nonzero_abs_dt <- abs(dt[apply(dt, 1, function(row) !all(row ==0)), ])
grays <- c("gray30", "gray15", "gray0")

nonzero_comb_matrix <- make_comb_mat(nonzero_abs_dt)
UpSet(nonzero_comb_matrix, pt_size = unit(5, "mm"),
      bg_col = "#F0F0F0", lwd = 3,
      column_title = "Overlapping differential genes",
      comb_col = grays[comb_degree(nonzero_comb_matrix)],
      right_annotation = upset_right_annotation(nonzero_comb_matrix, 
                                                gp = gpar(col = custom_colors, fill = custom_colors),
                                                annotation_name_side = "top",
                                                axis_param = list(side = "top")))

### Mean-difference plots (similar to volcano plots)

base.vs.mid <- topTreat(tfit, coef = 1, n = 20)
base.vs.mid_bylogFC <- topTreat(tfit, coef = 1, n = 20, sort.by = 'logFC')
base.vs.post <- topTreat(tfit, coef = 2, n = 20)
mid.vs.post <- topTreat(tfit, coef = 3, n = 20)

head(base.vs.mid)

plotMD(tfit, column = 1, status = dt[ ,1], main = colnames(tfit)[1], xlim = c(-5,17),
       hl.col = custom_colors)
plotMD(tfit, column = 2, status = dt[ ,2], main = colnames(tfit)[2], xlim = c(-5,17),
       hl.col = custom_colors)
plotMD(tfit, column = 3, status = dt[ ,3], main = colnames(tfit)[3], xlim = c(-5,17),
       hl.col = custom_colors)


### Heat map

base.vs.mid <- topTreat(tfit, coef = 1, n = Inf)
base.vs.post <- topTreat(tfit, coef = 2, n = Inf)
mid.vs.post <- topTreat(tfit, coef = 3, n = Inf)

topgenes_df <- rbind(base.vs.mid, base.vs.post, mid.vs.post)
topgenes_ordered_df <- topgenes_df[order(topgenes_df$adj.P.Val), ]
topgenes <- unique(c(topgenes_ordered_df$ENSEMBL))[1:20]

topgene_index <- which(v$genes$ENSEMBL %in% topgenes)

# Create a continuous color palette using colorRampPalette
heatmap_colors <- colorRampPalette(c("turquoise4", "turquoise3", "white", "salmon", "darkred"))(100)

heatmap.2(lcpm_filtered[topgene_index,], 
          scale = "row",
          labRow = v$genes$SYMBOL[topgene_index],
          labCol = factor(groups, labels = c("0h", "4h", "24h")),
          col = heatmap_colors, trace = "none", density.info = "none",
          dendrogram = "both")


### Strip plot of individual gene expression (leptin)

leptin_index <- which(genes$SYMBOL == 'Lep')

leptin_data <- data.frame(time_point = factor(groups, labels = c("0h", "4h", "24h")), LCPM = lcpm[leptin_index, ])

ggplot(leptin_data, aes(x = time_point, color = time_point, y = LCPM)) +
  geom_point(size = 3, alpha = 0.5) +
  stat_summary(fun = 'median', geom = 'point', shape = 5, size = 5) +
  theme(legend.position = 'none') +
  scale_color_manual(values = custom_colors) +
  labs(x = 'Time point') +
  ggtitle('Leptin expression')