Cell cycle effect

Assign each cell a score, based on its expression of G2/M and S phase markers as described in the Seurat workflow here.

s.genes <- cc.genes.updated.2019$s.genes
g2m.genes <- cc.genes.updated.2019$g2m.genes

seu <- CellCycleScoring(seu, s.features = s.genes, g2m.features = g2m.genes, 
                        set.ident = TRUE)

PCA of cell cycle genes.

DimPlot(seu, group.by = "Phase") -> p1

seu %>%
  RunPCA(features = c(s.genes, g2m.genes),
                      dims = 1:30, verbose = FALSE) %>%
  DimPlot(reduction = "pca") -> p2

(p2 / p1) + plot_layout(guides = "collect")

Past versions of unnamed-chunk-5-1.png

Distribution of cell cycle markers.

# Visualize the distribution of cell cycle markers across
RidgePlot(seu, features = c("PCNA", "TOP2A", "MCM6", "MKI67"), ncol = 2,
          log = TRUE)

Past versions of unnamed-chunk-6-1.png

Using the Seurat Alternate Workflow from here, calculate the difference between the G2M and S phase scores so that signals separating non-cycling cells and cycling cells will be maintained, but differences in cell cycle phase among proliferating cells (which are often uninteresting), can be regressed out of the data.

seu$CC.Difference <- seu$S.Score - seu$G2M.Score

Integrate RNA data

Split by batch for integration. Normalise with SCTransform. Increase the strength of alignment by increasing k.anchor parameter to 20 as recommended in Seurat Fast integration with RPCA vignette.

First, integrate the RNA data.

out <- here("data",
            "C133_Neeland_merged",
            glue("C133_Neeland_full_clean{ambient}_integrated_t_cells.SEU.rds"))

if(!file.exists(out)){
  DefaultAssay(seu) <- "RNA"
  VariableFeatures(seu) <- NULL
  seu[["pca"]] <- NULL
  seu[["umap"]] <- NULL
  
  seuLst <- SplitObject(seu, split.by = "Batch")
  rm(seu)
  gc()
  
  # normalise with SCTransform and regress out cell cycle score difference
  seuLst <- lapply(X = seuLst, FUN = SCTransform, method = "glmGamPoi",
                   vars.to.regress = "CC.Difference")
  # integrate RNA data
  features <- SelectIntegrationFeatures(object.list = seuLst,
                                        nfeatures = 3000)
  seuLst <- PrepSCTIntegration(object.list = seuLst, anchor.features = features)
  seuLst <- lapply(X = seuLst, FUN = RunPCA, features = features)
  anchors <- FindIntegrationAnchors(object.list = seuLst,
                                    normalization.method = "SCT",
                                    anchor.features = features,
                                    k.anchor = 20,
                                    dims = 1:30, reduction = "rpca")
  seu <- IntegrateData(anchorset = anchors, 
                       k.weight = min(100, min(sapply(seuLst, ncol)) - 5),
                       normalization.method = "SCT",
                       dims = 1:30)
  
  DefaultAssay(seu) <- "integrated"
  seu <- RunPCA(seu, dims = 1:30, verbose = FALSE) %>%
    RunUMAP(dims = 1:30, verbose = FALSE)

  saveRDS(seu, file = out)
  fs::file_chmod(out, "664")
  if(any(str_detect(fs::group_ids()$group_name, 
                    "oshlack_lab"))) fs::file_chown(out, 
                                                    group_id = "oshlack_lab")
  
} else {
  seu <- readRDS(file = out)
  
}

Integrate ADT data

out <- here("data",
            "C133_Neeland_merged",
            glue("C133_Neeland_full_clean{ambient}_integrated_t_cells.ADT.SEU.rds"))

# get ADT meta data
read.csv(file = here("data",
                     "C133_Neeland_batch1",
                     "data",
                     "sample_sheets",
                     "ADT_features.csv")) -> adt_data 
# cleanup ADT meta data
pattern <- "anti-human/mouse |anti-human/mouse/rat |anti-mouse/human "
adt_data$name <- gsub(pattern, "", adt_data$name)
# change ADT rownames to antibody names
DefaultAssay(seu) <- "ADT"
if(all(rownames(seu[["ADT"]]@counts) == adt_data$id)){
  adt_counts <- seu[["ADT"]]@counts
  rownames(adt_counts) <- adt_data$name
  seu[["ADT"]] <- CreateAssayObject(counts = adt_counts)
}

