Inflammation of Paediatric Pulmonary Diseases

DGE analysis of CF status in macrophages

workflowr

• Summary
• Checks
• Past versions

Last updated: 2025-02-13

Checks: 7 0

Knit directory: paed-inflammation-CITEseq/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

---

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20240216)

The command set.seed(20240216) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 20f857e

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 20f857e. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:

Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/obsolete/
    Ignored:    data/.DS_Store
    Ignored:    data/C133_Neeland_batch1/
    Ignored:    data/C133_Neeland_merged/
    Ignored:    output/dge_analysis/obsolete/
    Ignored:    renv/library/
    Ignored:    renv/staging/

Untracked files:
    Untracked:  analysis/cellxgene_submission.Rmd
    Untracked:  analysis/cellxgene_submission.nb.html
    Untracked:  analysis/figure/
    Untracked:  broad_markers_seurat.csv
    Untracked:  code/background_job.R
    Untracked:  code/reverse_modifier_severity_comparisons.sh
    Untracked:  data/.Rprofile
    Untracked:  data/.gitattributes
    Untracked:  data/.gitignore
    Untracked:  data/.renvignore
    Untracked:  data/Homo_sapiens.gene_info.gz
    Untracked:  data/analysis/
    Untracked:  data/cellxgene_cell_ontologies_ann_level_1.xlsx
    Untracked:  data/cellxgene_cell_ontologies_ann_level_3 (1).xlsx
    Untracked:  data/cellxgene_cell_ontologies_ann_level_3.xlsx
    Untracked:  data/code/
    Untracked:  data/data/
    Untracked:  data/homologs.rds
    Untracked:  data/intermediate_objects/CD4 T cells.CF_samples.fit.rds
    Untracked:  data/intermediate_objects/CD4 T cells.all_samples.fit.rds
    Untracked:  data/intermediate_objects/CD8 T cells.CF_samples.fit.rds
    Untracked:  data/intermediate_objects/CD8 T cells.all_samples.fit.rds
    Untracked:  data/intermediate_objects/DC cells.CF_samples.fit.rds
    Untracked:  data/intermediate_objects/DC cells.all_samples.fit.rds

Unstaged changes:
    Modified:   .gitignore
    Modified:   analysis/06.0_azimuth_annotation.Rmd
    Modified:   analysis/09.0_integrate_cluster_macro_cells.Rmd
    Modified:   analysis/13.1_DGE_analysis_macro-alveolar.Rmd
    Modified:   analysis/13.7_DGE_analysis_macro-proliferating.Rmd
    Modified:   analysis/15.0_proportions_analysis_ann_level_1.Rmd
    Modified:   analysis/99.0_Figure_4.Rmd
    Modified:   analysis/99.0_Figure_5.Rmd
    Modified:   analysis/99.0_Supplementary_Figure_3-5.Rmd
    Modified:   code/utility.R

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/13.0_DGE_analysis_macrophages.Rmd) and HTML (docs/13.0_DGE_analysis_macrophages.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
html	8a037c6	Jovana Maksimovic	2025-01-03	Build site.
Rmd	384a134	Jovana Maksimovic	2025-01-03	wflow_publish("analysis/13.0_DGE_analysis_macrophages.Rmd")
html	b27ca47	Jovana Maksimovic	2024-12-06	Build site.
Rmd	c03dd52	Jovana Maksimovic	2024-12-06	wflow_publish("analysis/13.0_DGE_analysis_macrophages.Rmd")
html	2c0c8ab	Jovana Maksimovic	2024-12-05	Build site.
Rmd	e782782	Jovana Maksimovic	2024-12-05	wflow_publish("analysis/13.0_DGE_analysis_macrophages.Rmd")
html	fafe1fb	Jovana Maksimovic	2024-12-05	Build site.
Rmd	6ee2579	Jovana Maksimovic	2024-12-05	wflow_publish("analysis/13.0_DGE_analysis_macrophages.Rmd")
html	d424655	Jovana Maksimovic	2024-12-01	Build site.
Rmd	5d08da3	Jovana Maksimovic	2024-12-01	wflow_publish("analysis/13.0_DGE_analysis_macrophages.Rmd")
html	a6f7d42	Jovana Maksimovic	2024-12-01	Build site.
Rmd	26a4a6d	Jovana Maksimovic	2024-12-01	wflow_publish("analysis/13.0_DGE_analysis_macrophages.Rmd")

---

Load libraries

suppressPackageStartupMessages({
  library(BiocStyle)
  library(tidyverse)
  library(here)
  library(glue)
  library(Seurat)
  library(patchwork)
  library(paletteer)
  library(limma)
  library(edgeR)
  library(RUVSeq)
  library(scMerge)
  library(SingleCellExperiment)
  library(scater)
  library(tidyHeatmap)
  library(org.Hs.eg.db)
  library(TxDb.Hsapiens.UCSC.hg38.knownGene)
  library(missMethyl)
  library(ComplexHeatmap)
})

source(here("code/utility.R"))

Load Data

ambient <- ""
file <- here("data",
            "C133_Neeland_merged",
            glue("C133_Neeland_full_clean{ambient}_macrophages_annotated_diet.SEU.rds"))

seu <- readRDS(file)
seu

An object of class Seurat 
21568 features across 165209 samples within 1 assay 
Active assay: RNA (21568 features, 0 variable features)

Prepare data

Create pseudobulk samples

Use cell type and sample as our two factors; each column of the output corresponds to one unique combination of these two factors.

cell <- "macrophages"

out <- here("data",
            "C133_Neeland_merged",
            glue("C133_Neeland_full_clean{ambient}_macrophages_all_pseudobulk.rds"))

sce <- SingleCellExperiment(list(counts = seu[["RNA"]]@counts),
                            colData = seu@meta.data)
sce <- sce[, !sce$ann_level_2 %in% c("macro-T", "macro-proliferating")]

if(!file.exists(out)){
  pseudoBulk <- aggregateAcrossCells(sce, 
                                 id = colData(sce)[, "sample.id"])
  saveRDS(pseudoBulk, file = out)
  
} else {
  pseudoBulk <- readRDS(file = out)
  
}

pseudoBulk

class: SingleCellExperiment 
dim: 21568 45 
metadata(0):
assays(1): counts
rownames(21568): A1BG A1BG-AS1 ... ZNRD2 ZRANB2-AS2
rowData names(0):
colnames(45): sample_1.1 sample_15.1 ... sample_6.1 sample_7.1
colData names(71): nCount_RNA nFeature_RNA ... ids ncells
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

Code micro information

Create a factor that identifies individuals that were infected with the top 4 clinically important pathogens at time of sample collection i.e. Pseudomonas aeruginosa, Staphylococcus aureus, Haemophilus influenzae, and Aspergillus.

important_micro <- c("Pseudomonas aeruginosa", "Staphylococcus aureus",
                     "Haemophilus influenzae", "Aspergillus", "S. aureus",
                     "Staph Aureus (Methicillin Resistant)", "MRSA")

pseudoBulk$Micro_code <- sapply(strsplit(pseudoBulk$Bacteria_type, ","), function(bacteria){
  any(tolower(str_trim(bacteria)) %in% tolower(important_micro))
})

table(pseudoBulk$Micro_code)

FALSE  TRUE 
   26    19 

Filter samples

Make a DGElist object from pseudobulk data.

yPB <- DGEList(counts = counts(pseudoBulk),
               samples = colData(pseudoBulk) %>% data.frame)
dim(yPB)

[1] 21568    45

Remove genes with zero counts in all samples.

keep <- rowSums(yPB$counts) > 0 
yFlt <- yPB[keep, ]
dim(yFlt)

[1] 21557    45

Identify any samples that have too few cells for downstream statistical analysis. Examine number of cells per sample. Identify outliers and cross-reference with MDS plot. Determine a threshold for minimum number of cells per sample.

 yFlt$samples %>%
  data.frame %>%
  arrange(Group) %>%
  ggplot(aes(x = fct_inorder(sample.id), 
             y = ncells, fill = Group)) +
  geom_col() + 
  scale_fill_brewer(palette = "Set2") +
  scale_y_log10() +
  labs(x = "Sample",
       y = "Log10 No. cells") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5,
                                   size = 8),
        legend.position = "bottom") +
  geom_hline(yintercept = 500, linetype = "dashed") +
  geom_hline(yintercept = 100, linetype = "dotted") +
  geom_hline(yintercept = 50, linetype = "dashed") +
  geom_hline(yintercept = 25, linetype = "dotted")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Examine MDS plot for outlier samples.

