First, we need to select k for use with RUVseq. Examine the structure of the raw pseudobulk data.

x1 <- as.factor(ySub$samples$Batch)
cols1 <- RColorBrewer::brewer.pal(7, "Set2")

par(mfrow = c(1,3))
EDASeq::plotRLE(edgeR::cpm(ySub$counts), 
                col = cols1[x1], ylim = c(-0.5, 0.5),
                main = "Raw RLE by batch", las = 2)
EDASeq::plotPCA(edgeR::cpm(ySub$counts), 
                col = cols1[x1], labels = FALSE,
                pch = 19, main = "Raw PCA by batch")
x2 <- as.factor(ySub$samples$Group)
cols2 <- RColorBrewer::brewer.pal(4, "Set1")
EDASeq::plotPCA(edgeR::cpm(ySub$counts), 
                col = cols2[x2], labels = FALSE,
                pch = 19, main = "Raw PCA by disease")

Past versions of unnamed-chunk-24-1.png

Select the value for the k parameter i.e. the number of columns of the W matrix that will be included in the modelling.

# define the sample groups
group <- factor(ySub$samples$Group)
#micro <- factor(ySub$samples$Micro_code)
sex <- factor(ySub$samples$Sex)
age <- log2(ySub$samples$Age)

for(k in 1:6){
  adj <- RUVg(ySub$counts, ctl, k = k)
  W <- adj$W
  
  # create the design matrix
  design <- model.matrix(~0 + group + W + sex + age)
  colnames(design)[1:length(levels(group))] <- levels(group)
  
  # add the factors for the replicate samples
  dups <- unique(ySub$samples$Participant[duplicated(ySub$samples$Participant)])
  dups <- sapply(dups, function(d){
    ifelse(ySub$samples$Participant == d, 1, 0)  
  }, USE.NAMES = TRUE)
  
  contr <- makeContrasts(CF.NO_MODvNON_CF.CTRL = CF.NO_MOD - NON_CF.CTRL,
                         CF.IVAvNON_CF.CTRL = CF.IVA - NON_CF.CTRL,
                         CF.LUMA_IVAvNON_CF.CTRL = CF.LUMA_IVA - NON_CF.CTRL,
                         levels = design)
  
  y <- DGEList(counts = ySub$counts)
  y <- calcNormFactors(y)
  y <- estimateGLMCommonDisp(y, design)
  y <- estimateGLMTagwiseDisp(y, design)
  fit <- glmFit(y, design)
  
  x1 <- as.factor(ySub$samples$Batch)
  cols1 <- RColorBrewer::brewer.pal(7, "Set2")
  
  par(mfrow = c(2,3))
  EDASeq::plotRLE(edgeR::cpm(adj$normalizedCounts), 
                  col = cols1[x1], ylim = c(-0.5, 0.5),
                  main = paste0("K = ", k, " RLE by batch"))
  EDASeq::plotPCA(edgeR::cpm(adj$normalizedCounts), 
                  col = cols1[x1], labels = FALSE,
                  pch = 19,
                  main = paste0("K = ", k, " PCA by batch"))
  
  x2 <- as.factor(ySub$samples$Group)
  cols2 <- RColorBrewer::brewer.pal(5, "Set1")
  EDASeq::plotPCA(edgeR::cpm(adj$normalizedCounts), 
                  col = cols2[x2], labels = FALSE,
                  pch = 19,
                  main = paste0("K = ", k, " PCA by disease"))
  
  lrt <- glmLRT(fit, contrast = contr[, 1])
  hist(lrt$table$PValue, main = paste0("K = ", k, " ", colnames(contr)[1]),
       cex.main = 0.8)
  lrt <- glmLRT(fit, contrast = contr[, 2])
  hist(lrt$table$PValue, main = paste0("K = ", k, " ", colnames(contr)[2]),
       cex.main = 0.8)
  lrt <- glmLRT(fit, contrast = contr[, 3])
  hist(lrt$table$PValue, main = paste0("K = ", k, " ", colnames(contr)[3]),
       cex.main = 0.8)

}

Past versions of unnamed-chunk-25-1.png

Past versions of unnamed-chunk-25-2.png

Past versions of unnamed-chunk-25-3.png

Past versions of unnamed-chunk-25-4.png

Past versions of unnamed-chunk-25-5.png

