ProjecTILs系列教程（三）：不同器官的T细胞异质性

说在前面

在慢性感染和肿瘤生长过程中，长期的抗原刺激导致T细胞逐渐失去效应功能，这一过程通常被称为T细胞的衰竭。淋巴细胞脉络丛脑膜炎病毒(LCMV)模型是研究得最好的小鼠病毒感染模型系统之一，在阐明T细胞衰竭生物学方面发挥了重要作用。在慢性LCMV感染的背景下，大多数研究集中在来自脾脏的病毒特异性T细胞。

在Sandu等人最近的研究中，应用单细胞测序技术研究了六种不同组织(脾脏、血液、骨髓、淋巴结、肝脏和肺)中的CD8 T细胞多样性，以确定组织微环境如何影响和塑造T细胞表型。他们发现T细胞采用组织特异性的转录组学谱，在不同器官中特定基因的表达存在差异，例如与抗原接触和TCR激活相关的基因。

在本期的案例演示中，Immugent将通过实操使用ProjecTILs重新分析Sandu等人的单细胞数据，并在LCMV参考图谱的背景下研究组织特异性T细胞的异质性。

---

代码实现

library(ggplot2)
library(reshape2)
library(gridExtra)
library(ProjecTILs)

ddir <- "input/P14_tissues"

tissues <- c("Blood", "Spleen", "BM", "LN", "Lung", "Liver")
folders <- c("1_blood", "2_spleen", "3_BM", "4_LN", "5_Lung", "6_Liver")

query.list <- list()

for (i in seq_along(folders)) {

    mpath <- sprintf("%s/%s", ddir, folders[i])
    query.list[[i]] <- read.sc.query(mpath, type = "10x", min.cells = 3, min.features = 50)
    query.list[[i]] <- RenameCells(query.list[[i]], add.cell.id = paste0("S", i))
    query.list[[i]]$Tissue <- tissues[i]
}

query.merged <- Reduce(merge, query.list)

percent.ribo.dv <- PercentageFeatureSet(query.merged, pattern = "^Rp[ls]")
percent.mito.dv <- PercentageFeatureSet(query.merged, pattern = "^mt-")

query.merged <- AddMetaData(query.merged, metadata = percent.ribo.dv, col.name = "percent.ribo")
query.merged <- AddMetaData(query.merged, metadata = percent.mito.dv, col.name = "percent.mito")

Idents(query.merged) <- "Tissue"
pl <- VlnPlot(query.merged, features = c("nFeature_RNA", "nCount_RNA", "percent.ribo",
    "percent.mito"), ncol = 2, pt.size = 0.001, combine = FALSE, split.by = "Tissue")

pll <- list()
for (i in seq_along(pl)) {
    pll[[i]] = pl[[i]] + theme(axis.text.x = element_blank())
}

plot(do.call("arrangeGrob", c(pll, ncol = 2)))

image.png

query.merged <- subset(query.merged, subset = nFeature_RNA > 500 & nFeature_RNA <
    4000 & nCount_RNA > 1000 & nCount_RNA < 15000 & percent.ribo < 50 & percent.mito <
    5)
dim(query.merged)
#[1] 17505 23902

ref <- load.reference.map("ref_LCMV_Atlas_mouse_v1.rds")

# reproduce cluster colors
library(scales)
functional.cluster.colors <- hue_pal()(7)
functional.cluster.colors <- functional.cluster.colors[c(4, 3, 2, 5, 1, 6, 7)]
names(functional.cluster.colors) <- levels(ref$functional.cluster)


DimPlot(ref, reduction = "umap", label = TRUE, pt.size = 0.5, group.by = "functional.cluster",
    dims = c(2, 1), cols = functional.cluster.colors) + NoLegend() + theme(aspect.ratio = 1) +
    scale_x_reverse() + scale_y_reverse() + ggtitle("Reference CD8 LCMV atlas")

image.png

# For example, to allow 3GB of RAM, and run on 6 computing cores:
mem_in_mb <- 3000
ncores <- 6

