Differential expression between groups using pseudobulk
Sarah Williams
Last updated: 2024-07-05
Checks: 7 0
Knit directory: spatialsnippets/
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.
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.
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.
The command set.seed(20231017) 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.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
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 8580064. 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/glossary.nb.html
Ignored: renv/library/
Ignored: renv/staging/
Unstaged changes:
Modified: analysis/d_cosmxIBD.Rmd
Modified: analysis/index.Rmd
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/e_DEPseudobulk_insitu.Rmd)
and HTML (docs/e_DEPseudobulk_insitu.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 | 8580064 | swbioinf | 2024-07-05 | Build site. |
| Rmd | 1706732 | swbioinf | 2024-07-05 | wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd") |
| html | b2fc2b0 | swbioinf | 2024-06-12 | Build site. |
| Rmd | 74dbb4e | swbioinf | 2024-06-12 | wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd") |
| html | 9cd7910 | swbioinf | 2024-06-11 | Build site. |
| Rmd | 561981e | swbioinf | 2024-06-11 | wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd") |
| html | 61d4499 | swbioinf | 2024-05-31 | Build site. |
| Rmd | cf7833a | swbioinf | 2024-05-31 | wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd") |
| html | 6978154 | swbioinf | 2024-05-30 | Build site. |
| html | 97393ec | swbioinf | 2024-05-28 | Build site. |
| Rmd | dd0f089 | swbioinf | 2024-05-28 | wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd") |
| html | d6d841a | swbioinf | 2024-05-22 | Build site. |
| Rmd | 9c53965 | swbioinf | 2024-05-22 | wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd") |
| html | 58ee2ec | swbioinf | 2024-05-22 | Build site. |
| Rmd | 36dd228 | swbioinf | 2024-05-22 | wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd") |
| html | 5a9d7e9 | swbioinf | 2024-05-16 | Build site. |
| Rmd | c67ef86 | swbioinf | 2024-05-16 | wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd") |
| html | c77f76c | swbioinf | 2024-05-07 | Build site. |
| Rmd | 7b0be93 | swbioinf | 2024-05-07 | wflow_publish(c("analysis/index.Rmd", "analysis/e_DEPseudobulk_insitu.Rmd", |
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6096346/
Overview
Once we have identified cell types present in the samples, its common to test how gene expression changes between experimental conditions, within each different cell type.
Some cell types may be dramatically affected by the experimental conditions, while others are not. Likewise some genes may change only in a specific cell type, whereas others show a more general difference.
This document describes how to apply a pseudobulk approach to test for differences between groups. In a pseudobulk approach counts are obtained by pooling together groups of cells; in this case cells from the of the same celltype from the same fov. These pooled counts can then be analysed more like a bulk RNAseq experiment.
This is very similar to how a non-spatial single cell experiment may be analysed.
However, unlike a non-spatial single cell experiment, this experimental design has multiple FOVs from the same tissue sample; these pseudoreplicates are informative but cannot be considered true biological replicates. They can be incorporated into the statistics with the correlate duplicates approach (Smyth, Michaud, and Scott 2005).
Note that there are many other approaches to calculate differential expression in this kind of data - including those that make use of individual cells; see review (Soneson and Robinson 2018).
This test requires:
- Biological replicates for each group
- Assigned cell types
- [Optionally] Multiple fovs measured per sample
For example:
- What genes are differentially expressed in epithelial cells in Crohn’s disease vs healthy individuals?
- How do genes change with treatment in each different cell type in my sample?
Steps:
- Calculate pseudobulk
- Filter to testable pseudobulk groups (enough cells to pool)
- Filter to testable genes (enough expression to see changes)
- Test for changes in gene expression
- Plot DE results and individual genes.
Worked example
How does gene expression change within each cell type between Ulcerative colitis or Crohn’s disease, and Healthy controls?
Using data from Macrophage and neutrophil heterogeneity at single-cell spatial resolution in human inflammatory bowel disease from (Garrido-Trigo et al. 2023);
The study included 9 cosmx slides of colonic biopsies
- 3x HC - Healthy controls
- 3x UC - Ulcerative colitis
- 3x CD - Chrones’s disease
Load libraries and data
library(Seurat)
library(speckle)
library(tidyverse)
library(limma)
library(DT)
library(edgeR)data_dir <- file.path("~/projects/spatialsnippets/datasets/GSE234713_IBDcosmx_GarridoTrigo2023/processed_data")
seurat_file_01_loaded <- file.path(data_dir, "GSE234713_CosMx_IBD_seurat_01_loaded.RDS")so <- readRDS(seurat_file_01_loaded)Experimental design
There are three individuals per condition, one tissue sample from each individual. 9 slides in total. Each tissue sample has multiple FOVs captured; these are our pseudoreplicates. Each ‘sample’ is listed below.
NB: ‘FOV’ Feild Of View: In the seurat package, an number of functions include an ‘fov’ parameter. This corresponds to the slide on which one or more samples are present on. In the cosMX output an ‘fov’ refers to the rectangular regions on the slide that are measured; there are multiple fovs per slide. This document uses the latter.
sample_table <- select(as_tibble(so@meta.data), condition, individual_code, fov_name) %>%
unique() %>%
group_by(condition, individual_code) %>%
summarise(n_fovs= n(), item = str_c(fov_name, collapse = ", "))
DT::datatable(sample_table)Cell filtering and counting
We need to define the grouping that will make one pseudobulk sample for analysis. In this experiment, we can consider all FOVs within a tissue sample as somewhat equivalent. Therefore that is all cells of a particular celltype within an FOV region.
If we only had one fov per sample, celltype + sample would be directly equivalent.
If we had a tiled fovs across a tissue sample with distinct regions it would be necessary to group them by subregion fov + celltype. Alternatively pool the entire sample if each one is equivalent - it depends on where transcriptional differences are expected.
Once the pseudobulk replicate grouping is decided, we need to check there will be enough cells to perform our analyses. Note there are several levels of filtering here!
- Need at least x reads in a cell to include it
- Need at least x cells of a celltype within an fov to include a sample
- Can only test where we have at least 2 samples on each side of a contrast.
What these thresholds should be needs to be determined for each experiment.
Filter counts per cell
To start, check out the distribution of reads per cell. Here, we choose and apply a filter of 200 reads. This is low compared to what you might see in a single scRNAseq experiment, but we are working with only 1000 genes.
min_reads_per_cell <- 200
ggplot(so@meta.data, aes(x=nCount_RNA)) +
geom_density() +
geom_vline(xintercept = min_reads_per_cell, lty=3) +
scale_x_log10() +
theme_bw()+
ggtitle("How many reads per cell?")
| Version | Author | Date |
|---|---|---|
| c77f76c | swbioinf | 2024-05-07 |
so<- so[,so$nCount_RNA >= min_reads_per_cell]Filter cells per group
Next we pool each celltype within each fov (naming those groupings fov_cluster). But there needs to be a certain number of cells for that to work - less than a certain number of cells and a pseudobulk pool will be excluded. The table below shows the cells per grouping.
Note there are much fewer t-cells overall, but given that we have a high number of samples, there should still be enough to include. Its typical that some of the less common cell types are difficult or impossible to reliably test.
min_cells_per_fovcluster <- 20
so$fov_cluster <- paste0(so$fov_name,"_", so$celltype_subset)
celltype_summary_table <- so@meta.data %>%
group_by(condition, group, individual_code, fov_name, celltype_subset, fov_cluster) %>%
summarise(cells=n(), .groups = 'drop')
DT::datatable(celltype_summary_table)ggplot(celltype_summary_table, aes(x=cells, col=celltype_subset)) +
geom_density() +
geom_vline(xintercept=min_cells_per_fovcluster, lty=3) +
geom_rug() +
scale_x_log10() +
theme_bw() +
ggtitle("How many cells per fov-cluster?")
Record the names of those fov_clusters that contain enough cells to be used. Will use this later to filter.
celltype_summary_table.passed <- celltype_summary_table[celltype_summary_table$cells >= min_cells_per_fovcluster,]
passed_fov_clusters <- celltype_summary_table.passed$fov_clusterSamples per contrast.
Ignoring pseodureplicate fovs from a single sample, how many bioloical samples are represented in our filtered table?
In this experiment, there are enough samples to run the contrasts for every cell type.
celltype_summary_table.passed %>% select(condition, individual_code, celltype_subset) %>%
unique() %>%
group_by(condition, celltype_subset) %>%
summarise(n=n()) %>%
pivot_wider( names_from=condition, values_from = n)# A tibble: 5 × 4
celltype_subset `Crohn's disease` `Ulcerative colitis` `Healthy controls`
<fct> <int> <int> <int>
1 epi 3 3 3
2 myeloids 3 3 3
3 plasmas 3 3 3
4 stroma 3 3 3
5 tcells 3 2 2NB: What if you do need to skip constrasts? Often its easiest to count your biological replicates at the differential expression step and skip them there. You might need to run contrasts between groups one at a time to do this.
Calculate pseudobulk
Now use the PseudobulkExpression() function to sum up each gene’s gene expression across each fov_cluster.
pseudobulk_counts <- PseudobulkExpression(so, assays = "RNA", layer="counts", method = 'aggregate', group.by = 'fov_cluster')
pseudobulk_counts_matrix <- pseudobulk_counts[["RNA"]]
# Change - back to _. Ideally we'd have neither, but - will cause problems later
colnames(pseudobulk_counts_matrix)<-gsub("-","_",colnames(pseudobulk_counts_matrix))Now instead of counts for each individual cell, we have our pseodubulk matrix with the pooled sum of counts for each celltype within each fov region. So the numbers are much higher, with fewer zeros. But note that each pool is a different size, so we can’t compare counts directly. We will address this later with normalisation.
pseudobulk_counts_matrix[1:10,1:4]10 x 4 sparse Matrix of class "dgCMatrix"
CD_a_001_epi CD_a_001_myeloids CD_a_001_plasmas CD_a_001_stroma
AATK 38 4 16 13
ABL1 65 1 6 18
ABL2 33 2 11 10
ACE 28 . 4 6
ACE2 59 . 17 18
ACKR1 16 1 13 21
ACKR3 37 4 12 15
ACKR4 36 2 9 7
ACTA2 23 1 7 45
ACTG2 32 2 15 27Next, filter to only those passed fov_clusters we saw earlier - this filtered table is what we’ll use for calculating differential expression.
We also need an annotation table to tell us what is in each fov_cluster (sample, celltype, condition…). Build this by filtering the celltype summary table made earlier.
For convenience we filter it to the same samples, in the same order as the pseudobulk matrix, so we can then use those columns of data to build our differential expression model in the next step. We need to be careful the order remains the same.
# Filter the pseudobulk table
pseudobulk_counts_matrix <- pseudobulk_counts_matrix[,passed_fov_clusters]
# And filter the celltype summary into an annotation table for only those fov_clusters
# pull in relevant annotation in a matched order
pseudobulk_anno_table <- celltype_summary_table
match_order <- match(passed_fov_clusters, pseudobulk_anno_table$fov_cluster)
pseudobulk_anno_table <- pseudobulk_anno_table[match_order,]
# Double check that the order of samples in the annotation table matches the pseudobulk table.
# If this is wrong then the results will be nonsense!
stopifnot(all(colnames(pseudobulk_counts_matrix) == pseudobulk_anno_table$fov_cluster ))Calculate Differential Expression
We now have a counts matrix, and an annotation table that describes its samples.
This looks very much like a bulk RNAseq experiment. Except - instead of X samples across Y conditions, we have X samples of Y conditions for each of Z celltypes.
pseudobulk_counts_matrix[1:10,1:4]10 x 4 sparse Matrix of class "dgCMatrix"
CD_a_001_epi CD_a_001_plasmas CD_a_001_stroma CD_a_002_myeloids
AATK 38 16 13 8
ABL1 65 6 18 3
ABL2 33 11 10 11
ACE 28 4 6 3
ACE2 59 17 18 4
ACKR1 16 13 21 4
ACKR3 37 12 15 2
ACKR4 36 9 7 3
ACTA2 23 7 45 6
ACTG2 32 15 27 3head(pseudobulk_anno_table)# A tibble: 6 × 7
condition group individual_code fov_name celltype_subset fov_cluster cells
<fct> <fct> <chr> <fct> <fct> <chr> <int>
1 Crohn's dise… CD CD_a CD_a_001 epi CD_a_001_e… 279
2 Crohn's dise… CD CD_a CD_a_001 plasmas CD_a_001_p… 99
3 Crohn's dise… CD CD_a CD_a_001 stroma CD_a_001_s… 78
4 Crohn's dise… CD CD_a CD_a_002 myeloids CD_a_002_m… 28
5 Crohn's dise… CD CD_a CD_a_002 stroma CD_a_002_s… 230
6 Crohn's dise… CD CD_a CD_a_002 tcells CD_a_002_t… 21We will use a fairly standard limma differential expression analysis.
We will process each celltype one at a time as follows:
- Subset the counts matrix and annotation table to the celltype.
- Build and fit the model. In our case, a simple one that looks for various by disease group, blocked on individual.
- Run the test.
How to build the model and run the test will need to be customised for your experimental design. Online resources that describe approaches for bulk RNAseq analyses can be applied to pseudobulk analyses - see the ‘More Information’ section of this document for suggestions.
In this case the line model.matrix( ~0 + group) builds a
model where gene expression is expected to vary by group (UC/CD/HC). The
~0 part indicates a ‘intercept’, which is convenient because it means we
can include all 3 groups in our contrasts, rather than one of them being
treated as the baseline. ( NB: Mathematically, you could use ~group
alone without intercept, and define contrasts accordingly to get the
same result. )
But, we also need to take into account our individuals. We do this by blocking on individual, and calculating a intra duplicate correlation with duplicateCorrelation()(Smyth, Michaud, and Scott 2005). That’s a method from microarray times, see discussion for how it is applied to bulk RNAseq data.
If instead we didn’t have multiple fovs per replicate, can skip blocking or use of duplicateCorrelation(). That is more typical of a typical of a bulk rnaseq or single cell RNAseq pseuobulk analysis.
Some example models:
- ~0 + group, (blocked on individual with using duplicate Correlation) : Used here, takes into account pseodureplicate fovs from the same samples
- ~0 + group, (no blocking or duplicate correlation) : If there was only one pseudoreplicate per individual. Either from a single fov, or from pooling fovs.
- ~0 + group + individual (no blocking or duplicate correlation) : If this was a treatment / time point experiment with paired samples from individuals. But only one pseodureplicate per individual. If there were multiple fovs per tissue sample, you could add blocking and duplicate correlation on tissue ID.
- ~0 + group + celltype (with or without blocking + duplicate correlation) : Instead of subsetting to celltype, includes cell type in the model. This pulls more data into the stataticacs, but makes building the contrasts challenging! Not covered here.
# Empty list to collect results
de_result_list <- list()
# celltype_subset is a matrix
for (the_celltype in levels(so$celltype_subset)) {
# Subset counts andn annotation to one cell type.
# Ensure order remains identical!
print(the_celltype)
anno_table.this <- pseudobulk_anno_table[pseudobulk_anno_table$celltype_subset == the_celltype,]
count_matrix.this <- pseudobulk_counts_matrix[,anno_table.this$fov_cluster]
## Check for sufficient replicates ##
# To do any calculations, we need at least 2 pseudobulk groups per contrast.
# there are plenty in this experiemnt, but with less replicates and rare cell types
# it may be neccesary to check and skip certain contrasts. Here woudl be a good
# if (not enouch samples to run test ) {next}
# skip clusters with no samples after filtering
if( nrow(anno_table.this) < 1 ) {next}
# Setup objects for limma
dge <- DGEList(count_matrix.this)
dge <- calcNormFactors(dge)
# Build model
group <- anno_table.this$group
individual_code <- anno_table.this$individual_code
# Model design
design <- model.matrix( ~0 + group)
# Run Voom
vm <- voom(dge, design = design, plot = FALSE)
# Adding dupliate correlation to use individual fovs, rather than pooled per biosample
corrfit <- duplicateCorrelation(vm, design, block=individual_code)
fit <- lmFit(vm, design, correlation = corrfit$consensus, block=individual_code)
# Define and fit contrasts and run ebayes
contrasts <- makeContrasts(UCvHC = groupUC - groupHC,
CDvHC = groupCD - groupHC,
levels=coef(fit))
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)
## Look through each contrast, and extract a results table.
for ( the_coef in colnames(contrasts) ) {
de_result.this <- topTable(fit, n = Inf, adjust.method = "BH", coef = the_coef) %>%
rownames_to_column("target") %>%
mutate(contrast=the_coef,
celltype=the_celltype) %>%
select(celltype,contrast,target,everything()) %>%
arrange(P.Value)
de_result_list[[paste(the_celltype, the_coef, sep="_")]] <- de_result.this
}
}[1] "epi"
[1] "myeloids"
[1] "plasmas"
[1] "stroma"
[1] "tcells"# Join together results for all celltypes, and pull out those with a singificant adjusted p-value
de_results_all <- bind_rows(de_result_list)
de_results_sig <- filter(de_results_all, adj.P.Val < 0.01)Table of significant results.
DT::datatable(mutate(de_results_sig, across(is.numeric, signif, digits = 3)))DE plots
The below plots show the logFC calculate for each gene versus its average expression across all samples. This is a useful diagnostic plot to evaluate your differential expression results.
- Only big changes are significant at lower expressions, simply becuase higher expression means more statistical confidence.
- A lopsided plot might indicate some oddness around the normalisation of your data (e.g. extremely different cell counts).
[[ADD pvalue plot?]]
library(ggrepel) # gg_repel, For non-overlapping gene labels
make_ma_style_plot <- function(res_table, pval_threshold = 0.01, n_genes_to_label = 10) {
p <- ggplot(res_table, aes(x=AveExpr, y=logFC, col=adj.P.Val < pval_threshold) ) +
geom_hline(yintercept = c(0), col='grey80') +
geom_point(pch=3) +
geom_text_repel(data = head(arrange(filter(res_table , adj.P.Val < pval_threshold ), P.Value), n=5),
mapping = aes(label=target), col="red" ) +
theme_bw() +
geom_hline(yintercept = c(-1,1), lty=3) +
scale_colour_manual(values = c('FALSE'="black", 'TRUE'="red")) +
theme(legend.position = 'none')
return(p)
}#res_table.UCvHC.epi <- filter(de_results_all, contrast == "UCvHC", celltype=="epi")
p1 <- make_ma_style_plot(res_table = filter(de_results_all, contrast == "UCvHC", celltype=="epi")) + ggtitle("UC vs HC - epi")
p2 <- make_ma_style_plot(res_table = filter(de_results_all, contrast == "UCvHC", celltype=="tcells"))+ ggtitle("UC vs HC - T-cells")
p3 <- make_ma_style_plot(res_table = filter(de_results_all, contrast == "UCvHC", celltype=="stroma")) + ggtitle("UC vs HC - stroma")
p1 + p2 + p3
Check some examples
Its always worth visualising how the expression of your differentially expressed genes really looks, with respect to your experimental design. How best to do this depends on your experiment.
The results suggests that TNFRSF18 was significantly DE between individuals with Ulcerative Colitis and Healthy Controls in plasma cells. As a first step, there’s some very convenient seurat plots below;
p1 <- VlnPlot(subset(so, celltype_subset == "plasmas"), features = "TNFRSF18", group.by = 'group', alpha = 0.1)
p2 <- FeaturePlot(so, "TNFRSF18", split.by = "group")
p1 / p2
But it gets difficult to summarise data at the single cell level. We can also use the the normalised pseudobulk expression to see how gene expression varies within each fov,individual,celltype and condition - The plot below shows an overview of normalised TNFRSF18 expression across the entire experiment.
# Get tmm normalised coutns for all pseudobulk
# WHen we did the DE we calculated this a celltype at a time, so values might differ slightly!
dge <- DGEList(pseudobulk_counts_matrix)
dge <- calcNormFactors(dge)
norm_pseudobulk <- cpm(dge , log=TRUE) # uses tmm normalisation
# Plot expression for TNFRSF18
plottable <- cbind(pseudobulk_anno_table, expression = norm_pseudobulk["TNFRSF18",])
ggplot(plottable, aes(x=individual_code, y=expression, col=condition )) +
geom_boxplot(outlier.shape = NA) +
geom_point() +
theme_bw() +
theme(axis.text.x=element_text(angle = -90, hjust = 0)) +
facet_wrap(~celltype_subset) 
Compared with IGHG1, which was flagged at differentially expressed across multiple cell types.
plottable <- cbind(pseudobulk_anno_table, expression = norm_pseudobulk["IGHG1",])
ggplot(plottable, aes(x=individual_code, y=expression, col=condition )) +
geom_boxplot(outlier.shape = NA) +
geom_point() +
theme_bw() +
theme(axis.text.x=element_text(angle = -90, hjust = 0)) +
facet_wrap(~celltype_subset) 
| Version | Author | Date |
|---|---|---|
| 8580064 | swbioinf | 2024-07-05 |
Code Snippet
library(Seurat)
library(edgeR)
library(limma)
# Set threhoehsolds
min_reads_per_cell <- 200
min_cells_per_fovcluster <- 20
# Remove cells with too few counts
so <- so[,so$nCount_RNA >= min_reads_per_cell]
# Define pseodoreplicate groups, with all relevant sample annotation
# remove those with too few cells.
so$fov_cluster <- paste0(so$fov_name,"_", so$celltype_subset)
celltype_summary_table <- so@meta.data %>%
group_by(condition, group, individual_code, fov_name, celltype_subset, fov_cluster) %>%
summarise(cells=n(), .groups = 'drop')
## Calculate pseudobulk
pseudobulk_counts <- PseudobulkExpression(so, assays = "RNA", layer="counts", method = 'aggregate', group.by = 'fov_cluster')
pseudobulk_counts_matrix <- pseudobulk_counts[["RNA"]]
# Change - back to _. Ideally have neither and skip this step
colnames(pseudobulk_counts_matrix)<-gsub("-","_",colnames(pseudobulk_counts_matrix))
# Determine fov_clusters with enough cells
# Filter both pseudobulk matrix and pseudobulk annotation
passed_fov_clusters <- celltype_summary_table$fov_cluster[celltype_summary_table$cells >= min_cells_per_fovcluster]
pseudobulk_counts_matrix <- pseudobulk_counts_matrix[,passed_fov_clusters]
pseudobulk_anno_table <- celltype_summary_table[passed_fov_clusters,]
# Calculate DE across every celltype
# Empty list to collect results
de_result_list <- list()
# celltype_subset is a matrix
for (the_celltype in levels(so$celltype_subset)) {
# Subset counts andn annotation to one cell type.
# Ensure order remains identical!
print(the_celltype)
anno_table.this <- pseudobulk_anno_table[pseudobulk_anno_table$celltype_subset == the_celltype,]
count_matrix.this <- pseudobulk_counts_matrix[,anno_table.this$fov_cluster]
## Check for sufficient replicates ##
# To do any calculations, we need at least 2 pseudobulk groups per contrast.
# there are plenty in this experiemnt, but with less replicates and rare cell types
# it may be neccesary to check and skip certain contrasts. Here woudl be a good
# if (not enouch samples to run test ) {next}
# skip clusters with no samples after filtering
if( nrow(anno_table.this) < 1 ) {next}
# Setup objects for limma
dge <- DGEList(count_matrix.this)
dge <- calcNormFactors(dge)
# Build model
group <- anno_table.this$group
individual_code <- anno_table.this$individual_code
# Model design
design <- model.matrix( ~0 + group)
# Run Voom
vm <- voom(dge, design = design, plot = FALSE)
# Adding dupliate correlation to use individual fovs, rather than pooled per biosample
corrfit <- duplicateCorrelation(vm, design, block=individual_code)
fit <- lmFit(vm, design, correlation = corrfit$consensus, block=individual_code)
# Define and fit contrasts and run ebayes
contrasts <- makeContrasts(UCvHC = groupUC - groupHC,
CDvHC = groupCD - groupHC,
levels=coef(fit))
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)
## Look through each contrast, and extract a results table.
for ( the_coef in colnames(contrasts) ) {
de_result.this <- topTable(fit, n = Inf, adjust.method = "BH", coef = the_coef) %>%
rownames_to_column("target") %>%
mutate(contrast=the_coef,
celltype=the_celltype) %>%
select(celltype,contrast,target,everything()) %>%
arrange(P.Value)
de_result_list[[paste(the_celltype, the_coef, sep="_")]] <- de_result.this
}
}
# Join together results for all celltypes, and pull out those with a singificant adjusted p-value
de_results_all <- bind_rows(de_result_list)
de_results_sig <- filter(de_results_all, adj.P.Val < 0.01)Results
DT::datatable(mutate(head(de_results_sig), across(is.numeric, signif, digits = 3)))This table is the typical output of limma tests; With a couple of extra columns added by our code.
- celltype: The celltype being tested (Added by example code)
- contrast: The contrast being tested (Added by example code)
- target : The gene name (Added by example code, is the rowname in limma output)
- rownames : The tested cell types
- logFC : Log 2 fold change between tested groups.
For a test of Test-Con;
- At logFC +1, A is doubled B.
- At logFC -1, A is half of B.
- A logFC 0 indicates no change.
- AveExpr : Average expression of a gene across all replicates.
- t : Moderated T-statistic. See Limma documentation.
- P.Value : P.value
- adj.P.Val : A multiple-hypothesis corrected p-value
- B : B statistic (rarely used). See Limma documentation.
More Information
- ‘Ochestrating single cell analysis with bioconductor’ book chapter ‘DE analyses between conditions’ : An explanation of the ‘pseudobulk’ approch to single cell differential expression calculation.
- Pseudobulk and differential expression (glmGamPoi documentation): Part of the documentation for the glmGamPoi R package (not used here), that features another clear explanation of pseudobulking.
- Bias, robustness and scalability in single-cell differential expression analysis (Soneson and Robinson 2018) : A review of single cell differential expression calculation methods.
- Seurat Differntial expression Vignette : How to run differential expression analyses using Seurat.
- limma documentation (Ritchie et al. 2015): The complete manual to limma.
- Use of within-array replicate spots for assessing differential expression in microarray experiments(Smyth, Michaud, and Scott 2005) : Paper describing the duplicate correlation approach for pseudoreplicates (in the original context of micorarray analysis).
- Bioconductor support: Further clarification on when not to use duplicateCorrelation with technical replicates (RNA-seq) : A very helpful discussion about when to use correlateDuplicates for pseuodreplicates.
- Differential Expression with Limma-Voom UC davis bioinformatics training : A more accessible explanation of bulk RNAseq analyses using limma.
- Interactions and contrasts : An excellent visual explanation of how to build linear models for more complex multi-factor experimental designs (e.g. treatment and genotype). Part of a larger Data Analysis for Genomics class resource.
References
sessionInfo()R version 4.3.2 (2023-10-31)
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] stats graphics grDevices datasets utils methods base
other attached packages:
[1] ggrepel_0.9.5 edgeR_4.0.16 DT_0.33 limma_3.58.1
[5] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[9] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[13] ggplot2_3.5.0 tidyverse_2.0.0 speckle_1.2.0 Seurat_5.1.0
[17] SeuratObject_5.0.2 sp_2.1-3 workflowr_1.7.1
loaded via a namespace (and not attached):
[1] RcppAnnoy_0.0.22 splines_4.3.2
[3] later_1.3.2 bitops_1.0-7
[5] polyclip_1.10-6 fastDummies_1.7.3
[7] lifecycle_1.0.4 rprojroot_2.0.4
[9] globals_0.16.3 processx_3.8.4
[11] lattice_0.22-6 MASS_7.3-60.0.1
[13] crosstalk_1.2.1 magrittr_2.0.3
[15] plotly_4.10.4 sass_0.4.9
[17] rmarkdown_2.26 jquerylib_0.1.4
[19] yaml_2.3.8 httpuv_1.6.15
[21] sctransform_0.4.1 spam_2.10-0
[23] spatstat.sparse_3.0-3 reticulate_1.35.0
[25] cowplot_1.1.3 pbapply_1.7-2
[27] RColorBrewer_1.1-3 abind_1.4-5
[29] zlibbioc_1.48.2 Rtsne_0.17
[31] GenomicRanges_1.54.1 BiocGenerics_0.48.1
[33] RCurl_1.98-1.14 git2r_0.33.0
[35] GenomeInfoDbData_1.2.11 IRanges_2.36.0
[37] S4Vectors_0.40.2 irlba_2.3.5.1
[39] listenv_0.9.1 spatstat.utils_3.0-4
[41] goftest_1.2-3 RSpectra_0.16-1
[43] spatstat.random_3.2-3 fitdistrplus_1.1-11
[45] parallelly_1.37.1 leiden_0.4.3.1
[47] codetools_0.2-20 DelayedArray_0.28.0
[49] tidyselect_1.2.1 farver_2.1.1
[51] matrixStats_1.2.0 stats4_4.3.2
[53] spatstat.explore_3.2-7 jsonlite_1.8.8
[55] progressr_0.14.0 ggridges_0.5.6
[57] survival_3.5-8 tools_4.3.2
[59] ica_1.0-3 Rcpp_1.0.12
[61] glue_1.7.0 gridExtra_2.3
[63] SparseArray_1.2.4 xfun_0.43
[65] MatrixGenerics_1.14.0 GenomeInfoDb_1.38.8
[67] withr_3.0.0 BiocManager_1.30.22
[69] fastmap_1.1.1 fansi_1.0.6
[71] callr_3.7.6 digest_0.6.35
[73] timechange_0.3.0 R6_2.5.1
[75] mime_0.12 colorspace_2.1-0
[77] scattermore_1.2 tensor_1.5
[79] spatstat.data_3.0-4 utf8_1.2.4
[81] generics_0.1.3 renv_1.0.5
[83] data.table_1.15.4 httr_1.4.7
[85] htmlwidgets_1.6.4 S4Arrays_1.2.1
[87] whisker_0.4.1 uwot_0.1.16
[89] pkgconfig_2.0.3 gtable_0.3.4
[91] lmtest_0.9-40 SingleCellExperiment_1.24.0
[93] XVector_0.42.0 htmltools_0.5.8.1
[95] dotCall64_1.1-1 scales_1.3.0
[97] Biobase_2.62.0 png_0.1-8
[99] knitr_1.45 rstudioapi_0.16.0
[101] tzdb_0.4.0 reshape2_1.4.4
[103] nlme_3.1-164 cachem_1.0.8
[105] zoo_1.8-12 KernSmooth_2.23-22
[107] parallel_4.3.2 miniUI_0.1.1.1
[109] pillar_1.9.0 grid_4.3.2
[111] vctrs_0.6.5 RANN_2.6.1
[113] promises_1.2.1 xtable_1.8-4
[115] cluster_2.1.6 evaluate_0.23
[117] cli_3.6.2 locfit_1.5-9.9
[119] compiler_4.3.2 rlang_1.1.3
[121] crayon_1.5.2 future.apply_1.11.2
[123] labeling_0.4.3 ps_1.7.6
[125] getPass_0.2-4 plyr_1.8.9
[127] fs_1.6.3 stringi_1.8.3
[129] viridisLite_0.4.2 deldir_2.0-4
[131] munsell_0.5.1 lazyeval_0.2.2
[133] spatstat.geom_3.2-9 Matrix_1.6-5
[135] RcppHNSW_0.6.0 hms_1.1.3
[137] patchwork_1.2.0 future_1.33.2
[139] statmod_1.5.0 shiny_1.8.1.1
[141] highr_0.10 SummarizedExperiment_1.32.0
[143] ROCR_1.0-11 igraph_2.0.3
[145] bslib_0.7.0