Update group labels

seuInt@meta.data %>%
  data.frame %>%
  mutate(Status = ifelse(str_detect(Treatment, "ivacaftor"),
                         "CF ivacaftor",
                         ifelse(str_detect(Treatment, "orkambi"),
                                "CF lumacaftor-ivacaftor",
                                ifelse(Treatment == "untreated",
                                       "CF no-modulator",
                                       "non-CF control"))),
         Status_sub = ifelse(str_detect(Treatment, "ivacaftor"),
                         "CF.IVA",
                         ifelse(str_detect(Treatment, "orkambi"),
                                "CF.LUMA_IVA",
                                ifelse(Treatment == "untreated",
                                       "CF.NO_MOD",
                                       "NON_CF.CTRL"))),
         Group = ifelse(!Status_sub %in% "NON_CF.CTRL", 
                        paste(Status_sub, 
                              toupper(substr(Severity, 1, 1)),
                              sep = "."), 
                        Status_sub),
         Severity = tolower(Severity),
         Participant = strsplit2(sample.id, ".", fixed = TRUE)[,1]) -> seuInt@meta.data

Load manual annotations

labels <- read_excel(here("data",
                          "cluster_annotations",
                          "macrophages_ambientRNAremoval_01.05.24.xlsx"))

# set selected cluster resolution
grp <- "integrated_snn_res.0.6"
seuInt@meta.data %>%
  rownames_to_column(var = "cell") %>%
  left_join(labels %>%
              mutate(Cluster = as.factor(Cluster),
                     Annotation = as.factor(Annotation),
                     Broad = as.factor(Broad)),
            by = c("integrated_snn_res.0.6" = "Cluster")) %>%
  column_to_rownames(var = "cell") -> seuInt@meta.data

seuInt$Annotation <- fct_drop(seuInt$Annotation)
seuInt$Broad <- fct_drop(seuInt$Broad)
seuInt

An object of class Seurat 
41688 features across 141121 samples within 5 assays 
Active assay: integrated (3000 features, 3000 variable features)
 4 other assays present: RNA, ADT, ADT.dsb, SCT
 2 dimensional reductions calculated: pca, umap

Visualise annotations

options(ggrepel.max.overlaps = Inf)
DimPlot(seuInt, reduction = 'umap', label = TRUE, repel = TRUE, 
        label.size = 3, group.by = "integrated_snn_res.0.6") + 
  NoLegend() -> p1

