微科盟_10X单细胞转录组_分析报告


1. 概述

单细胞转录组测序是一项可以在单细胞水平测定细胞的转录组的技术,利用该技术可以回答样本中每个细胞基因表达情况与细胞间异质性问题,这让解析细胞行为,机制及相互关系成为可能。

10X Genomics 单细胞转录组测序平台利用微流控、油包水以及单细胞标签等技术实现了高通量的细胞捕获。能够一次性标记500-10000个细胞。并结合转录本3’端获取的转录组信息从而实现研究单个细胞基因表达情况的目的。与其他单细胞转录组技术相比,该技术具有通量高,成本低,周期短等显著优势。

该技术可用于细胞分型与细胞标记分子的鉴定并进一步实现对细胞群体的划分与细胞群体间基因表达差异的检测,此外还可以利用该平台产生的数据预测细胞分化与发育轨迹。在发育、免疫、疾病领域的研究中发挥越来越重要的作用。

2. 单细胞建库测序

从组织样本单细胞分离到最终数据获得,样品检测、单细胞分离、单细胞建库、测序等每一环节都会直接影响数据的数量和质量,从而影响后续信息分析的结果。为从源头保证测序数据准确可靠,承诺在数据的所有生产环节都严格把关,从根源上确保高质量数据的产出。建库测序的流程图如下:

图2.1 单细胞建库流程

2.1 单细胞悬液样本制备

根据不同样本的特性将样本制备成单细胞悬浮液,对悬液中的细胞质检包括细胞计数和细胞活性测定,最终悬液中的细胞活性应高于90%,细胞浓度为700~1200个细胞/μl。

2.2 形成与收集GEMs

10X单细胞平台采用基于微滴条形码分配系统,该系统为样本提供了百万级的携带唯一DNA条形码的凝胶珠,再结合微流控技术,将带有条形码和引物的凝胶珠与单细胞悬液中的细胞包裹在油包水的液滴(GEMs)中。

该系统共有约有75万种带有标签的凝胶珠,每个珠子上有40-80万探针。探针的引物序列包含四个部分:

(1) Illumina TruSeq Read 1测序引物。

(2) 16 bp的细胞识别条形码(Barcode)。

(3) 12 bp的RNA独立分子标签(UMI)。

(4) 30 bp的Poly(dT)反转录引物。

Illumina TruSeq Read 1测序引物,用于后续的上机测序。16 nt的细胞识别条形码有400万种细胞识别条形码,一个凝胶珠仅含有一种细胞识别条形码,通过这400万种条形码,从而将其对应的细胞区分开。12 bp的RNA独立分子标签是由随机碱基组成的一段序列,捕获到的每一个核酸分子都将分配有特异的UMI序列,能够区分哪些reads是来自于一个原始RNA分子,从而区分基因片段是否为PCR的扩增重复及区分是真实的单核苷酸多态性(SNP)位点还是PCR产生的突变。30bp的Poly(dT)反转录引物,是含有30个T碱基的同聚DNA片段,用于捕获有polyA尾的转录本。

图2.2 凝胶珠结构

图2.3 GEMs生成过程示意图

2.3 文库构建与质检

2.3.1 文库构建

1.油滴形成后,凝胶珠在油滴内溶解并释放大量探针,细胞裂解释放mRNA,mRNA与油滴内的探针结合并在逆转录酶的作用下产生用于测序的带有细胞标签和UMI信息的cDNA一链,再以SMART扩增方法完成第二链合成。

图2.4 标签添加与一链合成示意图

2.油滴破碎,使用磁珠纯化cDNA一链,然后进行PCR扩增。

图2.5 二链合成示意图

3.cDNA扩增完成后,利用化学方法将cDNA打断成200-300bp左右的片段,并构建含有P5和P7接头的cDNA文库。

图2.6 测序文库构建示意图

2.3 文库构建与质检

文库构建完成后,先使用Qubit2.0 Fluorometer进行初步定量,稀释文库至1.5ng/ul,随后使用Agilent 2100 bioanalyzer对文库的insert size进行检测,insert size符合预期后,qRT-PCR对文库有效浓度进行准确定量(文库有效浓度高于2nM),以保证文库质量。