mds_by_factor <- function(data, factor, lab){
  dims <- list(c(1,2), c(2:3), c(3,4), c(4,5))
  p <- vector("list", length(dims))
  
  for(i in 1:length(dims)){
    
    mds <- limma::plotMDS(edgeR::cpm(data, 
                                     log = TRUE), 
                          gene.selection = "common",
                          plot = FALSE, dim.plot = dims[[i]])
    
    data.frame(x = mds$x, 
               y = mds$y,
               sample = rownames(mds$distance.matrix.squared)) %>%
      left_join(rownames_to_column(data$samples, var = "sample")) -> dat
    
    p[[i]] <- ggplot(dat, aes(x = x, y = y, 
                              colour = eval(parse(text=(factor))))) +
      geom_point(size = 3) +
      ggrepel::geom_text_repel(aes(label = sample.id),
                               size = 2) +
      labs(x = glue("Principal Component {dims[[i]][1]}"),
           y = glue("Principal Component {dims[[i]][2]}"),
           colour = lab) +
      theme(legend.direction = "horizontal",
            legend.text = element_text(size = 8),
            legend.title = element_text(size = 9),
            axis.text = element_text(size = 8),
            axis.title = element_text(size = 9)) -> p[[i]]
  }
  
  wrap_plots(p, ncol = 2) + 
    plot_layout(guides = "collect") &
    theme(legend.position = "bottom")
}

mds_by_factor(yFlt, "as.factor(Batch)", "Batch") & scale_color_brewer(palette = "Set1")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "as.factor(Sex)", "Sex") & scale_color_brewer(palette = "Set2")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "log2(Age)", "Log2 Age") & scale_colour_viridis_c(option = "magma") 

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "as.factor(Group)", "Group") & scale_color_brewer(palette = "Dark2")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "as.factor(Severity)", "Severity") & scale_color_brewer(palette = "Accent")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "as.factor(Micro_code)", "Infection") & scale_color_brewer(palette = "Pastel1")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Filter out samples with less than previously determined minimum number of cells.

minCells <- 500 
yFlt <- yFlt[, yFlt$samples$ncells > minCells]
dim(yFlt)

[1] 21557    44

Re-examine MDS plots.

mds_by_factor(yFlt, "as.factor(Batch)", "Batch") & scale_color_brewer(palette = "Set1")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "as.factor(Sex)", "Sex") & scale_color_brewer(palette = "Set2")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "log2(Age)", "Log2 Age") & scale_colour_viridis_c(option = "magma") 

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "as.factor(Group)", "Group") & scale_color_brewer(palette = "Dark2")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "as.factor(Severity)", "Severity") & scale_color_brewer(palette = "Accent")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(yFlt, "as.factor(Micro_code)", "Infection") & scale_color_brewer(palette = "Pastel1")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Analyse data subsets

CF vs. non-CF controls

Prepare data

Filter genes

Filter out genes with no ENTREZ IDs and very low median expression.

gns <- AnnotationDbi::mapIds(org.Hs.eg.db,
                             keys = rownames(yFlt),
                             column = c("ENTREZID"),
                             keytype = "SYMBOL",
                             multiVals = "first")
keep <- !is.na(gns)
ySub <- yFlt[keep,]

thresh <- 0
m <- rowMedians(edgeR::cpm(ySub$counts, log = TRUE))
plot(density(m))
abline(v = thresh, lty = 2)

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

# filter out genes with low median expression
keep <- m > thresh
table(keep)

keep
FALSE  TRUE 
 5178 11274 

ySub <- ySub[keep, ]
dim(ySub)

[1] 11274    44

Examine covariates

Principal components analysis (PCA) allows us to mathematically determine the sources of variation in the data. We can then investigate whether these correlate with any of the specifed covariates.

Prepare the data.

PCs <- prcomp(t(edgeR::cpm(ySub$counts, log = TRUE)), 
              center = TRUE, retx = TRUE)
loadings = PCs$x # pc loadings


nGenes = nrow(ySub)
nSamples = ncol(ySub)

datTraits <- ySub$samples %>% dplyr::select(Batch, Disease, Micro_code,
                                            Severity, Age, Sex, ncells) %>%
  mutate(Batch = factor(Batch),
         Disease = factor(Disease, 
                            labels = 1:length(unique(Disease))),
         Sex = factor(Sex, labels = length(unique(Sex))),
         Severity = factor(Severity, labels = length(unique(Severity)))) %>%
  mutate(across(everything(), as.numeric))

moduleTraitCor <- suppressWarnings(cor(loadings[, 1:min(10, nSamples)], 
                                       datTraits, use = "p"))
moduleTraitPvalue <- WGCNA::corPvalueStudent(moduleTraitCor, (nSamples-2))

textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", 
                    signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) <- dim(moduleTraitCor)

Output results.

par(mfrow = c(2, 1))
plot(PCs, type="lines", main = cell) # scree plot

## Display the correlation values within a heatmap plot
par(cex=0.75, mar = c(3, 5, 2, 1))
WGCNA::labeledHeatmap(Matrix = t(moduleTraitCor),
                            xLabels = colnames(loadings)[1:min(10, nSamples)],
                            yLabels = names(datTraits),
                            colorLabels = FALSE,
                            colors = WGCNA::blueWhiteRed(6),
                            textMatrix = t(textMatrix),
                            setStdMargins = FALSE,
                            cex.text = 1,
                            zlim = c(-1,1),
                            main = paste0("PCA-trait relationships: Top ", 
                                          min(10, nSamples), 
                                          " PCs"))

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

RUVseq analysis

Select negative control genes

Use house-keeping genes (HKG) identified from human single-cell RNAseq experiments.

data("segList", package = "scMerge")

HKGs <- segList$human$bulkRNAseqHK
ctl <- rownames(ySub) %in% HKGs
table(ctl)

ctl
FALSE  TRUE 
 7815  3459 

Plot HKG expression profiles across all the samples.

edgeR::cpm(ySub$counts, log = TRUE) %>% 
  data.frame %>%
  rownames_to_column(var = "gene") %>%
  pivot_longer(-gene, names_to = "sample") %>%
  left_join(rownames_to_column(ySub$samples, 
                               var = "sample")) %>%
  dplyr::filter(gene %in% HKGs) %>%
  mutate(Batch = as.factor(Batch)) -> dat

