单细胞入门-听一篇Cell报告
这篇文章是 Landscape of Infiltrating T Cells in Liver Cancer Revealed by Single-Cell Sequencing
第一个单细胞尺度的肝癌浸润性T细胞图谱
年前一作的郑春红博士在VIPaper有个在线讲座，我前一晚啃paper时不幸直接睡到开讲，听得云里雾里~
In Brief
对肝癌患者的T细胞群体作了单细胞测序分析，看到了浸润性淋巴细胞的独特亚型和克隆扩增。
Highlights
· 单细胞RNA-seq揭示肝癌T细胞的独特功能组份
· 浸润调节T细胞和耗竭性CD8 T细胞在肝癌的克隆富集
· LAYN基因与肿瘤调节T细胞以及耗竭性CD8 T细胞的抑制功能有关
· 结合表达和基于TCR的分析显示T细胞子集间的连通性
SUMMARY
在癌症里面，对肿瘤浸润淋巴细胞（TILs）的系统性问诊是发展免疫疗法和预测临床反应的关键。文章作者从六个肝癌患者的外周血、肿瘤、癌旁组织分离出5063个T细胞，做了深度单细胞RNA测序。这些单个细胞的转录表达谱结合有T细胞受体（TCR）序列，根据分子和功能性质鉴定出11个T细胞子集，描绘了它们的发展轨迹。特定的子集比如耗竭性CD8^+ T细胞和调节T细胞在肝癌中优先富集，可能存在克隆扩增。他们鉴定了每个子集的信号基因，其中一个基因，layilin（LAYN），在激活的CD8^+ T细胞和调节T细胞中上调，抑制CD8^+ T细胞在体外的功能。这些转录组数据的概况为理解癌症的免疫环境提供了有价值的见解和丰富的资源。
研究背景
癌症免疫疗法在过去数十年里极大促进了肿瘤的治疗。但是这些疗法对于不同的癌症患者以及癌症类型，效率并不一致。目前最重要的一点是找到预测治疗反应的稳定生物学标记，这需要深入理解肿瘤T细胞，特别是导致CD8^+ T细胞耗竭和调节T细胞积累的机制和通路，才能为协调免疫系统消灭癌症提供更好的策略。
浸润淋巴细胞的单细胞分析使得在高度复杂的肿瘤微环境中了解到这些细胞的细节成为可能。scRNA-seq是一个很强大的工具，可以明确定义出每个细胞的TCR序列，从而识别肿瘤特异的新抗原。尽管TCR库相当庞大，鉴定TCR序列可以解释T细胞的克隆扩张模式和相关谱系。
肝癌是致死人数最高的几种癌症之一，在中国，感染HBV是引发肝癌的主要原因。目前缺乏免疫治疗的成功案例，治疗选择有限。肝癌肿瘤导致TILs水平的显著提高，经推测这些TILs已经丧失了杀死肿瘤细胞的能力。因此，肝癌肿瘤的研究可以提供一个描述TILs调节异常的绝佳模型。
实验&分析
Study Overview
实验方法
对于6个未经治疗的肝癌患者，首先做了外显子组DNAseq和大量细胞RNAseq，寻找CNV和体细胞变异，发现基因模式跟已有的肝癌研究是吻合的。并且进一步确认了存在不同的T细胞亚型。
然后使用FACS分选出感兴趣的细胞群体，CD8+和CD4+，做单细胞RNA-seq。
单细胞收集
外周血单个核细胞使用HISTOPAQUE-1077解决方案。
新鲜的肿瘤和邻近正常组织样本先切成1立方毫米的小块，经研磨过滤悬浮去红重悬浮得到。
我以前处理过很多外周血和组织，去除红细胞的步骤是很重要很重要滴
单细胞分类、反转录和扩增
使用FACS将细胞分选到有裂解液的96孔板里面，一个板子上都是同一类的细胞，裂解液里面包含了dNTP，OligodT引物以及加了RNA酶抑制剂的TritonX-100。分选后的板子冻存在-80℃，有两个病人还加了外源RNA质控。单细胞转录组扩增是按照Smart-seq2的流程走的，cDNA纯化用Beckman的XP磁珠。在第一轮磁珠纯化后，每个细胞的cDNA做CD3D qPCR和片段分析。
鉴定出高质量的单细胞，进一步纯化去除小片段，Qubit做浓度质检，使用Vazyme的试剂构建文库。患者P0205、P0508 和 P0322使用HiSeq2500 测序（100bp读长），患者P0407 和P1116使用Hiseq 4000测序（150bp读长）。而对于患者P1202，每个单细胞是人工用吸管挑到孔里并且转录组扩增用的是Tang2010方法，因为明显的3'端偏倚，这些细胞的TCR序列不会组装到单细胞转录组数据上，该患者没有做大量细胞测序。
测序的实验设计很周到，这种算珍贵样本，提前考虑一些可能出现的异常数据
大量细胞的DNA/RNA分离和测序
使用QIAGEN的Mini Kit提取外周血和组织样本的基因组DNA，Qubit测浓度，琼脂糖电泳质检。使用Agilent的试剂构建外显子文库，Hiseq4000测序。
对于RNA，肿瘤和邻近正常组织经手术切除放在冰上，先用 RNAlater RNA stabilization reagent 保存。用QIAGEN的RNeasy Mini Kit释放出RNA，NanoDrop测浓度，AATI片段分析质检。然后使用NEB的试剂构建文库，Hiseq4000测序。
免疫组化
人体组织标本经过相关审核批准由北京世纪坛医院提供。肿瘤切除后福尔马林处理48小时，之后30分钟内收集，脱水包埋使用常规方法。
离体T细胞激活和FACS分析
培养外周血单个核细胞（PBMCs），在指定的时间点采集细胞并染色，使用LSR-II分析仪和Flowjo软件分析每一个存活/单线态亚群的LAYN基因和T细胞激活/耗竭标记的表面表达。所有数据来自至少三组独立实验和6个不同的PBMCs源。
LAYN在原代CD8^+ T细胞的过表达
把人layilin基因克隆到逆转录病毒载体，pAmpho包装，感染从人PBMCs分离的T细胞，免疫磁珠激活。使用FACS分离GFP+的人layilin过表达CD8^+ T细胞。
定量和统计分析
测序深度
5063个单细胞，平均每个细胞129w单一匹配读长（uniquely mapped reads）。
Smart-seq和Smart-seq2的原始文章建议在300w以上，不同的项目需求不同，一般超过100w就比较稳定了
为了说明在这种测序深度下的检测结果是可信赖的，做了一个饱和分析，主要是看重要的基因类别，特别是细胞激素和转录因子。
单细胞RNAseq数据处理
从Hiseq2500/4000下机的数据，首先过滤掉低质量的reads：1）有10%读长的碱基无法识别；2）50%读长的碱基质量小于5；3）含有接头序列。过滤后通常能剩下90%的原始reads，先使用GSNAP工具把这些reads比对到RFam下载的核糖体RNA序列。
没有匹配上的reads（clean reads）再拿去比对参考序列，包括UCSC下载的hg19和spike-in的质控RNA。另外，从UCSC下载了一个基因模型文件knownGene.txt，用可调参数“–novelsplicing 1 -n 10 -i 1 -M 2”来定义已知的exon-exon junctions，其他参数使用默认设置。使用R的函数findOverlaps检索跟已定义基因重合的reads，得到每个基因的读长计数。
使用DESeq方法对读长计数作进一步的标准化。同一个患者，首先计算每一个基因在所有细胞的读长计数的几何平均数，得到读长计数与其几何平均数的比值（只计算非零项）。然后得到比值的中位数作为每个细胞的大小因数（size factor），接着把标准化的读长计数（$RC_{normalized}$ ）定义为RC/size-factor。$log_2 (RC_{normalized} +1)$的值用于下游的分析，比如PCA和聚类。
还使用下面的步骤过滤掉低质量的细胞：1）去掉单一匹配读长少于25w的细胞；2）去掉CD3D基因的RPKM小于1.0的细胞；3）去掉FACS分选的CD8^+ T细胞里CD8A基因的RPKM小于1.0或者CD4基因的RPKM大于10的细胞；4）去掉FACS分选的CD4^+ T细胞里CD4基因的RPKM小于1.0或者CD8A基因的RPKM大于10的细胞；5）去掉离群点分析(outlier analysis)失败的细胞。
3)和4)是不信任FACS么，还是别的原因？
在离群点分析中，PCA分析使用了$log_2 (RC_{normalized} +1)$，拟合分析使用了方差平方系数（CV^2 ）、每个FACS-分选的T细胞亚型的$RC_{normalized}$ 的平均值，手动检查结果以确定异常细胞。典型的异常细胞定位在其他细胞很远的地方，通常对主成分有很高的贡献。此外，异常细胞也表现出很低的大小因数，表示细胞里的RNA含量很少，通过拟合曲线可以更好地移除这些细胞。过滤掉异常细胞之后，6个患者总共剩下4128个细胞，其中5个HBV阳性患者共有4070个细胞。
大量细胞的转录组和外显子组数据处理
大量细胞的RNA-seq数据处理跟单细胞的RNA-seq数据使用同一套GSNAP流程。对于外显子测序数据，首先过滤掉低质量的reads。高质量的reads使用bwa-0.7.12回贴到人的参考序列（从Broad下载，加了诱饵序列的b37版本）。根据建议的最佳做法使用GATK作比对后处理，包括bam文件分类、duplicate reads标记，候选INDEl附近的reads本地重排，以及碱基质量重新校正。使用Strelka的默认参数分析体细胞突变，包括SNV和INDEL。那些所谓的突变特别是低频突变，使用IGV工具作进一步的手动检查。用ADTex得到体细胞拷贝数变化和杂合缺失（LOH）。根据与已有的驱动基因（driver gene）名单的比较，进一步得到这些患者的候选驱动突变（driver mutations）。一些已知的突变热点：TP53 G245C (P0205), PIK3CA H1047R(P0508), PTEN C136R (P1116), TP53 R249S P0322) ，TP53 I195F (P0427)。后期患者 P0322有两个不同的TP53突变。
TCR分析
使用TraCeR给单细胞RNA-seq数据的每一个T细胞组装上TCR序列，进而鉴定出CDR3序列、重排的TCR基因以及它们的表达丰度（TPM）。首先，去掉没有明显TCR组成的细胞。然后按照后面的步骤分别整理TCR的α和β链，第一条TCR的α（β）链定义为：1）所有单个T细胞中，保证只能存在一条有用的α（β）链。2）如果一个T细胞中不只一条α（β）链，保留只有一个α（β）链的主要形态被检测到的细胞。通常来讲，一条α（β）链是有用的而其他的都没用，或者说其中一条的表达水平远远高于替代等位基因（alternative allele），以此可以识别出来。接下来，根据表达分布的直方图分析
过滤掉那些第二条TCR的α链TPM小于10或者β链TPM小于15的细胞。从成功组装TCR序列的4032个细胞里面，找到3792个细胞有TCR的α/β对。
非监督聚类
每个患者的给定基因的相对表达量通过Y-mean（Y）计算，Y是标准化的表达量$log_2 (RC_{normalized} +1)$，mean(Y)是这个患者所有细胞的Y的平均值。这个过程移除了由患者差异引起的不想要的变异。这种集中数据随后结合所有的患者。前n个高标准差的基因被认为是高变异基因，试验过不同的n值：1500、2000、2500和3000。所有基因产生的集中表达数据用于改进的SC3聚类流程。首先创建基于相关分析（Spearman correlation）的距离矩阵，然后计算拉普拉斯算子的特征向量。接着在最初的d特征向量（d取不同的值，从细胞总数的4%到7%）多次运用k-means算法。平均不同次运行的结果得到一个一致性矩阵（consensus matrix），从中可以看出两个细胞在不同次运行中聚到一起的频率是多少。最后在一致性矩阵执行完整的层次聚类（hierarchical clustering with complete agglomeration），得到k群（k clusters）。用于k-means和层次聚类的SC3的参数k，从2迭代到10。每一次运行SC3，都要计算剪影（silhouette），绘制一致性矩阵，鉴定集群特异基因（cluster specific genes）。所有这三个方面可以帮助研究者根据经验得到最优的k和n。一旦确定出稳定的集群，使用的程序将被迭代到每一个集群，用以揭示每个集群中细胞间的高变异基因，然后就可以使用这些变异基因来识别亚群。
中间数据要用眼睛去看，要根据经验去选，跑一个分析流程远不是那么简单
得到11个稳定的集群，包括5个CD8^+ 集群和6个CD4^+ 集群，每一个集群都有独特的标记基因，CD8^+ 在不同患者间可以比较。
得到聚类结果之后，用ANOVA的R函数aov识别差异表达基因（differential expressed genes），再用R函数TukeyHSD的杜凯氏方法（Tukey's range test）实现每个集群对的差异检验。每个集群的差异表达基因符合这些标准：1）FDR调整后的F检验P值 < 0.01；2）研究的集群与其他集群间的log2倍数变化绝对值>1；3）比较 研究的集群与其他集群 的杜凯氏检验的P值<0.01。t-SNE方法用作降低的2维空间里细胞距离的可视化。
此外，细胞周期也被认为是细胞群体分析的重要混杂因素。用ScLVM检查了细胞周期的贡献度，发现对这里的聚类结果影响非常小，细胞周期因素影响的方差只占生物学因素的4.5-13.6%（患者P0508, P0322, P0407, P1116）和37.7%（患者P0205），后面这个也远低于之前报道的85%。聚类结果确实包含了一些细胞周期相关的基因，像TUBB4B、TUBA1B、MCM7、STMN1和MKI67，但这些只有62个细胞，不是主要的亚群，也没有进行tSNE可视化以及差异基因表达分析。
基于综合 表达和TCR克隆 的CD8^+ 和CD4^+ T细胞的状态转换分析
CD8^+ /CD4^+ T细胞在二维状态空间的拟时间（pseudotime）排列使用Monocle2定义。细胞的顺序由CD8^+ /CD4^+ T细胞群（无MAIT）的最分散基因的表达推断而来。每个点代表一个细胞，每种颜色代表一个集群。
CD8^+ /CD4^+ T细胞集群的细胞状态转换由共有的TCR推断而来。只有在不同细胞集群共享TCR的才会着色，每种颜色代表一种独特的TCR克隆类型，线条连接基于不同的TCR共享程度，线条厚度代表共享的TCR数目。插入图片显示两个集群的共享TCR数目，条形中的两种颜色代表特别的两个集群共享的TCR。
T细胞耗竭状态分析
为了揭示耗竭性T细胞的特征，首先移除了一些肿瘤浸润T细胞，它们与正常组织的T细胞带有相同的TCR序列。MAIT细胞因为特殊的特性也被移除掉。根据SC3的聚类分析，C4_CD8-LAYN 集群中的肿瘤浸润CD8^+ T细胞被定义为耗竭性T细胞，其他的则是非耗竭性。耗竭性与非耗竭性T细胞间的差异表达基因检测使用R包limma的线性模型和经验贝叶斯方法实现，使用严格的显著性阈值，调整后的P值（BH多重检验校正）<0.01，倍数变化≥4。设计严格的阈值是为了降低假阳性找到可信的耗竭性标记。使用limma是因为它对高度平行的基因组数据的杠杆作用，可以借用gene-wise模型间的信息。计算每个患者每个DEG的相对表达（log2 fold-change），计算每个患者所有DEG的相对表达均值，从而得到每个患者的耗竭性分数（per-patient exhaustion score）。耗竭性分数随着临床分期靠后而增加。肿瘤调节T细胞特异的基因使用同样的方法检测，但是采用不同的阈值（adjusted p value<0.01，fold change≥2）。
TCGA数据分析
肿瘤基因组图谱（TCGA）是美国NCI和NHGRI的联合项目。
TCGA数据用来检测选取的基因与患者存活率的关联，存活分析使用R包survival实现。
文章的raw data 和 processed data分别是EGA: EGAS00001002072 和 GEO:GSE98638。
http://hcc.cancer-pku.cn 有这些数据的分析结果。
读了这篇文章的收获是知道了具体的单细胞项目可以是这样的。单细胞测序也不是完全推新的技术，是现有技术的改进和细化。很多分析的套路都是延用的，另需要照顾单细胞的数据量和异质性。怎么设计质控，丢掉哪些数据，使用哪些数据，如何整合出一个可用的模型呢？
菜:bird:表示还有n多基础的东西要刷~
点击可以加入单细胞数据处理学习交流小组