Introduction

肿瘤内异质性是决定肿瘤生物学，治疗反应以及病人存活情况的一个关键性因素。因此，对肿瘤中出现的各种类型细胞的表型以及相互之间的联系进行一个全面的描述是很有必要的。传统的针对 bulk-tissue 的分析会掩盖肿瘤内不同细胞群体的特征，虽然有研究通过一些 marker 来对细胞进行分类，但这种方法需要我们有足够的先验知识，并且依赖特异性、高亲和力的抗体。随着单细胞测序技术的出现以及流行，我们能够对肿瘤中细胞多样性进行全面的，无偏好性的分析。不过，在单细胞转录组领域，还是存在诸多因素导致难以准确区分不同类型的细胞。本文献中，作者提出了一种被称为 reference component analysis (RCA) 的新算法用于聚类分析。
结直肠癌 (colorectal cancer, CRC) 是全世界第三常见的肿瘤。本文献中，作者针对从11个原发结直肠癌及对应正常样本中分理出590个（经过质控）细胞进行 scRNA-seq 分析，并利用 RCA 进行聚类。此外，作者对不同类型的细胞进行了表达特征的描述，利用他们对病人存活情况进行分层。

RCA与现有聚类方法的比较

现有针对 scRNA-seq 数据进行聚类的方法的准确性没办法进行全面的量化，因为我们并不知道真实的细胞类型，一般只能通过某一类中的细胞表达了预期的细胞类型标记做一个定性的验证。但是，这种基于 marker 的方法有时候可能会误导我们，因为 single-cell 测序往往难以发现那些表达量较低的基因，以及 marker 基因的异位表达等情境也可能干扰我们的判定。
为了对不同的聚类方法做一个比较，作者从7个细胞系中经质控挑出561个细胞的单细胞转录组数据。这些细胞有两个不同来源：GM12787(类淋巴母细胞)和H1胚胎干细胞，可以评估批次效应带来的影响。利用现有的8大 scRNA-seq 聚类方法对这些数据进行了分析。这些方法包括：

hierarchical clustering using all expressed genes (All-HC)
hierarchical clustering using principal-component analysis (PCA)-based feature selection (HiLoadG-HC)
BackSPIN10
RaceID2
Seurat1
VarG-HC
VarG-PCAproj-HC
VarG-tSNEproj-HC

其中后三种方法是基于那些表达差异特别大的基因进行聚类的。针对聚类的结果，利用 adjusted Rand Index (ARI) 进行量化：ARI 为0时，表示聚类是随机的；为1时，表示结果完美。各方法的 ARI 结果如下图所示:

从图中作者得出以下结论：

考虑到现有聚类方法的准确性有限，作者认为 scRNA-seq 数据的系统方差主要是由于技术误差而非生物学差异组成。作者也确实发现有一种已知的技术性协变量，即每个细胞中发现的基因数量（number of detected gene, NODG），能够解释第一主成分中几乎所有的变异
现有的聚类算法效果都不理想

于是，作者效仿最近针对低覆盖度全基因组测序数据进行祖源分析所发展的 reference projection strategy，希望能将嘈杂的，不完整的数据集映射到由一个独立高质量地参考 panel 中的变量所定义的低维空间上。为了构建这样一个参考性 panel，作者整合了来自不同组织，不同类型细胞混池样本的转录组数据。将 scRNA-seq 数据映射到该 global reference panel 上，然后将标准化后的空间坐标提供给传统的层次聚类程序进行聚类，这种方法被称为 RCA。从上图我们可以看出，RCA 的表现明显好于其他的聚类方法

利用 RCA 确定 CRC 肿瘤和正常粘膜中的细胞类型

作者从11个 CRC 病人切除的原发瘤中取了969个细胞，并从这些肿瘤附近的正常粘膜出取了622个细胞。经过严格的质控，剩下375个肿瘤细胞和215个正常粘膜细胞。利用 RCA 对这些细胞进行聚类后，得到的结果如下图所示。
图中，最上面的系统进化树表示不同细胞类型簇。上部热图中，每一列代表一个细胞，每一行对应 reference panel 中的成分，颜色表示映射分值；中部热图中，每一行代表某种类型细胞的一个 marker gene，颜色深浅对应 log10(FPKM)；下部热图中，每一行代表一个病人，每个 block 代表细胞取自哪位病人。
从左图a中，我们可以看到结直肠整肠粘膜中包含7种不同类型的细胞，且每个细胞簇中都包含来自不同病人的细胞，表明 RCA 并不受细胞来源带来的批次效应的影响。这些类型的细胞在右图b的 CRC 肿瘤细胞中同样存在。