dat %>%
  heatmap(gene, sample, value,
          scale = "row",
          show_row_names = FALSE,
          show_column_names = FALSE) %>%
  add_tile(Group) %>%
  add_tile(Severity) %>%
  add_tile(Batch) %>%
  add_tile(Participant) %>%
  add_tile(Age) %>%
  add_tile(Sex)

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

Version	Author	Date
d424655	Jovana Maksimovic	2024-12-01
a6f7d42	Jovana Maksimovic	2024-12-01

MDS plots based only on variablity captured by HKGs.

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Batch)", "Batch") & scale_color_brewer(palette = "Set1")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Sex)", "Sex") & scale_color_brewer(palette = "Set2")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "log2(Age)", "Log2 Age") & scale_colour_viridis_c(option = "magma") 

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Group)", "Group") & scale_color_brewer(palette = "Dark2")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Severity)", "Severity") & 
  scale_color_brewer(palette = "Accent")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Micro_code)", "Infection") & scale_color_brewer(palette = "Pastel1")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Investigate whether HKG PCAs correlate with any known covariates. Prepare the data.

PCs <- prcomp(t(edgeR::cpm(ySub$counts[ctl, ], log = TRUE)), 
              center = TRUE, retx = TRUE)
loadings = PCs$x # pc loadings


nGenes = nrow(ySub)
nSamples = ncol(ySub)

datTraits <- ySub$samples %>% dplyr::select(Batch, Disease, 
                                            Severity, Age, Sex, ncells, Micro_code) %>%
  mutate(Batch = factor(Batch),
         Disease = factor(Disease, 
                            labels = 1:length(unique(Disease))),
         Sex = factor(Sex, labels = length(unique(Sex))),
         Severity = factor(Severity, labels = length(unique(Severity)))) %>%
  mutate(across(everything(), as.numeric))

moduleTraitCor <- suppressWarnings(cor(loadings[, 1:min(10, nSamples)], 
                                       datTraits, use = "p"))
moduleTraitPvalue <- WGCNA::corPvalueStudent(moduleTraitCor, (nSamples-2))

textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", 
                    signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) <- dim(moduleTraitCor)

Output results.

par(mfrow = c(2, 1))
plot(PCs, type="lines", main = cell) # scree plot

## Display the correlation values within a heatmap plot
par(cex=0.75, mar = c(3, 5, 2, 1))
WGCNA::labeledHeatmap(Matrix = t(moduleTraitCor),
                            xLabels = colnames(loadings)[1:min(10, nSamples)],
                            yLabels = names(datTraits),
                            colorLabels = FALSE,
                            colors = WGCNA::blueWhiteRed(6),
                            textMatrix = t(textMatrix),
                            setStdMargins = FALSE,
                            cex.text = 1,
                            zlim = c(-1,1),
                            main = paste0("PCA-trait relationships: Top ", 
                                          min(10, nSamples), 
                                          " PCs"))

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Select k value

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-21-1.png

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Select the value for the k parameter i.e. the number of columns of the W matrix that will be included in the modelling based on RLE and PCA plots and p-value histograms.

# define the sample groups
group <- factor(ySub$samples$Group_severity)
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 = 0.5*(CF.NO_MOD.M + CF.NO_MOD.S) - NON_CF.CTRL,
                         CF.IVAvNON_CF.CTRL = 0.5*(CF.IVA.M + CF.IVA.S) - NON_CF.CTRL,
                         CF.LUMA_IVAvNON_CF.CTRL = 0.5*(CF.LUMA_IVA.M + CF.LUMA_IVA.S) - 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-22-1.png

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
a6f7d42	Jovana Maksimovic	2024-12-01

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 = 5)
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()

CF.IVA.M	CF.IVA.S	CF.LUMA_IVA.M	CF.LUMA_IVA.S	CF.NO_MOD.M	CF.NO_MOD.S	NON_CF.CTRL	WW_1	WW_2	WW_3	WW_4	WW_5	sexM	age	sample_34	sample_35	sample_36	sample_37	sample_38	sample_39
0	0	0	0	0	0	1	-0.2450061	-0.1265758	0.0101429	0.0111911	0.1386761	1	-0.2590872	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.1431836	-0.0324202	-0.0033971	0.1181469	0.1068408	1	-0.0939001	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.1040383	0.0068883	0.0074602	0.1336005	0.1051636	0	-0.1151479	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.0328332	0.0397275	-0.0005200	-0.0422153	-0.0780554	0	-0.0441471	0	0	0	0	0	0
0	0	0	0	1	0	0	0.0906955	0.1672770	0.0046210	0.1969342	-0.0235249	1	0.1428834	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.0888784	-0.0141316	-0.0886558	-0.2263283	0.1965863	0	-0.0729608	0	0	0	0	0	0
0	0	0	0	0	0	1	-0.2663606	-0.1165286	0.0836698	0.3708690	0.1688643	1	0.1464588	0	0	0	0	0	0
0	0	0	0	0	1	0	0.1409506	0.1900579	-0.1223152	-0.1484214	0.0435024	1	0.5597097	0	0	0	0	0	0
0	0	0	0	0	1	0	-0.0256575	0.0699757	-0.0201385	0.0863657	0.0172435	0	1.5743836	0	0	0	0	0	0
1	0	0	0	0	0	0	0.0021913	0.1083388	0.0263644	0.3517771	0.0203369	1	1.5993830	0	0	0	0	0	0
1	0	0	0	0	0	0	0.1419716	0.2165477	-0.1321409	-0.0090838	-0.0059728	1	2.3883594	0	0	0	0	0	0
0	0	0	0	0	1	0	0.1937529	0.1213905	0.3616959	0.0400948	0.0433166	0	2.2957230	0	0	0	0	0	0
0	0	0	0	1	0	0	0.0410352	-0.0050045	0.2822219	-0.1125862	0.1188874	1	2.3360877	0	0	0	0	0	0
1	0	0	0	0	0	0	0.1266179	0.0683927	0.2899064	-0.0756459	0.0360264	1	2.2980155	0	0	0	0	0	0
0	0	0	0	1	0	0	0.0788375	0.1328774	-0.1921674	-0.2585038	0.0372584	0	2.5790214	0	0	0	0	0	0
0	0	0	0	0	1	0	-0.1025824	-0.1129317	0.3031475	-0.0626907	-0.1251431	0	2.5823250	0	0	0	0	0	0
0	0	0	0	0	0	1	0.0122669	0.0841661	-0.1009722	-0.1413771	0.0865024	1	0.1321035	0	0	0	0	0	0
0	0	0	0	0	1	0	0.0619390	0.0177711	0.2864460	-0.0960927	-0.1758058	0	2.5889097	0	0	0	0	0	0
0	0	0	0	1	0	0	0.0579710	0.0187786	0.3135267	0.0106381	0.0743106	0	2.5583683	0	0	0	0	0	0
0	0	0	0	1	0	0	0.1628802	0.1035554	0.2783385	-0.0487939	-0.0066800	0	2.5670653	0	0	0	0	0	0
0	1	0	0	0	0	0	-0.1997791	-0.2030963	0.2201163	-0.3661194	-0.3477940	1	2.5730557	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.2319633	-0.0938591	-0.1241579	0.0625620	-0.4719827	0	-0.9343238	1	0	0	0	0	0
0	0	0	0	1	0	0	-0.2342364	-0.1021734	-0.0699669	0.0100297	-0.0131395	0	0.0918737	1	0	0	0	0	0
0	0	0	0	1	0	0	0.0088398	0.1142844	-0.1338937	0.0603029	-0.0737981	0	1.0409164	1	0	0	0	0	0
0	0	0	0	1	0	0	0.0120871	0.1185286	-0.1371865	0.0521637	-0.1111260	1	0.0807044	0	1	0	0	0	0
0	0	0	0	1	0	0	0.0361703	0.1400663	-0.1279626	0.0899382	-0.0463356	1	0.9940589	0	1	0	0	0	0
0	0	0	0	0	1	0	-0.1138454	-0.0103438	-0.1431217	-0.0771837	-0.0096871	0	-0.0564254	0	0	1	0	0	0
0	0	0	1	0	0	0	-0.0787942	0.0290972	-0.1187562	-0.0138534	-0.1567066	0	1.1764977	0	0	1	0	0	0
0	0	0	0	1	0	0	-0.0805890	0.0114450	-0.0130459	0.0611262	0.0774385	0	1.5597097	0	0	0	1	0	0
0	0	1	0	0	0	0	0.0315536	0.0882093	-0.0528393	-0.0540565	0.1251916	0	2.1930156	0	0	0	1	0	0
0	0	1	0	0	0	0	-0.1342613	-0.0428337	-0.0078508	0.0469788	0.0477800	0	2.2980155	0	0	0	1	0	0
1	0	0	0	0	0	0	-0.0265792	0.0205819	-0.1189050	-0.3244326	0.2354583	1	1.5703964	0	0	0	0	1	0
1	0	0	0	0	0	0	-0.0481662	0.0205167	-0.0725035	-0.1607994	0.2496796	1	2.0206033	0	0	0	0	1	0
1	0	0	0	0	0	0	-0.1326300	-0.0509487	-0.0709927	-0.1772589	0.2669947	1	2.3485584	0	0	0	0	1	0
0	0	0	0	1	0	0	-0.0338131	0.0425317	-0.0220971	-0.0227872	-0.1481780	0	1.9730702	0	0	0	0	0	1
0	0	1	0	0	0	0	0.0028655	0.0839202	0.0032780	0.1142677	-0.1877387	0	2.6297159	0	0	0	0	0	1
0	0	0	0	0	0	1	0.0545678	0.1344350	-0.1728295	-0.1011463	-0.3237410	1	0.2923784	0	0	0	0	0	0
0	0	0	0	0	1	0	0.1265987	-0.4291945	-0.1293918	0.0626068	-0.0631813	1	1.5801455	0	0	0	0	0	0
0	0	0	0	1	0	0	0.2710167	-0.3627479	-0.0752238	-0.0120605	-0.0418061	1	1.5801455	0	0	0	0	0	0
0	1	0	0	0	0	0	0.3905018	-0.2990957	-0.0967870	0.1471434	0.0930224	1	1.5993178	0	0	0	0	0	0
0	0	0	0	0	0	1	0.3526026	-0.3608052	-0.1022576	0.0198213	0.0029629	1	1.5849625	0	0	0	0	0	0
0	0	0	0	0	0	1	-0.2641826	-0.1488508	0.0405886	0.1528308	0.1551274	0	3.0699187	0	0	0	0	0	0
0	0	0	0	0	0	1	0.1079356	0.1882191	-0.0317602	0.1186608	0.0024396	1	2.4204621	0	0	0	0	0	0
0	0	0	0	0	0	1	0.0815308	0.1739614	-0.0296872	0.2133871	-0.0392139	0	2.2356012	0	0	0	0	0	0

