CN102952854B - 单细胞分类和筛选方法及其装置 - Google Patents

Abstract

本发明提出单细胞分类和筛选的方法，以及用于所述方法的装置。将样本测序得到的reads与参考基因组进行比对，并将比对结果进行数据过滤；根据过滤后的数据确定每个单细胞样本的一致基因型，并将所有单细胞样本的一致基因型保存为SNP数据集；从已保存的SNP数据集提取与参考基因组SNP数据集位置对应的位点的基因型文件；挑选细胞突变SNP位点，并根据细胞突变SNP位点的基因型文件，对细胞进行分类及功能基因筛选。本发明避免对细胞进行标记的操作，解决了传统的单细胞分类方法中某些细胞亚群无对应的特异性标记物而无法分类的问题；另一方面，全面完整地分析单细胞基因组的遗传变异信息，大大提高了细胞亚群分类的准确性。

Description

单细胞分类和筛选方法及其装置

技术领域

本发明涉及生物信息学，尤其涉及单细胞分类和筛选方法以及用于所述方法的装置。

背景技术

不同个体之间，个体的不同组织之间，甚至同一组织的不同部位在基因表达、拷贝数变异、表观遗传等方面都存在显著差异。细胞之间也存在异质性，即使是体外培养遗传背景完全相同的细胞群体。对于干细胞或前体细胞，因为任何状态改变都是可遗传的，细胞异质性尤为明显。为了更好地研究细胞生物学，揭示细胞异质性的规律，非常需要开发应用于单个细胞研究的技术方法，因此有学者提出“单细胞分析(SCA)”概念，从“组学(Omics)”角度进行阐述。单细胞分类和筛选为单细胞分析提供了重要基础。

单细胞分类可以有效应用于各种干细胞分化过程的研究中，如肿瘤干细胞、胚胎干细胞的定向分化、造血干细胞的研究中，需要筛选不同分化阶段的干细胞，进行各种干细胞的检测。在耐药性研究中，需要对给药不同时期的细胞进行精确分类，从而进一步分析该细胞亚群的耐药性和耐药基因，例如可进行癌症病人的多药耐药性及多药耐药基因与药物滥用、药物耐受、药物依赖的关系的研究。同样地，在药物靶点基因的筛选中，由于药物与细胞，特别是敏感细胞相互作用，将引起细胞外部形态及内部正常代谢过程的一系列变化，因此筛选出敏感细胞是关键的第一步，为后期精确定位药物靶点基因提供重要基础。单细胞分类和筛选应用于建立药效筛选模型，为药物设计、靶点的选择和用药方案的确定提供理论依据，同时使药物筛选有了更高的特异性。

目前，常用的筛选单细胞方法多为物理机械、化学或生物的方法，如流式细胞仪、磁性细胞分选仪等方法。一方面，这些技术采用表面活性剂、荧光染料、抗原抗体，细胞毒性大，只能对特异标记的或非特异标记的单细胞悬液进行分选，前期样本制备过程繁琐，且目前对众多荧光探针、单抗(包括细胞表面CD分子)的特异性争论较多，许多细胞亚群并无对应的特异性标记物/特异性抗原；另一方面，这些技术采用生物学、免疫学、化学方法，通过表型测定(包括细胞大小、细胞粒度、细胞表面积、核浆比例等)，进行统计学分析，对于亚群分类、筛选和检测的灵敏度低，缺乏有效的准确性评估。

发明内容

在本发明中，除另有说明，否则本文中使用的科学和技术术语具有本领域技术人员所通常理解的含义。同时，为了更好地理解本发明，下面提供相关术语的定义和解释。

术语“基因型的可能性文件”，是指利用SNP检测软件，设置先验概率参数利用贝叶斯公式计算出的样本目标区域可能的基因型的后验概率的数值集合；当利用的SNP检测软件是SOAPsnp时，生成的“基因型的可能性文件”即为CNS文件。

如本文使用的，“基因型文件”是指选择上述“基因型的可能性文件”中概率最大的基因型作为每个细胞的一致基因型后，根据参考基因组SNP数据集位置信息，提取每个细胞基因型的相应位点，获得的群体SNP在各细胞相应位点的基因型集合。

鉴于现有单细胞的分类和筛选方法存在的问题，本发明提出单细胞分类和筛选方法，以及所述方法的装置。

本发明提出单细胞分类方法，包括以下步骤：

将每个单细胞样本经测序得到的reads(读段)结果与参考基因组序列进行比对，并将比对结果进行数据过滤；

根据过滤后的数据确定每个单细胞样本的一致基因型(genotype)，并将所有单细胞样本的一致基因型保存为SNP数据集；

从已保存的SNP数据集提取与参考基因组SNP数据集位置对应的位点的基因型文件；

挑选细胞突变SNP位点，并根据细胞突变SNP位点的基因型文件，对细胞进行分类。

本发明还提出单细胞分类装置，包括：

数据过滤模块，将每个单细胞样本经测序得到的reads与参考基因组序列进行比对，并将比对结果进行数据过滤；

基因型确定模块，根据过滤后的数据确定每个单细胞样本的一致基因型，并将所有单细胞样本的一致基因型保存为SNP数据集；

基因型文件提取模块，从已保存的SNP数据集提取与参考基因组SNP数据集位置对应的位点的基因型文件；

分类模块，挑选细胞突变SNP位点，根据细胞群体突变SNP的基因型文件，对细胞进行分类。

本发明还提出单细胞筛选方法，包括以下步骤：

得到已预测基因组中基因的起止位置；

根据细胞分类结果获得已分类群体，计算每类群体中每个基因所有SNP位点的统计量，并累加统计量；

对获得的统计量作差异检验，获得检验值；

将已预测基因按统计量或检验值进行排序，筛选出统计量或检验值最高的基因。

本发明还提出单细胞筛选装置，包括：

获取单元，得到已预测基因组中基因的起止位置；

计算单元，根据细胞分类结果获得已分类群体，计算每类群体中每个基因所有SNP位点的统计量，并累加统计量；对获得的统计量作差异检验，获得检验值；

排序单元，耦合于获取单元和计算单元，将已预测基因按统计量或检验值进行排序，筛选出统计量或检验值最高的基因。

