Inflammation of Paediatric Pulmonary Diseases

DGE analysis of CF status in lipid macrophages

workflowr

• Summary
• Checks
• Past versions

Last updated: 2025-01-03

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: 6e66284

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 6e66284. 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/C133_Neeland_batch1/
    Ignored:    data/C133_Neeland_merged/
    Ignored:    output/dge_analysis/obsolete/
    Ignored:    renv/library/
    Ignored:    renv/staging/

Untracked files:
    Untracked:  analysis/14.0_DGE_analysis_T-cells.Rmd
    Untracked:  analysis/99.0_Figure_1.Rmd
    Untracked:  analysis/99.0_Figure_2.Rmd
    Untracked:  analysis/99.0_Figure_3.Rmd
    Untracked:  analysis/99.0_Figure_4.Rmd
    Untracked:  analysis/99.0_Figure_5.Rmd
    Untracked:  analysis/99.0_Supplementary_Figure_ADTs.Rmd
    Untracked:  broad_markers_seurat.csv
    Untracked:  code/background_job.R
    Untracked:  code/reverse_modifier_severity_comparisons.sh
    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
    Untracked:  data/intermediate_objects/T cells.CF_samples.fit.rds
    Untracked:  data/intermediate_objects/T cells.all_samples.fit.rds
    Untracked:  output/dge_analysis/T cells/

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
    Deleted:    analysis/14.0_proportions_analysis_ann_level_1.Rmd
    Deleted:    analysis/14.1_proportions_analysis_ann_level_3_non-macrophages.Rmd
    Deleted:    analysis/14.2_proportions_analysis_ann_level_3_macrophages.Rmd
    Deleted:    analysis/15.0_Figure_1.Rmd
    Deleted:    analysis/16.0_Figure_2.Rmd
    Deleted:    analysis/17.0_Supplementary_Figure_ADTs.Rmd
    Modified:   analysis/index.Rmd
    Modified:   code/utility.R
    Modified:   data/cluster_annotations/marker_proteins_TNK_supp.xlsx
    Modified:   data/cluster_annotations/marker_proteins_macrophages_supp.xlsx
    Modified:   data/cluster_annotations/marker_proteins_other_supp.xlsx
    Modified:   data/cluster_annotations/seurat_markers_all_cells.rds
    Modified:   data/intermediate_objects/macro-APOC2+.CF_samples.fit.rds
    Modified:   data/intermediate_objects/macro-APOC2+.all_samples.fit.rds
    Modified:   data/intermediate_objects/macro-CCL.CF_samples.fit.rds
    Modified:   data/intermediate_objects/macro-CCL.all_samples.fit.rds
    Modified:   data/intermediate_objects/macro-IFI27.CF_samples.fit.rds
    Modified:   data/intermediate_objects/macro-IFI27.all_samples.fit.rds
    Modified:   data/intermediate_objects/macro-alveolar.CF_samples.fit.rds
    Modified:   data/intermediate_objects/macro-alveolar.all_samples.fit.rds
    Modified:   data/intermediate_objects/macrophages.CF_samples.fit.rds
    Modified:   data/intermediate_objects/macrophages.all_samples.fit.rds
    Modified:   output/dge_analysis/macro-APOC2+/ORA.FIBROSIS.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.FIBROSIS.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.FIBROSIS.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.GO.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.GO.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.GO.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.HALLMARK.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.HALLMARK.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.HALLMARK.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.REACTOME.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.REACTOME.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.REACTOME.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.WP.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.WP.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-APOC2+/ORA.WP.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.FIBROSIS.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.FIBROSIS.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.FIBROSIS.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.FIBROSIS.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.GO.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.GO.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.GO.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.GO.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.HALLMARK.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.HALLMARK.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.HALLMARK.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.HALLMARK.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.REACTOME.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.REACTOME.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.REACTOME.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.REACTOME.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.WP.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.WP.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.WP.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-CCL/ORA.WP.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.FIBROSIS.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.FIBROSIS.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.FIBROSIS.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.FIBROSIS.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.GO.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.GO.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.GO.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.GO.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.HALLMARK.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.HALLMARK.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.HALLMARK.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.HALLMARK.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.REACTOME.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.REACTOME.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.REACTOME.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.REACTOME.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.WP.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.WP.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.WP.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-IFI27/ORA.WP.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.FIBROSIS.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.FIBROSIS.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.FIBROSIS.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.FIBROSIS.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.GO.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.GO.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.GO.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.GO.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.HALLMARK.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.HALLMARK.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.HALLMARK.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.HALLMARK.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.REACTOME.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.REACTOME.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.REACTOME.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.REACTOME.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.WP.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.WP.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.WP.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macro-alveolar/ORA.WP.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.FIBROSIS.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.FIBROSIS.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.FIBROSIS.CF.LUMA_IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.FIBROSIS.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macrophages/ORA.FIBROSIS.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.GO.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.GO.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.GO.CF.LUMA_IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.GO.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macrophages/ORA.GO.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.HALLMARK.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.HALLMARK.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.HALLMARK.CF.LUMA_IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.HALLMARK.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macrophages/ORA.HALLMARK.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.REACTOME.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.REACTOME.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.REACTOME.CF.LUMA_IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.REACTOME.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macrophages/ORA.REACTOME.CF.NO_MODvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.WP.CF.IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.WP.CF.IVAvNON_CF.CTRL.csv
    Modified:   output/dge_analysis/macrophages/ORA.WP.CF.LUMA_IVAvCF.NO_MOD.csv
    Modified:   output/dge_analysis/macrophages/ORA.WP.CF.NO_MOD.SvCF.NO_MOD.M.csv
    Modified:   output/dge_analysis/macrophages/ORA.WP.CF.NO_MODvNON_CF.CTRL.csv

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.5_DGE_analysis_macro-lipid.Rmd) and HTML (docs/13.5_DGE_analysis_macro-lipid.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
Rmd	6e66284	Jovana Maksimovic	2025-01-03	wflow_publish("analysis/13.5_DGE_analysis_macro-lipid.Rmd")
html	e7ac317	Jovana Maksimovic	2024-12-06	Build site.
Rmd	900a349	Jovana Maksimovic	2024-12-06	wflow_publish("analysis/13.5_DGE_analysis_macro-lipid.Rmd")
html	3e98f1b	Jovana Maksimovic	2024-12-05	Build site.
Rmd	20fa04b	Jovana Maksimovic	2024-12-05	wflow_publish("analysis/13.5_DGE_analysis_macro-lipid.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.

# select the cell type for pseudobulks
cell <- "macro-lipid"

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

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

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] 21056    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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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

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

[1] 21056    41

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
3e98f1b	Jovana Maksimovic	2024-12-05

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

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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 <- 1
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
3e98f1b	Jovana Maksimovic	2024-12-05

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

keep
FALSE  TRUE 
 5319 10753 

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

[1] 10753    41

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
3e98f1b	Jovana Maksimovic	2024-12-05

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 
 7301  3452 

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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 = 4)
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	sexM	age	sample_34	sample_35	sample_36	sample_37	sample_38	sample_39
0	0	0	0	0	0	1	-0.2869949	0.1079046	0.0211944	0.4329522	1	-0.2590872	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.1574322	0.0027149	0.0320608	-0.0019883	1	-0.0939001	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.2441917	0.0333598	-0.1007447	-0.4614589	0	-0.1151479	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.0300439	-0.0659044	-0.0227775	-0.0587969	0	-0.0441471	0	0	0	0	0	0
0	0	0	0	1	0	0	0.2376656	-0.1959905	0.0665754	-0.0185020	1	0.1428834	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.2372713	0.0586405	-0.0269206	0.1733991	0	-0.0729608	0	0	0	0	0	0
0	0	0	0	0	0	1	-0.0119293	-0.0888899	-0.0484010	-0.1320510	1	0.1464588	0	0	0	0	0	0
0	0	0	0	0	1	0	-0.1580288	0.0138009	0.0617520	0.0944636	1	0.5597097	0	0	0	0	0	0
0	0	0	0	0	1	0	0.0273550	-0.1069516	0.0411273	-0.0667953	0	1.5743836	0	0	0	0	0	0
1	0	0	0	0	0	0	0.1995619	-0.1901310	0.0365145	-0.1204073	1	1.5993830	0	0	0	0	0	0
1	0	0	0	0	0	0	0.1320080	-0.1490622	0.1597634	0.0307166	1	2.3883594	0	0	0	0	0	0
0	0	0	0	0	1	0	0.2239199	-0.0501041	-0.3713149	-0.0034687	0	2.2957230	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.0055859	0.0554640	-0.2961729	0.0917326	1	2.3360877	0	0	0	0	0	0
1	0	0	0	0	0	0	0.1328905	-0.0145603	-0.2891798	0.0761856	1	2.2980155	0	0	0	0	0	0
0	0	0	0	1	0	0	0.1550195	-0.1311763	0.2012728	0.1654682	0	2.5790214	0	0	0	0	0	0
0	0	0	0	0	1	0	0.0902573	0.0044708	-0.2922606	0.0849805	0	2.5823250	0	0	0	0	0	0
0	0	0	0	0	0	1	-0.2048293	0.0224167	-0.0533283	-0.0599834	1	0.1321035	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.0003229	0.0377773	-0.3689912	-0.0252052	0	2.5583683	0	0	0	0	0	0
0	0	0	0	1	0	0	0.1468431	-0.0226471	-0.3271941	-0.0019623	0	2.5670653	0	0	0	0	0	0
0	1	0	0	0	0	0	-0.0986092	0.0874085	-0.2417670	0.2801207	1	2.5730557	0	0	0	0	0	0
0	0	0	0	1	0	0	-0.0638331	-0.0636796	0.1377969	0.0002638	0	-0.9343238	1	0	0	0	0	0
0	0	0	0	1	0	0	0.0025773	-0.1143108	0.1001281	-0.0959259	0	1.0409164	1	0	0	0	0	0
0	0	0	0	1	0	0	-0.0685978	-0.0881100	0.0974809	-0.1410045	1	0.0807044	0	1	0	0	0	0
0	0	0	0	1	0	0	0.1043814	-0.1536691	0.1546282	-0.0257335	1	0.9940589	0	1	0	0	0	0
0	0	0	0	0	1	0	-0.1920921	0.0148940	0.1200105	0.0100200	0	-0.0564254	0	0	1	0	0	0
0	0	0	1	0	0	0	-0.1792726	0.0085209	0.0905782	-0.0530025	0	1.1764977	0	0	1	0	0	0
0	0	0	0	1	0	0	0.0669056	-0.1130240	0.0436639	-0.0209908	0	1.5597097	0	0	0	1	0	0
0	0	1	0	0	0	0	-0.1072536	-0.0307969	0.0094561	-0.0985189	0	2.1930156	0	0	0	1	0	0
0	0	1	0	0	0	0	0.0126531	-0.0728746	0.0546795	0.0508093	0	2.2980155	0	0	0	1	0	0
1	0	0	0	0	0	0	-0.1740434	0.0312695	0.0579085	0.1661045	1	1.5703964	0	0	0	0	1	0
1	0	0	0	0	0	0	-0.0933254	-0.0036180	0.0795214	0.1322104	1	2.0206033	0	0	0	0	1	0
1	0	0	0	0	0	0	-0.0151928	-0.0371099	0.1583696	0.2516214	1	2.3485584	0	0	0	0	1	0
0	0	0	0	1	0	0	0.1488699	-0.1271052	0.0955022	0.0918902	0	1.9730702	0	0	0	0	0	1
0	0	1	0	0	0	0	0.1029548	-0.1256285	0.0483399	0.0002565	0	2.6297159	0	0	0	0	0	1
0	0	0	0	0	0	1	-0.0931588	-0.0614826	0.0796403	-0.0484220	1	0.2923784	0	0	0	0	0	0
0	0	0	0	0	1	0	0.0264424	0.4191125	0.0971220	-0.1403524	1	1.5801455	0	0	0	0	0	0
0	0	0	0	1	0	0	0.2135632	0.3958041	0.1115306	-0.0009611	1	1.5801455	0	0	0	0	0	0
0	1	0	0	0	0	0	0.3582889	0.3616107	0.1226946	-0.1003081	1	1.5993178	0	0	0	0	0	0
0	0	0	0	0	0	1	0.1321226	0.4833822	0.1378039	-0.0026873	1	1.5849625	0	0	0	0	0	0
0	0	0	0	0	0	1	0.1554134	-0.1625755	0.0802817	-0.0149088	1	2.4204621	0	0	0	0	0	0
0	0	0	0	0	0	1	-0.2476845	0.0308503	-0.0583451	-0.4397603	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
3e98f1b	Jovana Maksimovic	2024-12-05

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
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
3e98f1b	Jovana Maksimovic	2024-12-05

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
e7ac317	Jovana Maksimovic	2024-12-06
3e98f1b	Jovana Maksimovic	2024-12-05