querydata <- SplitObject(query.merged, split.by = "Tissue")
# For serial processing, just set ncores=1
query.projected <- make.projection(querydata, ref = ref, ncores = ncores, future.maxSize = mem_in_mb)

pll <- list()
for (i in seq_along(query.projected)) {
    s <- names(query.projected)[i]
    query.projected[[s]] <- cellstate.predict(ref = ref, query = query.projected[[s]])
    pll[[i]] <- plot.projection(ref, query.projected[[s]], pointsize = 0.3, linesize = 0.3,
        cols = functional.cluster.colors) + theme_bw() + theme(aspect.ratio = 0.8) +
        scale_x_reverse() + scale_y_reverse() + coord_flip() + NoLegend() + ggtitle(s)
}
g <- do.call("arrangeGrob", c(pll, ncol = 3))
plot(g)

image.png

# Reorder colors (to resemble order in Sandu Figure 2J)
cols_use <- functional.cluster.colors[c(3, 1, 4, 2, 7, 6, 5)]
states_all <- levels(factor(names(functional.cluster.colors), levels = names(cols_use)))

m <- matrix(nrow = length(names(query.projected)), ncol = length(states_all))
rownames(m) <- names(query.projected)
colnames(m) <- states_all
for (i in seq_along(query.projected)) {
    tb <- table(factor(query.projected[[i]]$functional.cluster, levels = states_all))
    m[i, ] <- tb * 100/sum(tb)
}

melt <- melt(m)
colnames(melt) <- c("Tissue", "Cell_state", "Percent")

p <- ggplot(melt, aes(x = Tissue, y = Percent, fill = Cell_state)) + geom_bar(stat = "identity",
    position = "stack") + scale_fill_manual(values = functional.cluster.colors) +
    theme_light() + theme(legend.position = "right")
p

image.png

genes4radar = c("Cd8a", "Tcf7", "Ccr7", "Gzmb", "Slamf6", "Pdcd1", "Havcr2", "Tox",
    "Entpd1", "Cx3cr1", "Cxcr6", "Xcl1", "Mki67")

g <- plot.states.radar(ref, query = query.projected, min.cells = 200, genes4radar = genes4radar,
    return = T)
plot(g)

image.png

# BiocManager::install('EnhancedVolcano')
library(EnhancedVolcano)

deg <- find.discriminant.genes(ref, query.projected[["Blood"]], query.control = query.projected[["Spleen"]],
    state = "SLEC", min.pct = 0.1, logfc.threshold = 0.1, query.assay = "RNA")
EnhancedVolcano(deg, lab = rownames(deg), x = "avg_log2FC", y = "p_val", pCutoff = 1e-09,
    FCcutoff = 0.5, labSize = 5, legendPosition = "none", drawConnectors = F, title = "Blood vs. Spleen (SLEC)")

image.png

deg2 <- find.discriminant.genes(ref, query.projected[["Liver"]], query.control = query.projected[["Spleen"]],
    state = "Tex", min.pct = 0.1, logfc.threshold = 0.1)
EnhancedVolcano(deg2, lab = rownames(deg2), x = "avg_log2FC", y = "p_val", pCutoff = 1e-09,
    FCcutoff = 0.5, labSize = 5, legendPosition = "none", drawConnectors = F, title = "Liver vs. Spleen (Tex)")

image.png

---

小结

在Sandu等人的研究中，作者使用基于无监督聚类、分类和差异表达的传统方法，描述了CD8 T细胞在多个组织中的异质性。细胞群的定义，包括批效应与组织相关的生物学差异的考虑，有意义的细胞类型的注释，差异表达分析以确定亚型间差异和组织间差异，所有这些都需要大量的系统生物学背景和专业知识。

通过这个案例研究，小编展示了投射如何用最少的努力和领域专业知识得到非常相似的结果。通过分析结果我们可以发现，ProjecTILs预测的T细胞亚型的组织特异性组成，与原始研究以无监督的方式定义的亚型有很好的相关性，可以检测到与特定组织和T细胞亚型相关的特定基因和基因模块。