结直肠癌细胞异质性大揭秘！
文献详解栏目
每个人的时间精力有限，必须优先阅读相关文献，开设这个栏目也是希望为大家推荐高质量的单细胞相关文献。如果大家对单细胞转录组感兴趣可以关注一下，哪怕每天只学一点点，积土成山，积水成渊。
当然一个人的力量终归是小的，我也希望汇聚一群人，形成一个场，这里头最重要的生产力不是单个人多聪明，多厉害，而是每个人相互作用，形成的那个氛围。
希望大家能有所收获！
文章信息
在之前的单细胞转录组研究中，许多实验检测的肿瘤微环境分辨率不足，或者忽视了肿瘤与对应的正常样品的比较。并且许多单细胞转录组分析算法则受到高背景噪音、批次效应、技术误差的干扰（原因主要是临床样本采集的条件、批次不同），限制了分析的准确性。
开发了一种叫做RCA(reference component analysis)即参考成分分析的新型算法，相比旧算法可以显著提高聚类准确性，对细胞类型、肿瘤微环境的识别率提高。利用这个方法，发现了两种独特的成纤维细胞（CAFs）亚型，其中上皮间质转化（EMT）相关基因显著上调表达；将结直肠肿瘤细胞依据不同的细胞状态和生存概率分成不同的亚型，更准确的细胞分类对预后有一定的指导意义。
背景知识
◆ ◆ ◆ ◆ ◆
肿瘤异质性
肿瘤异质性包括肿瘤间异质性（不同肿瘤细胞之间的基因与表型不同）和肿瘤内异质性（相同肿瘤细胞以内的基因与表型也不同），其中肿瘤内异质性又粗略分为空间异质性（未扩增的细胞背景中有成簇扩增细胞；少量扩增背景中有未扩增的细胞；孤立的细胞扩增）与时间异质性（原初肿瘤与次生肿瘤）。
Nature一篇文章按照其来源分为患者异质性 (patient heterogeneity) 和瘤内异质性 (intratumoral heterogeneity)，其中患者异质性就是指同种类型肿瘤在不同患者中表现出差异，瘤内异质性指同一患者的同一部位肿瘤内存在的明显差异。
瘤内异质性一直是研究重点，异质性是产生抗药性的主要原因，而异质性产生的主要原因又与肿瘤干细胞的发生、遗传物质不稳定、细胞间竞争等相关。
肿瘤干细胞：是指肿瘤中存在的一小群具有无限自我更新能力的干细胞样细胞,可以演化成表型不同的细胞群,具有强致瘤性。受细胞分化、克隆演化和微环境等因素影响,导致表型与功能的异质性。这也是现在关于异质性来源最为普遍的的观点
遗传物质不稳定：肿瘤细胞是异倍体，一篇2013年Nature在结肠癌中发现3个与DNA复制相关的基因MCD4、MEX3C和ZNF516，基因的缺失导致DNA复制时出现应激=》异倍体出现=》异质性肿瘤细胞出现
◆ ◆ ◆ ◆ ◆
上皮间质转化
EMT（epithelial-mesenchymal transition）是Greenberg和Hay于1982年提出的，指的是可逆的上皮细胞转化为间质细胞 (2014 Trends Pharmacol Sci)的生物学过程。上皮细胞将经历持续的细胞活动, 包括失去顶-底极性结构, 细胞间连接受到破坏, 细胞骨架重构, 改变细胞形态并最终呈现出间质及侵袭表型, 从而增加细胞运动性并提高其降解细胞外基质的能力（2009 Cell）。
恶性肿瘤中绝大多数是上皮性肿瘤，发生EMT的恶性肿瘤细胞将获得增强的迁移及侵袭能力并侵入周围的细胞外基质, 最终向远处位点转移，在结肠癌、甲状腺癌及乳腺癌的侵袭表型中发挥重要作用。发生EMT的细胞具有抵抗细胞凋亡、衰老、化疗和免疫治疗的能力, EMT通过诱导免疫耐受使恶性肿瘤细胞逃避免疫监视。也有研究表明EMT与肿瘤干细胞的特性相关(2016，Sci Report；2016，Sci Report；2015，Cancer Letter)
方法
◆ ◆ ◆ ◆ ◆
样本选择
选择11例 II-IV期结直肠癌样本，由新加坡总医院的新加坡国立癌症中心提供。原发肿瘤手术切除后，置于无菌、低温的RPMI培养基中培养。
◆ ◆ ◆ ◆ ◆
选出单个细胞
分离得到的单个细胞生活力由Calcein-AM and Ethidium Homodimer 1 (Live/Dead Kit, Life Technologies) 评估，然后加到 RNA–seq IFC-C1芯片上进行细胞捕获，然后基于 bright-field and Calcein-AM imaging进行芯片成像，只有真正单个的细胞可以继续进行文库制备和测序。
◆ ◆ ◆ ◆ ◆
文库制备和测序
根据cDNA浓度和质量丢弃低质量的单细胞后，利用Illumina 的Nextera XT DNA样本制备试剂盒进行单独的文库制备，然后采用Hiseq2000平台101bp PE方案对总共1591个CRC分离的肿瘤细胞和正常组织配对细胞，以及7个细胞系的630个细胞进行测序。
◆ ◆ ◆ ◆ ◆
测序数据预处理
原始fq数据利用Tophat 2.1.0比对到hg19基因组，利用的是GENCODEv19注释文件(Tophat参数：--read-edit-dist设置为3，--read-realign-edit-dist设置为0)。max-multihits参数设为1，保证BAM文件中只有一个uniquely mapped reads。表达定量使用Cuffdiff-2.2.1，设置参数--frag-bias-correct 以及--library-norm-method 选择FPKM标准化。为了减少rRNAs, tRNAs以及小RNA（snRNAs and snoRNAs ）对FPKM的影响，需要利用Cuffdiff的--mask-file参数，最后导出原始的reads count矩阵。
◆ ◆ ◆ ◆ ◆
质控过滤
计算几个统计值：
FASTQR： the number of paired-end sequences in the FASTQ file
BAMMR（Mapped reads in BAM）：the number of uniquely mapped reads in the aligned BAM file
ER（Exonic reads ）：summing raw read counts across all genes
ROER (rate of exonic reads)：the ratio between ER and BAMMR
NODG(number of detected genes)： total number of genes with FPKM ≥1
加入几个看家基因作为参考：TFRC, ACTB,RPLP0, PGK1, GAPDH, LDHA, NONO, B2M, GUSB and PPIH
选择NODG>1000, ROER>5%, ER>0.1M的细胞，总共626个CRC细胞与正常组织配对细胞，561个来自细胞系的细胞。
The number of detected genes (NODG) 与the rate of exonic reads (ROER) 比值。 颜色表示看家基因表达量的log10 (FPKM).
◆ ◆ ◆ ◆ ◆
对细胞系scRNA-seq数据集的聚类方法评估
使用7种细胞系的原始630个细胞(其中561个通过了质控过滤)；为了评估批次效应的处理，其中包含了两个批次： GM12878 (lymphoblastoid) cells 和 H1 embryonic stem cells。
八种方法：All-HC、HiLoadG-HC、BackSPIN、RaceID2、Seurat、VarG-HC、VarG-PCAprojHC、VarG-tSNEproj-HC
All-HC：选择FPKM大于等于0.001且至少存在于2个细胞中的基因，用log10(FPKM)值来聚类，使用
‘average’-linkage neighbor joining算法HiLoadG-HC：选择PC1、PC2、或PC3的排名前100或后100基因，利用log10(FPKM)值进行层次聚类，使用
average linkage算法BackSPIN：使用参数
-f 2000 -v -d 4对细胞进行聚类，这样会比默认参数得到更精确地类群划分RaceID2：使用默认参数对原始FPKM进行聚类，不需要对数据取log值降低维度(改变默认参数结果没有明显改进)，使用
within-cluster-dispersion算法，聚类数k=7Seurat：使用
RegreeeOut函数的latent.vars参数设定为nUMI，使用model.use参数设定为linear(当model.use = negbinom是结果相似)VarG-HC：使用BackSPIN挑选基因的方法，选出前1000变化最大的基因，然后进行层次聚类，表达量值也是取log10以后的
VarG-PCAproj-HC：方法与VarG-HC相同，只是选择基因后利用PCA进行降维，保留表达矩阵中90%的方差
VarG-tSNEproj-HC：方法VarG-HC相同，只是后来使用t-SNE降维
评级聚类准确度的指标为ARI(adjusted Rand index)，0-1之前，数值越大表示聚类效果越好
◆ ◆ ◆ ◆ ◆
利用黑色素瘤的scRNA数据集对算法评估
数据集在GSE72056，其中有4645个黑色素瘤细胞的表达数据，主要研究6种细胞类型以及根据marker基因过滤得到的细胞：T细胞(CD2, CD3D, CD3E, CD3G)，B cells (CD19,CD79A, CD79B, BLK), macrophages (CD163, CD14, CSF1R), endothelial cells (PECAM1, VWF, CDH5), CAFs (FAP, THY1, DCN, COL1A1, COL1A2,COL6A1, COL6A2, COL6A3) and malignant melanoma cells (MIA, TYR,SLC45A2)。排除了没有检测到任何marker的细胞，其余细胞在marker表达量的基础上进行层次聚类，最后将4057个细胞划分成了6类
结果
◆ ◆ ◆ ◆ ◆
RCA比现有的聚类算法表现更优
颜色代表7种不同细胞系，形状表示H1 and GM12878数据的批次；上面是聚类树，下面是PCA(RaceID2和Seurat使用的是tSNE)
尽管All-HC可以将大多数细胞准确地分开，但还有约15%的细胞聚类不正确(ARI=0.66)，另外，H1和 GM12878 都按批次被分开，也就是说同一批次没被放在一起；
HiLoadG-HC的效果不好(ARI=0.53)，有大量细胞没有正确聚类，然后也是H1和 GM12878 都按批次被分开；
BackSPIN的ARI结果和All-HC相似；
RaceID2的准确度要低许多(ARI=0.15)，受批次效应影响更大；
Seurat相对于All-HC有一些改善(ARI=0.70)，但是也存在批次影响，大量细胞分类错误；
另外三种VarG方法(图片在附件中)的ARI从0.0005-0.13，效果比较差；
RCA得到了紧密的细胞群，而且几乎都由相同类型细胞组成，另外即使在不同批次中的相同类型细胞也会聚在一起(ARI=0.91)
另外利用以上算法分析了黑色素瘤的scRNA数据，结果也是RCA具有近乎完美的准确性，ARI远超其他算法
◆ ◆ ◆ ◆ ◆
利用RCA鉴定了CRC肿瘤细胞多种细胞类型
分析了11个CRC患者的969个细胞以及作为对照的7个患者正常组织的622个粘膜细胞。严格过滤得到375个肿瘤细胞和215个正常粘膜细胞
左图是RCA对正常粘膜细胞聚类；右图是对CRC肿瘤细胞聚类其中聚类图表示细胞类型聚类；上面的热图中，列表示细胞，行表示RCA中的reference panel，颜色表示数据投射到reference panel 的分值；中间的热图表示原始log10(FPKM)在各个marker中的表达量；最下面的热图中行为不同的患者，每行的黑色小块表示该患者的细胞
利用RCA global panel模式一共得到CRC与正常细胞的7种细胞类型：上皮细胞（ epithelial cells）、纤维原细胞（fibroblasts）、内皮细胞（endothelial cells,）、B细胞、T细胞、肥大细胞（mast cells）和髓系细胞(myeloid cells )。其中上皮细胞并没有继续分亚群(比如enterocytes, goblet cells and transit-amplifying (TA) cells)。
为了测试是否为RCA算法导致，又使用了RCA的self-projection模式，选择了上皮单细胞转录组数据作为reference panel，结果正常的上皮细胞分成9种亚型，没有发现批次效应，并且分出的亚型和上皮细胞的marker也是可以对应上的，说明了确实可以得到不同的细胞类型或细胞状态。
关于这两种RCA的模式的具体算法 ，可以在文中methods出查到
对肿瘤上皮细胞进行分群，结果得到三个亚群： stem/TA-like, enterocyte 2B–like and goblet-like。其中 stem/TA-like占到93%，而之前正常粘膜上皮细胞中只有30%，这也与CRC细胞的增殖特性一致。
因此，即使临床样本存在批次效应，但是正常粘膜与CRC细胞中依然可以较为准确地识别细胞类型。
◆ ◆ ◆ ◆ ◆
揭示CRC肿瘤差异表达基因
先进行一个假设：scRNA与bulk RNA对肿瘤-正常组织得到的差异基因存在显著差异。对scRNA的结果分析：从肿瘤样本中选择了5个最大的细胞群( stem/TA-like epithelial cells, fibroblasts, B cells, T cells and myeloid cells)，每个细胞群与正常粘膜细胞群比对，共得到129个差异基因(FDR<0.05)。其中 stem/TA-like cells 和fibroblasts的差异基因数量最多。所有的差异基因放在了：https://media.nature.com/original/nature-assets/ng/journal/v49/n5/extref/ng.3818-S1.pdf 的 Supplementary Note中，并且有许多在之前的CRC功能研究中被报道。总体来说，大部分在scRNA中检测的差异基因在bulk RNA中没有检测到，从另一个方面展示了scRNA的优越性
a-d分别是 stem/TA cells and tumor stem/TA-like cells 、normal fibroblasts and tumor fibroblasts 、normal and tumor myeloid cells、normal and tumor B cells的 bulk RNA和scRNA差异基因绿色是bulkRNA检测的DEGs，灰色是bulkRNA检测的non-DEGs，红三角是scRNA检测的上调DEGs，蓝三角是scRNA检测的下调DEGs
◆ ◆ ◆ ◆ ◆
通路的改变和CRC中的CAFs多样性
通路改变
利用IPA(Ingenuity Pathway Analysis)鉴定了四种细胞类型的上游调控元件，发现前3名的元件包含了2个小分子：放线酰胺素Actinonin和莫特沙芬钆 Motexafin gadolinium。Actinonin在多种癌症细胞（包括结肠直癌）中起到抑制生长的作用，Motexafin gadolinium主要在干细胞中发挥作用，诱导癌细胞凋亡
另外，TGFB1( transforming growth factor B1)是排名最高的内源性上游调控因子，在CAFs激活通路中活性增强，同时SMAD3（ TGF-β effector ）也增强。另外TGFB3和LTPB1(编码 TGF-β结合蛋白)在CAFs中表达比正常黏膜成纤维细胞表达更多。在正常黏膜成纤维细胞中，TGF-β响应基因如TGFBI, CTGF and BHLHE40 主要在pithelial cells 和CAFs 中表达，说明由CAFs分泌的TGF-β可能激活肿瘤上皮细胞中的TGF-β通路。
图b是TGF-β 信号通路；d：正常黏膜与CRC主要细胞类型的TGF-β信号通路中基因表达量热图(采用 log10 (FPKM)
CRC中的CAFs多样性
上图中红色的CAFs的一些TGF-β基因表现出了"双峰表达Bimodal expression"的情况(遗传相同表型多样)，然后用RCA的self-projection模式对CAFs进行聚类，分成了NMFs(normal mucosa fibroblasts)、CAF-A and CAF-B三类，之后CAF-B细胞基因差异分析得到marker基因ACTA2,TAGLN and PDGFA 上调表达，它们在CAF-A中下调；CAF-A中上调的有MMP2,DCN and COL1A2
因此，CRC肿瘤中的CAFs主要有两种亚型，具有不同的转录表达
◆ ◆ ◆ ◆ ◆
分子表型分析干细胞特性和EMT
干细胞特性(stemness)
利用SI(stemness index) 【由已知的42种结肠或肠道干细胞marker基因平均表达水平得到】在肿瘤和正常组织中的分布，结果发现肿瘤细胞的SI值更高，约5%的肿瘤细胞SI值高于正常组织。然后将干细胞的SI值与基因表达量结合，发现了5个基因与干细胞作用相关，其中有3个与WNT信号通路相关(GPX2，OLFM4 and RNF43)。
有趣的是，lncRNA编码基因XIST （可以介导女性X染色体失活）也被发现与干细胞性相关，不过在正常样本中XIST 与SI评分无关。
总的来说，肿瘤干细胞的特性是一个不断变化的过程，并与XIST和WNT介质的表达相关
上皮间质转化（EMT）
之前有一种假设："上皮细胞转化为间充质状态，这样就允许它脱离并转移到其它位置"，转录水平上，EMT的特征是上皮组织的marker CDH1 下调表达，间充质转录因子TWIST, SNAIL and ZEB家族的上调表达。为了检测CRC中的EMT发展情况，检测了以上基因及其它细胞类型的marker基因的表达
橙色为EMT相关的转录因子log10(FPKM)值，灰色部分是一些细胞类型marker基因的表达量
并没有观察到肿瘤和正常粘膜上皮细胞的CDH1的表达差异，而且肿瘤与正常组织上皮细胞都没检测到EMT转录因子。相反，发现一些转录因子在CAFs（图中红色区域）中上调，不过这个变化和EMT的发生没有直接关系，因为fibroblast(成纤维细胞是间充质细胞)。为了进一步研究，文章结合了上图中的9个EMT上调基因的表达量，定义了每个细胞中EMT的评分，确定了肿瘤上皮细胞缺乏EMT信号，而CAFs中EMT信号显著增强。
之后为了排除组织切除带来的人工干扰，又使用了单分子FISH(smFISH)对肿瘤mRNA分子进行可视化，结果发现EMT的marker SNAI1、ZEB1和TWIST1 的表达与成纤维细胞的marker SPARC 表达一致，并且在表达EPCAM的上皮细胞中检测不到。免疫组化与生信分析结果是一致的。
◆ ◆ ◆ ◆ ◆
根据不同患者的生存率将CRC分成亚组
目前已经利用常规bulk转录组根据治疗反应和患者生存率，将CRC分成不同亚型。但是这种分类会受到不同类型间质细胞转录组的影响。因此文章利用了scRNA得到的6这种细胞型，然后结合了6个独立的带有生存信息注释的常规转录组数据：TCGA、GSE14333、PRECOG37以及一项相关研究的三个数据集。将bulk数据投射到6个细胞型上，并在每个细胞型中都鉴定到了三个肿瘤组：S1、S2、S3（图a）。其中S1肿瘤上皮细胞特征较弱，成纤维细胞较强，髓样特征较强；S2所有特征处于中等水平；S3上皮细胞特征较强，其余较弱（图b）。
接着讨论了新定义的CRC亚型的预后意义。6个细胞型中，S3肿瘤生存率最高(图a)，这和之前研究成纤维细胞与患者预后的关系结果一致。来自GSE14333数据集的bulk转录组数据中 ‘enterocyte’ and ‘goblet-like’ CRC肿瘤亚型中包含了这里分析的S2、S3亚型（图c）。在 ‘enterocyte’ 和 ‘goblet-like’ 亚型中，S3肿瘤比S2的存活率更高。进一步强调了CAFs与CRC预后的潜在相关性 ，利用单细胞转录组推断的细胞类型对预后具有一定的价值。
往期精彩
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
单细胞天地欢迎你