if(!file.exists(out)){
  tmp <- DietSeurat(subset(seu, cells = which(seu$Batch != 0)), 
                    assays = "ADT")
  DefaultAssay(tmp) <- "ADT"
  
  seuLst <- SplitObject(tmp, split.by = "Batch")
  seuLst <- lapply(X = seuLst, FUN = function(x) {
    # set all ADT as variable features
    VariableFeatures(x) <- rownames(x)
    x <- NormalizeData(x, normalization.method = "CLR", margin = 2)
    x
  })
  features <- SelectIntegrationFeatures(object.list = seuLst)
  seuLst <- lapply(X = seuLst, FUN = function(x) {
    x <- ScaleData(x, features = features, verbose = FALSE) %>%
      RunPCA(features = features, verbose = FALSE)
    x
  })
  anchors <- FindIntegrationAnchors(object.list = seuLst, reduction = "rpca",
                                    dims = 1:30)
  tmp <- IntegrateData(anchorset = anchors, dims = 1:30)
  
  DefaultAssay(tmp) <- "integrated"
  tmp <- ScaleData(tmp) %>%
    RunPCA(dims = 1:30, verbose = FALSE) %>%
    RunUMAP(dims = 1:30, verbose = FALSE)

  # create combined object that only contains cells with RNA+ADT data
  seuADT <- subset(seu, cells = which(seu$Batch !=0))
  seuADT[["integrated.adt"]] <- tmp[["integrated"]]
  seuADT[["pca.adt"]] <- tmp[["pca"]]
  seuADT[["umap.adt"]] <- tmp[["umap"]]

  saveRDS(seuADT, file = out)
  fs::file_chmod(out, "664")
  if(any(str_detect(fs::group_ids()$group_name, 
                    "oshlack_lab"))) fs::file_chown(out, 
                                                    group_id = "oshlack_lab")
  
} else {
  seuADT <- readRDS(file = out)
  
}

View integrated data

DefaultAssay(seuADT) <- "integrated"

DimPlot(seu, group.by = "Batch", reduction = "umap") -> p1
DimPlot(seuADT, group.by = "Batch", reduction = "umap") -> p2
DimPlot(seuADT, group.by = "Batch", reduction = "umap.adt") -> p3

(p1 / ((p2 | p3) +
  plot_layout(guides = "collect"))) &
  theme(axis.title = element_text(size = 8),
        axis.text = element_text(size = 8)) 

Past versions of unnamed-chunk-10-1.png

DimPlot(seu, group.by = "Phase", reduction = "umap") -> p1
DimPlot(seuADT, group.by = "Phase", reduction = "umap") -> p2
DimPlot(seuADT, group.by = "Phase", reduction = "umap.adt") -> p3

(p1 / ((p2 | p3) +
  plot_layout(guides = "collect"))) &
  theme(axis.title = element_text(size = 8),
        axis.text = element_text(size = 8)) 

Past versions of unnamed-chunk-11-1.png

Dimensionality reduction (RNA)

Exclude any mitochondrial, ribosomal, immunoglobulin and HLA genes from variable genes list, to encourage clustering by cell type.

# remove HLA, immunoglobulin, RNA, MT, and RP genes from variable genes list
var_regex = '^HLA-|^IG[HJKL]|^RNA|^MT-|^RP' 
hvg <- grep(var_regex, VariableFeatures(seuADT), invert = TRUE, value = TRUE)
# assign edited variable gene list back to object 
VariableFeatures(seuADT) <- hvg

# redo PCA and UMAP 
seuADT <- RunPCA(seuADT, dims = 1:30, verbose = FALSE) %>%
    RunUMAP(dims = 1:30, verbose = FALSE)

DimHeatmap(seuADT, dims = 1:30, cells = 500, balanced = TRUE,
           reduction = "pca", assays = "integrated")

Past versions of unnamed-chunk-12-1.png

ElbowPlot(seuADT, ndims = 30, reduction = "pca")

Past versions of unnamed-chunk-13-1.png

Dimensionality reduction (ADT)

DimHeatmap(seuADT, dims = 1:30, cells = 500, balanced = TRUE,
           reduction = "pca.adt", assays = "integrated.adt")

Past versions of unnamed-chunk-14-1.png

ElbowPlot(seuADT, ndims = 30, reduction = "pca.adt")

Past versions of unnamed-chunk-15-1.png

Run WNN clustering

Perform clustering at a range of resolutions and visualise to see which is appropriate to proceed with.

