3.13 与蛋白丰度信息结合

刘小泽写于2020.7.17

1 前言

CITE-seq（Cellular indexing of transcriptomes and epitopes by sequencing）技术是纽约基因组中心和Satija实验室（即Seurat的开发团队）联合在2017年开发的，能够对数千个单细胞的MRNA进行测定，同时还能够对同一个单细胞中的表面蛋白标记物进行测序。

以往的研究方法是在将细胞放在平板上进行scRNA测序之前，通过流式细胞仪捕获单个细胞的蛋白质信息，但这样通量很低，并且受限于相对少量的蛋白质标记物。而CITE-seq可以不受抗体之间信号干扰的影响，大大增加了表面蛋白标记的数量，从而可以一次性获得多种由抗体蛋白指示的免疫表型的信息。

在这个方法中，重点是设计抗体-寡核苷酸结合物，并预估使用量。细胞首先标记上带有RNA标签的抗体，如果一个细胞预估的目标蛋白丰度越高，那么就应该带更多的抗体，也就应该带更多的antibody-derived tag (ADT)。之后利用液滴法将细胞分离，ADTs和内源转录本分别进行反转录、做成cDNA文库、测序。也就是说，蛋白丰度和基因表达量是受ADTs和内源转录本数量影响。

下面来自：https://en.wikipedia.org/wiki/CITE-Seq

• Step1——ADT preparation：It involves labeling an antibody directed against a cell surface protein of interest with oligonucleotides for barcoding the antibody.

• Step2——Bind：ADT labelled cells are encapsulated within a droplet as single cells with DNA-barcoded microbeads

• Step3——Release：Cell lysed to release both bound ADTs as well as mRNA and converted to cDNA. cDNA is prepared from both ADTs and cellular mRNAs.

• Step4——PCR：cDNA is PCR-amplified and ADT cDNA and mRNA cDNA are separated based on size（ADT-derived cDNAs are < 180bp and mRNA-derived cDNAs are > 300bp）

• Step5——Sequence： The independent libraries are pooled together and sequenced to obtain proteomics and transcriptomics data

图片来自：https://www.protocols.io/view/cite-seq-protocols-ngzdbx6

那么问题来了：怎么才能把ADT数据和常规的表达量数据整合到一起呢？

首先需要明确一下：

• 这两种数据的本质是不同的，因此也不太可能直接把ADT行也想表达量数据一样当做是一个feature信息；

• 数量差异巨大：大部分实验只选用了少量感兴趣的蛋白（不超过20种），它们是研究者提前指定好的实验设计，而转录组却是整个转录组，包含几万个基因。

• 测序资源不均等：ADTs的测序深度会比较深，因为它们是从转录本中分离出来单独测序，因此更容易使测序资源产生”二八分布“，即ADTs会占用更多的测序资源

数据准备

将会使用10X PBMC的数据，其中就会包括一些感兴趣的表面蛋白定量结果

library(BiocFileCache)
bfc <- BiocFileCache(ask=FALSE)
stuff <- bfcrpath(bfc, file.path("http://cf.10xgenomics.com",
    "samples/cell-exp/3.0.0/pbmc_10k_protein_v3",
    "pbmc_10k_protein_v3_filtered_feature_bc_matrix.tar.gz"))
untar(stuff, exdir=tempdir())

library(DropletUtils)
sce <- read10xCounts(file.path(tempdir(), "filtered_feature_bc_matrix"))
sce
## class: SingleCellExperiment 
## dim: 33555 7865 
## metadata(1): Samples
## assays(1): counts
## rownames(33555): ENSG00000243485 ENSG00000237613 ... IgG1 IgG2b
## rowData names(3): ID Symbol Type
## colnames: NULL
## colData names(2): Sample Barcode
## reducedDimNames(0):
## altExpNames(0):

> names(rowData(sce))
[1] "ID"     "Symbol" "Type"  

# 看到行中就包含了ADT的信息
> table(rowData(sce)$Type)

Antibody Capture  Gene Expression 
              17            33538

2 预处理

2.1 数据分离

拿到这个sce对象，它里面是整合了转录组和ADT表达信息的，我们首先把ADT分离出来（用splitAltExps()），并且之前介绍sce的时候说到，sce有一个模块是alternative Experiment，使用altExp()获取，可以存储其他类型的实验结果。于是我们可以将ADT的信息放进去

sce <- splitAltExps(sce, rowData(sce)$Type)
altExpNames(sce)
## [1] "Antibody Capture"

# 现在这里就存储着17个ADT的信息
altExp(sce)
## class: SingleCellExperiment 
## dim: 17 7865 
## metadata(1): Samples
## assays(1): counts
## rownames(17): CD3 CD4 ... IgG1 IgG2b
## rowData names(3): ID Symbol Type
## colnames: NULL
## colData names(0):
## reducedDimNames(0):
## altExpNames(0):