Plot expression level of sex genes between males and females for raw and adjusted counts to check that we are not over-adjusting the counts with RUV.

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-24-1.png

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
a6f7d42	Jovana Maksimovic	2024-12-01

Create the contrast matrix for the sample group comparisons.

contr <- makeContrasts(CF.NO_MODvNON_CF.CTRL = 0.5*(CF.NO_MOD.M + CF.NO_MOD.S) - NON_CF.CTRL,
                       CF.IVAvNON_CF.CTRL = 0.5*(CF.IVA.M + CF.IVA.S) - NON_CF.CTRL,
                       CF.LUMA_IVAvNON_CF.CTRL = 0.5*(CF.LUMA_IVA.M + CF.LUMA_IVA.S) - NON_CF.CTRL,
                       levels = design)
  
contr %>% knitr::kable()

	CF.NO_MODvNON_CF.CTRL	CF.IVAvNON_CF.CTRL	CF.LUMA_IVAvNON_CF.CTRL
CF.IVA.M	0.0	0.5	0.0
CF.IVA.S	0.0	0.5	0.0
CF.LUMA_IVA.M	0.0	0.0	0.5
CF.LUMA_IVA.S	0.0	0.0	0.5
CF.NO_MOD.M	0.5	0.0	0.0
CF.NO_MOD.S	0.5	0.0	0.0
NON_CF.CTRL	-1.0	-1.0	-1.0
WW_1	0.0	0.0	0.0
WW_2	0.0	0.0	0.0
WW_3	0.0	0.0	0.0
WW_4	0.0	0.0	0.0
WW_5	0.0	0.0	0.0
sexM	0.0	0.0	0.0
age	0.0	0.0	0.0
sample_34	0.0	0.0	0.0
sample_35	0.0	0.0	0.0
sample_36	0.0	0.0	0.0
sample_37	0.0	0.0	0.0
sample_38	0.0	0.0	0.0
sample_39	0.0	0.0	0.0

Fit the model.

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

DEG results

Overall summary

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-27-1.png

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
a6f7d42	Jovana Maksimovic	2024-12-01

Save the contrast matrix, edgeR fit object and RUVseq adjusted data as an RDS object for downstream use in plotting, etc.

# Save group in fit object
fit$samples$group <- group

# save LRT results
deg_results <- list(
  contr = contr,
  fit = fit,
  adj = adj)

saveRDS(deg_results, file = here("data",
                                 "intermediate_objects",
                                 glue("{cell}.all_samples.fit.rds")))

Detailed summary

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

fibrosis <- create_custom_gene_lists_from_file(here("data/fibrosis_gene_sets.csv"))

# add fibrosis sets from REACTOME and WIKIPATHWAYS
fibrosis <- c(lapply(fibrosis, function(l) l[!is.na(l)]),
              Hs.c2.all[str_detect(names(Hs.c2.all), "FIBROSIS")])

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")],
                       FIBROSIS = fibrosis) 

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 = group, 
                                                 sample = rownames(ySub$samples)),
                         layout,
                         pal,
                         severity = rep(FALSE, ncol(contr))) -> p
p

[[1]]

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

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
8a037c6	Jovana Maksimovic	2025-01-03
b27ca47	Jovana Maksimovic	2024-12-06
2c0c8ab	Jovana Maksimovic	2024-12-05
fafe1fb	Jovana Maksimovic	2024-12-05
d424655	Jovana Maksimovic	2024-12-01
a6f7d42	Jovana Maksimovic	2024-12-01

[[2]]

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

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
8a037c6	Jovana Maksimovic	2025-01-03
b27ca47	Jovana Maksimovic	2024-12-06
2c0c8ab	Jovana Maksimovic	2024-12-05
fafe1fb	Jovana Maksimovic	2024-12-05
d424655	Jovana Maksimovic	2024-12-01
a6f7d42	Jovana Maksimovic	2024-12-01