out <- here("data",
            "C133_Neeland_merged",
            glue("C133_Neeland_full_clean{ambient}_integrated_clustered_t_cells.ADT.SEU.rds"))

if(!file.exists(out)){
  DefaultAssay(seuADT) <- "integrated"
  seuADT <- FindMultiModalNeighbors(seuADT, reduction.list = list("pca", "pca.adt"), 
                                 dims.list = list(1:30, 1:10), 
                                 modality.weight.name = "RNA.weight")
  seuADT <- FindClusters(seuADT, algorithm = 3, 
                      resolution = seq(0.1, 1, by = 0.1),
                      graph.name = "wsnn")
  seuADT <- RunUMAP(seuADT, dims = 1:30, nn.name = "weighted.nn", 
                 reduction.name = "wnn.umap", reduction.key = "wnnUMAP_",
                 return.model = TRUE)
  saveRDS(seuADT, file = out)
  fs::file_chmod(out, "664")
  if(any(str_detect(fs::group_ids()$group_name, 
                    "oshlack_lab"))) fs::file_chown(out, 
                                                    group_id = "oshlack_lab")
  
} else {
  seuADT <- readRDS(file = out)
  
}

clustree::clustree(seuADT, prefix = "wsnn_res.")

Past versions of unnamed-chunk-16-1.png

View clusters

Choose most appropriate resolution based on clustree plot above.

grp <- "wsnn_res.0.6"
# change factor ordering
seuADT@meta.data[,grp] <- fct_inseq(seuADT@meta.data[,grp])

DimPlot(seuADT, group.by = grp, label = T) + 
  theme(legend.position = "bottom")

Past versions of unnamed-chunk-17-1.png

Weighting of RNA and ADT data per cluster.

 VlnPlot(seuADT, features = "integrated.weight", group.by = grp, sort = TRUE, 
         pt.size = 0.1) +
  NoLegend()

Past versions of unnamed-chunk-18-1.png

Map data

Find transfer anchors.

# use WNN clustered batches 1-6 as reference
reference <- seuADT
DefaultAssay(seu) <- "RNA"
# batch 0 RNA data is the query
query <- DietSeurat(subset(seu, cells = which(seu$Batch == 0)),
                    assays = "RNA")

DefaultAssay(reference) <- "integrated"

anchors <- FindTransferAnchors(
  reference = reference,
  query = query,
  normalization.method = "SCT",
  reference.reduction = "pca",
  dims = 1:50
)

Map batch 0 samples onto reference.

query <- MapQuery(
  anchorset = anchors,
  query = query,
  reference = reference,
  refdata = list(
    wsnn = grp,
    ADT = "ADT"
  ),
  reference.reduction = "pca", 
  reduction.model = "wnn.umap"
)

query

An object of class Seurat 
21753 features across 5080 samples within 3 assays 
Active assay: RNA (21568 features, 0 variable features)
 2 other assays present: prediction.score.wsnn, ADT
 2 dimensional reductions calculated: ref.pca, ref.umap

Visualise batch 0 samples on reference UMAP.

query$predicted.wsnn <- fct_inseq(query$predicted.wsnn)

DimPlot(query, reduction = "ref.umap", group.by = "predicted.wsnn", 
        label = TRUE, label.size = 3 ,repel = TRUE) + 
  theme(legend.position = "bottom")

Past versions of unnamed-chunk-21-1.png

Distribution of Azimuth prediction scores per WNN cluster.

ggplot(query@meta.data, aes(x = predicted.wsnn, 
                            y = predicted.wsnn.score,
                            fill = predicted.wsnn)) +
  geom_boxplot() + NoLegend()

Past versions of unnamed-chunk-22-1.png

Compute combined UMAP

Computing a new UMAP can help to identify any cell states present in the query but not reference.

# merge reference (integrated + WNN clustered) and query (RNA only samples)
reference$id <- 'reference'
query$id <- 'query'
DefaultAssay(reference) <- "integrated"
refquery <- merge(DietSeurat(reference,
                             assays = c("RNA","ADT","integrated","SCT"),
                             dimreducs = c("pca")),
                  DietSeurat(query,
                             assays = c("RNA"),
                             dimreducs = "ref.pca"))
refquery[["pca"]] <- merge(reference[["pca"]], query[["ref.pca"]])
refquery <- RunUMAP(refquery, reduction = 'pca', dims = 1:50, assay = "integrated")

View combined UMAP.