[[2]]

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

[[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
3e98f1b	Jovana Maksimovic	2024-12-05

[[2]]

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

[[3]]
NULL

CF modifiers and severity

Prepare data

Filter genes

Extract only the CF samples.

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

[1] 21056    34

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 <- 1
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
3e98f1b	Jovana Maksimovic	2024-12-05

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

keep
FALSE  TRUE 
 5359 10713 

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

[1] 10713    34

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
3e98f1b	Jovana Maksimovic	2024-12-05

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 
 7261  3452 

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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 = 4)
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	sexM	microTRUE	age	sample_34	sample_35	sample_36	sample_37	sample_38	sample_39
0	0	0	0	1	0	-0.2000315	0.0104253	-0.0318121	0.0600555	1	0	-0.0939001	0	0	0	0	0	0
0	0	0	0	1	0	-0.2999245	0.0728900	0.0982552	-0.7654051	0	0	-0.1151479	0	0	0	0	0	0
0	0	0	0	1	0	-0.0533688	-0.0604991	0.0315079	-0.0604420	0	0	-0.0441471	0	0	0	0	0	0
0	0	0	0	1	0	0.2546575	-0.2266961	-0.0413403	-0.0508797	1	0	0.1428834	0	0	0	0	0	0
0	0	0	0	1	0	-0.2920805	0.0840707	0.0251430	0.3090124	0	0	-0.0729608	0	0	0	0	0	0
0	0	0	0	0	1	-0.2007192	0.0190746	-0.0630893	0.0689074	1	1	0.5597097	0	0	0	0	0	0
0	0	0	0	0	1	0.0126649	-0.1176782	-0.0287696	-0.0986565	0	1	1.5743836	0	0	0	0	0	0
1	0	0	0	0	0	0.2108749	-0.2125253	-0.0137910	-0.1286831	1	1	1.5993830	0	0	0	0	0	0
1	0	0	0	0	0	0.1329665	-0.1844703	-0.1404402	0.0125352	1	0	2.3883594	0	0	0	0	0	0
0	0	0	0	0	1	0.2390551	0.0026236	0.3689327	-0.0102930	0	0	2.2957230	0	0	0	0	0	0
0	0	0	0	1	0	-0.0250941	0.1130877	0.2842331	0.0867431	1	1	2.3360877	0	0	0	0	0	0
1	0	0	0	0	0	0.1342905	0.0271406	0.2873718	0.0653480	1	0	2.2980155	0	0	0	0	0	0
0	0	0	0	1	0	0.1593843	-0.1724758	-0.1820020	0.1158291	0	1	2.5790214	0	0	0	0	0	0
0	0	0	0	0	1	0.0851947	0.0533591	0.2863650	0.0572022	0	1	2.5823250	0	0	0	0	0	0
0	0	0	0	1	0	-0.0189786	0.1050938	0.3571792	-0.0187148	0	0	2.5583683	0	0	0	0	0	0
0	0	0	0	1	0	0.1503548	0.0282193	0.3222842	-0.0154325	0	0	2.5670653	0	0	0	0	0	0
0	1	0	0	0	0	-0.1323436	0.1486962	0.2287088	0.2753585	1	1	2.5730557	0	0	0	0	0	0
0	0	0	0	1	0	-0.0924224	-0.0818700	-0.1279539	0.0708524	0	0	-0.9343238	1	0	0	0	0	0
0	0	0	0	1	0	-0.0159244	-0.1346512	-0.0838970	-0.0760771	0	0	1.0409164	1	0	0	0	0	0
0	0	0	0	1	0	-0.0977788	-0.1064554	-0.0817957	-0.1332542	1	0	0.0807044	0	1	0	0	0	0
0	0	0	0	1	0	0.1012129	-0.1902156	-0.1338976	-0.0181000	1	1	0.9940589	0	1	0	0	0	0
0	0	0	0	0	1	-0.2400505	0.0126871	-0.1211613	0.0500265	0	0	-0.0564254	0	0	1	0	0	0
0	0	0	1	0	0	-0.2252462	0.0105548	-0.0909726	0.0073038	0	1	1.1764977	0	0	1	0	0	0
0	0	0	0	1	0	0.0581223	-0.1265232	-0.0288963	-0.0698584	0	1	1.5597097	0	0	0	1	0	0
0	0	1	0	0	0	-0.1422923	-0.0247105	-0.0044763	-0.1512977	0	0	2.1930156	0	0	0	1	0	0
0	0	1	0	0	0	-0.0043513	-0.0819305	-0.0459170	0.0108362	0	1	2.2980155	0	0	0	1	0	0
1	0	0	0	0	0	-0.2192519	0.0387505	-0.0589824	0.1794108	1	0	1.5703964	0	0	0	0	1	0
1	0	0	0	0	0	-0.1263095	-0.0035827	-0.0807965	0.1102489	1	1	2.0206033	0	0	0	0	1	0
1	0	0	0	0	0	-0.0364526	-0.0602894	-0.1510529	0.2164292	1	0	2.3485584	0	0	0	0	1	0
0	0	0	0	1	0	0.1524148	-0.1534962	-0.0788919	0.0733766	0	1	1.9730702	0	0	0	0	0	1
0	0	1	0	0	0	0.0996306	-0.1416831	-0.0321267	-0.0094323	0	1	2.6297159	0	0	0	0	0	1
0	0	0	0	0	1	0.0115497	0.4948664	-0.2134203	-0.1088818	1	0	1.5801455	0	0	0	0	0	0
0	0	0	0	1	0	0.2268472	0.4523224	-0.2221178	0.0136405	1	0	1.5801455	0	0	0	0	0	0
0	1	0	0	0	0	0.3934001	0.4058906	-0.2323802	-0.0677081	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
3e98f1b	Jovana Maksimovic	2024-12-05

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
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
3e98f1b	Jovana Maksimovic	2024-12-05

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
e7ac317	Jovana Maksimovic	2024-12-06
3e98f1b	Jovana Maksimovic	2024-12-05

[[2]]
NULL

[[3]]

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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
3e98f1b	Jovana Maksimovic	2024-12-05

[[2]]
NULL

[[3]]

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

Version	Author	Date
3e98f1b	Jovana Maksimovic	2024-12-05

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