3) Composition & expression shifts and sample structure

This script uses the R package cacoa (https://github.com/kharchenkolab/cacoa). In this script analysis of compositional changes and expression shifts as well as inspection of sample structure is performed. Only 2 conditions can be compared at a time.

Setup

#Load helper functions
source("/people/nrq364/NeuSiC/scRNA_helper_new.R")

library(cacoa)
library(magrittr)
library(dplyr)
library(conos)
library(pagoda2)
library(qs)
library(ggplot2)
library(ggrastr)

1. Read in data

Read in conos object and annotation factor generated in notebook 2 clustering and annotation.

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

2. Generate cacoa object

condition <- setNames(c("control", "control", "control", "control","control", "control", "treated", "treated", "treated","treated", "treated", "treated"), names(con$samples))

cao <- Cacoa$new(con, cell.groups = anno, sample.groups=condition, n.cores = 10, ref.level = "control", target.level = "treated")

qsave(cao, "cao.qs", nthreads = 10)

3. Compositional changes

Determine which cell types exhibit (significant) changes in abundance in response to the treatment/condition.

cao$estimateCellLoadings()

3.1 Default plot of compositional changes

Cell types above the red line are significantly altered. Negative cell loading means decreased abundance of the cell type in the treated samples, while positive cell loading points towards increased abundance in treated condition.

cao$plotCellLoadings(show.pvals = F) 

3.2 Polished plot of compositional changes

p2 <- cao$plotCellLoadings(show.pvals = F)
p2 <- p2$data

library(data.table)
p2 <- data.table(p2)
lvl <- p2[, abs(median(values)), by = ind][order(-V1),ind] %>% as.character()
levels(p2$ind) <- rev(lvl)

p2$ind <- factor(p2$ind, levels = rev(lvl))

palcomp <- setNames(ifelse(p2[, (median(values)), by = ind][order(-V1),V1] > 0, "blue", "grey50"), p2[, (median(values)), by = ind][order(-V1),ind])

ord <- p2[, abs(median(values)), by = ind][order(V1)]$ind

#define function to calculate interquartile range (IQR) for a numeric vector z
 iqr = function(z, lower = 0.25, upper = 0.75) {
  data.frame(
    y = median(z),
    ymin = quantile(z, lower),
    ymax = quantile(z, upper)
  )
}

p <- p2 %>% ggplot( aes( x = values, y = factor(ind, levels = ord), color = ind)) +
  geom_violin(scale = "width", 
                #fill = "#ededed", 
                color = NA, aes(fill = ind), alpha = 0.2) +
  theme_light() + theme(legend.position = "none") + scale_fill_manual(values = palcomp) +
  geom_vline(xintercept = 0, linewidth = 0.1) +
  #new_scale_fill() + 
  stat_summary(fun.data = iqr, geom = "errorbar", show.legend = F, aes(color = ind), stroke = 0.6, alpha = 0.9, color = "grey30", width = 0.25) + 
  stat_summary(fun.data = iqr, geom = "point", show.legend = F,aes(fill = ind), shape = 23, stroke = 1, color = "grey30", alpha = 0.9, size = 2) +  scale_fill_manual(values = palcomp) + xlab("separation coefficient") + ylab("") + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) + ggtitle("  enriched in control               enriched in treated")
 # annotate("text", x = -0.5,  y = length(levels(p2$ind)) + 0.5, label = "Left Heading", hjust = 1, angle = 0) +
 # annotate("text", x = 0.5,  y = length(levels(p2$ind)) + 0.5, label = "Right Heading", hjust = 0, angle = 0) +
  theme(plot.margin = unit(c(1, 1, 1, 1), "cm"))

List of 1
 $ plot.margin: 'simpleUnit' num [1:4] 1cm 1cm 1cm 1cm
  ..- attr(*, "unit")= int 1
 - attr(*, "class")= chr [1:2] "theme" "gg"
 - attr(*, "complete")= logi FALSE
 - attr(*, "validate")= logi TRUE

p

#p+scale_y_discrete(labels= rev(c("Immature_PN_Tbr1", "PIC_immature_PN","Highly_prolif_RG_S", "CGE_progenitors", "Cajal_Retzius_cells", "Roof_plate", "Immature_PN","MGE_progenitors", "Dorsal_IPC","Microglia","Dorsal_RG","Highly_prolif_RG_G2/M","Anteriomedial_RG","Lateral_RG","Cortical_hem","Immature_GABAergic_N")))+theme(axis.text.y = element_text( vjust = 0.5)) +theme(axis.text.y =  element_text(size = 10))
#ggsave("Compositionshifts.pdf",width=7, height=4.5)

4. Expression shifts

Specific number of top DE genes between conditions are used to calculate expression shifts. top.n.genes parameter needs to be adjusted per dataset. In case of neurons use a lower number than the default (1000) because neurons are very similar. If you normalize among many similar genes, you will not get a meaningful result. Check sample positions in UMAP plot in expression.space (see 5.) when using different numbers of top.n.genes.

cao$estimateExpressionShiftMagnitudes(top.n.genes = 300) 

4.1 Default plot of expression changes

y-axis shows magnitude of changes, while asterisks on top of bars show their significance. Expression shifts are calculated per cell-type on sample-wise pseudobulk gene expression profiles using only the top x DEG When using the default way to calculate the pairwise expression shifts ( shift = mean(A-B) - 1/2*(mean(A-A)+mean(B-B)) ), negative values mean that the average distance within groups (could be control (A) and/or condition (B)) is bigger than the average distance across groups.