# combine cluster annotations from reference and query
refquery@meta.data[, grp] <- ifelse(is.na(refquery@meta.data[,grp]),
                          refquery$predicted.wsnn,
                          refquery@meta.data[,grp])
# change factor ordering
refquery@meta.data[,grp] <- fct_inseq(refquery@meta.data[,grp])

DimPlot(refquery, reduction = "umap", group.by = grp, 
             label = TRUE, label.size = 3) + NoLegend() 

Past versions of unnamed-chunk-24-1.png

DimPlot(refquery, reduction = "umap", group.by = "Phase", 
             label = FALSE, label.size = 3)

Past versions of unnamed-chunk-24-2.png

DimPlot(refquery, reduction = "umap", group.by = "Disease", 
             label = FALSE, label.size = 3)

Past versions of unnamed-chunk-24-3.png

Save results.

out <- here("data",
            "C133_Neeland_merged",
            glue("C133_Neeland_full_clean{ambient}_integrated_clustered_mapped_t_cells.ADT.SEU.rds"))
if(!file.exists(out)){
  saveRDS(refquery, file = out)
  fs::file_chmod(out, "664")
  if(any(str_detect(fs::group_ids()$group_name, 
                    "oshlack_lab"))) fs::file_chown(out, 
                                                    group_id = "oshlack_lab")
}

Examine combined clusters

Number of cells per cluster.