如果现在看表达量的话，它是一个稀疏矩阵

> counts(altExp(sce))[1:3,1:3]
3 x 3 sparse Matrix of class "dgCMatrix"

CD3   18  30  18
CD4  138 119 207
CD8a  13  19  10

需要先把它变成常规的矩阵，因为它本身的数据并不稀疏，只是因为使用了sce对象，才默认使用稀疏矩阵节省空间和计算资源，而下游分析与稀疏矩阵的兼容性并不好

counts(altExp(sce)) <- as.matrix(counts(altExp(sce)))
> counts(altExp(sce))[1:3,1:3]
     [,1] [,2] [,3]
CD3    18   30   18
CD4   138  119  207
CD8a   13   19   10

2.2 质控

大多数情况下，质控还是根据转录组数据，去掉空的液滴和低质量细胞。因为细胞质量不好，又怎么能体现真实的转录本和ADT的数量呢？这一部分，低质量细胞的衡量标准还要依靠线粒体水平，而它还是转录组的范畴。因此ADT数据对于质控这一步的贡献不大

我们这里拿到的数据其实是经过CellRanger过滤了的，已经除去了空的液滴，因此只需要关注线粒体水平去掉低质量细胞即可。

library(scater)
mito <- grep("^MT-", rowData(sce)$Symbol)
df <- perCellQCMetrics(sce, subsets=list(Mito=mito))
mito.discard <- isOutlier(df$subsets_Mito_percent, type="higher")
summary(mito.discard)
##    Mode   FALSE    TRUE 
## logical    7569     296

如果只根据转录组过滤感觉不全面的话，也是可以用一下ADT的。基于ADT的特点，它会占用大量的测序资源，它的测序深度一般都很高，因此如果发现某个细胞的ADT表达量很低，那么可以它没有成功加上ADT信息，可以去除

isOutlier使用的是MAD指标（绝对中位差来估计方差,先计算出数据与它们的中位数之间的偏差，然后这些偏差的绝对值的中位数就是MAD，median absolute deviation）。函数认为超过中位数3倍MAD的值就是离群值

# 由于ADT表达量一般很大， 因此需要log2转换
# min_diff含义: minimum difference from the median to consider as an outlier
ab.discard <- isOutlier(df$`altexps_Antibody Capture_detected`,
    log=TRUE, type="lower", min_diff=1)
# 勉强过滤掉一个细胞
summary(ab.discard)
##    Mode   FALSE    TRUE 
## logical    7864       1

最后整合一下过滤信息

discard <- ab.discard | mito.discard
> table(discard)
discard
FALSE  TRUE 
 7568   297 

sce <- sce[,!discard]

2.3 归一化

这一部分，ADT就需要单独进行处理了

由于内源转录本和ADT的生物属性不同，导致它们的细胞捕获效率不同，因此这两种信息的细胞捕获相关的数据偏差也不太可能一致，因此需要各自处理。

另外，ADT的组成差异是很明显的：

• 蛋白的丰度变化更加复杂，蛋白之间是存在互作关系的（试想一下PPI蛋白互作网络），有一点”牵一发动全身“的意思。蛋白丰度的变化，也会直接导致ADT的数量变化

• 靶标蛋白都是提前选好的，因此细胞中也更容易出现丰度的差异

为了减少这些非生物因素差异，有几种方法可以对每个细胞计算size factor：

第一种：基于ADTs的文库大小的计算

首先回顾一下什么是文库大小： We define the library size as the total sum of counts across all genes for each cell, the expected value of which is assumed to scale with any cell-specific biases.

然后看看根据文库大小得到的size factor： The “library size factor” for each cell is then directly proportional to its library size where the proportionality constant is defined such that the mean size factor across all cells is equal to 1

这也是最直接最简单的方法，对分群是够用的，但它做了一个无偏估计，忽略了细胞间原本要上调或者下调的差异，对于分群结果的解释，我们却正需要差异分析的结果去判断每群的marker基因

sf.lib <- librarySizeFactors(altExp(sce))
summary(sf.lib)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.02668  0.52418  0.90485  1.00000  1.26625 22.81740

第二种：基于DESeq的计算

来自DESeq的名词解释： "ratio" uses the standard median ratio method introduced in DESeq. The size factor is the median ratio of the sample over a "pseudosample": for each gene, the geometric mean of all samples.

这里我们也需要指定一个pseudosample，其实最好是来自空液滴的barcode定量结果，但是这里的数据中空液滴被去除，所以只能用其他的均值来代替

