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$Status_sub)
cols2 <- RColorBrewer::brewer.pal(4, "Set1")
EDASeq::plotPCA(edgeR::cpm(ySub$counts), 
                col = cols2[x2], labels = FALSE,
                pch = 19, main = "Raw PCA by disease")

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$Status_sub)
#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$Status_sub)
  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)

}

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

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)

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.

if(!file.exists(here("data/Hs.c5.all.v7.1.entrez.rds")))
  download.file("https://bioinf.wehi.edu.au/MSigDB/v7.1/Hs.c5.all.v7.1.entrez.rds",
                here("data/Hs.c5.all.v7.1.entrez.rds"))
Hs.c5.all <- readRDS(here("data/Hs.c5.all.v7.1.entrez.rds"))

if(!file.exists(here("data/Hs.c2.cp.reactome.v7.1.entrez.rds")))
  download.file("https://bioinf.wehi.edu.au/MSigDB/v7.1/Hs.c2.cp.reactome.v7.1.entrez.rds",
                here("data/Hs.c2.cp.reactome.v7.1.entrez.rds"))
Hs.c2.reactome <- readRDS(here("data/Hs.c2.cp.reactome.v7.1.entrez.rds"))

if(!file.exists(here("data/Hs.h.all.v7.1.entrez.rds")))
  download.file("https://bioinf.wehi.edu.au/MSigDB/v7.1/Hs.h.all.v7.1.entrez.rds",
                here("data/Hs.h.all.v7.1.entrez.rds"))
Hs.h.all <- readRDS(here("data/Hs.h.all.v7.1.entrez.rds"))

Test two fibrosis gene sets using both over-represenationa analysis (ORA) and gene set enrichment analysis (GSEA). To adjust for gene length bias in ORA, setup the transcript length.

# translate gene SYMBOLS to entrez gene IDs 
gns <- AnnotationDbi::mapIds(org.Hs.eg.db, 
                             keys = rownames(fit), 
                             column = c("ENTREZID"),
                             keytype = "SYMBOL",
                             multiVals = "first")

# get transcript lengths
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
txlen <- transcriptLengths(txdb)
txlen %>%
  dplyr::select(gene_id, tx_len) %>%
  group_by(gene_id) %>%
  summarise(max_len = max(tx_len)) %>%
  dplyr::filter(gene_id %in% gns) -> supertxlen

m <- match(supertxlen$gene_id, gns)
all(supertxlen$gene_id == gns[m])

[1] TRUE

p <- vector("list", ncol(contr))

