5) Cluster-free analysis

This script performs cluster-free compositional analysis and expression shift analysis using the R package cacoa (v0.4.0). Since clustering limits our resolution at which changes can be observed, cluster-free analysis is an annotation-independent approach. Cluster-free analysis can reveal changes in subpopulations below the resolution of the utilized annotation.

More info about the package cacoa here: https://github.com/kharchenkolab/cacoa pre-print: https://www.biorxiv.org/content/10.1101/2022.03.15.484475v1.full.pdf Vignette: https://pklab.med.harvard.edu/viktor/cacoa/walkthrough_short.html

Setup

devtools::load_all("/people/nrq364/cacoa_dev") # need to double check if dev version is necessary here
#library(cacoa)
library(magrittr)
library(qs)

1. Read in data

Read in cacoa object generated in notebook 3.

cao <- qread("cao.qs", nthreads = 10)

2. Cluster-free compositional changes

Compute cluster-free compositional changes.

cao$estimateCellDensity(method='graph',n.cores = 50) #alternative: kde (identifies more very small subclusters as changed)

cao$estimateDiffCellDensity(type='wilcox', n.cores = 50) # type = method to calculate differential cell density, options: permutation, t.test, wilcox or subtract (target subtract ref density)

Plot annotated UMAP and for multiple comparison corrected z-scores on the UMAP next to each other. Positive Z adj. means those cells are enriched in the treated condition. Negative Z adj means those cells are more abundant in controls.

plot_grid(cao$plotEmbedding(color.by='cell.groups'), cao$plotDiffCellDensity(legend.position=c(0, 1)), ncol=2)

3. Cluster-free expression shifts

Estimate cluster-free expression shifts (takes a long time, choose many cores). Z score is based on p-value and results in this version of cacoa (v0.4.0) are very conservative.

cao$estimateClusterFreeExpressionShifts(n.top.genes=300, gene.selection = "expression", normalize.both = T, n.cores = 50) 
# normalize.both whether to normalize results relative to distances within both conditions (TRUE) or only to the control (FALSE), "z" means it uses the not adjusted z scores?
# Warning: Setting normalize.both=TRUE likely leads to wrong results for cluster-free shiftsWarning: Please run estimateClusterFreeDE() first to use gene.selection='z' or 'lfc'. Fall back to gene.selection='expression'

p <- cao$plotClusterFreeExpressionShifts(legend.position=c(0, 1), font.size=c(2,3), min.z = 0.5, build.panel = F) # min.z default=qnorm(0.9)=1.28; when I change min.z to 0.5 I see significant results in adj z-scores

p1 <- p[[1]]$data %>%
         arrange(Color)  %>% ggplot(aes(x = x, y = y, color = Color)) + geom_point(size = 0.1) + scale_color_gradient2(mid = "grey90", high = "firebrick",  name = "Z score") + theme_void() # + guides(color=guide_legend(title="Z score")) #midpoint=0.0075,low = "grey95",

p2 <- p[[2]]$data%>%
         arrange(Color) %>% ggplot(aes(x = x, y = y, color = Color)) + geom_point(size = 0.1) + scale_color_gradient2(mid = "grey90", high = "firebrick", name = "Adj. Z score") + theme_void() #+ guides(color=guide_legend(title="Adj.Z score"))  midpoint = 0.3, low = "grey95",

library(patchwork)
p1+p2