# pseudosample
ambient <- rowMeans(counts(altExp(sce)))
# DESeq-like size factors
sf.amb <- medianSizeFactors(altExp(sce), reference=ambient)
summary(sf.amb)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0030  0.5927  0.8282  1.0000  1.1439 41.7972

第三种：基于对照的计算

一些实验包含了对照的抗体（例如IgG），这些对照与真实抗体特性相似，但在细胞中缺乏特异性靶标。我们可以做出假设：对象的ADTs在细胞间的丰度无差异。如果发现其中的任何差异，都可以基于这些差异进行归一化。

controls <- grep("^Ig", rownames(altExp(sce)))
> controls
[1] 15 16 17
sf.control <- librarySizeFactors(altExp(sce), subset_row=controls) 
summary(sf.control)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.6854  0.8757  1.0000  1.1423 44.0155

这里推荐使用第二种基于DESeq的计算方法，最后使用logNormCounts会对转录本和ADTs分别计算size factor、进行归一化和log转换

sizeFactors(altExp(sce)) <- sf.amb
sce <- logNormCounts(sce, use_altexps=TRUE)

assayNames(sce)
## [1] "counts"    "logcounts"
assayNames(altExp(sce))
## [1] "counts"    "logcounts"

3 聚类

与转录本分析不同的是，ADT数据不需要挑选特征信息，因为：

• ADT的特征信息选择早在实验设计阶段就已经完成，我们是先选出感兴趣的蛋白，然后再做的实验得到的数据

• 由于ADTs的数据中行数很少，因此也没有必要去根据这些挑选HVGs或者选PCs

• 每个ADT数据是捕获不同的生物信号，因此这里面的外源噪音并不大，不需要降维也可以轻松去除

综上，对于ADT数据，可以直接对log归一化的结果进行分群和可视化

# 设置NA可以跳过PCA，直接分群
g.adt <- buildSNNGraph(altExp(sce), d=NA) 
clusters.adt <- igraph::cluster_walktrap(g.adt)$membership

# 可视化
set.seed(1010010)
altExp(sce) <- runTSNE(altExp(sce))
colLabels(altExp(sce)) <- factor(clusters.adt)
plotTSNE(altExp(sce), colour_by="label", text_by="label", text_col="red")

因为ADTs数量很少，所以对每个cluster的检查也变得简单，将每个ADT的log表达量均值画成热图即可

se.averaged <- sumCountsAcrossCells(altExp(sce), clusters.adt,
    exprs_values="logcounts", average=TRUE)
# 17个ADT，28个cluster
> dim(se.averaged)
[1] 17 28

> assay(se.averaged)[1:3,1:3]
            1        2        3
CD3  4.621559 9.556522 5.673089
CD4  4.348560 9.795961 9.078986
CD8a 4.797436 5.262174 5.306986

# 热图
library(pheatmap)
averaged <- assay(se.averaged)
pheatmap(averaged - rowMeans(averaged),
    breaks=seq(-3, 3, length.out=101))

4 与基因表达数据的结合

4.1 继续分亚群

上面我们对ADT进行了分群，现在根据这个分群结果+转录组数据再分亚群。

这个思路就很像scRNA数据分析后，再补充一个FACS实验，都是为了更好地分离检查细胞类型。既然得到了ADT的许多”干净“分群结果（因为ADT每群的表达量都很高并且代表的生物信号很强），那就相当于先得到了一些可靠的细胞类型，那么就可以结合转录组表达量去研究其中更细小的变化。

下面就对ADT得到的cluster进行循环分亚群

使用了quickSubCluster函数，它会对现在的sce对象中已有的分组，进一步分群

其中的sce表示对这个对象进行操作，clusters.adt指定分组，prepFUN指的是在分群之前的准备工作，clusterFUN指的是分群的操作

set.seed(101010)
all.sce <- quickSubCluster(sce, clusters.adt,
    prepFUN=function(x) {
        dec <- modelGeneVar(x)
        top <- getTopHVGs(dec, prop=0.1)
        x <- runPCA(x, subset_row=top, ncomponents=25)
    },
    clusterFUN=function(x) {
        g.trans <- buildSNNGraph(x, use.dimred="PCA")
        igraph::cluster_walktrap(g.trans)$membership
    }
)
# 结果是一个列表
> all.sce
List of length 28
names(28): 1 2 3 4 5 6 7 8 ... 21 22 23 24 25 26 27 28

如果对应上面的t-SNE图，看到第3群细胞数量还蛮多的，可以看一下

> all.sce[[3]]
class: SingleCellExperiment 
dim: 33538 1828 
metadata(1): Samples
assays(2): counts logcounts
rownames(33538): ENSG00000243485 ENSG00000237613
  ... ENSG00000277475 ENSG00000268674
rowData names(3): ID Symbol Type
colnames: NULL
colData names(4): Sample Barcode sizeFactor
  subcluster