在上述分析中，我们可以看到上皮细胞并未分出不同的亚簇，例如肠细胞，杯状细胞等。因此，作者对 RCA 聚类做了一个变形，将上皮细胞的单细胞转录组用做 reference panel，称为 self projection (自映射)，从而得到了9个亚簇，如下图所示。根据上皮细胞标记的表达模式，可以确定这些确实是不同的细胞类型及（或）状态。

（图中中部热图每行代表各细胞类型的标记）

对于 CRC 肿瘤细胞，作者用9个正常粘膜的上皮细胞转录组构建了一个 reference panel，称为 colon epithelium panel。利用该 panel，RCA 检测出三种类型的肿瘤上皮细胞。
综上，RCA 能够有效地鉴定正常粘膜和 CRC 肿瘤中的不同细胞类型，而不受临床样本批次效应带来的影响。

scRNA-seq 鉴定肿瘤中差异表达的基因

作者假设从单细胞测序得到的肿瘤-正常样本差异表达的基因与传统混池样本得到的结果会有很大不同。因此，他们优先选择了5个最大的细胞簇：stem/TA-like epithelial cells, fibroblasts, B cells, T cells and myeloid cells，鉴定出129个 FDR≤ 0.05 的差异表达基因。这些基因主要由 stem/TA-like epithelial cells 和 fibroblasts 贡献。将这些基因同一个大型 CRC 队列的混池样本得到的数据进行比较，发现两者之间其实有较明显的一致性，不少利用单细胞测序发现的差异表达基因之前也在混池测序中发现过（如下图所示）。不过，大部分差异表达的基因还是第一次发现。这说明两者之间有一定的一致性，但有限。

利用混池样本进行差异表达分析的可能会导致把某种类型细胞特异表达的基因误认为在肿瘤中调节异常的基因，例如肿瘤中可能包含大量 stem/TA-like cells，因此可能导致该类型细胞的标志基因相对正常样本呈现出差异表达的假象。为了研究这样的偏差性，作者从正常粘膜单细胞数据中鉴定出292个不同细胞类型或者上皮细胞亚型之间差异表达的基因，在剔除48个肿瘤-正常单细胞数据中也存在的差异表达基因后，对比混池样本得到的结果，发现这些基因确实在混池样本中富集。
综上，单细胞测序数据能更好地帮助我们找出肿瘤中表达异常的基因。

scRNA-seq 揭示 CRC 中 CAFs 的多样性和通路的改变

为了鉴定差异表达基因的上游调控者，作者进行了 Ingenuity Pathway Analysis (IPA)，将4种细胞 (stem/TA-like cells, fibroblasts, myeloid cells and B cells) 的结果混合并按照 P 值排序后，发现排名最靠前的三个调控者中包含两个小分子：actinonin 和 motexafin gadolinium。actinonin 能抑制多种肿瘤细胞系的生长，包括结直肠癌；而motexafin gadolinium，在 stem/TA-like epithelial cells 中高表达，能够促使癌症细胞凋亡。所以我们可以从下图中看到这两个分子的表达都是显著下降的。

此外，从图中我们可以看到 TGFB1 也是一个高排名的内源性上游调控者，z 值为正代表其在 cancer-associated fibraoblasts (CAFs) 通路活性的增加。与之相关的 TGF-β 效应因子 SMAD3 也被预测为活性增加。在 CAFs 中，编码 TGF-β 结合蛋白的基因 TGFB3 和 LTPB1，相较正常粘膜纤维母细胞，表达也上调。
如下图所示，在正常粘膜中，TGF-β 应答基因如 TGFB1，CTGF 和 BHLHE40 主要在骨髓细胞中表达；而在 CRC 中，这些基因主要在上皮细胞和 CAFs 中表达。以上结果表明 CAFs 分泌的 TGF-β 可能促使肿瘤上皮细胞 TGF-β 通路的活化。

