一种基于单细胞转录组聚类结果注释细胞身份的方法

摘要

本发明公开了一种基于单细胞转录组聚类结果注释细胞身份的方法，包括步骤：S1、提供存活细胞，利用10X genomics建库并测序，获得转录组测序数据；S2、对步骤S1所得测序数据进行过滤，然后用软件cellranger进行初始分析，输出原始数据；S3、分析步骤S2输出的原始数据；S4、细胞身份注释，包括步骤：S41，将Cell Marker数据库的标签基因按细胞类型进行分类；S42，将FindAllMarkers函数筛选出来的标签基因按P值排序；S43，将Cell Marker数据库中每类细胞的标签基因与FindAllMarkers函数筛选出来的细胞亚群标签基因取交集，根据交集基因，进行打分；S44，根据基因交集打分大小排序，将分值最大者对应的细胞类型注释为当前亚群的细胞身份。采用本发明的方法能准确快速地对细胞分类和身份注释。

说明书

技术领域

本发明涉及细胞分类技术领域，尤其是一种基于单细胞转录组聚类结果注释细胞身份的方法。

背景技术

10X genomics测序平台能够同时检测成千上万个单细胞的转录组，因此，根据它们的基因表达谱可以将它们聚成不同的类别。当前遇到的最大困难是如何注释这些不同类别细胞的身份。因为细胞的种类成百上千，我们无法完全掌握每种细胞表达的标签基因。10X genomics测序的主要目的又是要发现新的细胞类型。这就给我们造成一个极大的挑战。

10x Genomics测序平台(工作流程参见图1)利用油包水的微反应体系，每个油滴包含单个细胞。不同的油滴添加不同的人工序列标签，使得油滴内细胞的所有核酸标记相同的人工序列，而不同的细胞标记不同的人工序列，混在一起建库测序以后，很容易根据人工序列，鉴别来自不同细胞的转录组，从而实现数千甚至数万个单细胞群体分析，发现细胞的个体差异，鉴别新的细胞类型。例如，实体肿瘤组织中，既包含炎症细胞，也包含肿瘤细胞；而肿瘤细胞中，既包含普通肿瘤细胞，也包含肿瘤干细胞。因此，通过10x Genomics测序平台对实体肿瘤组织进行高通量单细胞测序，就有可能鉴别出其中的肿瘤干细胞。

但是，现实的困难是：需要发现的细胞，其基因表达特征往往是未知的。而且测序的群体细胞中往往掺杂多种细胞类型，其基因表达特征也是未知的。为了解决这个问题，必须事先尽可能了解清楚，测序的群体细胞中包含哪些细胞，其基因表达特征是什么。为此，必须查阅大量的文献，而这是一件艰巨的挑战，至少需要几个月的时间。

发明内容

基于上述问题，本发明的目的在于克服上述现有技术的不足之处而提供一种可基于单细胞转录组测序数据准确、快速地对细胞进行分类和身份注释的方法。

为实现上述目的，本发明采取的技术方案为：

一种基于单细胞转录组聚类结果注释细胞身份的方法，包括如下步骤：

S1、提供存活细胞，利用10X genomics建库并测序，获得转录组测序数据；

S2、对步骤S1所得测序数据进行过滤，然后用软件cellranger进行初始分析，输出原始数据；

S3、分析步骤S2输出的原始数据，包括：

1)选择至少在三个细胞表达的基因，基因表达数量在200～8000之间以及线粒体基因表达占比不超过0.4的细胞；

2)采用LogNormalize方法的全局数据标准化；

3)用ScaleData函数去除由于细胞总表达量不同和线粒体基因表达占比不同引起的基因表达波动；

4)PCA主成分分析降低数据的维度；

5)用FindClusters函数方法进行聚类分析，结果用tSNE呈现出来；

6)用FindAllMarkers函数筛选细胞亚群的标签基因(即表达上调的差异表达基因)；

S4、细胞身份注释，包括步骤：

