教程背景

這個(gè)例子演示了允許用戶使用Seurat分析和探索多模式數(shù)據(jù)的新特性。雖然這只是一個(gè)初始版本俏蛮，但我們很高興在未來為多模態(tài)數(shù)據(jù)集發(fā)布重要的新功能也祠。

在這里，我們分析了由CITE-seq產(chǎn)生的8,617個(gè)臍帶血單核細(xì)胞(CBMCs)的數(shù)據(jù)集琅束，我們同時(shí)測(cè)量了單細(xì)胞轉(zhuǎn)錄組和11個(gè)表面蛋白的表達(dá)，其水平通過DNA-barcode抗體來量化算谈。首先涩禀，我們加載兩個(gè)計(jì)數(shù)矩陣:一個(gè)用于RNA測(cè)量，另一個(gè)用于抗體衍生標(biāo)簽ADT然眼。你可以在這里下載ADT文件和RNA文件

1.加載RNA UMI矩陣

注意艾船，這個(gè)數(shù)據(jù)集還包含了大約5%的小鼠細(xì)胞，我們可以將其作為蛋白質(zhì)測(cè)量的陰性對(duì)照高每。因此屿岂，基因表達(dá)矩陣在每個(gè)基因的開頭都附加了HUMAN_或MOUSE_。

cbmc.rna <- as.sparse(read.csv(file = "Analysis/data/GSE100866_CBMC_8K_13AB_10X-RNA_umi.csv.gz", sep = ",", header = TRUE, row.names = 1))

為了讓生活更容易繼續(xù)下去鲸匿，我們將拋棄除前100個(gè)高度表達(dá)的老鼠基因外的所有基因爷怀，并從CITE-seq前綴中去除“HUMAN_”

# 這個(gè)函數(shù)有幾個(gè)默認(rèn)參數(shù)，ncontrols = 100带欢，prefix = "HUMAN_"运授，controls = "MOUSE_"
cbmc.rna <- CollapseSpeciesExpressionMatrix(cbmc.rna)

2.加載ADT UMI矩陣

cbmc.adt <- as.sparse(read.csv(file = "Analysis/data/GSE100866_CBMC_8K_13AB_10X-ADT_umi.csv.gz", sep = ",", header = TRUE, row.names = 1))

在向Seurat添加多模態(tài)數(shù)據(jù)時(shí)，可以使用重復(fù)的特性名稱乔煞。每一組模態(tài)數(shù)據(jù)(例如吁朦。RNA、ADT等)儲(chǔ)存在自己的Assay對(duì)象中渡贾。其中一個(gè)Assay對(duì)象被稱為“default assay”喇完，意思是它用于所有的分析和可視化。若要從非默認(rèn)Assay中提取數(shù)據(jù)，可以指定一個(gè)與Assay鏈接的鍵锦溪，用于提取特征不脯。要查看所有對(duì)象的所有鍵，請(qǐng)使用Key function刻诊。
最后防楷，我們觀察到CCR5、CCR7和CD10的低富集则涯，因此將它們從矩陣中去除(可選)复局。(關(guān)鍵是這三個(gè)蛋白得低富集怎么看出來的？)

cbmc.adt <- cbmc.adt[setdiff(rownames(x = cbmc.adt), c("CCR5", "CCR7", "CD10")), ]

3.設(shè)置一個(gè)Seurat對(duì)象粟判，并基于RNA表達(dá)對(duì)細(xì)胞進(jìn)行聚類

下面的步驟表示基于scRNA-seq數(shù)據(jù)對(duì)pbmc進(jìn)行快速聚類亿昏。有關(guān)單個(gè)步驟或更高級(jí)選項(xiàng)的詳細(xì)信息，請(qǐng)參見這里的PBMC集群指導(dǎo)教程档礁。

cbmc <- CreateSeuratObject(counts = cbmc.rna) # 創(chuàng)建seurat對(duì)象
cbmc <- NormalizeData(cbmc) # 標(biāo)準(zhǔn)化
cbmc <- FindVariableFeatures(cbmc) # 選擇高變基因
cbmc <- ScaleData(cbmc) # 歸一化
cbmc <- RunPCA(cbmc, verbose = FALSE) # 運(yùn)行PCA降維
ElbowPlot(cbmc, ndims = 50) # 確定取多少個(gè)PCs

Snipaste_2020-09-23_15-09-57.png

這里教程說選擇13個(gè)PCs角钩，但是看代碼選的是25個(gè)，蜜汁迷惑呻澜。

cbmc <- FindNeighbors(cbmc, dims = 1:25)
cbmc <- FindClusters(cbmc, resolution = 0.8)
cbmc <- RunTSNE(cbmc, dims = 1:25, method = "FIt-SNE")