Past versions of unnamed-chunk-25-6.png

Test for DGE using RUVSeq and edgeR. First, create design matrix to model the sample groups and take into account the unwanted variation, age, sex, severity and replicate samples from the same individual.

# use RUVSeq to identify the factors of unwanted variation
adj <- RUVg(ySub$counts, ctl, k = 3)
W <- adj$W
  
# create the design matrix
design <- model.matrix(~ 0 + group + W + sex + age)
colnames(design)[1:length(levels(group))] <- levels(group)

# add the factors for the replicate samples
dups <- unique(ySub$samples$Participant[duplicated(ySub$samples$Participant)])
dups <- sapply(dups, function(d){
  ifelse(ySub$samples$Participant == d, 1, 0)  
}, USE.NAMES = TRUE)

design <- cbind(design, dups)
design %>% knitr::kable()

edgeR::cpm(ySub$counts, log = TRUE) %>% 
      data.frame %>%
      rownames_to_column(var = "gene") %>%
      pivot_longer(-gene, 
                   names_to = "sample", 
                   values_to = "raw") %>%
      inner_join(edgeR::cpm(adj$normalizedCounts, log = TRUE) %>% 
                   data.frame %>%
                   rownames_to_column(var = "gene") %>%
                   pivot_longer(-gene, 
                                names_to = "sample", 
                                values_to = "norm")) %>%
      left_join(rownames_to_column(ySub$samples, 
                                   var = "sample")) %>%
      mutate(Batch = as.factor(Batch)) %>%
      dplyr::filter(gene %in% c("ZFY", "EIF1AY", "XIST")) %>%
      ggplot(aes(x = Sex,
                 y = norm,
                 colour = Sex)) +
      geom_boxplot(outlier.shape = NA, colour = "grey") +
      geom_jitter(stat = "identity",
                  width = 0.15,
                  size = 1.25) +
      geom_jitter(aes(x = Sex,
                      y = raw), stat = "identity",
                  width = 0.15,
                  size = 2, 
                  alpha = 0.2,
                  stroke = 0) +
     ggrepel::geom_text_repel(aes(label = sample.id),
                             size = 2) +
      theme_classic() +
      theme(axis.text.x = element_text(angle = 90,
                                       hjust = 1,
                                       vjust = 0.5),
            legend.position = "bottom",
            legend.direction = "horizontal",
            strip.text = element_text(size = 7),
            axis.text.y = element_text(size = 6)) +
      labs(x = "Group", y = "log2 CPM") +
      facet_wrap(~gene, scales = "free_y") + 
      scale_color_brewer(palette = "Set2") +
      ggtitle("Sex gene expression check") -> p2

p2

Past versions of unnamed-chunk-27-1.png

Create the contrast matrix for the sample group comparisons.

  contr <- makeContrasts(CF.NO_MODvNON_CF.CTRL = CF.NO_MOD - NON_CF.CTRL,
                         CF.IVAvNON_CF.CTRL = CF.IVA - NON_CF.CTRL,
                         CF.LUMA_IVAvNON_CF.CTRL = CF.LUMA_IVA - NON_CF.CTRL,
                         levels = design)

contr %>% knitr::kable()

Fit the model.

y <- DGEList(counts = ySub$counts)
y <- calcNormFactors(y)
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)

cutoff <- 0.05
  
dt <- lapply(1:ncol(contr), function(i){
  decideTests(glmLRT(fit, contrast = contr[,i]),
                            p.value = cutoff)
})

s <- sapply(dt, function(d){
  summary(d)
})
colnames(s) <- colnames(contr)
rownames(s) <- c("Down", "NotSig", "Up")

pal <- c(paletteer::paletteer_d("RColorBrewer::Set1")[2:1], "grey") 

s[-2,] %>% 
  data.frame %>%
  rownames_to_column(var = "Direction") %>%
  pivot_longer(-Direction) %>%
  ggplot(aes(x = name, y = value, fill = Direction)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = value), 
            position = position_dodge(width = 0.9),
            vjust = -0.5,
            size = 3) +
  labs(y = glue("No. DGE (FDR < {cutoff})"),
       x = "Contrast") +
      scale_fill_manual(values = pal) +
  theme(axis.text.x = element_text(angle = 45,
                                   hjust = 1,
                                   vjust = 1)) +
      scale_fill_manual(values = pal)