从上图中我们也可以看到，CAFs 中，部分 TGF-β 通路上的基因呈现出双峰表达模式。因此，作者利用 self-projection 模式的 RCA 对 CAFs 和正常粘膜纤维母细胞进行聚类，鉴定出三个不同的类别，如下图所示。可以看到，CAFs 分成了两个类别：

对差异表达的基因绘制热图，发现 CAF-B 细胞中表达的肌成纤维细胞标志基因如 ACAT2，TAGLN 和 PDGFA 等在 CAF-A 细胞中表达下调，反之主要表达 MMP2，DCN 和 COL1A2 等基因。这两种不同亚型在多位病人中均有所发现，排除了批次效应的干扰。
综上，CRC 肿瘤中的 CAFs 存在两种亚型，具有不同的转录表达谱。

单细胞的分子表型：细胞干性和 EMT

为了衡量细胞干性，作者基于42个已知结直肠干细胞标志基因定义了一个 stemness index (SI)。如下图所示，肿瘤和正常粘膜中 SI 值的都是连续分布的。肿瘤细胞明显具有更高的 SI 值，且其中5%的 stem/TA-like cells 的 SI 值要比所有正常粘膜细胞都高。

综合 stem/TA-like cells 的 SI 值和基因表达情况，进行相关性分析，作者鉴定出5个与 stemness 相关的基因，如下图所示。其中 GPX2，OLFM4 和 RNF43 均为 WNT 信号调节因子，与已知 WNT 信号在结直肠中所起的作用相符。此外，能够编码 lncRNA ，调控女性 X 染色体失活的 XIST 基因也被发现与细胞干性相关。当仅针对女性肿瘤进行分析时，这一结果得到复现；而在正常粘膜样本中，未发现这一关联。
综上，细胞干性在肿瘤 stem/TA-like cells 中呈现连续分布，且与 XIST 和 WNT 调控因子的表达相关。

肿瘤中上皮细胞能够转换成间叶细胞，实现分离并迁移至新的位点，从而促进肿瘤的转移。从转录层面看，EMT 被定义为上皮细胞标志基因 CDH1 的下调，而来自 TWIST，SNAIL 和 ZEB 家族的间叶细胞转录因子上调。为了研究 CRC 中 EMT 的广泛性，作者检测了上述基因以及其他一些细胞类型标记基因的表达情况，结果如下图：

从图中，我们可以看到，CDH1 在肿瘤和正常粘膜上皮细胞中并未呈现出差异表达，而且这些细胞中 EMT 转录因子的表达也没有很明显。相反，一些EMT 转录因子在 CAFs中表达上调，但是这种表达变化与 EMT 无关，因为纤维母细胞是间质充细胞。
为了量化这些观察，作者通过联合上图中9个 EMT 上调表达的基因定义了一个 EMT score。结果如下图所示，证明肿瘤上皮细胞中缺乏 EMT 特征，而在 CAFs 中，这一特征显著提升。

作者还通过对肿瘤细胞进行单分子 FISH 分析排除组织分离时引入的人为错误。该技术发现 EMT 标志基因 SNAI1，ZEB1 和 TWIST1 与纤维母细胞标志基因 SPARC 表达一致，且未在表达 EPCAM 的上皮细胞中出现。此外，组织免疫荧光和对已知上皮细胞和纤维母细胞的生信分析也进一步证明了之前的结论。

单细胞转录组谱将 CRC 肿瘤分成不同亚组，对应不同的患者存活率

过去，结直肠癌根据混池转录组数据分成不同亚组，表现出不同的治疗反应和患者存活率。这样的分类方法受基质细胞类型转录组变化的影响。因此，作者希望能够探索 CRC 肿瘤中主要几种类型细胞对混池转录组谱和患者存活率的影响。首选，作者将6种类型的肿瘤细胞的单细胞转录组取平均构建出对应的表达谱，然后汇编了来自6个相互独立的队列的混池转录组数据，并将这些数据映射到6种类型细胞的表达谱上，进行聚类，结果如下图。