# 找到定義每個(gè)cluster的標(biāo)記递礼，并使用它們對(duì)cluster進(jìn)行注釋，我們使用max.cells.per加快進(jìn)程
cbmc.rna.markers <- FindAllMarkers(cbmc, max.cells.per.ident = 100, min.diff.pct = 0.3, only.pos = TRUE)
head(cbmc.rna.markers)
             p_val avg_logFC pct.1 pct.2    p_val_adj cluster  gene
IL7R  4.280959e-21 1.1368139 0.872 0.274 8.776395e-17       0  IL7R
LTB   3.149773e-20 1.0246049 0.985 0.525 6.457351e-16       0   LTB
LDHB  8.557593e-20 0.9786093 0.964 0.527 1.754392e-15       0  LDHB
TRBC2 9.896619e-18 0.7739586 0.868 0.425 2.028906e-13       0 TRBC2
IL32  9.951232e-18 0.9599525 0.896 0.327 2.040102e-13       0  IL32
ITM2A 2.335897e-17 0.8866571 0.786 0.298 4.788822e-13       0 ITM2A

注意羹幸，為了簡(jiǎn)單起見脊髓，我們將兩個(gè)CD14+單核細(xì)胞簇(HLA-DR基因表達(dá)不同)和NK細(xì)胞簇(細(xì)胞周期階段不同)合并。
這里直接進(jìn)行了細(xì)胞類型注釋栅受，注釋得過程将硝？

new.cluster.ids <- c("Memory CD4 T", "CD14+ Mono", "Naive CD4 T", "NK", "CD14+ Mono", "Mouse", "B", 
    "CD8 T", "CD16+ Mono", "T/Mono doublets", "NK", "CD34+", "Multiplets", "Mouse", "Eryth", "Mk", 
    "Mouse", "DC", "pDCs")
names(new.cluster.ids) <- levels(cbmc)
cbmc <- RenameIdents(cbmc, new.cluster.ids)
DimPlot(cbmc, label = TRUE) + NoLegend()

Snipaste_2020-09-23_15-19-12.png

4.將蛋白質(zhì)表達(dá)水平添加到Seurat對(duì)象

Seurat v3.0允許您將來自多個(gè)分析的信息存儲(chǔ)在同一個(gè)對(duì)象中，只要數(shù)據(jù)是多模態(tài)的(在同一組細(xì)胞中收集)屏镊。您可以使用SetAssayData和GetAssayData訪問器函數(shù)來添加和從其他分析中獲取數(shù)據(jù)依疼。
我們將定義一個(gè)ADT試驗(yàn)，并為它存儲(chǔ)原始計(jì)數(shù)闸衫。
如果你對(duì)這些數(shù)據(jù)內(nèi)部如何存儲(chǔ)感興趣涛贯，您可以查看Assay類诽嘉，它在object . R中定義蔚出。注意，所有的單細(xì)胞表達(dá)數(shù)據(jù)虫腋，包括RNA數(shù)據(jù)骄酗，仍然存儲(chǔ)在Assay對(duì)象中，也可以使用GetAssayData訪問悦冀。

cbmc[["ADT"]] <- CreateAssayObject(counts = cbmc.adt)

現(xiàn)在我們可以重復(fù)我們通常在RNA上運(yùn)行的預(yù)處理(標(biāo)準(zhǔn)化和縮放)步驟趋翻，但是需要修改‘a(chǎn)ssay’的參數(shù)。對(duì)于CITE-seq數(shù)據(jù)盒蟆，我們不建議使用典型的對(duì)數(shù)標(biāo)準(zhǔn)化踏烙。相反师骗，我們使用中心對(duì)數(shù)比(CLR)歸一化，獨(dú)立計(jì)算每個(gè)特性讨惩。與原始版本相比辟癌，這是一個(gè)略微改進(jìn)的過程，我們將很快發(fā)布更高級(jí)的CITE-seq規(guī)范化版本荐捻。

# 這里的標(biāo)準(zhǔn)化方法是CLR
cbmc <- NormalizeData(cbmc, assay = "ADT", normalization.method = "CLR")
cbmc <- ScaleData(cbmc, assay = "ADT")

5.觀察RNA簇上的蛋白水平

您可以使用任何ADT標(biāo)記的名稱黍少。“adt_CD4”)处面，在FetchData厂置、FeaturePlot、RidgePlot魂角、FeatureScatter昵济、DoHeatmap或任何其他可視化特性中。
Q1：這里得ADTfeatures為什么突然多了三個(gè)字母adt_或颊？
Q2：這里得蛋白和RNA對(duì)應(yīng)得表達(dá)水平砸紊，從下圖可以看出來什么？