Past versions of unnamed-chunk-30-1.png

Explore results of statistical analysis for each contrast with significant DGEs. First, setup the output directories.

outDir <- here("output","dge_analysis")
if(!dir.exists(outDir)) dir.create(outDir)
cellDir <- file.path(outDir, cell)
if(!dir.exists(cellDir)) dir.create(cellDir)

Also, perform gene set enrichment analysis (GSEA) using the cameraPR method. cameraPR tests whether a set of genes is highly ranked relative to other genes in terms of differential expression, accounting for inter-gene correlation. Prepare the Broad MSigDB Gene Ontology, Hallmark gene sets and Reactome pathways.

Hs.c2.all <- convert_gmt_to_list(here("data/c2.all.v2024.1.Hs.entrez.gmt")) 
Hs.h.all <- convert_gmt_to_list(here("data/h.all.v2024.1.Hs.entrez.gmt")) 
Hs.c5.all <- convert_gmt_to_list(here("data/c5.all.v2024.1.Hs.entrez.gmt"))

gene_sets_list <- list(HALLMARK = Hs.h.all,
                       GO = Hs.c5.all,
                       REACTOME = Hs.c2.all[str_detect(names(Hs.c2.all), "REACTOME")],
                       WP = Hs.c2.all[str_detect(names(Hs.c2.all), "^WP")]) 

Plot a detailed summary of the results.

layout <- "
      AAAA
      AAAA
      AAAA
      BBBB
      BBBB
      BBBB
      BBBB
      EEEE
      EEEE
      EEEE
      EEEE"

plot_ruv_results_summary(contr, cutoff, cellDir, gene_sets_list, gns,
                         raw_counts = ySub$counts, 
                         norm_counts = adj$normalizedCounts, 
                         group_info = data.frame(Group = ySub$samples$Group, 
                                                 sample = rownames(ySub$samples)),
                         layout,
                         pal) -> p
p

[[1]]

[[2]]

[[3]]
NULL

Compare log fold changes and statistical significance between various contrasts.

lapply(1:ncol(contr), function(i) {
  lrt <- glmLRT(fit, contrast = contr[,i])
  topTags(lrt, n = Inf) %>%
    data.frame %>%
    rownames_to_column(var = "Symbol") %>%
    dplyr::arrange(Symbol) %>%
    dplyr::rename_with(~ paste0(.x, ".", i))
}) %>% bind_cols -> all_lrt

all_lrt %>%
  mutate(IVA = ifelse(FDR.1 < 0.05 & FDR.2 < 0.05, "red",
                      ifelse(FDR.1 < 0.05 & FDR.2 >= 0.05, "orange", 
                             ifelse(FDR.1 >= 0.05 & FDR.2 < 0.05, "green",
                                    "grey")))) -> all_lrt

ggplot(all_lrt, aes(x = logFC.1,
                    y = logFC.2)) +
  geom_point(data = subset(all_lrt, IVA %in% "grey"), aes(colour = "grey")) +
  geom_point(data = subset(all_lrt, IVA %in% "green"), aes(colour = "green")) +
  geom_point(data = subset(all_lrt, IVA %in% "orange"), aes(colour = "orange")) +
  geom_point(data = subset(all_lrt, IVA %in% "red"), aes(colour = "red")) +
  ggrepel::geom_text_repel(data = subset(all_lrt, (!IVA %in% "grey")), 
                           aes(x = logFC.1, y = logFC.2, 
                               label = Symbol.1), 
                           size = 2, colour = "black", max.overlaps = 10) +
  labs(x = "log2FC CF.NO_MODvNON_CF.CTRL",
       y = "log2FC CF.IVAvNON_CF.CTRL") +
  scale_colour_identity(guide = "legend",
                        breaks = c("red", "green", "orange","grey"),
                        labels = c("Sig. in both", 
                                   "Sig. CF.IVAvNON_CF.CTRL & N.S. CF.NO_MODvNON_CF.CTRL", 
                                   "Sig. CF.NO_MODvNON_CF.CTRL & N.S. CF.IVAvNON_CF.CTRL",
                                   "N.S. in both"),
                        name = "Statistical significance") +
  theme(legend.position = "bottom",
        legend.direction = "vertical")

Past versions of unnamed-chunk-34-1.png