In HectorRDB/condimentsPaper: Package to accompany the condiments Paper

library(knitr)
opts_chunk$set(
  fig.pos = "!h", out.extra = "", warning = F, 
  fig.align = "center"
)

Load data

libs <- c("here", "tidyr", "dplyr", "Seurat", "scater", "slingshot", "cowplot",
          "condiments", "ggplot2", "readr", "tradeSeq",  "pheatmap", "fgsea",
          "msigdbr", "scales")
suppressMessages(
  suppressWarnings(sapply(libs, require, character.only = TRUE))
)
rm(libs)
theme_set(theme_classic() +
            theme(rect = element_blank(),
                  plot.background = element_blank(),
                  strip.background = element_rect(fill = "transparent")
                  ))

fibrosis <- condimentsPaper::import_fibrosis()

data("fibrosis", package = "condimentsPaper")

EDA

Reduced dimension coordinates are obtained from the original publication.

df <- bind_cols(reducedDim(fibrosis, "UMAP") %>% as.data.frame(),
                colData(fibrosis)[, 1:14] %>% as.data.frame())
ggplot(df, aes(x = UMAP_1, y = UMAP_2, col = Status)) +
  geom_point(size = .5) +
  scale_color_brewer(palette = "Accent") +
  labs(col = "Disease") +
  theme(legend.position = c(.15, .4),
        legend.background = element_blank()) +
  guides(col = guide_legend(override.aes = list(size = 3)))
ggplot(df, aes(x = UMAP_1, y = UMAP_2, col = celltype)) +
  geom_point(size = .5) +
  scale_color_brewer(palette = "Dark2") +
  labs(col = "Cell Type") +
  theme(legend.position = c(.2, .4),
        legend.background = element_blank()) +
  guides(col = guide_legend(override.aes = list(size = 3)))

Imbalance score

scores <- condiments::imbalance_score(
  Object = df %>% select(UMAP_1, UMAP_2) %>% as.matrix(), 
  conditions = df$Status, k = 20, smooth = 30)
df$scores <- scores$scaled_scores
ggplot(df, aes(x = UMAP_1, y = UMAP_2, col = scores)) +
  geom_point(size = .7) +
  scale_color_viridis_c(option = "C", breaks = c(0, 2, 4, 6)) +
  labs(col = "Scores") +
  theme(legend.position = c(.15, .4),
        legend.background = element_blank())

Differential Topology

To estimate the trajectory, we use slingshot [@Street2018a].

Fit slingshot

fibrosis <- slingshot(fibrosis, reducedDim = 'UMAP',
                 clusterLabels = colData(fibrosis)$celltype,
                 start.clus = 'AT1', approx_points = 100)

Topology Test

set.seed(821)
topologyTest(SlingshotDataSet(fibrosis), fibrosis$Status)

fibrosis_SCGB3A2 <- slingshot(fibrosis[, fibrosis$celltype != "SCGB3A2+"],
                         reducedDim = 'UMAP',
            clusterLabels = colData(fibrosis[, fibrosis$celltype != "SCGB3A2+"])$celltype,
                 start.clus = 'AT1', approx_points = 100)
set.seed(821)
topologyTest(SlingshotDataSet(fibrosis_SCGB3A2), fibrosis_SCGB3A2$Status)
rm(fibrosis_SCGB3A2)

Separate trajectories

sdss <- slingshot_conditions(SlingshotDataSet(fibrosis), fibrosis$Status)
curves <- bind_rows(lapply(sdss, slingCurves, as.df = TRUE),
                    .id = "Status")

p4a <- ggplot(df, aes(x = UMAP_1, y = UMAP_2, col = Status)) +
  geom_point(size = .7, alpha = .2) +
  scale_color_brewer(palette = "Accent") +
  labs(col = "Disaes") +
  geom_path(data = curves %>% arrange(Status, Lineage, Order),
            aes(group = interaction(Lineage, Status)), size = 1.5) +
  annotate("text", x = -10, y = 6, label = "Lineage1", size = 5) +
  annotate("text", x = -7, y = -2.7, label = "Lineage2", size = 5) +
  theme(legend.position = c(.15, .35),
        legend.background = element_blank()) +
  NULL
p4a

Common trajectory

