单细胞混样技术（Cell hashing）_一种新奇的单细胞测序技术

Hi，大家好，我是晨曦
今天这期推文来自学习群里小伙伴的提问，他向晨曦提问有没有Cell hashing技术的分析教程，其实当我们产生了一种诉求以后，只要这个诉求是普遍的，那么就一定会有教程或者是guidance在网络上，而晨曦所做的只是去学习然后转化成自己可以理解的语言
所以当时也是跟这位小伙伴打赌，说要是坚持打卡学习一周晨曦就去学一下然后写成教程，结果这位小伙伴真的坚持了一周的学习，所以这期推文就是晨曦愿赌服输的产物
也是正好有机会和大家一起学习一下单细胞混样技术（Cell hashing）的分析流程
从同一个细胞同时测量多种数据类型的能力，被称为多模态分析，代表了单细胞基因组学的一个新的和令人兴奋的前

那么，我们开始吧

引言

这一部分我们主要讲解两部分内容：一部分是技术本身；一部分则是我们后面演示数据集来自的参考文献

那么，我们首先讲解第一部分：技术本身（单细胞混样技术（Cell hashing）），下面展示一下要点总结：

Cell Hashing是在CITE-seq的基础上改进的，CITE-seq全称cellular indexing of transcriptomes and epitopes by sequencing，是一种同时对细胞内RNA和细胞表面蛋白进行测序的技术，而Cell Hashing技术不同于CITE-seq技术的点在于，其是为了解决如何将不同样本的细胞混起来测序（便宜），测完了还能区分哪个细胞来源于哪个样本，这样做也减少了批次效应

Cell Hashing是在CITE-seq的基础上改进，是给需要混合的样品提前加上HTO (A distinct Hashtag oligonucleotide) 标签，这样即使混合后也可以提供不同的HTO标签进行区分

了解了技术本身，我们再来看一下我们后面代码演示数据集来自的参考文献：

中文题目：人类骨髓（BM）和脐血（UCB）细胞的单细胞转录组揭示了红细胞分化连续体的假定调节因子

摘要：

了解成年人类造血干细胞（HSC）参与红细胞生成的的精细分化轨迹对于理解人类红细胞生成的动态发育变化至关重要。利用直接从成人

骨髓(BM)和脐带血(UCB)中分离的原代人类末端红细胞(CD34−CD235a+)的单细胞RNA测序(scRNA-seq)，本研究以前所未有的分辨率记录

了最终分化的人类成红细胞的转录组。这些见解使我们能够区分处于发育早期和晚期的多色成红细胞（PolyEs），以及正色成色细胞

（OrthoEs）的不同发育阶段。本研究通过监测细胞分化和凋亡，进一步鉴定出一套终端红系分化调节剂，并在功能上验证了已鉴定的三

个基因AKAP8L，TERF2IP和RNF10。本研究记录了AKAP8L的敲低抑制了造血干细胞对红系谱系和细胞增殖的命运，并抑制了红系集落

（CFU-E）到成红细胞期（ProE）的分化。相比之下，敲低TERF2IP和RNF10延迟PolyE分化为OrthoE阶段

数据集：GSE150774

样本信息：从六个年轻的健康男性供体中收集了10毫升的BM样品，从六个健康的新生儿脐带中收集了20毫升的UCB样品。基于10x 平台进行测序后，最终获得了来自三个BM样品的8,668个细胞和来自三个UCB样品的17,692个细胞的数据。

代码实操（思路）

下面尝试进行分析思路的复现：

我们首先面对的第一个问题其实就是，那个95.9Mb的文件到底是什么？

因为我们理解了单细胞表达矩阵的本质，其实就是观测是Gene symbol，变量是Barcode ID，然后表达信息是稀疏矩阵，所以说这些一个文件的我们都可以把其整理成Seurat对象，但是第一个文件GSE4558614_UCB2UCB3_.....就很让人费解，因为为什么会把脐带血两个样本混在一起？

其实这个样本就是所谓的单细胞混样技术（Cell hashing）,那么下面我们通过代码来逐步的进行操作

#准备工作library(Seurat)library(ggplot2)library(tidyverse)library(stringr)library(harmony)

#读取输入数据（这里我们主要基于GSE4558614_UCB2UCB3这个样本，因为只有这个样本是基于Cell hashing）

HP <- Read10X("HP0_filtered_gene_bc_matrices/")

我们可以很清楚的看到，我们的环境变量中有一个名称为HP的列表，列表内则包含了两类信息，分别是表达矩阵以及