# in this plot, protein (ADT) levels are on top, and RNA levels are on the bottom
FeaturePlot(cbmc, features = c("adt_CD3", "adt_CD11c", "adt_CD8", "adt_CD16", "CD3E", "ITGAX", "CD8A", 
    "FCGR3A"), min.cutoff = "q05", max.cutoff = "q95", ncol = 4)

在這個(gè)圖中囱挑，蛋白質(zhì)(ADT)水平在上面醉顽，RNA水平在下面。

Snipaste_2020-09-23_15-31-28.png

山脊圖繪制

RidgePlot(cbmc, features = c("adt_CD3", "adt_CD11c", "adt_CD8", "adt_CD16"), ncol = 2)

Snipaste_2020-09-23_15-32-53.png

繪制ADT散點(diǎn)圖(如FACS的雙軸圖)平挑。注意游添，如果需要的話，你甚至可以使用HoverLocator和FeatureLocator來“gate”細(xì)胞通熄。

FeatureScatter(cbmc, feature1 = "adt_CD19", feature2 = "adt_CD3")

Snipaste_2020-09-23_15-35-09.png

查看蛋白質(zhì)和RNA之間的關(guān)系

FeatureScatter(cbmc, feature1 = "adt_CD3", feature2 = "CD3E")

image.png

查看T細(xì)胞中CD4和CD8的水平

tcells <- subset(cbmc, idents = c("Naive CD4 T", "Memory CD4 T", "CD8 T"))
FeatureScatter(tcells, feature1 = "adt_CD4", feature2 = "adt_CD8")

Snipaste_2020-09-23_15-39-57.png

讓我們看看原始的(非標(biāo)準(zhǔn)化的)ADT的count值唆涝。你可以看到這些值相當(dāng)高，特別是與RNA值相比唇辨。這是由于細(xì)胞中蛋白質(zhì)拷貝數(shù)顯著增加廊酣，這大大減少了ADT數(shù)據(jù)中的“drop-out”。
對(duì)于drop-out如何理解赏枚？

FeatureScatter(tcells, feature1 = "adt_CD4", feature2 = "adt_CD8", slot = "counts")

image.png

如果你仔細(xì)觀察亡驰，你會(huì)發(fā)現(xiàn)我們的CD8 T細(xì)胞群富含CD8 T細(xì)胞，但仍然包含許多CD4+ CD8- T細(xì)胞饿幅。這是因?yàn)槌跏糃D4和CD8 T細(xì)胞在轉(zhuǎn)錄上非常相似凡辱，CD4和CD8的RNA丟失水平相當(dāng)高。這說明了單獨(dú)從scRNA-seq數(shù)據(jù)中定義細(xì)微的免疫細(xì)胞差異的挑戰(zhàn)栗恩。

saveRDS(cbmc, file = "Analysis/cbmc.rds")

6.識(shí)別簇間差異表達(dá)的蛋白

對(duì)cluster進(jìn)行采樣透乾，每個(gè)cluster最多300個(gè)細(xì)胞(使小cluster的熱圖更容易看到)。

cbmc.small <- subset(cbmc, downsample = 300)

找到所有cluster的蛋白質(zhì)標(biāo)記，并繪制一個(gè)熱圖乳乌。

adt.markers <- FindAllMarkers(cbmc.small, assay = "ADT", only.pos = TRUE)
DoHeatmap(cbmc.small, features = unique(adt.markers$gene), assay = "ADT", angle = 90) + NoLegend()

image.png

你可以看到我們的未知細(xì)胞同時(shí)表達(dá)骨髓和淋巴標(biāo)記(在RNA水平上也是如此)捧韵。它們可能是應(yīng)該被丟棄的細(xì)胞團(tuán)塊(多胞胎)。我們現(xiàn)在也要移除老鼠細(xì)胞汉操。

cbmc <- subset(cbmc, idents = c("Multiplets", "Mouse"), invert = TRUE)

7.直接在蛋白質(zhì)水平上聚類

還可以直接在CITE-seq數(shù)據(jù)上運(yùn)行降維和基于圖的聚類纫版。

#因?yàn)槲覀儗V泛地使用ADT數(shù)據(jù)，所以我們將把默認(rèn)assay改為“ADT”assay客情。這將導(dǎo)致所有函數(shù)默認(rèn)使用ADT數(shù)據(jù)其弊，而不是要求我們每次都指定它
DefaultAssay(cbmc) <- "ADT"
cbmc <- RunPCA(cbmc, features = rownames(cbmc), reduction.name = "pca_adt", reduction.key = "pca_adt_", 
    verbose = FALSE)
DimPlot(cbmc, reduction = "pca_adt")

image.png