S41，将Cell Marker数据库的标签基因按细胞类型进行分类；

S42，将FindAllMarkers函数筛选出来的标签基因按P值排序；

S43，将Cell Marker数据库中每类细胞的标签基因与FindAllMarkers函数筛选出来的细胞亚群标签基因取交集，根据交集基因，进行打分；

S44，根据基因交集打分大小排序，将分值最大者对应的细胞类型注释为当前亚群的细胞身份。

优选地，所述步骤S2中过滤包括去掉接头和去掉reads首尾碱基质量小于3的碱基。

优选地，所述步骤S2中初始分析选择软件cellranger的默认参数。

优选地，所述步骤S3中，用R语言Seurat软件包读取cellranger软件输出的原始数据。

优选地，所述步骤S41中，计数标签基因的出现次数，将标签基因分成三类并赋予不同的加权分值，出现一次为一类，加权分值为1(表示初步可靠)；出现两次为二类，加权分值为1.5(表示基本可靠)；出现三次及三次以上为三类，加权分值为2(表示绝对可靠)。需要说明的是，标签基因在文献中重复出现，表明该基因被广泛证实，因此，重复出现的基因比单次出现的基因更可靠。

优选地，所述步骤S42中，以P值最小的基因代表当前类别细胞的特征。需要说明的是，P值即每类别细胞与所有其它类别细胞的Wilcoxon秩和检验的统计学分值，P值越小，说明该基因在当前类别细胞中的表达一致性更好，且与其它类别细胞的表达差异更显著，因此，P值越小的基因越能代表当前类别细胞的特征。

优选地，所述步骤S43中，打分公式为：其中n为交集基因数目，N当前类别细胞标签基因的数目，I_i为第j个交集基因在所有标签基因中的排序，W_j为第j个交集基因在Cell>

综上所述，本发明的有益效果为：

基于单细胞转录组测序数据，采用本发明的方法能准确快速地对细胞进行分类和身份注释；现有技术根据各个类别细胞的标签基因查阅文献，从而判定细胞的名称和类别，这个过程耗时耗力，至少需要几个月的时间；而采用本发明的注释方法只需几分钟的时间就可完成细胞分类和名称注释。

附图说明

图1为10x Genomics测序平台工作流程示意图；

图2为鼻咽癌组织10x genomics单细胞转录组测序聚类结果图；

图3为鼻咽癌组织10x genomics单细胞转录组测序聚类的注释结果图。

具体实施方式

最近哈尔滨医科大学发布了CellMarker数据库(http://biocc.hrbmu.edu.cn/CellMarker/)。该数据库通过人工检索超过10万篇文献后，从4124篇文章中摘取细胞标志物信息、组织类型、细胞类型、肿瘤信息和来源等，收集了人和鼠的各种已知细胞的标签基因。结果得到人各种已知细胞(476种)的标签基因12605个和鼠的各种已知细胞(389种)的标签基因9148个。基于该数据库，可以很好地注释群体单细胞测序的聚类结果。由此，本发明提出一种注释方法，可以注释10X genomics单细胞转录组的聚类结果。

为更好的说明本发明的目的、技术方案和优点，下面将结合附图和具体实施例对本发明作进一步说明。

实施例1

本发明中基于单细胞转录组聚类结果注释细胞身份的方法的一种实施例，包括如下步骤：

1)建库测序

从广州市第十二人民医院耳鼻喉科获得一例新鲜鼻咽癌组织(男性、42岁、病理诊断为鼻咽低分化鳞癌)。癌组织离体马上放入组织培养液中，冰上保存，并在两个小时之内制成细胞悬液。细胞计数仪检测活细胞比例超过90％之后，送给广州序科码生物技术有限责任公司进行10X genomics建库。建库完，送给北京安诺优达基因科技有限公司进行测序。测序平台为Illumina Hiseq X，测序深度为包lane测序。

2)数据分析

