In saeyslab/muscatWrapper: Wrapper functions around muscat for differential expression analysis of multi-sample multi-condition scRNAseq datasets

knitr::opts_chunk$set(
  collapse = TRUE,
  # comment = "#>",
  warning = FALSE,
  message = FALSE
)

In this vignette, you can learn how to perform a muscat differential state (DS) analysis. A DS analysis can be performed if you have multi-sample, multi-group single-cell data. For each cell type of interest, muscat will compare the sample-wise expression of all genes between groups of interest. Therefore, the absolute minimum of meta data you need to have, are following columns indicating for each cell: the group, sample and cell type.

In addition to these factors, in this vignette we will also perform the DE analysis while correcting for covariates/batch effects.

As example expression data, we will use data from Puram et al. of the tumor microenvironment in head and neck squamous cell carcinoma (HNSCC) [See @puram_single-cell_2017] . The groups we have here are tumors scoring high for a partial epithelial-mesenschymal transition (p-EMT) program vs low-scoring tumors.

The different steps of the MultiNicheNet analysis are the following:

• 1. Preparation of the analysis: load packages, read in the single-cell expression data, and define the main settings of the muscat analysis
• 1. Check cell type abundance for the cell types of interest
• 1. Perform genome-wide differential expression (DS) analysis
• 1. Downstream analysis of the DS output, including visualization
• 1. Exploring the effect of the batch effect correction

In this vignette, we will demonstrate all these steps in detail.

Step 0: Preparation of the analysis: load packages, read in the single-cell expression data, and define the main settings of the muscat analysis

IMPORTANT: The current implementation of the muscat wrapper starts from a SingleCellExperiment object, therefore we will need to load the SingleCellExperiment library. If you start from a Seurat object, you can convert it easily to a SingleCellExperiment via sce = Seurat::as.SingleCellExperiment(seurat_obj, assay = "RNA").

library(SingleCellExperiment)
library(dplyr)
library(ggplot2)
library(muscatWrapper)

In this case study, we want to study differences in expression between pEMT-high and pEMT-low tumors. The meta data columns that indicate the pEMT status of tumors are 'pEMT' and 'pEMT_fine', cell type is indicated in the 'celltype' column, and the sample is indicated by the 'tumor' column.

sce = readRDS(url("https://zenodo.org/record/5196144/files/sce_hnscc.rds"))
scater::plotReducedDim(sce, dimred = "UMAP", colour_by = "celltype")
scater::plotReducedDim(sce, dimred = "UMAP", colour_by = "tumor")
scater::plotReducedDim(sce, dimred = "UMAP", colour_by = "pEMT")
scater::plotReducedDim(sce, dimred = "UMAP", colour_by = "pEMT_fine")

We will also look at at a potential covariate/batch effect of interest. For this dataset, no batch effect was known. Therefore we will simulate a mock batch, just to demosntrate the code. (In the future, we will use another example dataset in this vignette.)

set.seed(1919)
sample_id = "tumor"
extra_metadata = SummarizedExperiment::colData(sce)  %>% tibble::as_tibble()  %>% dplyr::select(all_of(sample_id)) %>% dplyr::distinct() %>% mutate(batch = sample(c("A","B"),nrow(.), replace = TRUE))
new_metadata = SummarizedExperiment::colData(sce)  %>% tibble::as_tibble()  %>% inner_join(extra_metadata)
new_metadata = new_metadata %>% data.frame()
rownames(new_metadata) = new_metadata$cell

sce = SingleCellExperiment::SingleCellExperiment(list(counts=SingleCellExperiment::counts(sce), logcounts=SingleCellExperiment::logcounts(sce)),
                                                 reducedDims = SingleCellExperiment::reducedDims(sce),
                                                 colData=new_metadata,
                                                 rowData=SingleCellExperiment::rowData(sce),
                                                 metadata=sce@metadata
)
scater::plotReducedDim(sce, dimred = "UMAP", colour_by = "batch")

Now we will define in which metadata columns we can find the group, sample, cell type and covariates IDs

For the group_id in this vignette, we choose the 'pEMT' column instead of 'pEMT_fine'.

sample_id = "tumor"
group_id = "pEMT"
celltype_id = "celltype"
covariates = "batch"

Now we will go to the first real step of the muscat analysis

Step 1: Check cell type abundance

Step 1a: Number of cells per celltype-sample combination

We will now check the number of cells per cell type condition combination, and the number of patients per condition. This is important because muscat performs pseudobulking to infer group differences at the sample level for each cell type. This means that we will group the information of all cells of a cell type in a sample together to get 1 sample-celltype estimate. The more cells we have, the more accurate this aggregated expression measure will be.