[[3]]
NULL

DEG heatmaps

Heatmaps of up to the top 50 significant DGEs.

p <- lapply(1:ncol(contr), function(i){
    lrt <- glmLRT(fit, contrast = contr[,i])
    top <- topTags(lrt, p.value = cutoff, n = Inf) %>% data.frame
    top_deg_heatmap(top = top,
                    comparison = lrt$comparison,
                    counts = adj$normalizedCounts,
                    sample_data = ySub$samples)
})

p

[[1]]

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

[[2]]

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

[[3]]
NULL

CF modifiers and severity

Prepare data

Filter genes

Extract only the CF samples.

ySub <- yFlt[, yFlt$samples$Disease != "Healthy"]
dim(ySub)

[1] 21557    36

Filter out genes with no ENTREZ IDs and very low expression.

gns <- AnnotationDbi::mapIds(org.Hs.eg.db,
                             keys = rownames(ySub),
                             column = c("ENTREZID"),
                             keytype = "SYMBOL",
                             multiVals = "first")
keep <- !is.na(gns)
ySub <- ySub[keep,]

thresh <- 0
m <- rowMedians(edgeR::cpm(ySub$counts, log = TRUE))
plot(density(m))
abline(v = thresh, lty = 2)

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

# filter out genes with low median expression
keep <- m > thresh
table(keep)

keep
FALSE  TRUE 
 5178 11274 

ySub <- ySub[keep, ]
dim(ySub)

[1] 11274    36

Examine covariates

Principal components analysis (PCA) allows us to mathematically determine the sources of variation in the data. We can then investigate whether these correlate with any of the specifed covariates.

Prepare the data.

PCs <- prcomp(t(edgeR::cpm(ySub$counts, log = TRUE)), 
              center = TRUE, retx = TRUE)
loadings = PCs$x # pc loadings


nGenes = nrow(ySub)
nSamples = ncol(ySub)

datTraits <- ySub$samples %>% dplyr::select(Batch, Treatment, Micro_code,
                                            Severity, Age, Sex, ncells) %>%
  mutate(Batch = factor(Batch),
         Treatment = factor(Treatment, 
                            labels = 1:length(unique(Treatment))),
         Sex = factor(Sex, labels = length(unique(Sex))),
         Severity = factor(Severity, labels = length(unique(Severity)))) %>%
  mutate(across(everything(), as.numeric))

moduleTraitCor <- suppressWarnings(cor(loadings[, 1:min(10, nSamples)], 
                                       datTraits, use = "p"))
moduleTraitPvalue <- WGCNA::corPvalueStudent(moduleTraitCor, (nSamples-2))

textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", 
                    signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) <- dim(moduleTraitCor)

Output results.

par(mfrow = c(2, 1))
plot(PCs, type="lines", main = cell) # scree plot

## Display the correlation values within a heatmap plot
par(cex=0.75, mar = c(3, 5, 2, 1))
WGCNA::labeledHeatmap(Matrix = t(moduleTraitCor),
                            xLabels = colnames(loadings)[1:min(10, nSamples)],
                            yLabels = names(datTraits),
                            colorLabels = FALSE,
                            colors = WGCNA::blueWhiteRed(6),
                            textMatrix = t(textMatrix),
                            setStdMargins = FALSE,
                            cex.text = 1,
                            zlim = c(-1,1),
                            main = paste0("PCA-trait relationships: Top ", 
                                          min(10, nSamples), 
                                          " PCs"))

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

RUVseq analysis

Negative control genes

Use house-keeping genes (HKG) identified from human single-cell RNAseq experiments.

data("segList", package = "scMerge")

HKGs <- segList$human$bulkRNAseqHK
ctl <- rownames(ySub) %in% HKGs
table(ctl)

ctl
FALSE  TRUE 
 7814  3460 

Plot HKG expression profiles across all the samples.

edgeR::cpm(ySub$counts, log = TRUE) %>% 
  data.frame %>%
  rownames_to_column(var = "gene") %>%
  pivot_longer(-gene, names_to = "sample") %>%
  left_join(rownames_to_column(ySub$samples, 
                               var = "sample")) %>%
  dplyr::filter(gene %in% HKGs) %>%
mutate(Batch = as.factor(Batch)) -> dat

dat %>%
  heatmap(gene, sample, value,
          scale = "row",
          show_row_names = FALSE,
          show_column_names = FALSE) %>%
  add_tile(Group) %>%
  add_tile(Severity) %>%
  add_tile(Batch) %>%
  add_tile(Participant) %>%
  add_tile(Age) %>%
  add_tile(Sex)

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

Version	Author	Date
d424655	Jovana Maksimovic	2024-12-01
a6f7d42	Jovana Maksimovic	2024-12-01

MDS plots based only on variablity captured by HKGs.

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Batch)", "Batch") & scale_color_brewer(palette = "Set1")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Sex)", "Sex") & scale_color_brewer(palette = "Set2")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "log2(Age)", "Log2 Age") & scale_colour_viridis_c(option = "magma") 

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Group)", "Group") & scale_color_brewer(palette = "Dark2")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Severity)", "Severity") & 
  scale_color_brewer(palette = "Accent")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

mds_by_factor(ySub[rownames(ySub) %in% HKGs,], "as.factor(Micro_code)", "Infection") & scale_color_brewer(palette = "Pastel1")

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Investigate whether HKG PCAs correlate with any known covariates. Prepare the data.

PCs <- prcomp(t(edgeR::cpm(ySub$counts[ctl, ], log = TRUE)), 
              center = TRUE, retx = TRUE)
loadings = PCs$x # pc loadings


nGenes = nrow(ySub)
nSamples = ncol(ySub)

datTraits <- ySub$samples %>% dplyr::select(Batch, Treatment, 
                                            Severity, Age, Sex, ncells, Micro_code) %>%
  mutate(Batch = factor(Batch),
         Treatment = factor(Treatment, 
                            labels = 1:length(unique(Treatment))),
         Sex = factor(Sex, labels = length(unique(Sex))),
         Severity = factor(Severity, labels = length(unique(Severity)))) %>%
  mutate(across(everything(), as.numeric))

moduleTraitCor <- suppressWarnings(cor(loadings[, 1:min(10, nSamples)], 
                                       datTraits, use = "p"))
moduleTraitPvalue <- WGCNA::corPvalueStudent(moduleTraitCor, (nSamples-2))

textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", 
                    signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) <- dim(moduleTraitCor)

Output results.

par(mfrow = c(2, 1))
plot(PCs, type="lines", main = cell) # scree plot

## Display the correlation values within a heatmap plot
par(cex=0.75, mar = c(3, 5, 2, 1))
WGCNA::labeledHeatmap(Matrix = t(moduleTraitCor),
                            xLabels = colnames(loadings)[1:min(10, nSamples)],
                            yLabels = names(datTraits),
                            colorLabels = FALSE,
                            colors = WGCNA::blueWhiteRed(6),
                            textMatrix = t(textMatrix),
                            setStdMargins = FALSE,
                            cex.text = 1,
                            zlim = c(-1,1),
                            main = paste0("PCA-trait relationships: Top ", 
                                          min(10, nSamples), 
                                          " PCs"))

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Select k value

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-43-1.png

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

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_severity)
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 + micro + 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.IVAvCF.NO_MOD = 0.5*(CF.IVA.S + CF.IVA.M) - 0.5*(CF.NO_MOD.S + CF.NO_MOD.M),
                         CF.LUMA_IVAvCF.NO_MOD = 0.5*(CF.LUMA_IVA.S + CF.LUMA_IVA.M) - 0.5*(CF.NO_MOD.S + CF.NO_MOD.M),
                         CF.NO_MOD.SvCF.NO_MOD.M = CF.NO_MOD.S - CF.NO_MOD.M,
                         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-44-1.png

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Test for differences

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. Also include a factor for presence of top 4 clinically important organisms as we are only comparing CF samples which have all been tested for the presence of various microorganisms.