首先将测序的原始数据进行过滤，包括：去掉接头，去掉reads首尾碱基质量小于3的碱基；在处理过程中以4个碱基为一组，若组的平均碱基质量小于15，则去掉该组碱基序列。然后，用10x genomics官网的分析软件cellranger进行初始分析(选择默认参数)。

然后，用R语言Seurat软件包读取cellranger软件包输出的原始数据，进行以下分析：(1)选择至少在三个细胞表达的基因，基因表达数量在200～8000之间以及线粒体基因表达占比不超过0.4的细胞；(2)采用LogNormalize方法的全局数据标准化；(3)用ScaleData函数去除由于细胞总表达量不同和线粒体基因表达占比不同引起的基因表达波动；(4)PCA主成分分析降低数据的维度，即把表达线性相关的基因合并成一个新的变量，从而大大减少变量的数量；(5)用FindClusters函数方法进行聚类分析(参数设置为：reduction.type＝"pca",dims.use＝1:20,resolution＝0.8；即用主成分1:20进行聚类分析，聚类精度为0.8)，结果用tSNE呈现出来；(6)用FindAllMarkers函数筛选细胞亚群的标签基因(即表达上调的差异表达基因)。

3)细胞身份注释

首先，将Cell Marker数据库的标签基因按细胞类型进行分类。计数标签基因的出现次数。标签基因在文献中重复出现，表明该基因被广泛证实。因此，重复出现的基因比单次出现的基因更可靠。根据出现次数，把标签基因分成三类并赋予不同的加权分值：出现一次为一类(初步可靠，加权分值为1)，两次为二类(基本可靠，加权分值为1.5)，三次及三次以上为三类(绝对可靠，加权分值为2)。

其次，将FindAllMarkers函数筛选出来的标签基因按P值(即每类别细胞与所有其它类别细胞的Wilcoxon秩和检验的统计学分值)排序。P值越小，说明该基因在当前类别细胞中的表达一致性更好，且与其它类别细胞的表达差异更显著。因此，P值越小的基因越能代表当前类别细胞的特征。

第三，将Cell Marker数据库中每类细胞的标签基因与FindAllMarkers函数筛选出来的细胞亚群标签基因取交集。然后根据交集基因，进行打分。打分公式为：其中n为交集基因数目，N当前类别细胞标签基因的数目，I为第j个交集基因在所有标签基因中的排序，W_j为第j个交集基因在Cell>

最后，根据基因交集打分大小排序，将分值最大者对应的细胞类型注释为当前亚群的细胞身份，结果如图2和3所示。

结果分析：

该鼻咽癌组织样品制成细胞悬液，经10x genomics测序平台建库测序，再经cellranger和seurat软件分析后，得到2727个细胞，共表达19166个基因。全基因组表达谱聚类分析得到13个细胞亚群(图2)。FindAllMarkers函数筛选细胞亚群的标签基因数目如表1所示。用这些标签基因分别与Cell Marker数据库各个细胞类型的标签基因取交集，再用打分系统进行打分排序后，注释结果如图3所示。

由于人的鼻咽部具有丰富的淋巴组织，且鼻咽癌是慢性鼻咽炎经过几十年演化的结果，所以鼻咽癌组织里富含炎症细胞和免疫细胞。图3的注释结果证实了这一点。该样品的鼻咽癌组织富含T细胞、B细胞、浆细胞、中性粒细胞、自然杀伤细胞、单核巨噬细胞等炎症细胞和免疫细胞。同时，其中上皮细胞(“Epithelial cell”)则代表鼻咽癌细胞亚群，而其中的肝祖细胞(“Liver progenitor cell”)则代表鼻咽癌肿瘤干细胞亚群。这些都证明本发明注释结果的可靠性。有趣的是，发明人从注释结果中发现了一群新的细胞“Neuralprogenitor cell”。这群细胞可能是免疫细胞的祖细胞，也可能是肿瘤细胞的干细胞，值得深入探讨与研究。

表1不同细胞亚群的标签基因数目

最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。