本发明采用新一代测序技术(NGS)，通过生物信息学方法，对单细胞基因组进行分析和研究，收集细胞亚群(或微粒)进行深入的后续研究。一方面，避免对细胞进行标记的操作，有效解决了传统的单细胞分类方法中某些细胞亚群无对应的特异性标记物而无法进行分类的问题；另一方面，全面完整地分析单细胞基因组的遗传变异信息，并设置多个对照样本，大大提高了细胞亚群分类的准确性。

本发明还提出单细胞筛选方法，能够筛选出细胞亚群(或微粒)，提高了细胞筛选的准确性。

附图说明

图1为现有技术的重复片段(Duplication Reads)示意图；

图2为现有技术的唯一比对上参考基因组的片段(Unique mappedreads)的示意图；

图3为本发明单细胞分类和筛选的方法流程图；

图4为本发明肾癌外显子组测序的N-J关系树；

图5为本发明肾癌外显子组测序的最大似然关系树；

图6为本发明肾癌外显子组测序PCA结果图，横坐标表示第一主成分向量，纵坐标表示第二主成分向量；

图7为本发明肾癌外显子组测序PCA结果图，横坐标表示第一主成分向量，纵坐标表示第三主成分向量；

图8为本发明肾癌外显子组测序PCA结果图，横坐标表示第一主成分向量，纵坐标表示第四主成分向量；

图9为本发明肾癌外显子组测序Structure结果图，其中“系列1”表示癌细胞群体，“系列2”表示癌旁细胞群体；

图10为本发明53个癌细胞与8个正常细胞分类关系示意图；

图11为本发明癌细胞与正常细胞聚类示意图，横坐标表示第一主成分向量，纵坐标表示第二主成分向量；

图12为本发明单细胞分类装置示意图；

图13为本发明单细胞分类装置中筛选模块示意图。

具体实施方式

本发明采用新一代测序技术(NGS)，通过生物信息学方法，对单细胞基因组进行分析和研究，筛选和收集细胞亚群(或微粒)进行深入的后续研究。以更高效、方便地应用于临床诊断和治疗(如产前诊断、胚胎植入前遗传诊断、个体化治疗、多点图谱制作、精子和卵子的分型、遗传病诊断、肿瘤(如淋巴瘤、白血病)分型等)、医学研究(如自闭症、神经系统疾病和自体免疫性疾病的研究、基因组变异率研究、干细胞研究、耐药性研究、药物靶点基因的筛选、转录组检测、细胞模型研究、种群鉴定等)、考古学研究、法医学检测中。

本发明中涉及的单细胞样本包括核酸(基因组DNA或RNA，如非编码RNA、mRNA)；单细胞来自生物体，采用常规方法制备。特别的，DNA或RNA可由细菌、原生动物、真菌、病毒及高等生物/高等动植物，如哺乳动物，特别是人类的单细胞提取或扩增得到。单细胞可经过体外培养，或者从临床样本(包括血浆、血清、脊液、骨髓、淋巴液、腹水、胸腔积液、口腔液体、皮肤组织，呼吸道、消化道、生殖道、泌尿道，眼泪、唾液、血细胞、干细胞、肿瘤)中直接分离得到，胎儿细胞可来自胚胎(如一个或几个拟胚/胚胎)或母体血液，可来自活体或者死亡生物体。样本包括单细胞悬液、石蜡包埋组织切片、穿刺活检组织。

样本可反映细胞特定状态，如细胞增殖、细胞分化、细胞凋亡/死亡、疾病状态、外界刺激状态、发育阶段。

特别地，单细胞样本来自于哺乳动物，包括植入前胚胎、干细胞、疑似癌细胞、病源生物，甚至来自犯罪现场。例如分析人类卵裂球细胞(八细胞阶段胚胎及以后)，可检测胎儿基因组是否出现遗传缺陷。

下面结合图3，对本发明单细胞分类方法的具体实现过程进行详细说明。其中，图3示出了从(7)开始的流程。

(1)分离单细胞：采用物理机械、化学、生物的方法，如微流控、口吸分离、梯度稀释、低熔点琼脂糖固定等方法，分离得到包含完整基因组的单个细胞。

(2)裂解细胞：对分离得到的单个细胞，采用去污剂法、煮沸法、碱变性法、溶菌酶法、有机溶剂法等方法，裂解细胞核，得到完整的细胞基因组DNA。

(3)单细胞全基因组扩增(WGA)：

目前全基因组扩增有2种策略：即基于PCR的扩增，如DOP-PCR、PEP-PCR、T-PCR，和线性DNA扩增，如OmniPlex WGA、多重置换扩增(MDA)。进行单细胞全基因组扩增，以达到新一代测序技术所要求的DNA起始量。

(4)全基因组扩增产物定量：采用凝胶电泳检测、Agilent 2100Bioanalyzer检测、Quant-iT^TM dsDNA BR Assay Kit检测等方法对单细胞全基因组的扩增产物进行定量，结果显示为无降解、符合新一代测序技术所要求的DNA起始量的样本才可以继续DNA文库构建以及上机测序。

(5)全基因组扩增产物检测：采用STR检测、Housekeeping Gene检测等方法，对单细胞全基因组的扩增产物进行检测，结果显示为扩增产物在相应物种的染色体上均匀分布的样本才可以继续DNA文库构建以及上机测序。

(6)DNA文库构建及上机测序：采用常规的全基因组DNA文库构建或外显子组(Exome)序列捕获技术进行DNA文库构建，质检合格后，采用新一代测序技术进行单细胞基因组测序，如Illumina HiSeq 2000测序系统、Illumina Genome Analyzer II测序系统、AB SOLiD^TM 4.0测序系统、Roche GS FLX Titanium System等。

(7)Reads定位

每个单细胞样本经测序得到的Reads结果用短序列比对软件(如SOAPaligner、BWA、Bowtie)与参考基因组序列(如人类基因组HG18、HG19)进行比对，并根据具体数据设置最优参数，以对Reads进行准确定位。

(8)基本数据统计

根据上述比对结果，计算出每个样本(单细胞或组织)相对于参考基因组序列的测序深度和覆盖度等结果。

测序深度(sequencing depth)指一个基因组被测到的平均深度，计算方法为测序的总碱基数除以基因组的大小。