table(SummarizedExperiment::colData(sce)$celltype, SummarizedExperiment::colData(sce)$tumor) # cell types vs samples
table(SummarizedExperiment::colData(sce)$celltype, SummarizedExperiment::colData(sce)$pEMT) # cell types vs conditions
table(SummarizedExperiment::colData(sce)$tumor, SummarizedExperiment::colData(sce)$pEMT) # samples vs conditions
table(SummarizedExperiment::colData(sce)$tumor, SummarizedExperiment::colData(sce)$batch) # samples vs batches

As you can see in the upper table, some Celltype-Sample combinations have 0 cells. It is possible that during DE analysis, some cell types will be removed from the analysis if there is not enough information to do a DE analysis. (More info later)

We can define the minimum number of cells we require per celltype-sample combination. It is recommened to have at least 10 (and preferably more) cells in each sample-celltype combination. Therefore we will set the min_cells parameter in the analysis to 10. Celltype-sample combinations with less cells will not be considered during the muscat DS analysis.

min_cells = 10

We will now calculate the cell type abundance table and make the visualizations:

abundance_output = get_abundance_info(sce, sample_id, group_id, celltype_id, min_cells, covariates = covariates)
head(abundance_output$abundance_data)

To visually see which celltype-sample combinations won't be considered, you can run the following code:

abundance_output$abund_plot_sample

Note: each group is here split into its batches as well!

Celltype-sample combinations that won't be considered are indicated in red (because they have less cells than the min_cells threshold indicated by the red dashed line)

If too many celltype-sample combinations don't pass this threshold, we recommend to define your cell types in a more general way it this would still be possible and make sense biologically (--> use one level higher of the cell type ontology hierarchy; eg TH17 CD4T cells --> CD4T cells | but not myeloid + T.cell together).

We can see here that quite many sample-celltype combinations are left out. For Endothelial, Myeloid, and T cells, we don't even have two or more samples in each group that have enough cells of those cell types. When we don't have two or more samples per group left, we cannot do a group comparison (we need at least 2 replicates per group for a statistical analysis). Therefore, those cell types will be removed before the DE analysis.

As stated before when seeing this, we would recommend to use a higher-level cell type annotation if possible. But the annotation here is already high-level, and grouping Endothelial cells, T cells and Myeloid cells eg would not make sense biologically. That we won't be able to include these cell types in our analysis is a limitation of the muscat approach compared to classic cell-level-based approaches (like Seurat::FindMarkers). On the contrary, those cell-level-based approaches don't reveal the lack of cells in many samples, and might lead to biased results.

Step 1b: Differential cell type abundance between the groups of interest

In another visualization, we can compare the cell type abundances between the groups of interest. This can be interesting because too strong abundance differences might have an effect on the DS analysis. Downstream results of these cell types should then be considered with some caution.

abundance_output$abund_plot_group

Note: each group is here split into its batches as well!

Differential abundance looks quite OK for the cell types kept for the DE analysis (i.e. CAF, Malignant and myofibroblast)

We can also look at cell type proportions per sample, and compare this between the different groups

abundance_output$abund_barplot

Note: each group is here split into its batches as well!

Conclusion of this step:

Important: Based on the cell type abundance diagnostics, we recommend users to change their analysis settings if required, before proceeding with the rest of the analysis.

Step 2: Perform genome-wide differential expression analysis

Now we will go over to the multi-group, multi-sample differential expression (DE) analysis (also called 'differential state' analysis by the developers of Muscat).

Define the contrasts and covariates of interest for the DE analysis.

Here, we want to compare the p-EMT-high vs the p-EMT-low group and find genes that are differentially expressed in high vs low pEMT. We don't have other covariates to correct for in this dataset. If you would have covariates you can correct for, we recommend doing this.

about contrasts and how to set them:

Note the format to indicate the contrasts! (This formatting should be adhered to very strictly, and white spaces are not allowed)

contrasts_oi = c("'High-Low','Low-High'")
contrast_tbl = tibble(contrast = 
                        c("High-Low","Low-High"), 
                      group = c("High","Low"))

Perform the DE analysis for each cell type.

muscat_output = muscat_analysis(
     sce = sce,
     celltype_id = celltype_id,
     sample_id = sample_id,
     group_id = group_id,
     covariates = covariates,
     contrasts_oi = contrasts_oi,
     contrast_tbl = contrast_tbl)

Check DE results

The muscat_output object contains both the default output as given by the muscat::pbDS() function, and a cleaner output table (through muscat::resDS). It also contains a table with per gene the average expression, the fraction of cells in a celltype expressing it, and the pseudobulk expression: and this grouped per celltype per sample.

Table with this latter cell type info

muscat_output$celltype_info %>% lapply(head)

Table with logFC and p-values for each gene-celltype-contrast:

muscat_output$celltype_de$celltype_de$de_output_tidy %>% arrange(p_adj) %>% head()