2.4 上机测序

库检合格后,把不同文库按照有效浓度及目标下机数据量的需求pooling后进行Illumina测序。测序的基本原理是边合成边测序(Sequencing by Synthesis)。在测序的flow cell中加入四种荧光标记的dNTP、DNA聚合酶以及接头引物进行扩增,在每一个测序簇延伸互补链时,每加入一个被荧光标记的dNTP就能释放出相对应的荧光,测序仪通过捕获荧光信号,并通过计算机软件将光信号转化为测序峰,从而获得待测片段的序列信息。测序过程如下图所示:

图2.7 Illumina测序原理示意图

3 信息分析流程

对于单细胞转录组的分析,首先使用CellRanager对测序数据进行质控,比对后根据细胞标签与分子标签产生单细胞基因表达矩阵并根据每个细胞的表达谱特征去除低质量的细胞。使用Seurat软件识别与解释单细胞原始的表达矩阵来中异质性的来源,该步骤主要包括细胞与基因的质量控制,降维与聚类,标记基因鉴定与差异分析等。在得到不同类的标记基因后使用Singler软件对不同的细胞类群进行细胞类型注释。在完成标记基因鉴定与差异分析后使用ClusterProfiler软件对挑选出的基因集进行GO分析与通路分析。最后使用monocle软件根据单细胞基因表达的相似性进行拟时序分析来说明细胞类群之间的关系。信息分析流程图如下。

图3.1 分析流程图

3.1 原始数据质量控制

3.1.1 测序数据说明

测序片段被高通量测序仪测得的图像数据经CASAVA碱基识别转化为序列数据(reads),文件为fastq格式,其中主要包含测序片段的序列信息以及其对应的测序质量信息。

单个样本的10X单细胞测序数据分成I1,R1,R2三个测序文件。 R1文件每条序列包含26个碱基细胞识别条形码,10个碱基UMI。I1文件包含8碱基的样本识别标签。R2文件包含转录本测序信息。

fastq格式文件中每个read由四行描述信息组成,如下所示:

上述文件中第一行以“@”开头,随后为Illumina测序标识符(Squence Identifiers)和描述文字;第二行是测序片段的碱基序列;第三行以“+”开头,随后为Illumina测序标识符(也可为空);第四行是测序片段每个碱基相对应的测序质量值,该行每个字符对应的ASCII值减去33,即为该碱基的测序质量值。

3.1.2 测序错误率分布

测序过程本身存在机器错误的可能性,测序错误率分布检查可以反映测序数据的质量,序列信息中每个碱基的测序质量值保存在fastq文件中。如果测序错误率用e表示,Illumina的碱基质量值用Qphred表示,则有:Qphred=-10log10(e)。Illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系见下表。

表3.1 Illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系

单细胞测序中对应转录本信息的只有R2单端的序列,我们仅对R2端序列进行测序质量检测

在该部分分析中,若样品80%的测序序列错误率在0.1%以下即为合格。下图为其中一个样本作为示例:

3.3 测序错误率分布

3.1.3 GC含量分布

核苷酸序列中鸟嘌呤(G)和胞嘧啶(C)所占的比例称为GC含量。GC含量在物种间存在一定特异性,前几位碱基在核苷酸组成上有一定偏好性,产生正常波动,随后则趋于稳定。下图为其中一个样本作为示例:

图3.4 样本GC含量分布

3.2 Cellranger分析

CellRanger 是 10X genomics 官方提供的一套针对单细胞 RNA 测序输出结果进行比对、定量、聚类及基因表达分析的分析流程。

CellRanger对单细胞测序数据进行以下5步处理:

(1) 对测序数据进行修剪,去除poly-A 尾与TSO 序列结构。

(2) 使用 STAR 对测序数据进行比对并判断reads 是比对到了外显子、内含子还是基因间区上。

(3) 在定量前,对UMI的测序错误进行矫正。

(4) 依据UMI进行定量。

(5) 过滤低质量细胞。

