Seurat4.0系列教程21：结合Cell Hashing分析双细胞 by 单细胞天地

分享是一种态度

Cell Hashing 由NYGC 技术创新小组与Satija实验室合作开发，使用寡核苷酸标记的抗体标记细胞表面表达的蛋白质，在每个单细胞上放置一个"样本条形码"，使不同的样品能够一起多路复用，并在单次实验中运行。欲了解更多信息，请参阅此文^[1]

此教程简要演示如何处理 Seurat 中与Cell Hashing一起生成的数据。应用于两个数据集，我们可以成功地将细胞分离到它们原始的来源，并识别跨样本的双细胞。

The demultiplexing 函数 HTODemux()执行了以下程序：

在标准化的 HTO 值上执行 k-medoid 聚类，该值最初将细胞分离为 K（样本的# ）+1 群。
计算 HTO 的"negative"分布。对于每个 HTO，我们使用平均值最低的群作为negative组。
对于每个 HTO，我们选合适的负二元分布到negative组。我们使用此分布的 0.99 分位作为阈值。
根据这些阈值，每个细胞根据 HTO被归类为 positive 或negative。
多于1个 HTO positive 的细胞被注释为双细胞。

来自人类 PBMC 的 8-HTO 数据集

数据集描述：

数据来自八个不同捐赠者的外周血单核细胞（PBMC）。
每个捐赠者的细胞都具有独特的标签，使用CD45作为hashing抗体。
样品随后被混样，并在 10X v2 系统的单个lane上运行。
您可以在此处^[2]下载RNA和HTO的计数矩阵，或来自GEO^[3]的 FASTQ 文件

基本设置

加载包

library(Seurat)

读取数据

# Load in the UMI matrix
pbmc.umis <- readRDS("../data/pbmc_umi_mtx.rds")

# For generating a hashtag count matrix from FASTQ files, please refer to
# https://github.com/Hoohm/CITE-seq-Count.  Load in the HTO count matrix
pbmc.htos <- readRDS("../data/pbmc_hto_mtx.rds")

# Select cell barcodes detected by both RNA and HTO In the example datasets we have already
# filtered the cells for you, but perform this step for clarity.
joint.bcs <- intersect(colnames(pbmc.umis), colnames(pbmc.htos))

# Subset RNA and HTO counts by joint cell barcodes
pbmc.umis <- pbmc.umis[, joint.bcs]
pbmc.htos <- as.matrix(pbmc.htos[, joint.bcs])

# Confirm that the HTO have the correct names
rownames(pbmc.htos)

## [1] "HTO_A" "HTO_B" "HTO_C" "HTO_D" "HTO_E" "HTO_F" "HTO_G" "HTO_H"

设置Seurat对象并添加 HTO 数据

# Setup Seurat object
pbmc.hashtag <- CreateSeuratObject(counts = pbmc.umis)

# Normalize RNA data with log normalization
pbmc.hashtag <- NormalizeData(pbmc.hashtag)
# Find and scale variable features
pbmc.hashtag <- FindVariableFeatures(pbmc.hashtag, selection.method = "mean.var.plot")
pbmc.hashtag <- ScaleData(pbmc.hashtag, features = VariableFeatures(pbmc.hashtag))

添加 HTO 数据作为独立assay

可以在此处阅读更多有关使用多模式数据的信息^[4]

# Add HTO data as a new assay independent from RNA
pbmc.hashtag[["HTO"]] <- CreateAssayObject(counts = pbmc.htos)
# Normalize HTO data, here we use centered log-ratio (CLR) transformation
pbmc.hashtag <- NormalizeData(pbmc.hashtag, assay = "HTO", normalization.method = "CLR")

基于 HTO 富集Demultiplex细胞

在这里，我们使用 Seurat 函数HTODemux()将单个细胞分配回其样本来源。

# If you have a very large dataset we suggest using k_function = 'clara'. This is a k-medoid
# clustering function for large applications You can also play with additional parameters (see
# documentation for HTODemux()) to adjust the threshold for classification Here we are using the
# default settings
pbmc.hashtag <- HTODemux(pbmc.hashtag, assay = "HTO", positive.quantile = 0.99)

可视化demultiplexing 的结果

运行HTODemux()的输出保存在metadata中。我们可以可视化多少细胞被归类为单细胞、双细胞和负细胞/模糊细胞。

# Global classification results
table(pbmc.hashtag$HTO_classification.global)

## 
##  Doublet Negative  Singlet 
##     2598      346    13972

山脊图可视化选定 HTO的富集

# Group cells based on the max HTO signal
Idents(pbmc.hashtag) <- "HTO_maxID"
RidgePlot(pbmc.hashtag, assay = "HTO", features = rownames(pbmc.hashtag[["HTO"]])[1:2], ncol = 2)

可视化HTO 信号对，确认单个细胞中的相互排他性

FeatureScatter(pbmc.hashtag, feature1 = "hto_HTO-A", feature2 = "hto_HTO-B")

比较单细胞、双细胞和negative 细胞的 UMI 数量

Idents(pbmc.hashtag) <- "HTO_classification.global"
VlnPlot(pbmc.hashtag, features = "nCount_RNA", pt.size = 0.1, log = TRUE)