reducedDimNames(1): PCA
altExpNames(1): Antibody Capture

# 它又可以分8个亚群
> table(all.sce[[3]]$subcluster)

3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 
 90 772  25 178 695  44  19   5

看一下每个ADT分群结果各自对应的亚群

# 每个ADT群下面有多少细胞
ncells <- sapply(all.sce, ncol,simplify = TRUE)
> ncells
   1    2    3    4    5    6    7    8    9   10   11 
 888   92 1828 1186  512   55   71   85  999  571   70 
  12   13   14   15   16   17   18   19   20   21   22 
 148   73   39  409   21   39   69   74  126   32   50 
  23   24   25   26   27   28 
  13   17   29   38   22   12 

# 每个ADT群下面又分几个亚群
nsubclusters <- sapply(all.sce, function(x) length(unique(x$subcluster)),
                       simplify = TRUE)
> nsubclusters
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 
21  3  8  7  3  3  2  3  9  8  3  2  3  1 12  1  1  2 
19 20 21 22 23 24 25 26 27 28 
 3  3  1  2  1  1  1  1  1  1

画个图看一下

plot(ncells, nsubclusters, xlab="Number of cells", type="n",
    ylab="Number of subclusters", log="xy")
text(ncells, nsubclusters, names(all.sce))

继续分亚群的一个好处是：我们可以基于之前ADT的分群结果有一个大体的预估，每个cluster大体是什么性质的。假如我们知道了ADT的一个clusterX中包含T细胞，那么就没有必须再对X.1、X.2之类的亚群再进行判断是不是T细胞了。而是可以把更多的精力放在findMarkers()后续分析上。比如已经发现cluster12中包含CD8+ T细胞，而且它又分成了两个亚群，就可以继续根据颗粒酶表达量判断亚群

of.interest <- "12"
# 左边是GZMH基因，右边是GZHK
plotExpression(all.sce[[of.interest]], x="subcluster",
    features=c("ENSG00000100450", "ENSG00000113088"))

注意

我们这里的亚群都是根据之前的分群结果继续向下走的，如果说之前的某个分群不对怎么办？那不就误导了后续的分析？

因此可以补充一些检查，看看每个亚群是不是和它们的上一级分群结果有相似的蛋白丰度

of.interest <- "12"
sce.cd8 <- all.sce[[of.interest]]
plotExpression(altExp(sce.cd8), x=I(sce.cd8$subcluster),
    features=c("CD3", "CD8a"))

4.2 利用转录组的分群找蛋白相关的marker基因

我们选用ADT看蛋白丰度的目的应该不是为了看细胞类型，因为这个事情使用转录组也能完成，并且完成的效果还不错。那么为什么使用蛋白数据呢？因为蛋白可以更直接地反映生物过程活性，可以减少后续手动注释细胞类型的麻烦和不准确性。因此，结合转录组的分群结果，这样会不会更快地帮助判断某群细胞的生物功能呢？

对转录组数据进行快速分群

sce <- logNormCounts(sce)
dec <- modelGeneVar(sce)
top <- getTopHVGs(dec, prop=0.1)

set.seed(1001010)
sce <- runPCA(sce, subset_row=top, ncomponents=25)

g <- buildSNNGraph(sce, use.dimred="PCA")
clusters <- igraph::cluster_walktrap(g)$membership
colLabels(sce) <- factor(clusters)

set.seed(1000010)
sce <- runTSNE(sce, dimred="PCA")
plotTSNE(sce, colour_by="label", text_by="label")

ADT数据+ 分群结果 =》找marker基因

markers <- findMarkers(altExp(sce), colLabels(sce))
of.interest <- markers[[16]]
pheatmap(getMarkerEffects(of.interest), breaks=seq(-3, 3, length.out=101))

图中可以看到，cluster16中的PD-1表达量升高，如果再把它和某个感兴趣的表型结合起来（例如T细胞衰竭），就可以对数据多一些认识

用于标记衰竭程度的一个分子标记是PD-1，以“阻断T细胞杀伤能力”而闻名。PD-1的沉默会通过抑制T细胞的增殖而损害其抗肿瘤功能。

来自：PD-1免疫疗法的那些事儿 Naive T细胞在TCR信号通路激活的时候，会飞快地表达PD-1。而当T细胞被连续不断地刺激时，T细胞会变成exhausted T，会持续性地表达很高的PD-1

整个过程不需要任何的预判，也不需要对ADT有任何的了解，所有的结果都是基于分群+表达量差异。既然后面是按照ADT的数据推断的marker基因，那么如果有直接针对ADT数据的聚类分群，结果就会更准确，只不过目前数据量还达不到像转录组一样的水平，做降维聚类效果不太好。