# use RUVSeq to identify the factors of unwanted variation
adj <- RUVg(ySub$counts, ctl, k = 5)
W <- adj$W

# create the design matrix
design <- model.matrix(~ 0 + group + W + sex + micro + 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()

CF.IVA.M	CF.IVA.S	CF.LUMA_IVA.M	CF.LUMA_IVA.S	CF.NO_MOD.M	CF.NO_MOD.S	WW_1	WW_2	WW_3	WW_4	WW_5	sexM	microTRUE	age	sample_34	sample_35	sample_36	sample_37	sample_38	sample_39
0	0	0	0	1	0	-0.1826918	-0.0421014	0.0178968	0.1651916	-0.2021682	1	0	-0.0939001	0	0	0	0	0	0
0	0	0	0	1	0	-0.1342832	0.0053763	0.0096333	0.1639781	-0.1533793	0	0	-0.1151479	0	0	0	0	0	0
0	0	0	0	1	0	-0.0464556	0.0368160	-0.0143949	-0.0087972	0.1234221	0	0	-0.0441471	0	0	0	0	0	0
0	0	0	0	1	0	0.1065643	0.1809380	-0.0793958	0.2236498	-0.0864980	1	0	0.1428834	0	0	0	0	0	0
0	0	0	0	1	0	-0.1160083	-0.0441283	-0.0736710	-0.2551975	0.0534093	0	0	-0.0729608	0	0	0	0	0	0
0	0	0	0	0	1	0.1680658	0.1669624	-0.2108829	-0.2039734	0.1804611	1	1	0.5597097	0	0	0	0	0	0
0	0	0	0	0	1	-0.0375059	0.0637373	-0.0498214	0.1046546	-0.1347553	0	1	1.5743836	0	0	0	0	0	0
1	0	0	0	0	0	-0.0026673	0.1212506	-0.0272326	0.4041595	-0.2119412	1	1	1.5993830	0	0	0	0	0	0
1	0	0	0	0	0	0.1693700	0.1905856	-0.2339872	-0.0388730	-0.0081535	1	0	2.3883594	0	0	0	0	0	0
0	0	0	0	0	1	0.2348835	0.2314087	0.2825858	0.0476305	-0.0232675	0	0	2.2957230	0	0	0	0	0	0
0	0	0	0	1	0	0.0457799	0.0684965	0.2718871	-0.0987946	-0.0440518	1	1	2.3360877	0	0	0	0	0	0
1	0	0	0	0	0	0.1515844	0.1516500	0.2417190	-0.0693998	0.0365935	1	0	2.2980155	0	0	0	0	0	0
0	0	0	0	1	0	0.0910096	0.0797471	-0.2487636	-0.3199122	-0.0793758	0	1	2.5790214	0	0	0	0	0	0
0	0	0	0	0	1	-0.1318280	-0.0486203	0.3527620	0.0232739	-0.0571623	0	1	2.5823250	0	0	0	0	0	0
0	0	0	0	0	1	0.0714860	0.0921738	0.2696563	-0.0561743	0.2316445	0	0	2.5889097	0	0	0	0	0	0
0	0	0	0	1	0	0.0668300	0.1053596	0.2892008	0.0369661	-0.0861322	0	0	2.5583683	0	0	0	0	0	0
0	0	0	0	1	0	0.1964004	0.1880283	0.2133152	-0.0436351	0.0830792	0	0	2.5670653	0	0	0	0	0	0
0	1	0	0	0	0	-0.2525909	-0.1759350	0.3260070	-0.2813416	0.2342432	1	1	2.5730557	0	0	0	0	0	0
0	0	0	0	1	0	-0.2931118	-0.1471316	-0.0621076	0.2076569	0.4225603	0	0	-0.9343238	1	0	0	0	0	0
0	0	0	0	1	0	-0.2956079	-0.1359512	-0.0101738	0.0928318	0.1543568	0	0	0.0918737	1	0	0	0	0	0
0	0	0	0	1	0	0.0047586	0.0815965	-0.1835836	0.0842886	0.1362957	0	0	1.0409164	1	0	0	0	0	0
0	0	0	0	1	0	0.0087770	0.0860591	-0.1871385	0.0777956	0.2040882	1	0	0.0807044	0	1	0	0	0	0
0	0	0	0	1	0	0.0386566	0.1123423	-0.1919839	0.1015817	0.1153893	1	1	0.9940589	0	1	0	0	0	0
0	0	0	0	0	1	-0.1469942	-0.0585424	-0.1279904	-0.0458979	0.0783640	0	0	-0.0564254	0	0	1	0	0	0
0	0	0	1	0	0	-0.1036406	-0.0090533	-0.1213031	0.0338177	0.2052745	0	1	1.1764977	0	0	1	0	0	0
0	0	0	0	1	0	-0.1054066	0.0028403	-0.0153273	0.0863246	-0.2355324	0	1	1.5597097	0	0	0	1	0	0
0	0	1	0	0	0	0.0330547	0.0751584	-0.0935320	-0.0714013	-0.1104228	0	0	2.1930156	0	0	0	1	0	0
0	0	1	0	0	0	-0.1717769	-0.0564827	0.0172470	0.1022984	-0.2104867	0	1	2.2980155	0	0	0	1	0	0
1	0	0	0	0	0	-0.0391679	-0.0188990	-0.1211888	-0.4029631	-0.1652483	1	0	1.5703964	0	0	0	0	1	0
1	0	0	0	0	0	-0.0656091	-0.0041107	-0.0779579	-0.2176521	-0.2568057	1	1	2.0206033	0	0	0	0	1	0
1	0	0	0	0	0	-0.1700388	-0.0808587	-0.0403258	-0.2160066	-0.2803701	1	0	2.3485584	0	0	0	0	1	0
0	0	0	0	1	0	-0.0476931	0.0329069	-0.0386872	0.0224894	0.0860196	0	1	1.9730702	0	0	0	0	0	1
0	0	1	0	0	0	-0.0021918	0.0847070	-0.0346133	0.1663061	0.0424025	0	1	2.6297159	0	0	0	0	0	1
0	0	0	0	0	1	0.1509259	-0.5317037	-0.0001366	0.0864766	-0.2130449	1	0	1.5801455	0	0	0	0	0	0
0	0	0	0	1	0	0.3296316	-0.4381351	0.0107532	-0.0195189	0.1056561	1	0	1.5801455	0	0	0	0	0	0
0	1	0	0	0	0	0.4774913	-0.3664876	-0.0584640	0.1181672	0.0655361	1	0	1.5993178	0	0	0	0	0	0

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-46-1.png

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Create the contrast matrix for the sample group comparisons.

contr <- makeContrasts(CF.IVAvCF.NO_MOD = 0.5*(CF.IVA.S + CF.IVA.M) - 0.5*(CF.NO_MOD.S + CF.NO_MOD.M),
                       CF.LUMA_IVAvCF.NO_MOD = 0.5*(CF.LUMA_IVA.S + CF.LUMA_IVA.M) - 0.5*(CF.NO_MOD.S + CF.NO_MOD.M),
                       CF.NO_MOD.SvCF.NO_MOD.M = CF.NO_MOD.S - CF.NO_MOD.M,
                       levels = design)

contr %>% knitr::kable()

	CF.IVAvCF.NO_MOD	CF.LUMA_IVAvCF.NO_MOD	CF.NO_MOD.SvCF.NO_MOD.M
CF.IVA.M	0.5	0.0	0
CF.IVA.S	0.5	0.0	0
CF.LUMA_IVA.M	0.0	0.5	0
CF.LUMA_IVA.S	0.0	0.5	0
CF.NO_MOD.M	-0.5	-0.5	-1
CF.NO_MOD.S	-0.5	-0.5	1
WW_1	0.0	0.0	0
WW_2	0.0	0.0	0
WW_3	0.0	0.0	0
WW_4	0.0	0.0	0
WW_5	0.0	0.0	0
sexM	0.0	0.0	0
microTRUE	0.0	0.0	0
age	0.0	0.0	0
sample_34	0.0	0.0	0
sample_35	0.0	0.0	0
sample_36	0.0	0.0	0
sample_37	0.0	0.0	0
sample_38	0.0	0.0	0
sample_39	0.0	0.0	0

Fit the model.

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

DEG results

Overall summary

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-49-1.png

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Save the contrast matrix, edgeR fit object and RUVseq adjusted data as an RDS object for downstream use in plotting, etc.

# Save group in fit object
fit$samples$group <- group

# save LRT results
deg_results <- list(
  contr = contr,
  fit = fit,
  adj = adj)

saveRDS(deg_results, file = here("data",
                                 "intermediate_objects",
                                 glue("{cell}.CF_samples.fit.rds")))

Detailed summary

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

fibrosis <- create_custom_gene_lists_from_file(here("data/fibrosis_gene_sets.csv"))
# add fibrosis sets from REACTOME and WIKIPATHWAYS
fibrosis <- c(lapply(fibrosis, function(l) l[!is.na(l)]),
              Hs.c2.all[str_detect(names(Hs.c2.all), "FIBROSIS")])

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")],
                       FIBROSIS = fibrosis)

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 = group,
                                                 sample = rownames(ySub$samples)),
                         layout,
                         pal,
                         severity = c(FALSE, FALSE, TRUE)) -> p