for(i in 1:(ncol(contr))){
  lrt <- glmLRT(fit, contrast = contr[,i])
  topTags(lrt, n = Inf) %>%
    data.frame -> top
  
  if(sum(top$FDR < cutoff) > 0){
    # top DGE results
    write.table(top, 
                file = file.path(cellDir, glue("{colnames(contr)[i]}.csv")),
                sep = ",", quote = F, col.names = NA)
    
    # test DGEs for over-representation of fibrosis gene sets
    # read in fibrosis gene sets
    read.csv2(file = here("data/fibrosis_gene_sets.csv"),
              header = TRUE, sep = ",") %>% 
      pivot_longer(cols = everything()) %>%
      mutate(value = str_trim(value),
             entrez = unname(unlist(AnnotationDbi::mapIds(org.Hs.eg.db, 
                                                          keys = value, 
                                                          column = c("ENTREZID"),
                                                          keytype = "SYMBOL",
                                                          multiVals = "first"))[value])) %>%
      dplyr::rename("symbol" = "value",
                    "source" = "name") -> fibrosis_table
      
      fibrosis_list <- split(fibrosis_table$entrez, fibrosis_table$source)
    
    deg <- rownames(top)[top$FDR < cutoff]
    gsa <- topGSA(gsaseq(unname(gns[deg]),
                         universe = unname(gns),
                         collection = fibrosis_list,
                         plot.bias = FALSE,
                         gene.length = supertxlen$max_len[match(unname(gns), 
                                                                supertxlen$gene_id)]))
    
    # get gene indices for all fibrosis gene set genes
    f.id <- ids2indices(fibrosis_list, unname(gns[rownames(lrt)]))
    
    # run camera competitive gene set test  
    # use signed likelihood ratio test statistic as recommended by GS here: https://support.bioconductor.org/p/112937/
    statistic <- sign(lrt$table$logFC) * sqrt(lrt$table$LR)
    f.cam <- cameraPR(statistic, f.id)
    
    # write fibrosis over-representation results to file
    write.table(gsa %>%
                  data.frame %>%
                  rownames_to_column(var = "Set"), 
                file = file.path(cellDir, glue("fibrosis.{colnames(contr)[i]}.csv")),
                sep = ",", quote = F, col.names = NA)    
    
    # write fibrosis GSEA results to file
    write.table(f.cam, 
                file = file.path(cellDir, glue("fibrosis.GSEA.{colnames(contr)[i]}.csv")),
                sep = ",", quote = F, col.names = NA)
    
    # get gene indices for all gene set genes
    h.id <- ids2indices(Hs.h.all, unname(gns[rownames(lrt)]))
    c2.id <- ids2indices(Hs.c2.reactome, unname(gns[rownames(lrt)]))
    c5.id <- ids2indices(Hs.c5.all, unname(gns[rownames(lrt)]))
    
    # run camera competitive gene set test  
    # use signed likelihood ratio test statistic as recommended by GS here: https://support.bioconductor.org/p/112937/
    h.cam <- cameraPR(statistic, h.id)
    c2.cam <- cameraPR(statistic, c2.id)
    c5.cam <- cameraPR(statistic, c5.id)
    
    # write top 50 GSEA results to file
    write.table(h.cam[1:50,], 
                file = file.path(cellDir, glue("h.{colnames(contr)[i]}.csv")),
                sep = ",", quote = F, col.names = NA)
    write.table(c2.cam[1:50,], 
                file = file.path(cellDir, glue("c2.{colnames(contr)[i]}.csv")),
                sep = ",", quote = F, col.names = NA)
    write.table(c5.cam[1:50,], 
                file = file.path(cellDir, glue("c5.{colnames(contr)[i]}.csv")),
                sep = ",", quote = F, col.names = NA)
    
    top %>% 
      mutate(sig = ifelse(FDR <= cutoff, glue("<= {cutoff}"), 
                          glue("> {cutoff}"))) %>%
      rownames_to_column(var = "SYMBOL") %>%
      left_join(dt[[i]][,1] %>% 
                  data.frame %>%
                  rownames_to_column(var = "SYMBOL") %>%
                  dplyr::rename(status = 2)) %>%
      mutate(status = ifelse(status == 1, "Up",
                             ifelse(status == -1, "Down",
                                    "NotSig"))) %>%
      mutate(status = as.factor(status)) %>%
      mutate(status = fct_relevel(status, "NotSig", after = Inf)) %>%
      ggplot(aes(x = logFC, y = -log10(PValue), color = status)) +
      geom_point(alpha = 0.75) +
      ggrepel::geom_text_repel(data = function(x) subset(x, FDR < cutoff), 
                               aes(x = logFC, y = -log10(PValue), 
                                   label = SYMBOL), 
                               size = 2, colour = "black", max.overlaps = 15) +
      labs(x = expression(~log[2]~"(Fold Change)"), 
           y = expression(~-log[10]~"(P-value)"),
           colour = glue("FDR < {cutoff}")) +
      scale_colour_manual(values = pal) +
      theme_classic() +
      theme(legend.position = "bottom") -> p1
    
    # plot up to top 9 DGE
    grps <- strsplit2(colnames(contr)[i], "v")[1,]
    
    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(Status_sub %in% names(contr[,i])[abs(contr[,i]) > 0],
                    gene %in% rownames(top)[1:min(9, max(which(top$FDR < cutoff)))]) %>%
      mutate(Group = ifelse(str_detect(Group, str_remove(grps[1], "CF.")),
                            grps[1], 
                            grps[2])) %>%
      ggplot(aes(x = Group,
                 y = norm,
                 colour = Group)) +
      geom_boxplot(outlier.shape = NA, colour = "grey") +
      geom_jitter(stat = "identity",
                  width = 0.15,
                  size = 1.25) +
      geom_jitter(aes(x = Group,
                      y = raw), stat = "identity",
                  width = 0.15,
                  size = 1.25, 
                  alpha = 0.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(colnames(contr)[i]) -> p2
    
    # plot fibrosis gene set over-representation results
    gsa %>% 
      data.frame %>%
      rownames_to_column(var = "Set") %>%
      mutate(Rank = 1:n()) %>%
      ggplot(aes(x = -log10(FDR), y = fct_reorder(Set, -Rank), colour = DE/N*100)) +
      geom_point(aes(size = N)) +
      geom_vline(xintercept = -log10(0.05),
                 linetype = "dashed")  +
      scale_colour_viridis_c(option = "plasma") +
      labs(y = "Gene set",
           colour = "% DEGs",
           size = "Set size") +
      ggtitle("Fibrosis over-represenation") +
      theme_classic(base_size = 8) +
      theme(legend.key.height = unit(0.3, "cm"),
            legend.key.width = unit(0.3, "cm"),
            legend.text = element_text(size = 7),
            title = element_text(size = 8),
            legend.position = "bottom",
            legend.direction = "vertical") -> p3
    
    # plot fibrosis GSEA results
    f.cam %>% 
      data.frame %>%
      rownames_to_column(var = "Set") %>%
      mutate(Rank = 1:n()) %>%
      ggplot(aes(x = -log10(FDR), y = fct_reorder(Set, -Rank),
                 colour = Direction)) +
      geom_point(aes(size = NGenes)) +
      geom_vline(xintercept = -log10(0.05),
                 linetype = "dashed")  +
      scale_colour_manual(values = pal) +
      labs(y = "Gene set") +
      ggtitle("Fibrosis GSEA") +
      theme_classic(base_size = 8) +
      theme(legend.key.height = unit(0.3, "cm"),
            legend.key.width = unit(0.3, "cm"),
            legend.text = element_text(size = 7),
            title = element_text(size = 8),
            legend.position = "bottom",
            legend.direction = "vertical") -> p4
    
    # GSEA top 10 plots
    h.cam %>%
      dplyr::slice(1:10) %>%
      rownames_to_column(var = "Set") %>%
      mutate(Type = "HALLMARK",
             Rank = 1:10) %>%
      bind_rows(c2.cam %>%
                  dplyr::slice(1:10) %>%
                  rownames_to_column(var = "Set") %>%
                  mutate(Type = "REACTOME",
                         Rank = 1:10)) %>%
      bind_rows(c5.cam %>%
                  dplyr::slice(1:10) %>%
                  rownames_to_column(var = "Set") %>%
                  mutate(Type = "GO",
                         Rank = 1:10)) %>%
      mutate(Set = str_wrap(str_replace_all(Set, "_", " "), width = 75),
             Set = str_remove_all(Set, "GO |REACTOME |HALLMARK ")) %>%
      ggplot(aes(x = -log10(FDR), y = fct_reorder(Set, -Rank),
                 colour = Direction)) +
      geom_point(aes(size = NGenes)) +
      facet_wrap(~Type, ncol = 1, scales = "free_y") +
      geom_vline(xintercept = -log10(0.05),
                 linetype = "dashed")  +
      scale_colour_manual(values = pal) +
      labs(y = "Gene set") +
      theme_classic(base_size = 10) -> p5

    layout <- "
      AAAA
      AAAA
      AAAA
      BBBB
      BBBB
      EEEE
      EEEE
      EEEE
      EEEE"
    
    p[[i]] <- wrap_elements(p1 + p2) + 
      wrap_elements(p3 + p4) + 
      wrap_elements(p5) +
      plot_layout(design = layout)
  }
}

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")