为HTOs生成一个二维tsNE图。在这里，我们将细胞按单细胞和双细胞分组，以实现简化。

# First, we will remove negative cells from the object
pbmc.hashtag.subset <- subset(pbmc.hashtag, idents = "Negative", invert = TRUE)

# Calculate a tSNE embedding of the HTO data
DefaultAssay(pbmc.hashtag.subset) <- "HTO"
pbmc.hashtag.subset <- ScaleData(pbmc.hashtag.subset, features = rownames(pbmc.hashtag.subset), 
    verbose = FALSE)
pbmc.hashtag.subset <- RunPCA(pbmc.hashtag.subset, features = rownames(pbmc.hashtag.subset), approx = FALSE)
pbmc.hashtag.subset <- RunTSNE(pbmc.hashtag.subset, dims = 1:8, perplexity = 100)
DimPlot(pbmc.hashtag.subset)

# You can also visualize the more detailed classification result by running Idents(object) <-
# 'HTO_classification' before plotting. Here, you can see that each of the small clouds on the
# tSNE plot corresponds to one of the 28 possible doublet combinations.

根据Cell Hashing文章中图 1C 创建 HTO 热图。

# To increase the efficiency of plotting, you can subsample cells using the num.cells argument
HTOHeatmap(pbmc.hashtag, assay = "HTO", ncells = 5000)

使用常规的 scRNA-seq 工作流对细胞进行聚类和可视化，并检查潜在的批次效应。

# Extract the singlets
pbmc.singlet <- subset(pbmc.hashtag, idents = "Singlet")

# Select the top 1000 most variable features
pbmc.singlet <- FindVariableFeatures(pbmc.singlet, selection.method = "mean.var.plot")

# Scaling RNA data, we only scale the variable features here for efficiency
pbmc.singlet <- ScaleData(pbmc.singlet, features = VariableFeatures(pbmc.singlet))

# Run PCA
pbmc.singlet <- RunPCA(pbmc.singlet, features = VariableFeatures(pbmc.singlet))

# We select the top 10 PCs for clustering and tSNE based on PCElbowPlot
pbmc.singlet <- FindNeighbors(pbmc.singlet, reduction = "pca", dims = 1:10)
pbmc.singlet <- FindClusters(pbmc.singlet, resolution = 0.6, verbose = FALSE)
pbmc.singlet <- RunTSNE(pbmc.singlet, reduction = "pca", dims = 1:10)

# Projecting singlet identities on TSNE visualization
DimPlot(pbmc.singlet, group.by = "HTO_classification")

来自四个人类细胞系的 12个HTO 数据集

数据集描述：

数据来源于从四个细胞系HEK、K562、KG1 和 THP1收集的单细胞：
每个细胞系被进一步分成三个样本（总共12个样本）。
每个样品都标有 hashing抗体混合物（CD29和CD45），汇集在一起，在10X的单lane上运行。
基于此设计，我们应该能够检测跨细胞类型和细胞类型内的双细胞
您可以在此处^[5]下载RNA和HTO的计数矩阵，并可在GEO上找到^[6]

创建Seurat对象，添加 HTO 数据并执行标准化

# Read in UMI count matrix for RNA
hto12.umis <- readRDS("../data/hto12_umi_mtx.rds")

# Read in HTO count matrix
hto12.htos <- readRDS("../data/hto12_hto_mtx.rds")

# Select cell barcodes detected in both RNA and HTO
cells.use <- intersect(rownames(hto12.htos), colnames(hto12.umis))

# Create Seurat object and add HTO data
hto12 <- CreateSeuratObject(counts = hto12.umis[, cells.use], min.features = 300)
hto12[["HTO"]] <- CreateAssayObject(counts = t(x = hto12.htos[colnames(hto12), 1:12]))

# Normalize data
hto12 <- NormalizeData(hto12)
hto12 <- NormalizeData(hto12, assay = "HTO", normalization.method = "CLR")

Demultiplex data

hto12 <- HTODemux(hto12, assay = "HTO", positive.quantile = 0.99)

可视化demultiplexing的结果

选定 HTO 的分布按照分类分组，用山脊图展示

RidgePlot(hto12, assay = "HTO", features = c("HEK-A", "K562-B", "KG1-A", "THP1-C"), ncol = 2)

在热图中可视化 HTO 信号

HTOHeatmap(hto12, assay = "HTO")

可视化RNA聚类

下面，我们使用标准的 scRNA-seq 工作流对细胞进行聚类。正如预期的那样，我们看到四个主要群，对应于细胞系
此外，我们看到中间有小群，表示正确注释为双细胞的混合转录组。
我们还看到细胞亚型内的双细胞，它们与同一细胞类型的单细胞混合在一起

# Remove the negative cells
hto12 <- subset(hto12, idents = "Negative", invert = TRUE)

# Run PCA on most variable features
hto12 <- FindVariableFeatures(hto12, selection.method = "mean.var.plot")
hto12 <- ScaleData(hto12, features = VariableFeatures(hto12))
hto12 <- RunPCA(hto12)
hto12 <- RunTSNE(hto12, dims = 1:5, perplexity = 100)
DimPlot(hto12) + NoLegend()