df <- bind_cols(
  as.data.frame(reducedDim(fibrosis, "UMAP")),
  slingPseudotime(fibrosis) %>% as.data.frame() %>%
    dplyr::rename_with(paste0, "_pst", .cols = everything()),
  slingCurveWeights(fibrosis) %>% as.data.frame(),
  ) %>%
  mutate(Lineage1_pst = if_else(is.na(Lineage1_pst), 0, Lineage1_pst),
         Lineage2_pst = if_else(is.na(Lineage2_pst), 0, Lineage2_pst),
         pst = if_else(Lineage1 > Lineage2, Lineage1_pst, Lineage2_pst),
         pst = max(pst) - pst)
curves <- slingCurves(fibrosis, as.df = TRUE)
ggplot(df, aes(x = UMAP_1, y = UMAP_2)) +
  geom_point(size = .7, aes(col = pst)) +
  scale_color_viridis_c() +
  labs(col = "Pseudotime") +
  geom_path(data = curves %>% arrange(Order),
            aes(group = Lineage), col = "black", size = 1.5) +
  annotate("text", x = -10, y = 6, label = "Lineage1", size = 5) +
  annotate("text", x = -7, y = -2.7, label = "Lineage2", size = 5) +
  theme(legend.position = c(.15, .35),
        legend.background = element_blank()) 

Differential Progression

Test

progressionTest(fibrosis, conditions = fibrosis$Status, lineages = TRUE)

Plot

df <-  slingPseudotime(fibrosis) %>% as.data.frame() 
df$Status <- fibrosis$Status
df <- df %>% 
  pivot_longer(-Status, names_to = "Lineage",
               values_to = "pst") %>%
  filter(!is.na(pst))
ggplot(df, aes(x = pst)) +
  geom_density(alpha = .8, aes(fill = Status), col = "transparent") +
  geom_density(aes(col = Status), fill = "transparent", size = 1.5) +
  labs(x = "Pseudotime", fill = "Status") +
  facet_wrap(~Lineage, scales = "free_x") +
  guides(col = "none", fill = guide_legend(
    override.aes = list(size = 1.5, col = c("#7FC97F", "#BEAED4"))
  )) +
  scale_fill_brewer(palette = "Accent") +
  scale_color_brewer(palette = "Accent") +
  theme(legend.position = c(.8, .8), legend.direction = "horizontal")

Differential fate selection

Test

fateSelectionTest(fibrosis, conditions = fibrosis$Status)

Plot

df <- slingCurveWeights(SlingshotDataSet(fibrosis), as.probs = TRUE) %>% as.data.frame()
df$Status <- fibrosis$Status
df <- df %>%
  pivot_longer(-Status, names_to = "Lineage", values_to = "ws")
ggplot(df, aes(x = ws)) +
  geom_density(alpha = .8, aes(fill = Status), col = "transparent") +
  geom_density(aes(col = Status), fill = "transparent", size = 1.5) +
  labs(x = "Weights", fill = "Status") +
  facet_wrap(~Lineage, scales = "free_x") +
  guides(fill = "none", col = "none") +
  scale_fill_brewer(palette = "Accent") +
  scale_color_brewer(palette = "Accent")

Differential expression

We use tradeSeq [@VandenBerge2020].

Select number of knots

library(tradeSeq)
set.seed(3)
BPPARAM <- BiocParallel::bpparam()
BPPARAM$workers <- 3
icMat <- evaluateK(counts = fibrosis, sds = SlingshotDataSet(fibrosis), 
                   conditions = factor(fibrosis$Status),
                   nGenes = 300, parallel = TRUE, BPPARAM = BPPARAM, k = 3:7)

Fit GAM

Next, we fit the NB-GAMs using 6 knots, based on the pseudotime and cell-level weights estimated by Slingshot. We use the conditions argument to fit separate smoothers for each condition.

library(tradeSeq)
set.seed(3)
fibrosis <- fitGAM(counts = fibrosis, conditions = factor(fibrosis$Status), nknots = 6)

Differential expression between conditions

condRes <- conditionTest(fibrosis, l2fc = log2(2), lineages = TRUE)
condRes$padj <- p.adjust(condRes$pvalue_lineage1, "fdr")
mean(condRes$padj <= 0.05, na.rm = TRUE)
sum(condRes$padj <= 0.05, na.rm = TRUE)
conditionGenes <- rownames(condRes)[condRes$padj <= 0.05]
conditionGenes <- conditionGenes[!is.na(conditionGenes)]

Differential expression between lineages