在我們?cè)贏DT級(jí)別上重新聚類數(shù)據(jù)之前，我們將稍后隱藏RNA集群id膀斋。

cbmc[["rnaClusterID"]] <- Idents(cbmc)

現(xiàn)在梭伐，我們僅在ADT(蛋白質(zhì))水平上使用PCA重新計(jì)算tSNE。

cbmc <- RunTSNE(cbmc, dims = 1:10, reduction = "pca_adt", reduction.key = "adtTSNE_", reduction.name = "tsne_adt")
cbmc <- FindNeighbors(cbmc, features = rownames(cbmc), dims = NULL)
cbmc <- FindClusters(cbmc, resolution = 0.2, graph.name = "ADT_snn")

我們可以比較RNA和蛋白質(zhì)的聚類仰担，并使用它來注釋蛋白質(zhì)聚類糊识。
(當(dāng)然，我們也可以使用findmarker)

clustering.table <- table(Idents(cbmc), cbmc$rnaClusterID)
clustering.table
     Memory CD4 T CD14+ Mono Naive CD4 T   NK    B CD8 T CD16+ Mono T/Mono doublets CD34+ Eryth   Mk   DC
  0          1754          0        1217   29    0    27          0               5     2     4   24    1
  1             0       2189           0    4    0     0         30               1     1     5   25   55
  2             3          0           2  890    3     1          0               0     1     3    7    2
  3             0          4           0    2  319     0          2               0     2     2    3    0
  4            24          0          18    4    1   243          0               0     0     1    2    0
  5             1         27           4  157    2     2         10              56     0     9   16    6
  6             4          5           0    1    0     0          0               1   113    81   16    5
  7             4         59           4    0    0     0          9             117     0     0    2    0
  8             0          9           0    2    0     0        179               0     0     0    1    0
  9             0          0           1    0    0     0          0               0     0     0    0    1
  10            1          0           2    0   25     0          0               2     0     0    0    0
    
     pDCs
  0     2
  1     0
  2     1
  3     0
  4     0
  5     2
  6     0
  7     1
  8     0
  9    43
  10    0

根據(jù)上面的結(jié)果對(duì)蛋白質(zhì)水平的細(xì)胞進(jìn)行細(xì)胞類型鑒定

new.cluster.ids <- c("CD4 T", "CD14+ Mono", "NK", "B", "CD8 T", "NK", "CD34+", "T/Mono doublets", 
    "CD16+ Mono", "pDCs", "B")
names(new.cluster.ids) <- levels(cbmc)
cbmc <- RenameIdents(cbmc, new.cluster.ids)

tsne_rnaClusters <- DimPlot(cbmc, reduction = "tsne_adt", group.by = "rnaClusterID") + NoLegend()
tsne_rnaClusters <- tsne_rnaClusters + ggtitle("Clustering based on scRNA-seq") + theme(plot.title = element_text(hjust = 0.5))
tsne_rnaClusters <- LabelClusters(plot = tsne_rnaClusters, id = "rnaClusterID", size = 4)

tsne_adtClusters <- DimPlot(cbmc, reduction = "tsne_adt", pt.size = 0.5) + NoLegend()
tsne_adtClusters <- tsne_adtClusters + ggtitle("Clustering based on ADT signal") + theme(plot.title = element_text(hjust = 0.5))
tsne_adtClusters <- LabelClusters(plot = tsne_adtClusters, id = "ident", size = 4)

注意:為了進(jìn)行比較摔蓝，RNA和蛋白的聚類都是在使用ADT距離矩陣生成的tSNE上顯示的赂苗。

wrap_plots(list(tsne_rnaClusters, tsne_adtClusters), ncol = 2)

image.png

基于adt的聚類產(chǎn)生了類似的結(jié)果，但有一些區(qū)別：

• 基于CD4贮尉、CD8拌滋、CD14和CD45RA的強(qiáng)大ADT數(shù)據(jù)，CD4/CD8 T細(xì)胞群的聚類得到了改善
• 然而猜谚，一些ADT數(shù)據(jù)中不包含有很好的區(qū)別蛋白標(biāo)記的簇(如Mk/Ery/DC)會(huì)失去分離
• 也可以在RNA水平上使用findmarker驗(yàn)證這一點(diǎn)

tcells <- subset(cbmc, idents = c("CD4 T", "CD8 T"))
FeatureScatter(tcells, feature1 = "CD4", feature2 = "CD8")

image.png

RidgePlot(cbmc, features = c("adt_CD11c", "adt_CD8", "adt_CD16", "adt_CD4", "adt_CD19", "adt_CD14"), 
    ncol = 2)

image.png

保存結(jié)果

saveRDS(cbmc, file = "Analysis/cbmc_multimodal.rds")