p

[[1]]

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

Version	Author	Date
4fb95f3	Jovana Maksimovic	2025-02-13
8a037c6	Jovana Maksimovic	2025-01-03
b27ca47	Jovana Maksimovic	2024-12-06
2c0c8ab	Jovana Maksimovic	2024-12-05
fafe1fb	Jovana Maksimovic	2024-12-05
d424655	Jovana Maksimovic	2024-12-01
a6f7d42	Jovana Maksimovic	2024-12-01

[[2]]

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

Version	Author	Date
8a037c6	Jovana Maksimovic	2025-01-03
a6f7d42	Jovana Maksimovic	2024-12-01

[[3]]

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

Version	Author	Date
8a037c6	Jovana Maksimovic	2025-01-03
a6f7d42	Jovana Maksimovic	2024-12-01

DEG heatmaps

Heatmaps of up to the top 50 significant DGEs.

p <- lapply(1:ncol(contr), function(i){
    lrt <- glmLRT(fit, contrast = contr[,i])
    top <- topTags(lrt, p.value = cutoff, n = Inf) %>% data.frame
    top_deg_heatmap(top = top,
                    comparison = lrt$comparison,
                    counts = adj$normalizedCounts,
                    sample_data = ySub$samples)
})

p

[[1]]

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

[[2]]

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

[[3]]

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

Version	Author	Date
a6f7d42	Jovana Maksimovic	2024-12-01

Session info

Session information

sessionInfo()

R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_AU.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_AU.UTF-8        LC_COLLATE=en_AU.UTF-8    
 [5] LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8   
 [7] LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
 [1] grid      parallel  stats4    stats     graphics  grDevices datasets 
 [8] utils     methods   base     

other attached packages:
 [1] ComplexHeatmap_2.18.0                              
 [2] missMethyl_1.36.0                                  
 [3] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
 [4] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1 
 [5] minfi_1.48.0                                       
 [6] bumphunter_1.44.0                                  
 [7] locfit_1.5-9.8                                     
 [8] iterators_1.0.14                                   
 [9] foreach_1.5.2                                      
[10] TxDb.Hsapiens.UCSC.hg38.knownGene_3.18.0           
[11] GenomicFeatures_1.54.3                             
[12] org.Hs.eg.db_3.18.0                                
[13] AnnotationDbi_1.64.1                               
[14] tidyHeatmap_1.8.1                                  
[15] scater_1.30.1                                      
[16] scuttle_1.12.0                                     
[17] SingleCellExperiment_1.24.0                        
[18] scMerge_1.18.0                                     
[19] RUVSeq_1.36.0                                      
[20] EDASeq_2.36.0                                      
[21] ShortRead_1.60.0                                   
[22] GenomicAlignments_1.38.2                           
[23] SummarizedExperiment_1.32.0                        
[24] MatrixGenerics_1.14.0                              
[25] matrixStats_1.2.0                                  
[26] Rsamtools_2.18.0                                   
[27] GenomicRanges_1.54.1                               
[28] Biostrings_2.70.2                                  
[29] GenomeInfoDb_1.38.6                                
[30] XVector_0.42.0                                     
[31] IRanges_2.36.0                                     
[32] S4Vectors_0.40.2                                   
[33] BiocParallel_1.36.0                                
[34] Biobase_2.62.0                                     
[35] BiocGenerics_0.48.1                                
[36] edgeR_4.0.15                                       
[37] limma_3.58.1                                       
[38] paletteer_1.6.0                                    
[39] patchwork_1.3.0                                    
[40] SeuratObject_4.1.4                                 
[41] Seurat_4.4.0                                       
[42] glue_1.8.0                                         
[43] here_1.0.1                                         
[44] lubridate_1.9.3                                    
[45] forcats_1.0.0                                      
[46] stringr_1.5.1                                      
[47] dplyr_1.1.4                                        
[48] purrr_1.0.2                                        
[49] readr_2.1.5                                        
[50] tidyr_1.3.1                                        
[51] tibble_3.2.1                                       
[52] ggplot2_3.5.0                                      
[53] tidyverse_2.0.0                                    
[54] BiocStyle_2.30.0                                   
[55] workflowr_1.7.1                                    