3.2.1 CellRnager分析结果

CellRanger分析结果主要包括:

(1) CellRanger网页版分析报告。包括异常结果警报,细胞和基因数目的评估,比对比例统计及部分基础分析的结果。

(2) CellRanger定量单细胞原始表达矩阵与过滤后的表达矩阵。

3.2.1.1 CellRnager网页版分析报告

CellRanger网页版的内容包括以下几个部分,示例如下:

(1) 异常结果警告

如果数据中存在异常,在网页的头部会给黄色的警告框,显示检查到了某些错误。点击detail可以查看详情。

(2) 细胞与基因数目定量结果

细胞和表达的基因个数定量结果,具体包括样本测到的细胞数,每个细胞平均序列数,每个细胞基因数的中位数,总序列数,校正后的有效细胞标签数目,有效分析标签数,测序饱和度,细胞标签Q30,分子标签Q30,序列Q30。

(3) 比对结果统计

统计序列比对到基因组及基因组的基因间区,外显子,内含子的比例。

(4) 细胞数目及细胞过滤情况

带有标签的凝胶珠与细胞形成油包水的过程中,会存在油滴中不含细胞的情况,这时候油滴中只有细胞标签的序列及溶液中的核酸片段,而含有细胞的凝胶珠中序列数目一般空油滴的10倍以上,由此来确定哪些是真实的细胞,哪些是背景。

曲线图横轴是细胞标签,纵轴是对应的分子标签的数量,深蓝色代表细胞,灰色代表背景,

(5) 样本基本信息

(6) 细胞表达量分布tSNE图

每个细胞的分子标签对应的基因的序列数目就是转录本的表达量,采用tSNE降维算法, 对进行降维后可视化,从而将具有相似表达谱的细胞聚类在一起。对结果进行聚类从而分出不同类型的细胞。示意如下

(7) 差异分析

对不同类的基因进行差异分析,从而分成该类与其他类两种,然后进行差异分析。

(8) 饱和度分析

3.2.1.2 CellRanger定量单细胞原始表达矩阵与过滤后的表达矩阵

在结果目录,可以看到raw_feature_bc_matrix和filtered_gene_bc_matrices两个目录。分别包含有单细胞原始表达矩阵与过滤后的表达矩阵

3.3 单细胞转录组基因与细胞质量控制

3.3.1 单细胞转录组基因与细胞质量控制简介

质控的主要参数有2种:

(1) 每个细胞测到的基因数目与序列数目。

(2) 每个细胞检测到的线粒体与核糖体基因的比例,理论上细胞状态较好时线粒体基因与核糖体基因表达比例不会过高,因此过高线粒体比例与核糖体比例的细胞会被过滤。

图3.3.1 单细胞序列数,基因数,线粒体比例与核糖体比例质量控制图

过滤前(上图)与过滤后(下图)单细胞序列数,基因数,线粒体比例与核糖体比例比较。

3.3.2 细胞质量控制参数及关系

在过滤后我们会对过滤前后的细胞基因数,表达量,线粒体与核糖体比例进行比较,并利用点图对其关系进行描述。

图3.3.2 单细胞序列数,基因数关系点图

过滤前(左图)与过滤后(右图)单细胞序列数,基因数关系比较。

图3.3.3 单细胞线粒体比例与序列数,基因数关系点图

过滤前(左图)与过滤后(右图)单细胞线粒体比例与序列数,基因数关系比较。

图3.3.4 单细胞核糖体比例与序列数,基因数关系点图

过滤前(左图)与过滤后(右图)单细胞核糖体比例与序列数,基因数关系比较。

3.4 单细胞转录组降维与聚类

3.4.1 高度变化基因的鉴定与数据归一化

细胞为了保证自身的正常功能均会高度表达一些维持状态的基因,这些基因表达量很高但是并不能显示不同细胞的特点,因此需要鉴定出细胞间表达高度变化的基因,并在研究中基于这部分基因进行后续分析。我们计算每一个基因的均值和方差,并且获取其关系。从而鉴定出高度变化的基因。