Observe that the top gene here is different than in the analysis without batch effect correction Multi-sample Multi-condition Differential Expression Analysis via Muscat: HNSCC application:vignette("basic_analysis", package="muscatWrapper")

We can also show the distribution of the p-values:

muscat_output$celltype_de$hist_pvals

(Note: this p-value histograms are the same for High-Low and Low-High because we only have two groups and compare them to each other - a DE gene in one comparison will then also be DE in the other comparison, with just a reversed sign of the logFC)

In order to trust the p-values, the p-value distributions should be uniform distributions, with a peak allowed between 0 and 0.05 if there would be a clear biological effect in the data. This clear effect (=clear DE) seems to be present here in the Malignant cell type populations, although the histogram is not very uniformly distributed for p-values between 0.05 and 0.25. This might point to issues in the DE model definition. Most common issues occur when we did not add all important covariates to the model, or when there is substructure present: meaning that one of the groups is actually a combination of multiple groups.

Step 3: Downstream analysis and visualization

Visualize the expression of some DE genes: here as example: genes higher in the pEMT high tumors in the Malignant cell type

celltype_oi = "Malignant"
group_oi = "High"

DE_genes = muscat_output$celltype_de$celltype_de$de_output_tidy  %>% inner_join(contrast_tbl) %>% filter(group == group_oi) %>% filter(cluster_id == celltype_oi) %>% filter(p_adj <= 0.05 & logFC >= 1) %>% arrange(p_adj) %>% pull(gene) %>% unique()
DE_genes

(Note 1 : Due to the pseudoubulking, single-cell level information is lost and Muscat can be underpowered. Therefore it is possible that are sometimes no significant DE genes after multiple testing correction. In that case, using less stringent cutoffs is better)

(Note 2 : If having a few samples per group (<5), it is likely that some DE genes will be driven by an outlier sample. Therefore it is always necessary to visualize the expression of the DE genes in the violin and dotplots shown here)

First, make a violin plot

gene_oi = DE_genes[1]

violin_plot = make_DEgene_violin_plot(sce = sce, gene_oi = gene_oi, celltype_oi = celltype_oi, group_id = group_id, sample_id = sample_id, celltype_id = celltype_id, covariate_oi = covariates)
violin_plot

Then a Dotplot

dotplots = make_DEgene_dotplot_pseudobulk_covariate(genes_oi = DE_genes, celltype_info = muscat_output$celltype_info, abundance_data = abundance_output$abundance_data, celltype_oi = celltype_oi, covariate_oi = covariates)
dotplots$pseudobulk_plot 
dotplots$singlecell_plot

Step 4: Exploring the effect of the batch effect correction

Run the analysis without batch effect correction

This is needed to compare the results with the batch-effect-correction analysis performed before in this vignette

covariates_uncorrected = NA

abundance_output_uncorrected = get_abundance_info(sce, sample_id, group_id, celltype_id, min_cells, covariates = covariates_uncorrected)

muscat_output_uncorrected = muscat_analysis(
     sce = sce,
     celltype_id = celltype_id,
     sample_id = sample_id,
     group_id = group_id,
     covariates = covariates_uncorrected,
     contrasts_oi = contrasts_oi,
     contrast_tbl = contrast_tbl)

Find genes with strongest batch effect correction coefficient in the statistical model

coef_df = muscat_output$celltype_de$celltype_de$de_output$fit[[celltype_oi]]$coefficients %>% as.data.frame() %>% tibble::rownames_to_column("gene") %>% tibble::as_tibble()
head(coef_df)

Sort the genes based on the covariate effect

coef_df %>% arrange(-ei.batchB) %>% head()
coef_df %>% arrange(ei.batchB) %>% head()

In the case your covariate corresponds to a well-known effect, such as biological sex, you will see genes relevant for this process pop up in this list (eg Y chromosome genes for the sex covariate)

Interpretation of these values: if value of the ei.batchB coefficient is positive: higher expression in batch B, if negative: higher expression in batch A.

We will save the top10 genes with highest coefficient values in both directions

top_A_genes = coef_df %>% arrange(ei.batchB) %>% head(10) %>% pull(gene)
top_B_genes = coef_df %>% arrange(-ei.batchB) %>% head(10) %>% pull(gene)

Merge DE table of analysis with and without batch effect correction

And compare the "uncorrected logFC" between the groups of interest with the "corrected logFC"

muscat_output_joined = muscat_output_uncorrected$celltype_de$celltype_de$de_output_tidy %>% distinct(gene, cluster_id, logFC, p_adj, contrast) %>% rename(logFC_uncorrected = logFC, p_adj_uncorrected = p_adj) %>% inner_join(
  muscat_output$celltype_de$celltype_de$de_output_tidy %>% distinct(gene, cluster_id, logFC, p_adj, contrast) %>% rename(logFC_corrected = logFC, p_adj_corrected = p_adj)
) %>% mutate(diff_logFC = logFC_corrected - logFC_uncorrected)