cao$plotExpressionShiftMagnitudes(show.pvalues = "adjusted", show.jitter = F, notch = F)

4.2 Polished plot of expression changes

# define function to calculate interquartile range (25th and 75th percentile) and median
 iqr = function(z, lower = 0.25, upper = 0.75) {
  data.frame(
    y = median(z),
    ymin = quantile(z, lower),
    ymax = quantile(z, upper)
  )
}

p3 <- cao$plotExpressionShiftMagnitudes()
p <- p3$data %>% ggplot(aes( y = value, x = Type, color = Type)) +
  geom_violin(scale = "width", 
                #fill = "#ededed", 
                color = NA, aes(fill = Type), alpha = 0.2, 
              trim = F) + theme_light() + theme(legend.position = "none", axis.text.x = element_text(angle = 90, hjust = 1), panel.grid.minor = element_blank(), panel.grid.major = element_blank()) + scale_color_manual(values = cao$cell.groups.palette)+ scale_fill_manual(values = cao$cell.groups.palette) +
  geom_hline(yintercept = 0, linewidth = 0.1) +
  stat_summary(fun.data = iqr, geom = "errorbar", show.legend = F, aes(color = Type),stroke = 0.5, alpha = 0.9,color = "grey30", width = 0.25) + 
  stat_summary(fun.data = iqr, geom = "point", show.legend = F,aes(fill = Type), shape = 23, color = "grey30", alpha = 0.9)  + ylab(("normalized expression\ndistance")) + xlab("") + theme(plot.margin = margin(0.5,0,0,2, "cm"))
#+scale_x_discrete(labels= c("Immature_GABAergic_N", "Highly_prolif_RG_S", "PIC_immature_PN","Dorsal_RG","Highly_prolif_RG_G2/M","Lateral_RG", "Cortical_hem","MGE_progenitors", "Dorsal_IPC", "Microglia", "CGE_progenitors", "Roof_plate", "Anteriomedial_RG", "Immature_PN","Immature_PN_Tbr1", "Cajal_Retzius_cells"))+theme(axis.text.x = element_text( vjust = 0.5)) +theme(axis.text.x =  element_text(size = 10))
p

#ggsave("Expressionshifts_Vln.png",width=7, height=4.5)

5. Inspect sample structure

To inspect sample structure, expression shift and compositional analysis need to be run first.

5.1 Metadata seperation

= proportion of variance explained by covariates. Create metadata table per sample.

metadata <- data.frame(row.names = levels(cao$sample.per.cell), condition=c("control", "control", "control", "control","control", "control", "treated", "treated", "treated","treated", "treated", "treated"),version =c("3.0", "3.0", "3.0","3.1","3.1","3.1","3.0", "3.0", "3.0","3.0","3.1","3.1"), sequencing= c("batch1", "batch1","batch2","batch3", "batch3","batch3","batch1", "batch1","batch2","batch2","batch3", "batch3"), sex = c("male","mix","male","mix", "mix","mix","mix","male" ,"mix", "mix","mix","mix") )
metadata

This function shows a warning if expression shifts are significantly confounded by covariates as for example sex, sequencing batch or 10x version, instead of being only driven by condition:

pvals <- cao$estimateMetadataSeparation(sample.meta = metadata,  space = "expression.shifts")

  |                                                                                                          |   0%, ETA NA
  |==========================                                                                             |  25%, ETA 00:12
  |====================================================                                                   |  50%, ETA 00:04
  |=============================================================================                          |  75%, ETA 00:01
  |===================================================================================================| 100%, Elapsed 00:05

Warning: Significant separation by: condition

sort(pvals$padjust)

 condition        sex sequencing    version 
0.00559888 0.14917017 0.17889755 0.40951810 

This function shows a warning if compositional changes are significantly confounded by covariates as for example sex, sequencing batch or 10x version, instead of being only driven by condition:

pvals <- cao$estimateMetadataSeparation(sample.meta = metadata,  space = "coda")

  |                                                                                                          |   0%, ETA NA
  |==========================                                                                             |  25%, ETA 00:12
  |====================================================                                                   |  50%, ETA 00:04
  |=============================================================================                          |  75%, ETA 00:02
  |===================================================================================================| 100%, Elapsed 00:05

sort(pvals$padjust)

 condition sequencing        sex    version 
0.09224822 0.09224822 0.09224822 0.83383323 

5.2 expression structure of samples

Now we can plot sample expression structure in a UMAP colored by the metadata. Ideally, we expect a clustering of the samples by condition (like here) and not by other factors like sex or sequencing.

cao$plotSampleDistances(space='expression.shifts', font.size=4, show.sample.size=T, method="UMAP", sample.colors = setNames(metadata$condition, names(cao$sample.groups)), color.title = "condition") #change colour by specifying palette = ; method can be changed to "MDS"

Warning: n_components > number of columns in input data: 2 > 1, this may give poor or unexpected results

cao$plotSampleDistances(space='expression.shifts', font.size=4, show.sample.size=T, method="UMAP", sample.colors = setNames(metadata$sequencing, names(cao$sample.groups)), color.title = "sequencing") #change colour by specifying palette = ; method can be changed to "MDS"

Warning: n_components > number of columns in input data: 2 > 1, this may give poor or unexpected results