在鉴定出高度变化的基因后我们将数据归一化到均值为0方差为1,从而消除基因表达的绝对量的影响。

3.4.2 降维分析

使用前面鉴定表达变化大的基因对归一化后的数据进行主成分分析,根据标准差的变化程度选择出最能体现“差异”的主成分进行后续的分析。

图3.4.1 主成分变异图

PCA降维后不同主成分对变异的解释程度图,左图的PC显示具有低p值的特征的强烈丰富(虚线上方的点图)。右图基于每个PC解释的方差百分比对主成分进行排序。

图3.4.2 主成分变异来源基因点图与热图

PCA降维后不同主成分变异来源基因图,上图为变异来源基因点图,下图为变异来源基因热图。

使用tSNE和UMAP的方法对降维后的数据进行可视化。并将聚集在一起的细胞定义为一个类。

图3.4.3 单细胞PCA,tSNE,UMAP降维聚类图

上图为PCA,中图为tSNE,下图为UMAP降维聚类图,横纵坐标代表了降维前两位主成分值,颜色代表聚类结果。

图3.4.4 输入基因UMAP表达分布图

感兴趣的基因在不同分组中的表达情况,颜色深浅代表表达量的高低。

3.5 单细胞转录组标记基因鉴定

3.5.1 类标记基因的鉴定

通过比较某一个类中的细胞与剩余其他类中细胞的表达谱的差异,我们可以找到找到各个类中显著区别于其他细胞的基因,并将其定义为该类的生物学标记基因。我们利用热图对这些基因进行可视化,并绘制基因在所有类中表达分布图与映射图来详细阐述标记基因的表达情况。

图3.5.1 单细胞标记基因在不同类中表达量琴图与细胞分布图(示例Cluster0)。

上图为标记基因在不同类中表达量琴图,下图为标记基因在细胞分布图(示例Cluster0)中的表达量。

3.6 单细胞转录组细胞类型鉴定

3.6.1 SingleR原理简介

SingleR是对单细胞转录组数据进行细胞类型自动注释的R包。它通过已知类型的细胞样本作为参考数据集,对数据集中与参考集相似的细胞进行标记注释。对每个细胞进行如下处理

(1) 计算每个细胞的基因表达与参考细胞数据集表达谱之间的相关性。

(2) 对每种细胞类型的得分进行计算并取固定分位数作为细胞类型得分。

(3) 对所有的标签重复此操作,然后将得分最高的标签作为此细胞的注释

SingleR提供了针对人和小鼠多种组织的细胞类型的参考数据集,可以根据实际需要选择使用。

3.6.2 细胞类型注释结果

在输出的结果表格中,包含了每个细胞注释的预测结果以及相关得分。并使用热图对细胞注释可靠性结果进行展示。

图3.6.1 细胞类型分布与注释热图

上图为tSNE细胞类型分布图,下图为细胞类型可靠性热图。

3.7 单细胞转录组标记基因富集分析

对于不同的细胞类我们会得到很多的标记基因,这无疑增加了后续分析的复杂度。因此,可以将这些差异基因按照不同的功能类别进行分类,即对差异基因进行富集分析。通过富集分析,可以找到不同条件下的差异基因与哪些生物学功能或通路显著性相关,从而揭示和理解生物学过程的基本分子机制。

我们采用clusterProfiler软件对差异基因集进行GO功能富集分析,KEGG通路富集分析等。富集分析基于超几何分布原理,其中差异基因集为差异显著分析所得差异基因并注释到GO或KEGG数据库的基因集,背景基因集为所有进差异显著分析的基因并注释到GO或KEGG数据库的基因集。富集分析结果是对每个差异比较组合的所有差异基因集进行富集。

图3.7.1 GO富集分析原理图

3.7.1 GO功能富集

3.7.1.1 GO富集列表