DimPlot(seuInt, reduction = 'umap', label = FALSE, group.by = "Annotation") + 
  scale_color_paletteer_d("miscpalettes::pastel") + 
  theme(text = element_text(size = 9),
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  NoLegend() -> p2

DimPlot(seuInt, reduction = 'umap', label = FALSE, group.by = "Broad") + 
  NoLegend() +
  scale_color_paletteer_d("miscpalettes::pastel") + 
  theme(text = element_text(size = 9),
        axis.text = element_blank(),
        axis.ticks = element_blank()) -> p3

p1

LabelClusters(p2, id = "Annotation",  fontface = "bold", repel = TRUE, size = 3.5)

LabelClusters(p3, id = "Broad",  fontface = "bold", repel = TRUE, size = 3.5)

seuInt@meta.data %>%
  ggplot(aes(x = Annotation, fill = Annotation)) +
  geom_bar() +
  geom_text(aes(label = after_stat(count)), stat = "count",
            vjust = -0.5, colour = "black", size = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  NoLegend() +
  scale_fill_paletteer_d("miscpalettes::pastel")

seuInt@meta.data %>%
  ggplot(aes(x = Annotation, fill = Status_sub)) +
  geom_bar(position = "dodge") +
  geom_text(aes(label = ..count..), stat = "count",
            vjust = -0.5, colour = "black", size = 2, 
            position=position_dodge(width=0.9)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
        legend.position = "bottom")

seuInt@meta.data %>% 
  count(Annotation) %>% 
  mutate(perc = round(n/sum(n)*100, 1)) %>%
  dplyr::rename(`Cell Label` = "Annotation", 
                `No. Cells` = n,
                `% Cells` = perc) %>%
  knitr::kable()

Test for marker genes using limma

# limma-trend for DE
Idents(seuInt) <- "Annotation"

logcounts <- normCounts(DGEList(as.matrix(seuInt[["RNA"]]@counts)),
                        log = TRUE, prior.count = 0.5)
entrez <- AnnotationDbi::mapIds(org.Hs.eg.db,
                                keys = rownames(logcounts),
                                column = c("ENTREZID"),
                                keytype = "SYMBOL",
                                multiVals = "first")
# remove genes without entrez IDs as these are difficult to interpret biologically
logcounts <- logcounts[!is.na(entrez),]
# remove confounding genes from counts table e.g. mitochondrial, ribosomal etc.
# remove HLA, immunoglobulin, RNA, MT, and RP genes from marker gene analysis
var_regex = '^HLA-|^IG[HJKL]|^RNA|^MT-|^RP' 
logcounts <- logcounts[!str_detect(rownames(logcounts), var_regex),]

maxclust <- length(levels(Idents(seuInt))) - 1

clustgrp <- seuInt$Annotation
clustgrp <- factor(clustgrp)
donor <- factor(seuInt$sample.id)
batch <- factor(seuInt$Batch)

design <- model.matrix(~ 0 + clustgrp + donor)
colnames(design)[1:(length(levels(clustgrp)))] <- levels(clustgrp)

# Create contrast matrix
mycont <- matrix(NA, ncol = length(levels(clustgrp)),
                 nrow = length(levels(clustgrp)))
rownames(mycont) <- colnames(mycont) <- levels(clustgrp)
diag(mycont) <- 1
mycont[upper.tri(mycont)] <- -1/(length(levels(factor(clustgrp))) - 1)
mycont[lower.tri(mycont)] <- -1/(length(levels(factor(clustgrp))) - 1)

# Fill out remaining rows with 0s
zero.rows <- matrix(0, ncol = length(levels(clustgrp)),
                    nrow = (ncol(design) - length(levels(clustgrp))))
fullcont <- rbind(mycont, zero.rows)
rownames(fullcont) <- colnames(design)

fit <- lmFit(logcounts, design)

fit.cont <- contrasts.fit(fit, contrasts = fullcont)
fit.cont <- eBayes(fit.cont, trend = TRUE, robust = TRUE)

summary(decideTests(fit.cont))

       macro-alveolar macro-CCL macro-interstitial macro-lipid
Down             5429      5266               5365        5631
NotSig           6722      8237               6586        7614
Up               3471      2119               3671        2377
       macro-monocyte-derived macro-MT macro-NFKB macro-proliferating
Down                     4538     4128       5052                3151
NotSig                   7296     8553       9659                6165
Up                       3788     2941        911                6306
       macro-regulation macro-vesicle macro-viral
Down               3209          5778        2760
NotSig             9098          8273        8989
Up                 3315          1571        3873

Test relative to a threshold (TREAT).

tr <- treat(fit.cont, lfc = 0.5)
dt <- decideTests(tr)
summary(dt)

       macro-alveolar macro-CCL macro-interstitial macro-lipid
Down                6        13                167          14
NotSig          15539     15540              15301       15567
Up                 77        69                154          41
       macro-monocyte-derived macro-MT macro-NFKB macro-proliferating
Down                       34        4          1                  17
NotSig                  15551    15539      15596               15307
Up                         37       79         25                 298
       macro-regulation macro-vesicle macro-viral
Down                732             8           2
NotSig            14694         15569       15508
Up                  196            45         112

Mean-difference (MD) plots per cluster.

par(mfrow=c(4,3))
par(mar=c(2,3,1,2))

for(i in 1:ncol(mycont)){
  plotMD(tr, coef = i, status = dt[,i], hl.cex = 0.5)
  abline(h = 0, col = "lightgrey")
  lines(lowess(tr$Amean, tr$coefficients[,i]), lwd = 1.5, col = 4)
}

limma marker gene dotplot

DefaultAssay(seuInt) <- "RNA"
contnames <- colnames(mycont)
top_markers <- NULL
n_markers <- 5

for(i in 1:ncol(mycont)){
  top <- topTreat(tr, coef = i, n = Inf)
  top <- top[top$logFC > 0, ]
  top_markers <- c(top_markers, 
                   setNames(rownames(top)[1:n_markers], 
                            rep(contnames[i], n_markers)))
}

top_markers <- top_markers[!is.na(top_markers)]
top_markers <- top_markers[!duplicated(top_markers)]
cols <- paletteer_d("miscpalettes::pastel")[factor(names(top_markers))]

DotPlot(seuInt,    
        features = unname(top_markers),
        group.by = "Annotation",
        cols = c("azure1", "blueviolet"),
        dot.scale = 2.5, 
        assay = "SCT") +
    RotatedAxis() +
    FontSize(x.text = 8, y.text = 8) +
    labs(y = element_blank(), x = element_blank()) +
    theme(axis.text.x = element_text(color = cols,
                                   angle = 90,
                                   hjust = 1,
                                   vjust = 0.5,
                                   face = "bold"),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 10))

Test for marker genes using Seurat

DefaultAssay(seuInt) <- "RNA"
Idents(seuInt) <- "Annotation"

out <- here("data/cluster_annotations/seurat_markers_macrophages.rds")

if(!file.exists(out)){
  # restrict genes to same set as for limma analysis
  markers <- FindAllMarkers(seuInt, only.pos = TRUE, 
                            features = rownames(logcounts))
  saveRDS(markers, file = out)

} else {
  markers <- readRDS(out)
  
}

head(markers) %>% knitr::kable()

Seurat marker gene dotplot

DefaultAssay(seuInt) <- "RNA"

maxGenes <- 5
markers %>%
    group_by(cluster) %>%
    top_n(n = maxGenes, wt = avg_log2FC) -> top5

sig <- top5$gene
geneCols <- paletteer_d("miscpalettes::pastel")[top5$cluster][!duplicated(sig)]

pal <- paletteer::paletteer_d("vapoRwave::cool")
DotPlot(seuInt,
        features = sig[!duplicated(sig)],
        group.by = "Annotation",
        cols = c("azure1", "blueviolet"),
        dot.scale = 2.5,
        assay = "SCT") +
    FontSize(x.text = 8, y.text = 8) +
    labs(y = element_blank(), x = element_blank()) +
    theme(axis.text.x = element_text(color = cols,
                                   angle = 90,
                                   hjust = 1,
                                   vjust = 0.5,
                                   face = "bold"),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 10))

–> –>