patRes <- patternTest(fibrosis, l2fc = log2(2))
patRes$padj <- p.adjust(patRes$pvalue, "fdr")
mean(patRes$padj <= 0.05, na.rm = TRUE)
sum(patRes$padj <= 0.05, na.rm = TRUE)
patternGenes <- rownames(patRes)[patRes$padj <= 0.05]
patternGenes <- patternGenes[!is.na(patternGenes)]

Visualize select genes from the og paper

library(RColorBrewer)
scales <- c("#FCAE91", "#DE2D26", "#BDD7E7", "#3182BD")

# plot genes
oo <- order(patRes$waldStat, decreasing = TRUE)
genes <- c("SCGB3A2", "SFTPC", "ABCA3", "AGER")

# most significant gene
ps <- lapply(genes, function(gene){
  plotSmoothers(fibrosis, assays(fibrosis)$counts, gene = gene, curvesCol = scales,
                border = FALSE, sample = .2) +
    scale_x_reverse(breaks = 0:4/.4, labels = 4:0/.4) +
    labs(title = gene, col = "Lineage and Condition") +
    scale_color_manual(values = scales) +
    theme(plot.title = element_text(hjust = .5),
          rect = element_blank())
  }
)
legend_gene <- get_legend(
  ps[[1]] + theme(legend.position = "bottom") +
    guides(col = guide_legend(title.position = "top"))
)
ps <- lapply(ps, function(p) {p + guides(col = "none")})
p6 <- plot_grid(
  plot_grid(plotlist = ps, ncol = 2, scale = .95,
            labels = c('a)', 'b)', 'c)', 'd)')),
  plot_grid(NULL, legend_gene, NULL, rel_widths = c(1, 1, 1), ncol = 3), 
  nrow = 2, rel_heights = c(4, .5)
)
p6

Heatmaps of genes DE between lineages

Below we show heatmaps of the genes DE between lineages The DE genes in the heatmaps are ordered according to a hierarchical clustering on the first lineage.

### based on mean smoother
yhatSmooth <- predictSmooth(fibrosis, gene = patternGenes, nPoints = 50, tidy = TRUE) %>%
  mutate(yhat = log1p(yhat)) %>%
  group_by(gene) %>%
  mutate(yhat = scales::rescale(yhat)) %>%
  filter(condition == "ILD") %>%
  select(-condition) %>%
  ungroup()
heatSmooth_Lineage1 <- pheatmap(
  yhatSmooth %>%
    filter(lineage == 1) %>%
    select(-lineage) %>%
    arrange(-time) %>%
    pivot_wider(names_from = time, values_from = yhat) %>%
    select(-gene) %>%
    `rownames<-`(patternGenes),
  cluster_cols = FALSE, show_rownames = FALSE, show_colnames = FALSE,
  main = "Lineage 1", legend = FALSE, silent = TRUE
)

heatSmooth_Lineage2 <- pheatmap(
  (yhatSmooth %>%
    filter(lineage == 2) %>%
    select(-lineage) %>%
    arrange(-time) %>%
    pivot_wider(names_from = time, values_from = yhat) %>%
    select(-gene))[heatSmooth_Lineage1$tree_row$order, ],
  cluster_cols = FALSE, cluster_rows = FALSE,
  show_rownames = FALSE, show_colnames = FALSE, main = "Lineage 2",
  legend = FALSE, silent = TRUE 
)

p7 <- plot_grid(heatSmooth_Lineage1[[4]], heatSmooth_Lineage2[[4]], ncol = 2)
p7

Gene set enrichment analysis

This is done using the fgsea package [@Korotkevich2016].

geneSets <- msigdbr(species = "Homo sapiens", category = "C5", subcategory = "BP") %>%
  mutate(gene_symbol = toupper(gene_symbol)) %>%
  filter(gene_symbol %in% names(fibrosis)) 
m_list <- geneSets %>% split(x = .$gene_symbol, f = .$gs_name)
statsPat <- patRes$waldStat
names(statsPat) <- rownames(patRes)
eaRes <- fgsea(pathways = m_list, stats = statsPat, nperm = 5e4, minSize = 10)
eaRes <- eaRes %>% arrange(pval)
knitr::kable(head(eaRes[, 1:3], n = 20))

Session info

sessionInfo()

References

HectorRDB/condimentsPaper documentation built on Oct. 16, 2021, 7:02 p.m.