GO(Gene Ontology)是描述基因功能的综合性数据库(http://www.geneontology.org/),可分为生物过程(Biological Process)和细胞组成(Cellular Component)分子功能(Molecular Function)三个部分。GO 功能显著性富集分析给出与基因组背景相比,在差异表达基因中显著富集的 GO 功能 条目,从而给出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向 Gene Ontology 数据库的各个 term 映射,计算每个 term 的基因数目,然后找出与整个基因组背景相比,在差异表达基因中显著富集。GO富集分析方法为clusterProfiler的enricher函数,该方法使用超几何分布检验的方法获得显著富集的GO Term。统计被显著富集的各个GOterm中的基因数,以柱状图的形式展示。GO富集以padj小于0.05为显著富集。

3.7.1.2 GO富集柱状图

根据GO富集分析结果,在GO的三个分类中(cellular_component:细胞组分;biological_process:生物学过程;molecular_function:分子功能)分别统计每个分类上调基因、下调基因及所有基因的数目,以柱状图展示。该部分结果图片见文件夹:customer_files/08EnrichAnalysis/08.01Goenrichment。

图3.7.2 差异基因GO富集柱状图

注:横坐标为GO三个大类,纵坐标为注释到该类别的基因数目。

3.7.2 KEGG通路富集

3.7.2.1 KEGG通路富集列表

在生物体内,不同基因相互协调行使其生物学功能,通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。

KEGG(Kyoto Encyclopedia of Genes and Genomes)是有关Pathway的主要公共数据库,其中整合了基因组化学和系统功能信息等众多内容。信息分析时,我们将KEGG数据库按照动物,植物,真菌等进行分类,依据研究物种的种属选择相应的类别进行分析。Pathway显著性富集分析以KEGG Pathway为单位,应用超几何检验,找出差异基因相对于所有有注释的基因显著富集的pathway,KEGG通路富集同样以padj小于0.05作为显著性富集的阈值。

3.7.2.2 KEGG通路富集气泡图

差异基因KEGG富集气泡图是KEGG富集分析结果的图形化展示方式。在此图中,KEGG富集程度通过GeneRatio、p.adjust和富集到此通路上的基因个数来衡量。其中GeneRatio指差异表达的基因中位于该pathway条目的基因数目与所有有注释基因中位于该pathway条目的基因总数的比值。p.adjust是做过多重假设检验校正之后的Pvalue,p.adjust的取值范围为[0,1],越接近于零,表示富集越显著。我们挑选了富集最显著的20条pathway条目在该图中进行展示,若富集的pathway条目不足20条,则全部展示。

图3.7.3 KEGG通路富集柱状图与气泡图

注:KEGG pathway富集气泡图:纵轴表示pathway名称,横轴表示pathway对应的GeneRatio,p.adjust的大小用点的颜色来表示,p.adjust越小则颜色越接近红色,每个pathway下包含的差异基因的多少用点的大小来表示。

3.8 单细胞转录组差异分析

除了可以比较某一类细胞和其他类细胞之间的差别来寻找标记基因,也可以在类间进行比较获得感兴趣的类之间具有显著差异的基因。

3.9 单细胞转录组差异基因富集分析

3.9.1 差异基因的富集分析结果同上

4 部分分析软件列表及版本

5 参考文献

1. Zheng, Grace X.Y., Terry, Jessica M., [...] Bielas, Jason H. (2017). Massively parallel digital transcriptional profiling of single cells. Nature Communications. 8: 1-12, doi:10.1038/ncomms14049.

2. Yuhan Hao, Stephanie Hao, [...] Rahul Satija (2021). Integrated analysis of multimodal single-cell data,Cell,Volume 184, Issue 13,Pages 3573-3587.e29,https://doi.org/10.1016/j.cell.2021.04.048.

3. Aran D, Looney AP, [...] Bhattacharya M (2019).Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat. Immunol, 20,163-172. doi:10.1038/s41590-018-0276-y.

4. Yu G, Wang L, Han Y and He Q*. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.

5. Cao, J. et. al. (2019) The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496–502. https://doi.org/10.1038/s41586-019-0969-x

6. Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

7. Philip Ewels, Måns Magnusson, Sverker Lundin and Max Käller. (2016) MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics. doi: 10.1093/bioinformatics/btw354


6. 联系我们

地址:广东省深圳市南山区留仙大道1201号大学城创客小镇

邮编:518055