[

](https://imgtu.com/i/X7occD

#先查看一下HTO标志的数据yhto <- as.data.frame(HP$`Antibody Capture`)rownames(yhto)#[1] "B0251 anti-human Hashtag1" "B0252 anti-human Hashtag2"#这个数据存在的目的其实就是帮助我们区分不同样本的细胞

可以很清楚的看到一共是两个样本的数据

那么，到这里我们需要再次明确我们分析的目的

我们分析的目的其实就是因为我们的数据基于Cell hashing技术，这就导致我们读取的这个数据其实是包含两个样本的数据，而我们想做的就是把两个样本分开，因为都混在一起了，我们也不太好进行后续的操作

#我们来看一下表达矩阵数据

表达矩阵数据显然只有一个表达矩阵，所以我们可以有理由推断，Cell hashing技术是把多样本的表达信息融合在一个表达矩阵之中，然后通过标志矩阵（HTO tag矩阵）识别不同样本来源的细胞ID，进而允许我们将其分离

明确了我们的目的后，我们继续往下探索

#简单进行质控，减少下游分析的计算机负荷pbmc.htos <- yhto[colSums(yhto)>0]#在两个样本中都是0的barcode ID显然就不是我们研究的重点#然后再查看表达量矩阵pbmc.umis <- HP$`Gene Expression`

joint.bcs <- intersect(colnames(pbmc.umis), colnames(pbmc.htos))#因为我们上一步对表面蛋白进行了质控，所以这里我们需要取交集，获得交集信息

#提取最终表达矩阵pbmc.umis <- pbmc.umis[, joint.bcs]pbmc.htos <- as.matrix(pbmc.htos[, joint.bcs])

然后我们到这里就简单完成了数据的EDA以及简单的过滤筛选，下面我们就可以导入Seurat中进行标准的下游分析流程

#创建Seurat对象pbmc.hashtag <- CreateSeuratObject(counts = pbmc.umis)#标准化pbmc.hashtag <- NormalizeData(pbmc.hashtag)#寻找高可变pbmc.hashtag <- FindVariableFeatures(pbmc.hashtag, selection.method = "mean.var.plot")#标准化pbmc.hashtag <- ScaleData(pbmc.hashtag, features = VariableFeatures(pbmc.hashtag))pbmc.hashtag#An object of class Seurat #32738 features across 18777 samples within 1 assay #Active assay: RNA (32738 features, 1749 variable features)

我们可以清楚的看到，我们这个Seurat对象中包含一个数据矩阵信息，也就是说Cell hashing是把样本混合起来成一个大矩阵然后展现出来的

#添加HTO标签矩阵（用来识别出不同样本来源）pbmc.hashtag[["HTO"]] <- CreateAssayObject(counts = pbmc.htos)pbmc.hashtag#An object of class Seurat #32740 features across 18777 samples within 2 assays #Active assay: RNA (32738 features, 1749 variable features)# 1 other assay present: HTO

这样我们就可以看到，这个时候我们的Seurat对象中就包含了两个方面的信息

#针对HTO矩阵需要进行CLR转换pbmc.hashtag <- NormalizeData(pbmc.hashtag, assay = "HTO", normalization.method = "CLR")

上面的步骤我们肯定会存在疑问，为什么需要进行CLR转换：

其实到这里我们大概可以理解到，我们之所以可以把不同所属的细胞找到其实并不是我们一开始所想的，加上一个标签，然后测序后跟着标签把细胞一个个归类到样本中（Ps:一开始我真的是这么想的！！），但是通过对于后续代码的理解，我们是通过HTO标签的然后进行聚类，通过聚类的方法把样本划归到所属的样本上，所以我们需要进行标准化数据，之所以用这种方法，则是方法学文献中要求的

head([email protected])# orig.ident nCount_RNA nFeature_RNA nCount_HTO nFeature_HTO#AAACCCAAGAAGGTAG-1 SeuratProject 933 196 5 2#AAACCCAAGCTGTGCC-1 SeuratProject 4791 1636 11 1#AAACCCAAGGAAACGA-1 SeuratProject 7220 439 1 1#AAACCCAAGGGAGGAC-1 SeuratProject 8998 1424 3 2#AAACCCAAGGGTCAAC-1 SeuratProject 8745 667 2 1#AAACCCAAGTGCAGCA-1 SeuratProject 564 90 1 1

根据文献中提示，我们到最后是根据HTO表达量的情况来把我们的样本通过聚类的方式划归到每一个所属的原始样本中去的

#核心步骤pbmc.hashtag <- HTODemux(pbmc.hashtag, assay = "HTO", positive.quantile = 0.99)

这一步就是我们Cell hashing的核心，也就是HTODemux函数的使用

1. 我们对标准化的HTO值执行K -medoid聚类，它最初将细胞分离为K(# of samples )+1聚类

假设这里我们有两个样本，那我们自然就是设置三个聚类，其中两个聚类对特定的HTO表达呈现高度阳性，而第三个聚类则是对所有HTO表达呈现高度阴性的一个集合

2. 我们计算了HTO的一个“负”分布。对于每个HTO，我们使用平均值最低的群作为背景组

首先，我们鉴定了HTO平均表达量最高的k-HTO簇，并排除了这些细胞。

在进一步排除最高的0.5%的值作为潜在的异常值后

我们接下来将剩余的HTO值拟合为一个负的二项分布，我们计算了拟合分布的q=0.99分位数，并基于这个HTO特定值对数据集中的每个单元进行阈值筛选

3. 对于每个HTO，我们对负的聚类拟合一个负的二项分布。我们使用这个分布的0.99分位数作为阈值。

4. 根据这些阈值，每个细胞被划分为阳性或阴性的HTO

我们使用上面的程序来确定每个条形码的“HTO分类”。只有一个HTO呈阳性的条形码被归类为单类体。两个或两个以上hto呈阳性的条形码被归类为多重体，并根据其两个最高表达的HTO分配样本id。所有8个hto均为阴性的条形码被归类为“阴性”

所以，简单来说就是HTODemux函数是为了帮助我们完成上述的流程而诞生的

table(pbmc.hashtag$HTO_classification.global)# Doublet Negative Singlet # 152 14650 3975

但是很奇怪的是阴性细胞居然会这么多，按照常理来说应该是Singlet细胞类型最多，可能这就是样本的问题

pbmc.hashtag@assays#$RNA#Assay data with 32738 features for 18777 cells#Top 10 variable features:# HBE1, S100A9, CLC, HBD, PF4, PPBP, S100A8, HSPA1A, LYZ, NKG7 #$HTO#Assay data with 2 features for 18777 cells#First 2 features:# B0251 anti-human Hashtag1, B0252 anti-human Hashtag2

下面的代码运行后，可以很清楚的看到我们的样本被成功的划分了出来

table(pbmc.hashtag$hash.ID)Doublet Negative B0251 anti-human Hashtag1 B0252 anti-human Hashtag2 152 14650 2144 1831

既然细胞已经被我们成功的划分了出来，那么自然我们就可以开始进行后续的一系列标准化的操作，但是在这之前，我们需要剔除掉数据中我们不需要的部分

#剔除阴性以及双重pbmc.hashtag.subset <- subset(pbmc.hashtag, idents = c("Negative","Doublet"), invert = TRUE)table(Idents(pbmc.hashtag.subset)) #B0251 anti-human Hashtag1 B0252 anti-human Hashtag2 # 2144 1831 #标准流程pbmc.hashtag.subset <- RunPCA(pbmc.hashtag.subset)pbmc.hashtag.subset <- RunUMAP(pbmc.hashtag.subset,dims=1:10)DimPlot(pbmc.hashtag.subset) #提取特定样本的细胞B0251 <- subset(pbmc.hashtag, idents = c("B0251 anti-human Hashtag1"))B0252 <- subset(pbmc.hashtag, idents = c("B0252 anti-human Hashtag2"))#合并sce.all <- merge( B0251, y= B0252 ,add.cell.ids = c('B0251','B0252')) #获得样本与细胞ID之间的对应关系并赋值到meta信息中phe=as.data.frame(str_split(colnames(sce.all),'_',simplify = T))[email protected]$orig.ident = phe$V1table([email protected]$orig.ident )head([email protected])

经过上面的处理，我们的数据处理就已经结束了，我们下面就可以进行标准流程

sce=sce.all#赋值一下，方便后续检查对象#进行harmony的数据要求是需要标准流程走到PCA为止sce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 1e4) #标准化数据sce <- FindVariableFeatures(sce)#寻找高可变sce <- ScaleData(sce)#标准化sce <- RunPCA(sce, features = VariableFeatures(object = sce))#PCAseuratObj <- RunHarmony(sce, "orig.ident")#合并数据集names(seuratObj@reductions)seuratObj <- RunUMAP(seuratObj, dims = 1:15,

                    reduction = "harmony")#合并完数据集以后就可以开始走标准流程的后半部分了

DimPlot(seuratObj,reduction = "umap",label=F ) sce=seuratObjsce <- FindNeighbors(sce, reduction = "harmony", dims = 1:15) sce <- FindClusters(sce, resolution = 0.8)sce.all=sce DimPlot(sce.all)

然后到这里我们就完成了Cell hashing的分析，那么我们来看看我们最后的图表展示

#图表1DimPlot(sce.all,label = T)

#图表2pbmc.markers <- FindAllMarkers(sce.all, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)marker <- pbmc.markers %>% group_by(cluster) %>% slice_max(n = 1, order_by = avg_log2FC)DotPlot(sce.all, features = marker$gene)+RotatedAxis()+ scale_x_discrete("")+scale_y_discrete("")