覆盖度(sequencing coverage)指一个基因组中被测到的部分所占的大概比例，计算方法为基因组被覆盖到的位点除以基因组的有效长度。

测序深度和覆盖度用来评价用于生物信息分析的数据量是否足够，测序是否均一。

(9)数据过滤

当存在多对重复双末端reads时，在多对reads的序列完全一致时随机选择一对reads，将所述多对reads中其他重复双末端reads去除；和/或去除未唯一比对上参考基因组序列的reads。

根据数据特点挑选出每个DNA文库中重复双末端(pair-end)reads，比如，由于PCR扩增次数过多引起的重复双末端reads，当然，不限于PCR扩增，也可以是其他扩增方式导致具有多对重复的双末端reads。

当存在多对重复双末端reads，且所述重复双末端reads的序列完全一致时，从中随机选择一对，将其他重复双末端reads去除。

如图1所示，三对reads A、B、C序列完全一致并且比对到基因组上的起止位置也完全一致，起止位置完全一致即具有重复双末端。在这种情况下只随机保留其中一对reads，将其他重复的reads去除。

为了确认数据的准确性，还可以去除未唯一比对上参考基因组序列的reads。以对人类基因组进行外显子组测序为例，当然，并不限于此，比如对其他哺乳动物等进行测序，测序方式也不限于外显子组测序，比如全基因组测序等方式。考虑到人类外显子区域在基因组上并不可能有多个拷贝，即不可能是重复序列，所以外显子组测序得到的reads应大都唯一比对上人类参考基因组。为排除错误比对造成的影响，只挑选出唯一比对上参考基因组的reads进行分析(即hit数为1的reads)，这样在极大程度上减少了错误比对造成的影响。

如图2所示，Reads D比对上参考基因组的多个位置，而Reads E只比对上唯一一个位置，由于外显子组在基因组上并不是重复区域，因此直接去除Reads D。

(10)个体基因型(genotype)的判定

我们充分考虑并利用参考基因组上的已有信息，用基因型判定软件(如SOAPsnp、SAMtools/Pileup/Mpileup)，判定每个细胞样本目标区域的可能的基因型，以得到每个细胞样本的基因型的可能性文件。

本发明中所测的是外显子区域的数据，在该实施例中，目标区域是外显子所在的区域。一般会指明所需测序和生物信息分析的特定区域如：

chr1  20038   20358

chr1  58832   59992

chr1  357410  358570

……

(11)SNP数据集

由于每个细胞基因组存在一些低深度的区域，本发明综合所有细胞基因型的可能性文件，利用最大似然法(maximum likelihood approach)整合所有细胞的数据，产生包含所有样本每个位点的伪基因组(Pseudo-genome)。选择概率最大的基因型作为每个细胞的一致基因型，通过基因型和测序深度等信息检测出高质量的SNPs。得到样本的一致性序列后，结果保存为群体SNPs格式的SNP数据集。

(12)群体SNP的基因型

根据参考基因组SNP数据集位置信息，从基因型的可能性文件中提取每个细胞相应位点的基因型，得到群体SNP在各细胞相应位点的基因型文件。位点指SNP所在的位置(position)。

(13)挑选与细胞突变相关的SNP位点

本发明主要在于寻找各细胞间的差异位点，因此须挑选出与细胞突变有关的位点。

首先，去除非群体间的SNP位点。

非群体间的SNP位点的定义：所有个体碱基型都一样，且相对参考序列都是SNP，那么该位点是非群体间的SNP位点。

例如，参考序列是A，所有个体在该位置都是杂合的碱基型R，该位点即是非群体间的SNP位点。如

chr1 319660  R R R R R R R R R R R R R R R

其次，还可以去除存在杂合丢失的位点。由于单细胞在WGA扩增时，存在只扩增出一对染色体中的一条的情况，导致等位基因丢失(AlleleDropout)，所以检测出来的每个细胞在某些位点存在杂合丢失的现象。排除这一类位点的干扰。

最后，去除已经公布的SNP位点，比如，以人类为例，去除正常人的SNP位点，即去除人类基因组HG18的dbSNP、炎黄一号的SNP和千人基因组的SNP。

上述的三个操作顺序不分先后，在执行完这三个操作之后，得到的SNP位点是细胞突变SNP位点。

(14)群体结构分析

根据细胞群体突变SNP位点的基因型文件，分别运用生物信息群体分析中常用的方法，如邻接N-J方法构树、MEGA软件、主成分分析(Principal Components Analysis，PCA)、群体结构(Group Structure)等对细胞进行分类。在进行细胞分类时，可以采用以上方法中的至少一种。作为本发明的一个实施例，采用以上全部方法，当各种方法的分类结果一致时，确认为最终的细胞分类结果。

14-1、邻接N-J方法构树

由于不同类别的细胞受到选择的程度不同，即单碱基的突变率不同，因此进化中的类别差异也在SNP上表现出来。两个细胞的差异程度可由SNP数据计算得到。经理论推算，得到如下公式：

${\mathrm{Dis}}_{\mathrm{ij}}=\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{diff}}_{\mathrm{ij}}^k$

上式中Dis_ij表示细胞i和细胞j差异距离，其中n是SNP的总数目，而表示细胞i和细胞j在位点k的差异程度，其中定义