每个队列中均包含三种不同的肿瘤细胞组：S1，S2 和 S3。S1 细胞具有微弱的上皮细胞特征，明显的纤维母细胞特征以及升高的骨髓细胞特征；S2 所有特征水平均呈中间状态；S3 具有明显的上皮细胞特征以及微弱的纤维母细胞和骨髓细胞特征。
作者针对三种亚型进行了预后情况的分析，发现 S3 肿瘤的存活率最高，与之前具有纤维母细胞特征的病人的研究结果一致。
此外，作者发现之前根据混池转录组数据分类得到的 enterocyte (肠上皮) 和 goblet-like (类杯状) 亚型中均包含 S2 和 S3 两种亚型（如下图）。在类杯状亚型中，S3 肿瘤的存活率要明显优于 S2；肠上皮亚型中，这种差异也存在，只不过没有那么显著。

以上结果进一步强调了 CAFs 与 CRC 结果之间可能存在的相关性，证明根据单细胞转录组得到的细胞类型特异性图谱的预后价值。

Discussion

现有的算法未能准确地进行聚类，可能是因为它们倾向于将批次效应和技术性差异误判为生物学差异；而使用 RCA，能有效降低这些因素带来的干扰。但是我们也看到依赖于外部参考，可能导致无法发现细胞亚型之间的微小差异，所以对于寻找细胞亚型或者不同细胞状态，RCA 的 self-projection 变式可能更适合。不过随着参考数据集越来越大，越来越多样化，预计 RCA 的分辨率也将会不断提升。
利用 RCA，我们从正常粘膜和 CRC 肿瘤中均鉴定出7中主要细胞类型，此外针对上皮细胞还确定了9种亚型，其中7种是第一次从正常粘膜中鉴定出来。利用这9种亚型作参考来推断肿瘤上皮细胞的分化类型，发现 stem/TA-like cells 在肿瘤样本中高度富集。相信正常粘膜单细胞转录组数据集除了能揭示 CRC 功能状态，对未来人类结直肠生物学和相关病理学的研究也是一个很有价值的资源。
传统利用混池测序寻找差异表达基因可能会受到正常和患病样本中细胞组成情况不同带来的干扰，在本次实验中也确实发现这一问题。所以相较而言，单细胞测序是更好的选择。
单细胞测序的另一个优点在于能够研究肿瘤通路相关基因在不同类型细胞中特异性表达情况。通路富集分析和基因 panel 分析表明 CRC 中 TGF-β 信号的增强，可能是由于上皮细胞和基质细胞之间交互驱动的。单细胞转录组和 IHC 分析表明 CAFs 可能通过分泌 TGF-β3 配体和 TGF-β 激活因子 MMP2 促进 TGF-β 信号。CAFs 和肿瘤上皮细胞中 TGF-β 应答基因的表达上调表明肿瘤内可能同时存在自分泌和旁分泌应答。
RCA 聚类揭示 CAFs 中存在异质性，这是首次利用单细胞转录组数据鉴定出 CAFs 亚型。之前有很多临床试验直接或间接靶向针对表达 FAP (fibroblast activation protein α) 的细胞，因为这是一种仅在 CAFs 中表达的细胞膜丝氨酸蛋白水解酶，但是效益未被证明。本文献发现只有 CAF-A 亚型表达 FAP，表明 CAF 一致性可能对靶向 FAP 的治疗产生妨碍。
之前有两项研究表明混池结直肠肿瘤中 EMT 相关基因的表达上调是由纤维母细胞中表达模式的变化引起的，而非上皮细胞。本实验结果与之一致，肿瘤上皮细胞中 EMT 相关基因未表达上调，而 CAFs 中 EMT 标志基因表达显著上升，这可能是因为细胞转向活化状态的纤维母细胞导致。综上，结直肠肿瘤中 EMT 发生可能并不普遍，但在特定条件下或者肿瘤特定亚区中，仍可能存在 EMT。
基于混池转录组数据的 CRC 肿瘤分类已经建立，且能够用于预测患者存活率。本文献发现将混池肿瘤转录组数据映射到细胞类型特异性转录组上，现有分类系统可以进一步改善，表明 RCA 能够基于病人混池数据进行分类。相信随着单细胞转录组成本的不断下降，未来能够利用单细胞数据对 CRC 进行分类。

REF

Li H, et al. Reference component analysis of single-cell transcriptomes elucidates cellular heterogeneity in human colorectal tumors. Nature Genetics, 2017, 49:708-718. doi: 10.1038/ng.3818.