loaded via a namespace (and not attached):
  [1] igraph_2.0.1.1            ica_1.0-3                
  [3] plotly_4.10.4             Formula_1.2-5            
  [5] rematch2_2.1.2            zlibbioc_1.48.0          
  [7] tidyselect_1.2.1          bit_4.0.5                
  [9] doParallel_1.0.17         clue_0.3-65              
 [11] lattice_0.22-5            rjson_0.2.21             
 [13] nor1mix_1.3-3             M3Drop_1.28.0            
 [15] blob_1.2.4                rngtools_1.5.2           
 [17] S4Arrays_1.2.0            base64_2.0.1             
 [19] scrime_1.3.5              png_0.1-8                
 [21] ResidualMatrix_1.12.0     cli_3.6.3                
 [23] askpass_1.2.0             openssl_2.1.1            
 [25] multtest_2.58.0           goftest_1.2-3            
 [27] BiocIO_1.12.0             bluster_1.12.0           
 [29] BiocNeighbors_1.20.2      densEstBayes_1.0-2.2     
 [31] uwot_0.1.16               dendextend_1.17.1        
 [33] curl_5.2.0                mime_0.12                
 [35] evaluate_0.23             leiden_0.4.3.1           
 [37] V8_6.0.0                  stringi_1.8.3            
 [39] backports_1.4.1           XML_3.99-0.16.1          
 [41] httpuv_1.6.14             magrittr_2.0.3           
 [43] rappdirs_0.3.3            splines_4.3.3            
 [45] mclust_6.1                jpeg_0.1-10              
 [47] doRNG_1.8.6               sctransform_0.4.1        
 [49] ggbeeswarm_0.7.2          DBI_1.2.1                
 [51] HDF5Array_1.30.0          genefilter_1.84.0        
 [53] jquerylib_0.1.4           withr_3.0.0              
 [55] git2r_0.33.0              rprojroot_2.0.4          
 [57] lmtest_0.9-40             bdsmatrix_1.3-6          
 [59] rtracklayer_1.62.0        BiocManager_1.30.22      
 [61] htmlwidgets_1.6.4         fs_1.6.5                 
 [63] biomaRt_2.58.2            ggrepel_0.9.5            
 [65] labeling_0.4.3            SparseArray_1.2.4        
 [67] DEoptimR_1.1-3            annotate_1.80.0          
 [69] reticulate_1.35.0         zoo_1.8-12               
 [71] knitr_1.45                beanplot_1.3.1           
 [73] timechange_0.3.0          fansi_1.0.6              
 [75] caTools_1.18.2            data.table_1.15.0        
 [77] rhdf5_2.46.1              ruv_0.9.7.1              
 [79] R.oo_1.26.0               irlba_2.3.5.1            
 [81] ellipsis_0.3.2            aroma.light_3.32.0       
 [83] lazyeval_0.2.2            yaml_2.3.8               
 [85] survival_3.7-0            scattermore_1.2          
 [87] crayon_1.5.2              RcppAnnoy_0.0.22         
 [89] RColorBrewer_1.1-3        progressr_0.14.0         
 [91] later_1.3.2               ggridges_0.5.6           
 [93] codetools_0.2-19          base64enc_0.1-3          
 [95] GlobalOptions_0.1.2       KEGGREST_1.42.0          
 [97] bbmle_1.0.25.1            Rtsne_0.17               
 [99] shape_1.4.6               startupmsg_0.9.6.1       
[101] filelock_1.0.3            foreign_0.8-86           
[103] pkgconfig_2.0.3           xml2_1.3.6               
[105] getPass_0.2-4             sfsmisc_1.1-17           
[107] spatstat.sparse_3.0-3     viridisLite_0.4.2        
[109] xtable_1.8-4              interp_1.1-6             
[111] fastcluster_1.2.6         highr_0.10               
[113] hwriter_1.3.2.1           plyr_1.8.9               
[115] httr_1.4.7                tools_4.3.3              
[117] globals_0.16.2            pkgbuild_1.4.3           
[119] beeswarm_0.4.0            htmlTable_2.4.2          
[121] checkmate_2.3.1           nlme_3.1-164             
[123] loo_2.6.0                 dbplyr_2.4.0             
[125] digest_0.6.34             numDeriv_2016.8-1.1      
[127] Matrix_1.6-5              farver_2.1.1             
[129] tzdb_0.4.0                reshape2_1.4.4           
[131] viridis_0.6.5             cvTools_0.3.2            
[133] rpart_4.1.23              cachem_1.0.8             
[135] BiocFileCache_2.10.1      polyclip_1.10-6          
[137] WGCNA_1.72-5              Hmisc_5.1-1              
[139] generics_0.1.3            proxyC_0.3.4             
[141] dynamicTreeCut_1.63-1     mvtnorm_1.2-4            
[143] parallelly_1.37.0         statmod_1.5.0            
[145] impute_1.76.0             ScaledMatrix_1.10.0      
[147] GEOquery_2.70.0           pbapply_1.7-2            
[149] dqrng_0.3.2               utf8_1.2.4               
[151] siggenes_1.76.0           StanHeaders_2.32.5       
[153] gtools_3.9.5              preprocessCore_1.64.0    
[155] gridExtra_2.3             shiny_1.8.0              
[157] GenomeInfoDbData_1.2.11   R.utils_2.12.3           
[159] rhdf5filters_1.14.1       RCurl_1.98-1.14          
[161] memoise_2.0.1             rmarkdown_2.25           
[163] scales_1.3.0              R.methodsS3_1.8.2        
[165] future_1.33.1             reshape_0.8.9            
[167] RANN_2.6.1                renv_1.0.3               
[169] Cairo_1.6-2               illuminaio_0.44.0        
[171] spatstat.data_3.0-4       rstudioapi_0.15.0        
[173] cluster_2.1.6             QuickJSR_1.1.3           
[175] whisker_0.4.1             rstantools_2.4.0         
[177] spatstat.utils_3.0-4      hms_1.1.3                
[179] fitdistrplus_1.1-11       munsell_0.5.0            
[181] cowplot_1.1.3             colorspace_2.1-0         
[183] quadprog_1.5-8            rlang_1.1.4              
[185] DelayedMatrixStats_1.24.0 sparseMatrixStats_1.14.0 
[187] circlize_0.4.15           mgcv_1.9-1               
[189] xfun_0.42                 reldist_1.7-2            
[191] abind_1.4-5               rstan_2.32.5             
[193] Rhdf5lib_1.24.2           bitops_1.0-7             
[195] ps_1.7.6                  promises_1.2.1           
[197] inline_0.3.19             RSQLite_2.3.5            
[199] DelayedArray_0.28.0       GO.db_3.18.0             
[201] compiler_4.3.3            prettyunits_1.2.0        
[203] beachmat_2.18.1           listenv_0.9.1            
[205] Rcpp_1.0.12               BiocSingular_1.18.0      
[207] tensor_1.5                MASS_7.3-60.0.1          
[209] progress_1.2.3            spatstat.random_3.2-2    
[211] R6_2.5.1                  fastmap_1.1.1            
[213] vipor_0.4.7               distr_2.9.3              
[215] ROCR_1.0-11               rsvd_1.0.5               
[217] nnet_7.3-19               gtable_0.3.4             
[219] KernSmooth_2.23-24        latticeExtra_0.6-30      
[221] miniUI_0.1.1.1            deldir_2.0-2             
[223] htmltools_0.5.8.1         RcppParallel_5.1.7       
[225] bit64_4.0.5               spatstat.explore_3.2-6   
[227] lifecycle_1.0.4           processx_3.8.3           
[229] callr_3.7.3               restfulr_0.0.15          
[231] sass_0.4.9                vctrs_0.6.5              
[233] spatstat.geom_3.2-8       robustbase_0.99-2        
[235] scran_1.30.2              sp_2.1-3                 
[237] future.apply_1.11.1       bslib_0.6.1              
[239] pillar_1.9.0              batchelor_1.18.1         
[241] prismatic_1.1.1           gplots_3.1.3.1           
[243] metapod_1.10.1            jsonlite_1.8.8           
[245] GetoptLong_1.0.5