由于人类基因组是二倍体，A表示等位基因中的两个位点都是A，M是杂合位点，即是A和C组合。基于上述步骤(13)获得的细胞群体突变SNP位点的基因分型文件，统计所有单细胞样本两两对比的差别，得到一个目标区域两两对比的差异矩阵，将上述差异矩阵应用到Fneighbor程序(http://emboss.bioinformatics.nl/cgi-bin/emboss/help/fneighbor)，该程序通过邻接(N-J)的方法可得到其进化树。

14-2、MEGA软件

利用MEGA软件(http://www.megasoftware.net)，将每个细胞所有SNP位点的基因型文件组成序列，一个细胞对应一个序列，作为MEGA的输入文件，MEGA根据各细胞序列上的差异，其中该软件有三种方法(Maximum likelihood、Least Squares和Maximum parsimony)构成关系树。

14-3、PCA主成分分析

在统计学中，主成分分析(Principal Components Analysis，PCA)是一种简化数据集的技术，是一个线性变换。这个变换把数据变换到一个新的坐标系统中，使得任何数据投影的第一大变量数在第一个坐标(称为第一主成分)上、第二大变量数在第二坐标(第二主成分)上，依次类推。主成分分析经常用于减少数据集的维数，同时保留对数据集贡献最大的特征变量。通过保留低阶主成分，忽略高阶主成分来实现的。这是由于低阶成分往往能够保留数据集中最重要的方面。

根据参考文献A tutorial on Principal Components Analysis.Lindsay ISmith，2002-02和真实的SNP数据特点，首先将SNP数据转换成数字矩阵(与参考序列一致为0，相反为2，简并碱基为1)并作均一化。然后通过上述介绍的方法构建线性向量方程。

$y_i=a_{i0}+a_{i1}{x}_i+a_{i2}{x}_i^2+...+a_{i20}{x}_i^{21}$

其中i从1到21表示第i个样本。

应用R语言软件包强大的解方程能力，解得矩阵a，根据各细胞的数据特点提取前四个主成分向量，以向量作为坐标轴展示各细胞聚类情况。

14-4、Structure群体结构

Structure软件(http://pritch.bsd.uchicago.edu/software/structure2_1.html)基于SNP位点的基因分型数据，推断是否存在不同群体并判断每个个体所归属的群体。根据软件说明，将群体SNP的基因型文件转换格式，作为Structure输入文件并在混合模型中采用高达5万次模拟，假设多个群体存在时，计算每个细胞归属各类群体的概率。

经过以上的方法流程，实现了对单细胞的分类。在分类的基础上，还可以进一步筛选单细胞，其流程如下：

(15)群体结构分析结果

根据上述群体结构分析结果，实现了对单细胞的分类，提取每个细胞样本信息，剔除存在异议的细胞，比如分类不明确或明显离群样本。

(16)筛选相关基因

根据细胞群体的SNP，通过一系列统计和检验进行这些群体在基因组的比较，找出有明显差异的区域或者基因，即可筛选到相关系数较高的基因。

以人类基因组为例，具体做法如下：

从人类基因组数据库中下载HG18对应的注释文件，得到目前已预测的人类基因组中30000多个基因的起止位置。

根据分类结果获得已分类群体，计算每类群体中每个基因所有SNP位点的统计量，并累加统计量。这里所说的每个基因是指已预测的基因组中的基因。

其中主要采用的计算统计量π的公式如下，π是衡量一个群体多态性高低的指标，a，b指某群体中两种碱基的样本数，公式可以为：

$\mathrm{\π}=\frac{a*b}{C_{a+b}^2}$

还可以对获得的统计量作差异检验，获得检验值。所采用的检验值至少为如下之一：Lod、Fst、Pbs。作为本发明的一个实施例，可以采用以上三种检验值，当以上三种检验值一致时，作为最终的检验值结果。

将这30000多个基因按统计量和/或检验值进行排序，挑选出统计量和/或检验值最高的基因。也就是说，可以根据统计量进行排序、可以根据检验值进行排序、也可以根据统计量和检验值进行排序。作为本发明的一个实施例，可以采用最后一种方法，当根据统计量得到的排序结果与根据检验值得到的排序结果一致时，作为最终筛选出的基因。

(17)基因功能分析

查看筛选出的基因的功能，分别进行功能分析。判断这些基因是否在某些通路中受到影响，从而与某些疾病的发生机制相关。

如图12所示，为本发明单细胞分类装置示意图。该装置包括：

数据过滤模块，将每个单细胞样本经测序得到的reads结果与参考基因组序列进行比对，并将比对结果进行数据过滤；

基因型确定模块，与数据过滤模块耦合，根据过滤后的数据确定每个单细胞样本的一致基因型，并将所有单细胞样本的一致基因型保存为SNP数据集；

基因型文件提取模块，与基因型确定模块耦合，从已保存的SNP数据集提取与参考基因组SNP数据集位置对应的位点的基因型文件；

分类模块，与基因型文件提取模块耦合，挑选细胞突变SNP位点，根据细胞群体突变SNP的基因型文件，对细胞进行分类，采用的分类方法至少包括如下之一：邻接N-J方法构树、MEGA软件、PCA主成分分析、Structure群体结构。

在另一实施例中，还如图12所示，单细胞分类装置还包括：

筛选模块，与分类模块耦合，用于得到已预测基因组中基因的起止位置；根据分类结果获得已分类群体，计算每类群体中每个基因所有SNP位点的统计量，并累加统计量；对获得的统计量作差异检验，获得检验值；将已预测基因按统计量或检验值进行排序，筛选出统计量或检验值最高的基因。

筛选模块可以进一步包括以下单元，如图13所示，其中：

获取单元，得到已预测基因组中基因的起止位置；

计算单元，根据细胞分类结果获得已分类群体，计算每类群体中每个基因所有SNP位点的统计量，并累加统计量；对获得的统计量作差异检验，获得检验值；

排序单元，耦合于获取单元和计算单元，将已预测基因按统计量或检验值进行排序，筛选出统计量或检验值最高的基因。

本发明单细胞分类装置中各个模块执行的具体操作体现在上述方法流程中，也可以根据以下实施例看出各个模块的具体操作。

下面通过具体实施例，对本发明作进一步说明。

实施例一：肾癌单细胞的分类

1-1、Reads定位

每个单细胞样本经测序得到的reads结果用SOAPaligner比对软件(http://soap.genomics.org.cn/soapaligner.html)与参考基因组序列(人类基因组HG18)进行比对，由于人类SNP在千份之二并且Reads的读长为100bp，所以SOAP比对时设置每条Reads最多有3个错误比对(Mismacth)，不容Gap的参数，以确保可比对上的Reads的位置准确。

1-2、基本数据统计

根据上述比对结果，计算出每个样本(单细胞或组织)相对于参考基因组序列的测序深度和覆盖度等结果，经统计得到全基因组测序且Mean Depth为3×附近时，由于PCR扩增存在一定的偏向性(Bias)，所以样本的覆盖度在55-90％之间波动较大。

表1肾癌单细胞样本全基因组测序的覆盖度和深度数据

其中RC-1至RC-15表示肾癌单细胞，共有15个单细胞样本；RN-1至RN-5表示癌旁单细胞；RN-T表示正常组织直接提取DNA测序，以作为对照用于数据分析和评估。主要以癌旁单细胞作为对照样本。也存在以癌旁单细胞以及正常组织都同时作为对照样本，如在去除杂合丢失的位点时就用到上述两个对照样本。

在外显子组测序中增加测序深度，当外显子目标区域的平均深度(Mean Depth)为30×附近时，目标区域覆盖度达到80-96％。从统计学意义上看，一个位点有四条reads支持则可以判断该位点碱基的准确性达到99％，而统计得到测序深度为4的位点占外显子区域的比例达到60-90％，说明外显子组测序的数据比全基因组测序得到的数据更好。

表2肾癌单细胞样本外显子组测序的覆盖度和深度数据

将上述2个表格进行对比，可得到全基因组测序的深度低而无法进行后续分析，而外显子组测序的深度高。另外考虑到测序成本问题，因此下面主要基于外显子组测序得到的数据进行分析。

1-3、数据过滤

根据数据特点挑选出每个DNA文库中由于PCR扩增次数过多引起的重复双末端reads，在多对重复双末端reads的序列完全一致时随机挑选一对reads，去除其他reads。

比如，图1中三对reads A、B、C序列完全一致并且比对到基因组上的起止位置也完全一致。在这种情况下只随机保留其中一对reads。

为了确认数据的准确性，考虑到人类外显子区域在基因组上并不可能有多个拷贝，即不可能是重复序列，所以外显子组测序得到的reads应大都唯一比对上人类参考基因组。为排除错误比对造成的影响，只挑选出唯一比对上参考基因组的reads进行分析(即hit数为1的reads)，这样在极大程度上减少了错误比对造成的影响。

如图2，Reads D比对上参考基因组的多个位置，而Reads E只比对上唯一一个位置，由于外显子在基因组上并不是重复区域，因此直接去除Reads D。

1-4、个体基因型的判定

我们充分考虑并利用人类基因组(该实施例中为参考基因组)上的已有信息，从NCBI网站下载人类基因组(HG18)对应的dbsnp，作为SOAPsnp的先验概率，并依据目前确定的研究结果，设置杂合位点SNP的先验概率为0.0010，纯合位点SNP的先验概率为0.0005。

在设置以上参数后，将步骤1-3经过滤的数据输入到SOAPsnp软件中，由SOAPsnp软件将过滤后数据与参考基因组比对，得到比对结果为CNS文件。

1-5、SNP数据集

由于每个细胞基因组存在一些低深度的区域，本发明综合所有细胞基因型的可能性的文件，利用最大似然法整合所有细胞的数据，产生包含所有样本每个位点的伪基因组(Pseudo-genome)。选择概率最大的基因型作为每个细胞的一致基因型，通过基因型和深度等信息检测出高质量的SNPs。得到样本的一致性序列后，结果保存为群体SNPs格式。

1-6、SNP的基因型

根据参考基因组SNP数据集位置信息，从CNS文件中提取每个细胞基因型的相应位点，得到群体SNP在各细胞相应位点的基因型文件。文件格式如表3所示。

“SNP位置”表示该SNP位点在染色体上的位置，“碱基型”对应某一个细胞基因组在这一位点的碱基型，深度为0的位点用“-”表示(即数据丢失的位点)。“样本ID”对应21个单细胞或组织DNA样本。

表3群体SNP在各细胞相应位点的基因型文件格式示意

该21个单细胞或组织DNA样本相对于人类基因组HG18在目标区域的群体SNP数目为93957个。其中杂合位点的组合由以下字母表示：

“M”代表“A和C”，“R”代表“A和G”，“W”代表“A和T”，“Y”代表“C和T”，“S”代表“C和G”，“K”代表“G和T”。

1-7、挑选与细胞突变相关的SNP位点

本发明主要在于寻找各细胞间的差异位点，因此须挑选出与细胞突变有关的位点。

表4非群体间SNP位点示意

首先，去除非群体间的SNP位点，这些位点如表4所示。在同一位点所有单细胞基因组的碱基型都一致，即在该位点21个细胞构成的群体都是一致的。经过计算，共有504个这样的位点，去除之后剩下93453个SNP位点。

其次，由于单细胞在WGA扩增时，存在只扩增出一对染色体中的一条的情况，导致等位基因丢失(Allele Dropout)，即原本是杂合位点，在测序时只测到其中一个碱基型，所以检测l出来的每个细胞在某些位点存在杂合丢失的现象，如表3所示的位点：

chr19  10079332  R R R R R R R R G R R R R R R R R R R R R

chr19  10079408  R R R R R R R R G R R R R R R R R R R R R

即第9个单细胞样本RC-9在该位点很可能原本是杂合位点，但由于只扩增到一个碱基型，所以判断为纯合位点。

为了排除这一类位点的干扰，考虑到同一位点同时出现几个样本都发生杂合丢失的概率极低，因此本发明采取策略如下：

第一，正常组织直接提取DNA测序的RN-T(即最后一列)必须为杂合位点。因为这个样本不经过WGA扩增，不可能出现杂合丢失。

第二，其它20个单细胞样本中杂合位点和数据丢失的样本数大于等于18。即同一位点最多只允许两个单细胞样本发生杂合丢失。因为三个或三个以上样本在同一位点同时发生杂合丢失的概率极低。

上述两个条件都要满足，即正常组织直接提取DNA测序的最后一列为杂合位点，其它20个单细胞样本中杂合位点和数据丢失的样本数大于等于18。这样的位点才除去。

经过计算，共有3975个这样的位点，经过这一步过滤得到SNP数目为89478个。

最后，为了得到与肾癌单细胞突变有关的位点，须去除已经公布的正常人的SNP位点，即去除人类基因组HG18的dbSNP、炎黄一号的SNP和千人基因组的SNP，得到与各细胞突变相关的SNP位点50524个。

1-8、群体结构分析

根据细胞群体突变SNP位点的基因型文件，分别运用生物信息群体分析中常用的方法对细胞进行分类。分类是由进化树的分枝和聚类情况来确定的。如图4，RC和RN在进化树中明显聚集为分离的两部分，因此将其分为两类。

1-9.1邻接N-J方法构树

如图4所示，根据进化树可以对细胞进行分类。

1-9.2MEGA软件

图5是最大似然法构成的关系树，根据关系树对细胞进行分类。

1-9.3PCA主成分分析

如图6、图7和图8所示的肾癌外显子组测序PCA结果，根据聚类情况对细胞进行分类。

1-9.4Structure群体结构

如图9，横坐标表示样本名称，纵坐标表示每个样本所属群体的概率，根据概率对单细胞进行分类。由图9所示，20个单细胞大致可以分为两个群体。肾癌外显子组测序Structure结果如图9所示。

1-10、群体结构分析结果

根据上述群体结构分析结果，提取每个细胞样本信息，剔除存在异议的细胞(分类不明确或明显离群样本)。从上述各种群体结构分析结果来看，取样正常，分类合理。这20个单细胞样本大致可以分为2个群体，即癌细胞群体(15个RC)和癌旁细胞群体(5个RN)，其中RC-1、RC-6、RC-7为癌细胞的一个亚群。

细胞样本信息指进行分析的单细胞中哪些是癌细胞、哪些是癌旁细胞(由取样时确定)，细胞样本信息仅作为参考，需要与聚类结果一起分析。若取样时的细胞样本信息认为是癌细胞和癌旁细胞，而聚类时正好分成两个群体聚集，说明结果互相对应；若取样时的细胞样本信息与聚类结果不一致，以聚类结果为准。

在癌细胞聚类中另行聚类在一起，由此确认RC-1、RC-6、RC-7为癌细胞的亚群。

1-11、筛选与肾癌相关的基因

根据上述RC和RN两个细胞群体在外显子区域SNP，通过一系列统计和检验进行这两个群体在外显子区域的比较，找出有明显差异的区域或者基因，即可筛选到与该例肾癌相关系数较高的基因。具体做法如下：

1-11.1从人类基因组数据库中下载HG18对应的注释文件，得到目前已预测的人类基因组中30000多个基因的起止位置。

1-11.2根据分类结果，得到RC和RN两个群体，计算每类群体中每个基因所有SNP位点的统计量，并累加所述统计量。

其中主要采用的计算统计量π的公式如下，π是衡量一个群体多态性高低的指标，a，b指某群体中两种碱基的样本数，公式可以为：

$\mathrm{\π}=\frac{a*b}{C_{a+b}^2}$

如在上述15个RC样本中，共有30条染色体，对于下述两个位点：位点1只有1条染色体为C，而另外29条染色体为T(a＝1，b＝29)；位点2有15条染色体为C，而另外15条染色体为T(a＝15，b＝15)。代入公式，得到位点1的π值是0.06，位点2的π值是0.517，则这2个位点的多态性具有明显差别。

在统计一个基因的多态性时，累加该基因所有位点的π值，由于非SNP位点的π值为0(a＝0或b＝0时，π＝0)，即对某群体累加该基因所有SNP位点的π值。

Pi_Gene＝∑_SnpInGeneπ

1-11.3将这30000多个基因按统计量或检验值进行排序，挑选出统计量或检验值最高的基因，查看这些基因的功能。

所采用的检验值至少为如下之一：检验值为Lod、Fst、Pbs，该实施例采用以上三种检验值。下面具体说明各个检验值以及计算过程。

分别对RC和RN这两个群体代入数据，可以得到和由于需要对这两个群体的差异进行比较，因此定义Lod如下：

${\mathrm{Lod}}_{\mathrm{Gene}}=1-({\mathrm{Pi}}_{\mathrm{Gene}}^{\mathrm{RC}}/{\mathrm{Pi}}_{\mathrm{Gene}}^{\mathrm{PN}})$

若和相差很小，即该基因在这两个群体中没有太大差异，Lod_Gene为0。可见Lod_Gene明显偏离于0，则可初步认为该基因是造成这两个群体分化的重要基因。

如上所述，分别统计HG18中30000多个基因的值Lod_Gene，然后从大到小进行排序，筛选出排名较前的基因。

F_ST(Fixation index)主要用来评价群体间的基因组距离和种群的差异，是度量种群间分化程度的一个指标，由Sewall Wright在1922年应用F-检验的一种特殊情况发展而来。

F_ST的零假设是在群体没有分化时，多态性位点在群内和群间的次等位碱基的频率差别是不显著的。计算F_ST的方法很多，虽然具体计算方法不同，但基本理论是一致的，即由Hudson(1992)给出的定义：

$F_{\mathrm{ST}}=\frac{{\mathrm{\Π}}_{\mathrm{Between}}-{\mathrm{\Π}}_{\mathrm{Within}}}{{\mathrm{\Π}}_{\mathrm{Between}}}$

其中∏_Between表示从两个群体(Between)中分别抽取一个样本，组成一对，计算这对样本SNP基因型的差异，如此可以计算所有成对样本SNP基因型的差异，最后求平均值即为∏_Between。

∏_Within表示从一个群体(Within)中分别抽取2个样本，组成一对，计算这对样本SNP基因型的差异，如此可以计算所有成对样本SNP基因型的差异，最后求平均值即为∏_Within。

如果有两个群体，是两个群体分别先计算∏_Within，然后累加。

结合目前已有SNP集的数据结构，基于上述原理，推导公式如下：

$F_{\mathrm{ST}}=\frac{{\mathrm{\Π}}_{\mathrm{Between}}-{\mathrm{\Π}}_{\mathrm{Within}}}{{\mathrm{\Π}}_{\mathrm{Between}}}$

$=1-\frac{{\mathrm{\Π}}_{\mathrm{Within}}}{{\mathrm{\Π}}_{\mathrm{Between}}}=1-\frac{\left[\underset{j}{\mathrm{\Σ}}\left(\begin{array}{c}{n}_j\\ 2\end{array}\right)\underset{i}{\mathrm{\Σ}}2\frac{n_{\mathrm{ij}}}{n_{\mathrm{ij}}-1}{x}_{\mathrm{ij}}\left({1-x}_{\mathrm{ij}}\right)\right]/\underset{j}{\mathrm{\Σ}}\left(\begin{array}{c}{n}_j\\ 2\end{array}\right)}{\underset{i}{\mathrm{\Σ}}2\frac{n_i}{n_i-1}{x}_i(1-x_i)}$

上式中x_ij是SNP位点i在群体j中的次等位碱基(第二碱基)的频率；而n_ij是SNP位点i在群体j中染色体上的物理位置；n_j则是群体j用于分析的SNP位点个数的总和。

其中变量j取RC和RN，变量i以最后判定的SNP位置代入。以基因为单位计算每个基因的Fst_Gene值，然后对HG18中30000多个基因的Fst_Gene值进行排序，筛选出排名较前的基因。

由于数据丢失(Missing data)的情况下，对SNP位点频率估计不精确，使得F_ST不能灵敏反映出数据原本的属性。根据参考文献Sequencingof 50 Human Exomes Reveals Adaptation to High Altitude.Scien ce，2July2010，329，75-78)采用的方法，对F_ST取log，并引入第三个群体(本实施例引入千人基因组中的部分数据，北京人基因组数据记为B)，定义了Pbs，公式如下：

T＝-log(1-Fst)

即三个群体两两对比的Fst如下：

T_RC-RN＝-log(1-Fst_RC-RN)

T_RC-B＝-log(1-Fst_RC-B)

T_RB-B＝-log(1-Fst_RN-B)

此时Pbs的公式如下：

$\mathrm{Pbs}=\frac{T_{\mathrm{RC}-\mathrm{RN}}+T_{\mathrm{RC}-B}-T_{\mathrm{RN}-B}}{2}$

以基因为计算单位，分别统计HG18中30000多个基因的Pbs_Gene值，然后排序筛选出排名较前的基因。

1-12、基因功能分析

根据以上三种检验值中至少一个，该实施例根据Lod、Fst以及Pbs筛选得到重要基因，分别进行功能分析。判断这些基因是否在某些通路中受到影响，从而与肾癌的发生机制相关。

实施例二：白血病单细胞的分类和筛选

2-1、Reads定位

对每个癌症单细胞进行30×深度的外显子组测序，得到的reads结果用SOAPaligner2.0比对软件与参考基因组序列(人类基因组HG18)进行比对。由于人的SNP在千份之二并且Reads的读长约为100bp，所以SOAP比对时我们设置每条Reads最多有2个错误比对(mismacth)，不允许出现Gap，以确保比对到参考基因组上Reads的准确性。

2-2、基本数据统计

共进行53个癌细胞、8个口腔上皮细胞(正常细胞)测序。表5为各细胞样本外显子组测序的覆盖度和深度数值信息。

表5各细胞样本外显子组测序的覆盖度和深度

2-3、数据过滤

与实施例一相同

2-4、个体基因型的判定

与实施例一相同

2-5、SNP数据集

在进行SNP数据集的确定时，考虑到血癌细胞数目较多，每个单细胞基因组外显子的覆盖率不是很高，且是基于每个个体进行SNPs的确定，所以我们选用较为严格的标准对得到的数据进行筛选。

标准如下：

Soapsnp软件中每个位点的一致性基因型的质量值不小于20，秩检验p值不小于1％；对于杂合变异的SNPs：位点基因型与参考基因组不同，最好碱基型测序质量值不小于20、测序深度不小于6，次好碱基型测序质量值不小于20，测序深度不小于2，且两基因型测序深度的比值在0.2-5范围内。

质量值越大，基因分型越正确，一般大于20时，错误率在万分之一以下，可忽略不计。

利用如上标准筛到可靠的SNPs以后，根据参考基因组SNP数据集位置信息确定位点，提取每个细胞每个位点的基因分型数据生成基因型文件。文件格式如表3所示。

2-6、群体结构分析

根据细胞群体突变SNP的基因型文件，我们分别运用多种生物信息群体分析中常用的方法对各细胞分类。

2-6.1、邻接N-J方法构成聚类构树

如图10所示为本发明53个癌细胞与8个正常细胞分类关系示意图，其中，ET-T1表示癌组织，NC-T1表示正常组织。

2-6.2、PCA主成分分析

如图11所示为本发明癌细胞与正常细胞聚类示意图，其中，LC表示癌细胞，L N表示正常细胞。

根据上面的群体分析结果，提取细胞样本信息，剔除存在异议的细胞(分类不明确，或离群样本)。从上面的群体结构说明取品正常，分类合理。

2-6.3、亚群的分类

依据进化树的形状或者情况，在所有53个癌细胞中可以清晰地分出4类亚群，说明在癌细胞中存在真实的差异。利用单细胞分析方法可以将同一癌症组织中不同的细胞亚群分类出来。

2-7、挑选高可信度的体细胞突变

从基因型文件中筛高可信度的体细胞突变位点，标准如下：

正常细胞拥有一致的纯合基因型，癌细胞存在两个或两个以上的杂合突变或纯合突变，且不能出现第三种纯合基因型以及和两种纯合基因型不一致的杂合基因型。比如正常细胞基因型为A，或突变类型是A-＞C，则癌细胞中只能出现三种基因型，即A、C和M，且C和M的数量不小于2。这种位点我们称为高可信度体细胞位点(High confident somaticmutation，HCSM)。由于我们利用的是外显子组测序技术，所以过滤掉非外显子区域的位点，共得到2296个HCSMs，其中同义位点879个，非同义位点1417(包含错义突变和截断突变位点)，非同义/同义突变比值为1.61，表6所示。

表6高可信度的体细胞突变位点

2-8、基因功能和途径分析

是在进行细胞分类和筛选之后可进行的下游分析，将基因突变位点的位置信息以及每个基因中存在的非同义突变位点数目作为基因功能富集的标准，利用Webgestalt在线分析工具对突变影响基因的功能以及途径进行研究(http://bioinfo.vanderbilt.edu/webgestalt/option.php)，发现突变主要集中在如下8类功能的基因中。

表7突变位点基因功能分析结果

途径分析发现突变基因主要集中在10种途径中，其中绝大多数与癌症的发生机制相关：

Metabolic pathways(代谢途径)

ECM-receptor interaction(ECM受体相互作用)

Pathways in cancer(癌症涉及的通路)

Viral myocarditis(病毒性心肌炎)

Type I diabetes mellitus(I型糖尿病)

MAPK signaling pathway(MAPK信号通路)

Focal adhesion(黏附斑)

Pantothenate and CoA biosynthesis(泛酸和辅酶A的生物合成)

Cell adhesion molecules(CAMs，细胞黏附分子)

Allograft rejection(移植排斥)

2-9、突变位点基因功能的预测

我们将外显子区域中非同义突变位点挑出，利用SIFT(http://sift.jcvi.org/)软件对这些突变位点对应的基因功能进行预测。结果分为4种情况，如下表8所示：

表8突变位点基因功能预测结果

4种情况是：对基因功能存在危害；低可信度危害，变异对于基因的功能影响不大以及无法判断。我们挑选含有功能危害性的突变位点并在上述功能富集和途径中出现的基因进行后续实验的验证。

Claims (12)

1.一种非诊断目的的单细胞分类方法，包括以下步骤：

将每个单细胞样本经测序得到的reads结果与参考基因组序列进行比对，并将比对结果进行数据过滤；

根据过滤后的数据确定每个单细胞样本的一致基因型，并将所有单细胞样本的一致基因型保存为SNP数据集；

从已保存的SNP数据集提取与参考基因组SNP数据集位置对应的位点的基因型文件；

挑选细胞突变SNP位点，并根据细胞突变SNP位点的基因型文件，对细胞进行分类；

其中，挑选细胞突变SNP位点的操作，包括以下步骤：

去除非群体间的SNP位点、去除存在杂合丢失的位点和/或去除已经公布的SNP位点，得到细胞突变SNP位点，以寻找各细胞间的差异位点；

根据过滤后的数据确定每个单细胞样本的一致基因型的操作，包括以下步骤：

根据过滤后的数据计算每个单细胞样本目标区域基因型的可能性文件；根据基因型的可能性文件产生包含所有单细胞样本每个位点的伪基因组，从伪基因组中选择概率最大的基因型作为每个单细胞样本的一致基因型。

2.根据权利要求1所述单细胞分类方法，将比对结果进行数据过滤的操作，包括以下步骤：

当存在多对重复双末端reads时，在多对reads的序列完全一致时随机选择一对reads，将所述多对reads中其他重复双末端reads去除；和/或去除未唯一比对上参考基因组序列的reads。

3.根据权利要求1所述单细胞分类方法，去除存在杂合丢失的位点的操作中，去除满足以下条件的位点：

当单细胞样本数为N时，N为自然数，正常组织直接提取DNA测序的最后一列为杂合位点；

其它N-1个单细胞样本中，杂合位点和数据丢失的样本数大于等于N-3。

4.根据权利要求1或2所述单细胞分类方法，从已保存的SNP数据集提取与参考基因组SNP数据集位置对应的位点的基因型文件的操作，包括以下步骤：

每个位点的一致性基因型的质量值不小于20，秩检验p值不小于1％；对于杂合变异的SNPs：最好碱基型测序质量值不小于20、测序深度不小于6，次好碱基型测序质量值不小于20，测序深度不小于2，且两基因型测序深度的比值在0.2-5范围内；

根据参考基因组SNP数据集位置信息确定位点，提取每个细胞每个位点的基因分型数据生成基因型文件。

5.根据权利要求1或2所述单细胞分类方法，对细胞进行分类之后，还包括以下步骤：

提取每个细胞样本信息，剔除存在异议的细胞。

6.根据权利要求1所述单细胞分类方法，对细胞进行分类之后，还包括以下步骤：

得到已预测基因组中基因的起止位置；

根据分类结果获得已分类群体，计算每类群体中每个基因所有SNP位点的统计量，并累加统计量；

对获得的统计量作差异检验，获得检验值；

将已预测基因按统计量或检验值进行排序，筛选出统计量或检验值最高的基因。

7.单细胞分类装置，包括：

数据过滤模块，将每个单细胞样本经测序得到的reads结果与参考基因组序列进行比对，并将比对结果进行数据过滤；

基因型确定模块，根据过滤后的数据确定每个单细胞样本的一致基因型，并将所有单细胞样本的一致基因型保存为SNP数据集；

基因型文件提取模块，从已保存的SNP数据集提取与参考基因组SNP数据集位置对应的位点的基因型文件；

分类模块，挑选细胞突变SNP位点，根据细胞群体突变SNP的基因型文件，对细胞进行分类；

其中，分类模块去除非群体间的SNP位点；去除存在杂合丢失的位点；和/或去除已经公布的SNP位点，得到细胞突变SNP位点，以寻找各细胞间的差异位点；

基因型确定模块根据过滤后的数据计算每个单细胞样本目标区域基因型的可能性文件；根据基因型的可能性文件产生包含所有单细胞样本每个位点的伪基因组，从伪基因组中选择概率最大的基因型作为每个单细胞样本的一致基因型。

8.根据权利要求7所述单细胞分类装置，其中：

数据过滤模块在有多对重复双末端reads时，从序列完全一致的多对reads中随机选择一对reads，将所述多对reads中其他重复双末端reads去除；和/或去除未唯一比对上参考基因组序列的reads。

9.根据权利要求7所述单细胞分类装置，其中：

分类模块去除的杂合丢失位点满足以下条件：当单细胞样本数为N时，N为自然数，正常组织直接提取DNA测序的最后一列为杂合位点；其它N-1个单细胞样本中，杂合位点和数据丢失的样本数大于等于N-3。

10.根据权利要求7或8所述单细胞分类装置，其中：

基因型文件提取模块根据参考基因组SNP数据集位置信息确定位点，提取每个细胞每个位点的基因分型数据生成基因型文件；

每个位点的一致性基因型的质量值不小于20，秩检验p值不小于1％；对于杂合变异的SNPs：最好碱基型测序质量值不小于20、测序深度不小于6，次好碱基型测序质量值不小于20，测序深度不小于2，且两基因型测序深度的比值在0.2-5范围内。

11.根据权利要求7或8所述单细胞分类装置，其中：

分类模块提取每个细胞样本信息，剔除存在异议的细胞。

12.根据权利要求7所述单细胞分类装置，还包括：

筛选模块，得到已预测基因组中基因的起止位置；根据分类结果获得已分类群体，计算每类群体中每个基因所有SNP位点的统计量，并累加统计量；对获得的统计量作差异检验，获得检验值；将已预测基因按统计量或检验值进行排序，筛选出统计量或检验值最高的基因。