muscat_output_joined %>% arrange(-abs(diff_logFC)) %>% filter(cluster_id == celltype_oi)

Check now also the DE analysis results and difference for the top A/B genes

muscat_output_joined %>% arrange(-abs(diff_logFC)) %>% filter(cluster_id == celltype_oi) %>% filter(gene %in% top_A_genes)
muscat_output_joined %>% arrange(-abs(diff_logFC)) %>% filter(cluster_id == celltype_oi) %>% filter(gene %in% top_B_genes)

Check now also once the batch coefficients for the genes with strongest logFC differences

top_diff_genes = muscat_output_joined %>% arrange(-abs(diff_logFC)) %>% filter(cluster_id == celltype_oi) %>% pull(gene) %>% head(20) %>% unique()
coef_df %>% filter(gene %in% top_diff_genes) # we see indeed quite high Batch coefficient values

And define a list of genes that are more clearly DE without batch effect correction than with

DE_without_cor = muscat_output_joined %>% arrange(diff_logFC) %>% filter(cluster_id == celltype_oi) %>% filter(contrast == "High-Low")  %>% filter(logFC_corrected < logFC_uncorrected) %>% filter(logFC_uncorrected > 1) %>% pull(gene) %>% head(10) %>% unique()
coef_df %>% filter(gene %in% DE_without_cor) # we see indeed quite high Batch coefficient values

Visualization of the batch effect correction

Show expression of these top-batch-effect-correction effected genes: with non-corrected pseudobulkexpression values

dotplots_topA = make_DEgene_dotplot_pseudobulk_covariate(genes_oi = top_A_genes, celltype_info = muscat_output_uncorrected$celltype_info, abundance_data = abundance_output$abundance_data, celltype_oi = celltype_oi, covariate_oi = covariates)
dotplots_topA$pseudobulk_plot 

As can be expected if we show non-corrected values, higher expression of topA genes in A-samples can be observed. The same is true for the B genes:

dotplots_topB = make_DEgene_dotplot_pseudobulk_covariate(genes_oi = top_B_genes, celltype_info = muscat_output_uncorrected$celltype_info, abundance_data = abundance_output$abundance_data, celltype_oi = celltype_oi, covariate_oi = covariates)
dotplots_topB$pseudobulk_plot 

Now: visualize the expression of the genes that are more clearly DE (high in pEMT-high group) without correction than with

dotplots_DE_without_cor = make_DEgene_dotplot_pseudobulk_covariate(genes_oi = DE_without_cor, celltype_info = muscat_output_uncorrected$celltype_info, abundance_data = abundance_output$abundance_data, celltype_oi = celltype_oi, covariate_oi = covariates)
dotplots_DE_without_cor$pseudobulk_plot 

Show expression of these top-batch-effect-correction effected genes: with corrected pseudobulkexpression values (default plot of muscatWrapper)

dotplots_topA = make_DEgene_dotplot_pseudobulk_covariate(genes_oi = top_A_genes, celltype_info = muscat_output$celltype_info, abundance_data = abundance_output$abundance_data, celltype_oi = celltype_oi, covariate_oi = covariates)
dotplots_topA$pseudobulk_plot 

In the original expression space, these genes have clearly higher expression in batch A samples. However, we can see that after batch correction, this effect is not visible anymore. The batch effect correction had its effect.

This is less clear though for the B genes:

dotplots_topB = make_DEgene_dotplot_pseudobulk_covariate(genes_oi = top_B_genes, celltype_info = muscat_output$celltype_info, abundance_data = abundance_output$abundance_data, celltype_oi = celltype_oi, covariate_oi = covariates)
dotplots_topB$pseudobulk_plot 

Now: visualize the expression of the genes that are more clearly DE (high in pEMT-high group) without correction than with

dotplots_DE_without_cor = make_DEgene_dotplot_pseudobulk_covariate(genes_oi = DE_without_cor, celltype_info = muscat_output$celltype_info, abundance_data = abundance_output$abundance_data, celltype_oi = celltype_oi, covariate_oi = covariates)
dotplots_DE_without_cor$pseudobulk_plot 

In contrast to the visualization without correction, these genes are less clearly higher in the p-EMT group

References

Crowell, H.L., Soneson, C., Germain, PL. et al. muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nat Commun 11, 6077 (2020). https://doi.org/10.1038/s41467-020-19894-4

Puram, Sidharth V., Itay Tirosh, Anuraag S. Parikh, Anoop P. Patel, Keren Yizhak, Shawn Gillespie, Christopher Rodman, et al. 2017. “Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer.” Cell 171 (7): 1611–1624.e24. https://doi.org/10.1016/j.cell.2017.10.044.

saeyslab/muscatWrapper documentation built on March 11, 2023, 6:14 p.m.