refquery@meta.data %>%
  ggplot(aes(x = !!sym(grp), fill = !!sym(grp))) +
  geom_bar() +
  geom_text(aes(label = ..count..), stat = "count",
            vjust = -0.5, colour = "black", size = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  NoLegend()

Past versions of unnamed-chunk-26-1.png

Visualise quality metrics by cluster. Cluster 17 potentially contains low quality cells.

refquery@meta.data %>%
  ggplot(aes(x = !!sym(grp),
             y = nCount_RNA,
             fill = !!sym(grp))) +
  geom_violin(scale = "area") +
  scale_y_log10() +
  NoLegend() -> p2

refquery@meta.data %>%
  ggplot(aes(x = !!sym(grp),
             y = nFeature_RNA,
             fill = !!sym(grp))) +
  geom_violin(scale = "area") +
  scale_y_log10() +
  NoLegend() -> p3

(p2 / p3) & theme(text = element_text(size = 8))

Past versions of unnamed-chunk-27-1.png

Check the batch composition of each of the clusters. Cluster 17 does not contain any cells from batch 0; could be a quality issue or the cell type was not captured in batch 0?

dittoBarPlot(refquery,
             var = "Batch", 
             group.by = grp)

Past versions of unnamed-chunk-28-1.png

Check the sample compositions of combined clusters.

dittoBarPlot(refquery,
             var = "sample.id", 
             group.by = grp) + ggtitle("Samples") +
  theme(legend.position = "bottom")

Past versions of unnamed-chunk-29-1.png

Test for marker genes using limma

# limma-trend for DE
Idents(refquery) <- grp

logcounts <- normCounts(DGEList(as.matrix(refquery[["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.
logcounts <- logcounts[!str_detect(rownames(logcounts), var_regex),]

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

clustgrp <- paste0("c", Idents(refquery))
clustgrp <- factor(clustgrp, levels = paste0("c", 0:maxclust))
donor <- factor(seu$sample.id)
batch <- factor(seu$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))

          c0    c1    c2    c3    c4    c5    c6    c7    c8    c9   c10   c11
Down     907  1314  3531  1496   483  1286   437   425  1618   597   230   500
NotSig 13972 14333 12505 13949 15180 14540 15214 15077 13403 15271 14702 14487
Up      1386   618   229   820   602   439   614   763  1244   397  1333  1278
         c12   c13   c14   c15   c16   c17   c18   c19   c20   c21
Down     486   291   297    65    95    17   105   222    98    66
NotSig 15424 15264 15804 15886 13276 16126 15855 15042 15705 15346
Up       355   710   164   314  2894   122   305  1001   462   853

Test relative to a threshold (TREAT).

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

          c0    c1    c2    c3    c4    c5    c6    c7    c8    c9   c10   c11
Down      20     8    24    24     4    38     9    23    97    12     1    32
NotSig 16215 16225 16237 16217 16229 16205 16224 16210 16152 16233 16240 16166
Up        30    32     4    24    32    22    32    32    16    20    24    67
         c12   c13   c14   c15   c16   c17   c18   c19   c20   c21
Down       9    13     5     2    11     0     1    21    13     7
NotSig 16224 16217 16254 16200 15959 16265 16245 16200 16208 16203
Up        32    35     6    63   295     0    19    44    44    55

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)
}

Past versions of unnamed-chunk-32-1.png

Past versions of unnamed-chunk-32-2.png

limma marker gene dotplot

DefaultAssay(refquery) <- "RNA"
contnames <- colnames(mycont)
top_markers <- NULL
n_markers <- 10

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::paletteer_d("pals::glasbey")[factor(names(top_markers))]

DotPlot(refquery,    
        features = unname(top_markers),
        group.by = grp,
        cols = c("azure1", "blueviolet"),
        dot.scale = 3, assay = "SCT") +
    RotatedAxis() +
    FontSize(y.text = 8, x.text = 12) +
    labs(y = element_blank(), x = element_blank()) +
    coord_flip() +
  theme(axis.text.y = element_text(color = cols)) +
  ggtitle("Top 10 cluster marker genes (no duplicates)")

Past versions of unnamed-chunk-33-1.png

Save marker genes and pathways

The Broad MSigDB Reactome pathways are tested for each contrast using cameraPR from limma. The cameraPR method tests whether a set of genes is highly ranked relative to other genes in terms of differential expression, accounting for inter-gene correlation.

Prepare gene sets of interest.

if(!file.exists(here("data/Hs.c2.cp.reactome.v7.1.entrez.rds")))
  download.file("https://bioinf.wehi.edu.au/MSigDB/v7.1/Hs.c2.cp.reactome.v7.1.entrez.rds",
                here("data/Hs.c2.cp.reactome.v7.1.entrez.rds"))

Hs.c2.reactome <- readRDS(here("data/Hs.c2.cp.reactome.v7.1.entrez.rds"))

gns <- AnnotationDbi::mapIds(org.Hs.eg.db, 
                             keys = rownames(tr), 
                             column = c("ENTREZID"),
                             keytype = "SYMBOL",
                             multiVals = "first")

Run pathway analysis and save results to file.

options(scipen=-1, digits = 6)
contnames <- colnames(mycont)
dirName <- here("output",
                "cluster_markers",
                glue("RNA{ambient}"),
                "t_cells")

if(!dir.exists(dirName)) dir.create(dirName, recursive = TRUE)

for(c in colnames(tr)){
  top <- topTreat(tr, coef = c, n = Inf)
  top <- top[top$logFC > 0, ]

  write.csv(top[1:100, ] %>%
              rownames_to_column(var = "Symbol"),
            file = glue("{dirName}/up-cluster-limma-{c}.csv"),
            sep = ",",
            quote = FALSE,
            col.names = NA,
            row.names = TRUE)

  # get marker indices
  c2.id <- ids2indices(Hs.c2.reactome, unname(gns[rownames(tr)]))
  # gene set testing results
  cameraPR(tr$t[,glue("{c}")], c2.id) %>%
    rownames_to_column(var = "Pathway") %>%
    dplyr::filter(Direction == "Up") %>%
    slice_head(n = 50) %>%
    write.csv(file = here(glue("{dirName}/REACTOME-cluster-limma-{c}.csv")),
            sep = ",",
            quote = FALSE,
            col.names = NA,
            row.names = TRUE)
}

Find all marker ADT using limma

# identify isotype controls for DSB ADT normalisation
read_csv(file = here("data",
                     "C133_Neeland_batch1",
                     "data",
                     "sample_sheets",
                     "ADT_features.csv")) %>%
  dplyr::filter(grepl("[Ii]sotype", name)) %>%
  pull(name) -> isotype_controls

# normalise ADT using DSB normalisation
adt <- seuADT[["ADT"]]@counts
adt_dsb <- ModelNegativeADTnorm(cell_protein_matrix = adt,
                                denoise.counts = TRUE,
                                use.isotype.control = TRUE,
                                isotype.control.name.vec = isotype_controls)

[1] "fitting models to each cell for dsb technical component and removing cell to cell technical noise"

Running the limma analysis on the normalised counts.

# limma-trend for DE
Idents(seuADT) <- grp

logcounts <- adt_dsb
# remove isotype controls from marker analysis
logcounts <- logcounts[!rownames(logcounts) %in% isotype_controls,]
maxclust <- length(levels(Idents(seuADT))) - 1

clustgrp <- paste0("c", Idents(seuADT))
clustgrp <- factor(clustgrp, levels = paste0("c", 0:maxclust))
donor <- seuADT$sample.id

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))

        c0  c1  c2  c3  c4  c5  c6  c7  c8  c9 c10 c11 c12 c13 c14 c15 c16 c17
Down    28  35  27  41  31  41  17  19  38  16   5  21  34  34  10   3   5   4
NotSig  90  99 121 107 112 100 111 121 108 126  84 104 120 114 138 141 111 139
Up      36  20   6   6  11  13  26  14   8  12  65  29   0   6   6  10  38  11
       c18 c19 c20 c21
Down    22  23  22   0
NotSig 125 122 125 151
Up       7   9   7   3

Test relative to a threshold (TREAT).

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

        c0  c1  c2  c3  c4  c5  c6  c7  c8  c9 c10 c11 c12 c13 c14 c15 c16 c17
Down     4   2   1   9   4   6   0   5   4   0   0   7   1  14   0   0   0   0
NotSig 136 141 153 141 144 145 148 144 150 149 114 131 153 136 151 152 145 148
Up      14  11   0   4   6   3   6   5   0   5  40  16   0   4   3   2   9   6
       c18 c19 c20 c21
Down     6  14   8   0
NotSig 142 138 141 153
Up       6   2   5   1

ADT marker dot plot

Dot plot of the top 5 ADT markers per cluster without duplication.

contnames <- colnames(mycont)
top_markers <- NULL
n_markers <- 5

for (i in 1:length(contnames)){
  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::paletteer_d("pals::glasbey")[factor(names(top_markers))][!duplicated(top_markers)]

# add DSB normalised data to Seurat assay for plotting
seuADT[["ADT.dsb"]] <- CreateAssayObject(data = logcounts)
DotPlot(seuADT, 
        group.by = grp,
        features = unname(top_markers), 
        cols = c("azure1", "blueviolet"),
        assay = "ADT.dsb")  +
  RotatedAxis() + 
  FontSize(y.text = 8, x.text = 9) +
  labs(y = element_blank(), x = element_blank()) +
  theme(axis.text.y = element_text(color = cols),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10)) +
    coord_flip() +
  ggtitle("Top 5 cluster markers ADTs (no duplicates)")

Past versions of unnamed-chunk-39-1.png

ADT marker heatmap

Make data frame of proteins, clusters, expression levels.

cbind(seuADT@meta.data %>%
        dplyr::select(!!sym(grp)),
      as.data.frame(t(seuADT@assays$ADT.dsb@data))) %>%
  rownames_to_column(var = "cell") %>%
  pivot_longer(c(-!!sym(grp), -cell), 
               names_to = "ADT",
               values_to = "expression") %>%
  dplyr::group_by(!!sym(grp), ADT) %>%
  dplyr::summarize(Expression = mean(expression)) %>%
  ungroup() -> dat

# plot expression density to select heatmap colour scale range
plot(density(dat$Expression))

Past versions of unnamed-chunk-40-1.png

dat %>%
  dplyr::filter(ADT %in% top_markers) |>
  heatmap(
    .column = !!sym(grp),
    .row = ADT,
    .value = Expression,
    row_order = top_markers, 
    scale = "none",
    rect_gp = grid::gpar(col = "white", lwd = 1),
    show_row_names = TRUE,
    cluster_columns = FALSE,
    cluster_rows = FALSE,
    column_names_gp = grid::gpar(fontsize = 10),
    column_title_gp = grid::gpar(fontsize = 12),
    row_names_gp = grid::gpar(fontsize = 8, col = cols[order(top_markers)]),
    row_title_gp = grid::gpar(fontsize = 12),
    column_title_side = "top",
    palette_value = circlize::colorRamp2(seq(-0.5, 2.5, length.out = 256),
                                         viridis::magma(256)),
    heatmap_legend_param = list(direction = "vertical"))

Past versions of unnamed-chunk-41-1.png

Save ADT markers

options(scipen=-1, digits = 6)
contnames <- colnames(mycont)
dirName <- here("output",
                "cluster_markers",
                glue("ADT{ambient}"),
                "t_cells")
if(!dir.exists(dirName)) dir.create(dirName, recursive = TRUE)

for(c in contnames){
  top <- topTreat(tr, coef = c, n = Inf)
  top <- top[top$logFC > 0, ]
  write.csv(top,
            file = glue("{dirName}/up-cluster-limma-{c}.csv"),
            sep = ",",
            quote = FALSE,
            col.names = NA,
            row.names = TRUE)
}