CN115527610B - 一种单细胞组学数据的聚类分析方法 - Google Patents

一种单细胞组学数据的聚类分析方法 Download PDF

Info

Publication number
CN115527610B
CN115527610B CN202211396624.0A CN202211396624A CN115527610B CN 115527610 B CN115527610 B CN 115527610B CN 202211396624 A CN202211396624 A CN 202211396624A CN 115527610 B CN115527610 B CN 115527610B
Authority
CN
China
Prior art keywords
clustering
cell
data
stability
clustering result
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202211396624.0A
Other languages
English (en)
Other versions
CN115527610A (zh
Inventor
刘阿显
邵志峰
李华
郭妍
胡传圣
丹尼尔·马克·恰科夫
胡苗苗
邓永卓
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Suzhou Jingmai Biotechnology Co ltd
Shanghai Jiaotong University
Original Assignee
Suzhou Jingmai Biotechnology Co ltd
Shanghai Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Suzhou Jingmai Biotechnology Co ltd, Shanghai Jiaotong University filed Critical Suzhou Jingmai Biotechnology Co ltd
Priority to CN202211396624.0A priority Critical patent/CN115527610B/zh
Publication of CN115527610A publication Critical patent/CN115527610A/zh
Application granted granted Critical
Publication of CN115527610B publication Critical patent/CN115527610B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biophysics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本申请涉及一种单细胞组学数据的聚类分析方法,包括步骤:S1,将单细胞基因组学数据进行标准化预处理,获得预处理数据;S2,将S1的预处理数据进行多层次聚类,根据不断增加干扰条件来分析单细胞聚类的结构特征,再采用细胞聚类流向获得聚类结果;S3,对S2获得的聚类结果,采用罚函数计算聚类结果的稳定性分值,从而得到稳定性高的聚类结果;S4,根据S3获得的稳定性高的聚类结果的稳定性变化体系及对应的聚类数目出现的频次,得到最适聚类数目及细胞子群的类型和聚类结果。

Description

一种单细胞组学数据的聚类分析方法
技术领域
本发明涉及生物技术领域,更具体地,涉及一种单细胞组学数据的聚类分析方法。
背景技术
单细胞组学数据分析主要通过对单细胞转录组、单细胞核转录组或染色质等单细胞组学数据进行聚类等分析,并基于细胞子群的特异性表达基因等信息,对具有相似特性的单细胞归类,判断样本中存在的各种细胞表型以及它们的功能状态等。单细胞组学技术包括:single-cell RNA sequening(单细胞转录组测序或单细胞RNA测序,scRNA-seq)、Single Nuclei RNA Sequencing(单细胞核转录组测序或单细胞RNA测序,snRNA-Seq)、single-cell ATAC sequencing(单细胞ATAC测序,sc-ATAC-Seq)、single-cell Hi-C(单细胞Hi-C技术或单细胞高通量/分辨率染色体构象捕获技术),由于测序反应步骤中动力学以及效率的限制,数据缺失度相对较高,所获得的单个细胞的组学数据信噪比低,例如单细胞转录组技术所能获得的转录本不到每个细胞中的10%。由此可见,基于每一个细胞的分析不具有任何统计意义,无法对其生物学状态进行表征,单细胞组学数据只有通过聚类分析方法将具有相似特性的细胞通过聚类归类,将具有不同特性的单细胞归类生成细胞子群(cluster),并为每一个细胞子群赋予特定的细胞表型,才能获得具有生物学意义的、可以确定特定细胞表型的结果。这些结果可以和先验的知识进行比对以确定新的表型或生物学状态,还可以通过网络分析、不同细胞子群之间的关系等揭示细胞与细胞之间的应答、已知表型的变异等。
聚类分析方法作为单细胞组学数据分析最基本的步骤,已经成为单细胞组学技术的核心之一。虽然,诸多单细胞组学数据,尤其是单细胞转录组数据,其聚类分析方法已经被发展并获得了广泛应用,但是,这些聚类分析方法都面临如下的挑战:
1、一般的聚类分析方法,需要直接或者间接的人为指定(自定义)聚类数目,其聚类结果不是直接或完全根据数据情况来驱动的。例如,K-means聚类算法、分层聚类算法和Seurat聚类算法,其解释在很大程度上依赖于用户主观判断的聚类数目和先验知识。但是,主观选择的聚类数目难以解决所选聚类数目对应聚类结果的质量问题,如聚类结果是否聚类不足或过度聚类。尽管聚类方法根据细胞相似性划分细胞,但细胞类型的注释一般是用户人工手动完成,导致了聚类结果中新颖且从未用传统方法表征的细胞表型很大程度上可能是不正确分类得到的结果,新表型的错误识别将导致基于这种不正确分类对相关研究的误导。
2、最新的聚类分析方法,证明一个单细胞组学数据集,确实存在由数据本身决定的一个或数个稳定聚类结果(最优聚类),使得利用数据本身来获得最适聚类数目成为可能。最优聚类是指能够对诸如参数的轻微扰动保持稳定,这种稳定趋势及其变化则表征了单细胞数据的稳定性。其通过多次取样作为聚类扰动条件(MultiK),以相同聚类数目出现频次作为判断最适聚类数目的主要依据。但是,这种基于多次采样的方法,需要从原单细胞数据集中随机构建数千次或更多的亚集。第一方面,该方法每次抽样得到的亚集之间只有部分单细胞是相同的,无法完全保证各亚集之间聚类结果的一致性。第二方面,该方法认为获得可靠的聚类结果需要抽样至少400次并对每个抽样执行40个不同的分辨率参数对应的聚类,其时间复杂度随细胞数量和抽样次数的增长呈线性增长,需要庞大的计算量和较高配置的设备,尤其是含有数万数十万细胞的单细胞数据集,其无法满足目前单细胞数据规模不断增长的发展现状。第三方面,以聚类数目的出现频次作为最适聚类数目的判据依据也存在其内在的局限性,例如,对于某些特定的单细胞数据,即使在不同划分尺度下得到了相同的聚类数目,对应细胞子群中的细胞归属并不能保证是完全相同的;换句话说,聚类数目的稳定,并不能完全代表每一个细胞子群中的细胞也是稳定的。
有鉴于此,本申请提供一种单细胞组学数据的分析方法,对于任何一个单细胞组学数据,由数据本身获得准确的最适聚类数目及相关细胞子群之间的层级关系,其聚类结果一致性高、细胞子群中的细胞稳定,并且计算量小。
发明内容
本发明的目的在于,提供一种单细胞组学数据的分析方法,对于任何一个单细胞组学数据,由数据本身获得准确的最适聚类数目及相关细胞子群之间的层级关系,其聚类结果一致性高、细胞子群中的细胞稳定,并且计算量小。
一种单细胞组学数据的聚类分析方法,包括步骤:
S1,将单细胞基因组学数据进行标准化预处理,获得预处理数据;
S2,将S1的预处理数据进行多层次聚类,根据不断增加干扰条件来分析单细胞聚类的结构特征,再采用细胞聚类流向获得聚类结果;
S3,对S2获得的聚类结果,采用罚函数计算聚类结果的稳定性分值,从而得到稳定性高的聚类结果;
S4,根据S3获得的稳定性高的聚类结果的稳定性变化体系及对应的聚类数目出现的频次,得到最适聚类数目及细胞子群的类型和聚类结果。
在一些实施方式中,在步骤S1中,所述单细胞基因组学数据为能够用于聚类分析的单细胞数据,单细胞数据采用软件与人类参考基因组进行比对,数据质控后获得相应的可直接用于聚类分析的单细胞基因组学数据;单细胞数据包括:单细胞转录组测序scRNA-Seq数据、单细胞核转录组测序snRNA-Seq数据、单细胞ATAC测序sc-ATAC-Seq数据、单细胞Hi-C测序数据。
在一些实施方式中,在步骤S1中,将单细胞基因组学数据采用相应的软件进行标准化预处理。
进一步的,对于单细胞转录组及单细胞核转录组的基因组学数据的标准化预处理包括:滤除低质量细胞和信息重叠细胞、基因表达量的对数变换、高可变基因的归一化及降维以及构建k近邻稀疏矩阵和边缘权重;对于单细胞ATAC测序peaks矩阵的预处理包括:根据fragments数目和TSS enrichment score的阈值滤除低质量细胞和双细胞、peaks标准化、降维、构建k近邻稀疏矩阵和边缘权重;对于单细胞Hi-C测序数据的预处理包括:使用线性卷积和漫步模型对数据进行imputation、数据降维、构建k近邻稀疏矩阵和边缘权重。
在一些实施方式中,在步骤S2中,对预处理数据,使用聚类算法中表征划分尺度参数的增加为聚类稳定性的干扰条件进行多层次聚类,多层次聚类包括步骤:
(1)将S1中预处理数据作为输入,根据单细胞信息表征的细胞间权重,使用起始聚类参数完成初始聚类划分;使用模块度对初始聚类划分结果进行质量评估,并根据模块度得分情况,对处于聚类边缘的细胞进行重新调整,得到可靠的初始聚类结果;
(2)定义聚类划分干扰单位和层级聚类变化范围,以初始聚类结果为基准,首先增加一个干扰单位,完成新的聚类划分;使用模块度对当前聚类划分进行评分,并根据当前干扰单位下细胞间权重进行聚类结果的进一步调整,进而得到当前聚类“粒度”下的可靠的聚类结果;
(3)以当前聚类结果为基准,计算下一个干扰单位下的聚类结果;通过依次对同一基准下的聚类结果增加干扰,得到单细胞组学数据多层级聚类结果的聚类数目和不同层级聚类之间细胞的流向。
进一步优选的,经使用不同规模的单细胞数据集进行测试,定义0.1至1.5的“粒度”单位为层级聚类变化范围,定义0.1为干扰单位,能够提供足够的聚类干扰能力,同时足以捕捉数据中的结构和子结构。
在一些实施方式中,在步骤S2中,借助可视化工具显示不同参数下的聚类结果之间细胞的流向关系,并追踪每一个细胞在cluster数目逐步增加的过程中流向或定位,初步探究单细胞数据在聚类分析中体现的数据结构特征;在随聚类分辨率参数变化的过程中,部分细胞出现归属跳跃及分群现象,被认为是一种相对不稳定的聚类状态,聚类结果在随分辨率参数增加后,其细胞归属保持不变,对应的聚类结果最为稳定,从而获得聚类结果。
进一步的,单细胞数据的整体结构的稳定性及其变化,可以利用聚类结果的无间断表征,通过构建连续参数的聚类结果稳定性评分系统,将相邻划分单位下的聚类稳定性进行关联,评估相邻分辨率参数下聚类结果的一致性得分。
在一些实施方式中,在步骤S3中,根据不同层级下每一个细胞子群中细胞的聚类流向定义罚函数,根据罚函数量化不同层级下聚类结果的稳定性特征(即细胞归属随参数变化的特征)。
进一步的,不同层级下每一个细胞子群中细胞的聚类流向包含:经聚类干扰条件保持不变的细胞子群、出现独立分支的细胞子群以及出现汇入行为的细胞子群;其中,细胞归属保持不变的细胞子群表征了该单细胞数据稳定的细胞类群,计作V,出现分支及出现汇入行为的细胞子群表征了该单细胞数据中不稳定的细胞对数目,分别计作C1、C2,根据公式c1计算出相邻划分单位下的聚类结果不一致的细胞对数目UI,公式c1为:
进一步的,根据标准化原理,构建罚函数F,计算相邻参数下聚类结果的稳定性,罚函数F的计算公式c2为:
其中,UI为相邻划分单位下的聚类结果不一致的细胞对数目,E(UI)为邻划分单位下的聚类结果不一致细胞对数目的期望值,max(UI)为随机模型下相邻划分单位聚类结果不一致的理论最大值。
进一步的,使用罚函数F计算得到各层级下聚类结果的对应的稳定性得分scScore(i),计算公式为C3:
scScore(i)=1-f(i)
其中,F(i)为对应干扰条件下稳定性罚分,scScore(i)为对应干扰单位下的聚类结果的稳定性得分。
进一步的,根据各干扰单位下的稳定性得分,构建不同层级下聚类结果稳定性得分的变化体系;以整体得分的平均值作为稳定性阈值,过滤掉稳定性低于稳定性阈值的聚类结果,得到稳定性高的聚类结果。
在一些实施方式中,在步骤S4中,包括步骤:
(1)根据S3获得的稳定性高的聚类结果的稳定性得分变化趋势,确定出现连续稳定性得分的变化区间(由于聚类结果的内部结构连续稳定性通常呈现出从稳定态到混乱态再回到稳定态的趋势);
(2)统计连续稳定性得分的变化区间聚类数目的变化情况,选择稳定区间内具有一定一致性的聚类数目作为待选最适聚类数目;
(3)比较待选最适数目的稳定性得分,选择前两个稳定性区间中稳定性得分最高的聚类数目,作为最终的最适聚类数目。
进一步优选的,所述待选最适聚类数目为2-3个,通常第一个稳定点下的聚类结果对应了该数据的基础细胞类型,第二个稳定点则是包含了亚细胞类型或亚细胞状态的聚类结果。
与现有技术相比,本发明的有益效果在于:
(1)本申请的聚类分析方法,由数据本身获得准确的最适聚类数目及相关细胞子群之间的层级关系。通过使用聚类划分尺度参数作为干扰条件对单细胞数据进行多尺度下的聚类,能够更加直接地发现数据本身的稳定性,从而摆脱了MultiK需要从原单细胞数据集中随机构建数千次或更多的亚集的限制,并且不同尺度的聚类结果均基于原单细胞数据,克服了MultiK由于各亚集之间包含不同细胞导致的聚类结果不一致的缺陷。
(2)本申请采用多层次聚类以及细胞的流向关系得到聚类结构,能够以聚类结果的在不同程度的干扰条件下的稳定性整体变化。克服了以往大多数方法仅在单个参数下讨论聚类结果的局限性,且稳定性的变化体系能够直观地体现不同干扰水平下的聚类结果的从属关系,更有利于解析不同细胞表型之间的相关性和系统进化关系。对于发育、疾病发生过程等重要生物学问题,具有重要的实用意义。
(3)本申请使用表征的单细胞数据聚类分析稳定性变化体系,定义不同细胞的轨迹行为所表征的数据内部结构的稳定性并构建罚函数。创新地对单细胞在不同划分尺度下的稳定性进行量化,以分数变化体系的形式表征单细胞数据的整体稳定性变化。
(4)本申请判断最适聚类数目,统计连续稳定性得分的变化区间聚类数目的变化情况,选择稳定区间内具有一定一致性的聚类数目作为待选最适聚类数目;然后比较待选最适数目的稳定性得分,选择前两个稳定性区间中稳定性得分最高的聚类数目,作为最终的最适聚类数目。基于此判断标准,不仅考虑了由于单细胞数据包含不同亚型和亚细胞状态引起的可靠聚类结果多种可能性,并且考虑到了聚类结果随参数增加存在过度划分的可能性,一般情况下仅取前两个稳定状态并根据稳定性变化趋势过滤掉了中间不稳定状态,最终再根据稳定性最高得分实现了更加准确的判断。
(5)本申请的聚类分析方法,对于任何单细胞组学数据仅需运行15次不同划分单位下的聚类结果。相比MultiK每次至少运行4000次的聚类过程,克服了其在时间复杂度的缺陷。在实施例1中,使用相同设备MultiK需要运行约58000s,本申请则仅需41s。因此,本申请能够花费较少的时间和资源且对分析设备的配置要求相对较低,实现准确可靠的最适聚类数目的判断且能够快速地在本地化平台快速完成最适聚类数目的判断。
附图说明
结合以下附图一起阅读时,将会更加充分地描述本申请内容的上述和其他特征。可以理解,这些附图仅描绘了本申请内容的若干实施方式,因此不应认为是对本申请内容范围的限定。通过采用附图,本申请内容将会得到更加明确和详细地说明。
图1为本申请实施例1的三种细胞系混合单细胞数据多分辨率下聚类结果变化关系图。
图2为本申请实施例1的三种细胞系混合单细胞数据多分辨率下聚类稳定性变化图。
图3为本申请实施例1的三种细胞系混合单细胞数据在第一个稳定点下的聚类结果图。
图4为本申请实施例1的三种细胞系混合单细胞数据在第二个稳定点下的聚类结果图。
图5为本申请实施例2的FVB3乳腺单细胞数据多分辨率下聚类结果变化关系图。
图6为本申请实施例2的FVB3乳腺单细胞数据多分辨率下聚类稳定性变化图。
图7为本申请实施例2的FVB3乳腺单细胞数据在第一个稳定点下的聚类结果图。
图8为本申请实施例2的FVB3乳腺单细胞数据在第二个稳定点下的聚类结果图。
具体实施方式
描述以下实施例以辅助对本申请的理解,实施例不是也不应当以任何方式解释为限制本申请的保护范围。
下列实施例中未注明具体条件的实验方法,按照常规实验条件,例如Sambrook等人的分子克隆实验室手册(New York:Cold Spring Harbor Laboratory Press,1989)中所述的条件,或按照制造厂商所建议的条件。除非另外说明,否则百分比和份数按重量计算。除非有特别说明,否则实施例所用的材料均为市售产品。
实施例1:三种细胞系混合单细胞数据
步骤1:获得样本的单细胞基因组学数据
采用三种人类细胞系混合样本的公共数据(ID:GSE136148):为人乳腺癌高转移细胞(MDA-MB-438),人乳腺癌细胞(MCF7)和人真皮成纤维细胞(HF)按照6:3:1比例混合,作为单细胞转录组测序原始数据,该数据包含的细胞类型及亚细胞状态已知。
步骤2:获得单细胞转录组表达矩阵
对于上述单细胞转录组测序原始数据,使用Cell Ranger Single Cell SoftwareSuite 1.3版进行解测序文库拆分、barcode和UMI处理、和单细胞3'端基因计数,并将reads比对到人类参考基因组Hg19上。
使用Samtools对生成的BAM文件进行排序和索引,并使用Picard执行数据质量控制,使用Salmon进行转录本读取计数,获得单细胞转录组表达矩阵。
步骤3:标准化预处理
读取单细胞转录测序数据表达矩阵并创建Seurat对象。根据每个单细胞数据表达量总计数数或检测到的基因数或表达的线粒体基因比例进行表达矩阵的过滤处理。若表达量总计数数或检测到的基因数或表达的线粒体基因比例大于或小于预定义阈值(所有细胞的中间值±3×中值绝对偏差)则移除该细胞,从而过滤去除可能的低质量细胞或含量两个或多个细胞转录组信息的单细胞。
需要对单细胞表达矩阵进行标准化预处理,也既使用每个细胞中的表达基因的表达量,除以每个单细胞的的所有基因的表达量总数,乘以10000的比例因子,然后进行对数变换来实现。为了识别高度可变的基因,采用在Seurat软件包中应用了方差稳定转换(“vst”)方法,选择前2000个高度可变基因进行归一化,也既使得每个基因在细胞中的表达量平均值为0,方差为1。使用所选前2000个高度可变基因执行PCA降维,使用前30个主成分中构造一个K-最近邻图,然后,使用Jaccard相似性度量,基于其局部邻域中的共享重叠来细化边缘权重。
步骤:4:多层次聚类
将标准化预处理后的数据作为输入,根据单细胞信息表征的细胞间权重,使用起始聚类参数完成初始聚类划分;使用模块度对初始聚类划分结果进行质量评估,并根据模块度得分情况,对处于聚类边缘的细胞进行重新调整,得到可靠的初始聚类结果。定义聚类划分干扰单位和层级聚类变化范围。经使用不同规模的单细胞数据集进行测试,定义0.1至1.5的“粒度”单位为层级聚类变化范围,定义0.1为干扰单位,能够提供足够的聚类干扰能力,同时足以捕捉数据中的结构和子结构。以初始聚类结果为基准,首先增加一个干扰单位,完成新的聚类划分;使用模块度对当前聚类划分进行评分,并根据当前干扰单位下细胞间权重进行聚类结果的进一步调整,进而得到当前聚类“粒度”下的可靠的聚类结果。以当前聚类结果为基准,计算下一个干扰单位下的聚类结果;通过依次对同一基准下的聚类结果增加干扰,得到单细胞组学数据多层级聚类结果的聚类数目和不同层级聚类之间细胞的流向。
借助可视化工具显示不同参数下的聚类结果之间细胞的流向关系,并追踪每一个细胞在cluster数目逐步增加的过程中流向或定位,初步探究单细胞数据在聚类分析中体现的数据结构特征;在随聚类分辨率参数变化的过程中,部分细胞出现归属跳跃及分群现象,被认为是一种相对不稳定的聚类状态,聚类结果在随分辨率参数增加后,其细胞归属保持不变,对应的聚类结果最为稳定,从而获得聚类结果。单细胞数据的整体结构的稳定性及其变化,可以利用聚类结果的无间断表征,通过构建连续参数的聚类结果稳定性评分系统,将相邻划分单位下的聚类稳定性进行关联,评估相邻分辨率参数下聚类结果的一致性得分,获得聚类结果。得到三种细胞系混合单细胞数据多分辨率下聚类结果变化关系图,如图1所示。从图1可以看出,该数据在不同干扰尺度下的整体稳定性变化,在以分辨率为0.1的初始聚类结果中,整体结构首先被初步划分成3个主体,伴随着干扰尺度的不断增加,部分主体的细胞出现了跳跃和分割的不稳定性表现并在一定区间内形成相对稳定的聚类结果。
步骤5:判断稳定性分值,得到稳定性高的聚类结果
根据不同层级下每一个细胞子群中细胞的聚类流向定义罚函数,根据罚函数量化不同层级下聚类结果的稳定性特征(即细胞归属随参数变化的特征)。不同层级下每一个细胞子群中细胞的聚类流向包含:经聚类干扰条件保持不变的细胞子群、出现独立分支的细胞子群以及出现汇入行为的细胞子群;其中,细胞归属保持不变的细胞子群表征了该单细胞数据稳定的细胞类群,计作V,出现分支及出现汇入行为的细胞子群表征了该单细胞数据中不稳定的细胞对数目,分别计作C1、C2,根据公式c1计算出相邻划分单位下的聚类结果不一致的细胞对数目UI,公式c1为:
根据标准化原理,构建罚函数F,计算相邻参数下聚类结果的稳定性,罚函数F的计算公式c2为:
其中,UI为相邻划分单位下的聚类结果不一致的细胞对数目,E(UI)为邻划分单位下的聚类结果不一致细胞对数目的期望值,max(UI)为随机模型下相邻划分单位聚类结果不一致的理论最大值。
使用罚函数F计算得到各层级下聚类结果的对应的稳定性得分scScore(i),计算公式为C3:
scScore(i)=1-f(i)
其中,F(i)为对应干扰条件下稳定性罚分,scScore(i)为对应干扰单位下的聚类结果的稳定性得分。根据各干扰单位下的稳定性得分,构建不同层级下聚类结果稳定性得分的变化体系;以整体得分的平均值作为稳定性阈值,过滤掉稳定性低于稳定性阈值的聚类结果,得到稳定性高的聚类结果。
其三种细胞系混合单细胞数据多分辨率下聚类稳定性变化图,如图2所示,从图2可知,该数据聚类结果的整体稳定性变化从稳定区间进入不稳定状态后进入了另一个稳定区间。第一个稳定区间位于分辨率0.1-0.3之间,第二个稳定性区间位于0.3-0.8之间。根据两个区间的稳定性得分情况,前两个稳定性区间最高得分点分别为0.89和0.99,并且均处于整体平均稳定性得分以上,分别对应3个和7个细胞子群,其中聚类数目为7的结果形成一个相对更稳定的区间。
步骤6:获得最适聚类数目及细胞子群的类型和聚类结果
对于稳定性高的聚类结果的稳定性得分变化趋势,确定出现连续稳定性得分的变化区间(由于聚类结果的内部结构连续稳定性通常呈现出从稳定态到混乱态再回到稳定态的趋势)。统计连续稳定性得分的变化区间聚类数目的变化情况,选择稳定区间内具有一定一致性的聚类数目作为待选最适聚类数目,选择了2个待选最适聚类数目。确定第一个稳定点分辨率为0.1,第二个稳定点分辨率为0.5。
使用ARI(Adjusted Rand Index)对原分析的结果和本申请的结果进行比较,细胞子群数目及细胞归属与原聚类分析结果保持高度一致。第一个稳定点对应的聚类结果为包含三种细胞类型,如图3所示,其对应的聚类结果,分别是成纤维细胞(COL1A1,VIM)、管腔细胞(GATA3,AREG)和基底细胞(KRT7,KRT17)。第二个稳定点对应的聚类结果为包含三种细胞类型,如图4所示,第二个稳定点包含7个细胞子群,其中基底细胞进一步地被分成了4个细胞亚群,其中3个是反映细胞周期不同阶段的基础细胞系群,最后一个是典型的基础细胞系特征群。通过与该套数据经过验证的原分析结果比较,对应的聚类数目及细胞归属高度一致。因此,本申请能够准确地判断该数据聚类结果的稳定性表征的最适聚类数目以及所包含的真实细胞类型及亚细胞类型或状态。
实施例2:FVB3乳腺单细胞数据
步骤1:获得样本的单细胞基因组学数据
采用公共数据人体细胞系FVB3乳腺细胞(ID:GSE136148),作为单细胞转录组测序原始数据。
步骤2:获得单细胞转录组表达矩阵
对于上述单细胞转录组测序原始数据,使用Cell Ranger Single Cell SoftwareSuite 1.3版进行解测序文库拆分、barcode和UMI处理、和单细胞3'端基因计数,并将reads比对到人类参考基因组Hg19上。使用Samtools对生成的BAM文件进行排序和索引,并使用Picard执行数据质量控制,使用Salmon进行转录本读取计数,获得单细胞转录组表达矩阵。
步骤3:标准化预处理
读取单细胞转录测序数据表达矩阵并创建Seurat对象。根据每个单细胞数据表达量总计数数或检测到的基因数或表达的线粒体基因比例进行表达矩阵的过滤处理。若表达量总计数数或检测到的基因数或表达的线粒体基因比例大于或小于预定义阈值(所有细胞的中间值±3×中值绝对偏差)则移除该细胞,从而过滤去除可能的低质量细胞或含量两个或多个细胞转录组信息的单细胞。
需要对单细胞表达矩阵进行标准化预处理,也既使用每个细胞中的表达基因的表达量,除以每个单细胞的的所有基因的表达量总数,乘以10000的比例因子,然后进行对数变换来实现。为了识别高度可变的基因,采用在Seurat软件包中应用了方差稳定转换(“vst”)方法,选择前2000个高度可变基因进行归一化,也既使得每个基因在细胞中的表达量平均值为0,方差为1。使用所选2000个高度可变基因执行PCA降维,使用前30个主成分中构造一个K-最近邻图,然后,使用Jaccard相似性度量,基于其局部邻域中的共享重叠来细化边缘权重。
步骤:4:多层次聚类
将标准化预处理后的数据作为输入,根据单细胞信息表征的细胞间权重,使用起始聚类参数完成初始聚类划分;使用模块度对初始聚类划分结果进行质量评估,并根据模块度得分情况,对处于聚类边缘的细胞进行重新调整,得到可靠的初始聚类结果。定义聚类划分干扰单位和层级聚类变化范围。经使用不同规模的单细胞数据集进行测试,定义0.1至1.5的“粒度”单位为层级聚类变化范围,定义0.1为干扰单位,能够提供足够的聚类干扰能力,同时足以捕捉数据中的结构和子结构。以初始聚类结果为基准,首先增加一个干扰单位,完成新的聚类划分;使用模块度对当前聚类划分进行评分,并根据当前干扰单位下细胞间权重进行聚类结果的进一步调整,进而得到当前聚类“粒度”下的可靠的聚类结果。以当前聚类结果为基准,计算下一个干扰单位下的聚类结果;通过依次对同一基准下的聚类结果增加干扰,得到单细胞组学数据多层级聚类结果的聚类数目和不同层级聚类之间细胞的流向。
借助可视化工具显示不同参数下的聚类结果之间细胞的流向关系,并追踪每一个细胞在cluster数目逐步增加的过程中流向或定位,初步探究单细胞数据在聚类分析中体现的数据结构特征;在随聚类分辨率参数变化的过程中,部分细胞出现归属跳跃及分群现象,被认为是一种相对不稳定的聚类状态,聚类结果在随分辨率参数增加后,其细胞归属保持不变,对应的聚类结果最为稳定,从而获得聚类结果。单细胞数据的整体结构的稳定性及其变化,可以利用聚类结果的无间断表征,通过构建连续参数的聚类结果稳定性评分系统,将相邻划分单位下的聚类稳定性进行关联,评估相邻分辨率参数下聚类结果的一致性得分,获得聚类结果。得到FVB3乳腺单细胞数据多分辨率下聚类结果变化关系图,如图5所示。从图5可以看出,该数据在不同干扰尺度下的整体稳定性变化,在以分辨率为0.1的初始聚类结果中,整体结构首先被初步划分成7个主体,伴随着干扰尺度的不断增加,部分主体的细胞出现了跳跃和分割的不稳定性表现并在一定区间内形成相对稳定的聚类结果。
步骤5:判断稳定性分值,得到稳定性高的聚类结果
根据不同层级下每一个细胞子群中细胞的聚类流向定义罚函数,根据罚函数量化不同层级下聚类结果的稳定性特征(即细胞归属随参数变化的特征)。不同层级下每一个细胞子群中细胞的聚类流向包含:经聚类干扰条件保持不变的细胞子群、出现独立分支的细胞子群以及出现汇入行为的细胞子群;其中,细胞归属保持不变的细胞子群表征了该单细胞数据稳定的细胞类群,计作V,出现分支及出现汇入行为的细胞子群表征了该单细胞数据中不稳定的细胞对数目,分别计作C1、C2,根据公式c1计算出相邻划分单位下的聚类结果不一致的细胞对数目UI,公式c1为:
根据标准化原理,构建罚函数F,计算相邻参数下聚类结果的稳定性,罚函数F的计算公式c2为:
其中,UI为相邻划分单位下的聚类结果不一致的细胞对数目,E(UI)为邻划分单位下的聚类结果不一致细胞对数目的期望值,max(UI)为随机模型下相邻划分单位聚类结果不一致的理论最大值。
使用罚函数F计算得到各层级下聚类结果的对应的稳定性得分scScore(i),计算公式为C3:
scScore(i)=1-f(i)
其中,F(i)为对应干扰条件下稳定性罚分,scScore(i)为对应干扰单位下的聚类结果的稳定性得分。根据各干扰单位下的稳定性得分,构建不同层级下聚类结果稳定性得分的变化体系;以整体得分的平均值作为稳定性阈值,过滤掉稳定性低于稳定性阈值的聚类结果,得到稳定性高的聚类结果。
其FVB3乳腺单细胞数据多分辨率下聚类稳定性变化图,如图6所示,从图6可知,该数据聚类结果的整体稳定性变化从稳定区间进入不稳定状态后进入了另一个稳定区间。第一个稳定区间位于分辨率0.1-0.3之间,第二个稳定性区间位于0.5-1.2之间。根据两个区间的稳定性得分情况,前两个稳定性区间最高得分点分别为0.98和0.99,并且均处于整体平均稳定性得分以上,分别对应8个和17个细胞子群,其中聚类数目为17的结果形成一个相对更稳定的区间。
步骤6:获得最适聚类数目及细胞子群的类型和聚类结果
对于稳定性高的聚类结果的稳定性得分变化趋势,确定出现连续稳定性得分的变化区间(由于聚类结果的内部结构连续稳定性通常呈现出从稳定态到混乱态再回到稳定态的趋势)。统计连续稳定性得分的变化区间聚类数目的变化情况,选择稳定区间内具有一定一致性的聚类数目作为待选最适聚类数目,选择了2个待选最适聚类数目。确定第一个稳定点分辨率为0.2,第二个稳定点分辨率为0.1。
使用ARI(Adjusted Rand Index)对原分析的结果和本申请的结果进行比较,细胞子群数目及细胞归属与原聚类分析结果保持高度一致。其中第一个稳定点对应的聚类结果为包含8种基础细胞类型,如图7所示。第二个稳定点包含17个细胞类型及亚型/亚状态,如图8所示,经表达特异性表达基因注释,准确地得到对应的细胞组成。
上述基于单细胞组学数据在多层次聚类下细胞轨迹定义的罚函数表征了数据内在结构的连续稳定性特征。根据“最合理”聚类的细胞归属随聚类参数变化保持稳定的判定标准,用于指导判断单细胞数据的最适聚类数目,对于基于聚类结果进行确定细胞表型并指导后续的相关分析和探究的单细胞分析具有广泛的应用价值。对于目前单细胞组学数据分析领域探究聚类数目问题提供了新的便捷可靠的方法思路。对于单细胞组学数据确定新的表型或生物学状态、通过网络分析、不同细胞子群之间的关系等揭示细胞与细胞之间的应答、已知表型的变异等,以及为更进一步的实验提供基础或提供机理性的分析具有重要参考意义和实用价值。
尽管本申请已公开了多个方面和实施方式,但是其它方面和实施方式对本领域技术人员而言将是显而易见的,在不脱离本申请构思的前提下,还可以做出若干变形和改进,这些都属于本申请的保护范围。本申请公开的多个方面和实施方式仅用于举例说明,其并非旨在限制本申请,本申请的实际保护范围以权利要求为准。

Claims (8)

1.一种单细胞组学数据的聚类分析方法,其特征在于,包括步骤:
S1,将单细胞基因组学数据进行标准化预处理,获得预处理数据;
S2,将S1的预处理数据进行多层次聚类,根据不断增加干扰条件来分析单细胞聚类的结构特征,再采用细胞聚类流向获得聚类结果,
其中,对预处理数据,使用聚类算法中表征划分尺度参数的增加为聚类稳定性的干扰条件进行多层次聚类,多层次聚类包括步骤:
(1)将S1中预处理数据作为输入,根据单细胞信息表征的细胞间权重,使用起始聚类参数完成初始聚类划分;使用模块度对初始聚类划分结果进行质量评估,并根据模块度得分情况,对处于聚类边缘的细胞进行重新调整,得到可靠的初始聚类结果;
(2)定义聚类划分干扰单位和层级聚类变化范围,经使用不同规模的单细胞数据集进行测试,定义0.1至1.5的“粒度”单位为层级聚类变化范围,定义0.1为干扰单位,能够提供足够的聚类干扰能力,同时足以捕捉数据中的结构和子结构;以初始聚类结果为基准,首先增加一个干扰单位,完成新的聚类划分;使用模块度对当前聚类划分进行评分,并根据当前干扰单位下细胞间权重进行聚类结果的进一步调整,进而得到当前聚类“粒度”下的可靠的聚类结果;
(3)以当前聚类结果为基准,计算下一个干扰单位下的聚类结果;通过依次对同一基准下的聚类结果增加干扰,得到单细胞组学数据多层级聚类结果的聚类数目和不同层级聚类之间细胞的流向;
S3,对S2获得的聚类结果,采用罚函数计算聚类结果的稳定性分值,从而得到稳定性高的聚类结果;
S4,根据S3获得的稳定性高的聚类结果的稳定性变化体系及对应的聚类数目出现的频次,得到最适聚类数目及细胞子群的类型和聚类结果。
2.如权利要求1所述的单细胞组学数据的聚类分析方法,其特征在于,在步骤S1中,所述单细胞基因组学数据为能够用于聚类分析的单细胞数据,单细胞数据采用软件与人类参考基因组进行比对,数据质控后获得相应的可直接用于聚类分析的单细胞基因组学数据;单细胞数据包括:单细胞转录组测序scRNA-Seq数据、单细胞核转录组测序snRNA-Seq数据、单细胞ATAC测序sc-ATAC-Seq数据、单细胞Hi-C测序数据;将单细胞基因组学数据采用相应的软件进行标准化预处理。
3.如权利要求1所述的单细胞组学数据的聚类分析方法,其特征在于,在步骤S2中,借助可视化工具显示不同参数下的聚类结果之间细胞的流向关系,并追踪每一个细胞在cluster数目逐步增加的过程中流向或定位,初步探究单细胞数据在聚类分析中体现的数据结构特征;在随聚类分辨率参数变化的过程中,部分细胞出现归属跳跃及分群现象,被认为是一种相对不稳定的聚类状态,聚类结果在随分辨率参数增加后,其细胞归属保持不变,对应的聚类结果最为稳定,从而获得聚类结果。
4.如权利要求3所述的单细胞组学数据的聚类分析方法,其特征在于,单细胞数据的整体结构的稳定性及其变化,可以利用聚类结果的无间断表征,通过构建连续参数的聚类结果稳定性评分系统,将相邻划分单位下的聚类稳定性进行关联,评估相邻分辨率参数下聚类结果的一致性得分。
5.如权利要求1所述的单细胞组学数据的聚类分析方法,其特征在于,在步骤S3中,根据不同层级下每一个细胞子群中细胞的聚类流向定义罚函数,根据罚函数量化不同层级下聚类结果的稳定性特征。
6.如权利要求5所述的单细胞组学数据的聚类分析方法,其特征在于,包括选自下组的一个或多个特征:
(1)不同层级下每一个细胞子群中细胞的聚类流向包含:经聚类干扰条件保持不变的细胞子群、出现独立分支的细胞子群以及出现汇入行为的细胞子群;其中,细胞归属保持不变的细胞子群表征了该单细胞数据稳定的细胞类群,计作V,出现分支及出现汇入行为的细胞子群表征了该单细胞数据中不稳定的细胞对数目,分别计作C1、C2,根据公式c1计算出相邻划分单位下的聚类结果不一致的细胞对数目UI,公式c1为:
(2)根据标准化原理,构建罚函数F,计算相邻参数下聚类结果的稳定性,罚函数F的计算公式c2为:
其中,UI为相邻划分单位下的聚类结果不一致的细胞对数目,E(UI)为相邻划分单位下的聚类结果不一致细胞对数目的期望值,max(UI)为随机模型下相邻划分单位聚类结果不一致的理论最大值;
(3)使用罚函数F计算得到各层级下聚类结果的对应的稳定性得分scScore(i),计算公式为C3:
scScore(i)=1-F(i)
其中,F(i)为对应干扰条件下稳定性罚分,scScore(i)为对应干扰单位下的聚类结果的稳定性得分;
(4)根据各干扰单位下的稳定性得分,构建不同层级下聚类结果稳定性得分的变化体系;以整体得分的平均值作为稳定性阈值,过滤掉稳定性低于稳定性阈值的聚类结果,得到稳定性高的聚类结果。
7.如权利要求1所述的单细胞组学数据的聚类分析方法,其特征在于,在步骤S4中,包括步骤:
(1)根据S3获得的稳定性高的聚类结果的稳定性得分变化趋势,确定出现连续稳定性得分的变化区间;
(2)统计连续稳定性得分的变化区间聚类数目的变化情况,选择稳定区间内具有一定一致性的聚类数目作为待选最适聚类数目;
(3)比较待选最适数目的稳定性得分,选择前两个稳定性区间中稳定性得分最高的聚类数目,作为最终的最适聚类数目。
8.如权利要求7所述的单细胞组学数据的聚类分析方法,其特征在于,所述待选最适聚类数目为2-3个,通常第一个稳定点下的聚类结果对应了该数据的基础细胞类型,第二个稳定点则是包含了亚细胞类型或亚细胞状态的聚类结果。
CN202211396624.0A 2022-11-09 2022-11-09 一种单细胞组学数据的聚类分析方法 Active CN115527610B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211396624.0A CN115527610B (zh) 2022-11-09 2022-11-09 一种单细胞组学数据的聚类分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211396624.0A CN115527610B (zh) 2022-11-09 2022-11-09 一种单细胞组学数据的聚类分析方法

Publications (2)

Publication Number Publication Date
CN115527610A CN115527610A (zh) 2022-12-27
CN115527610B true CN115527610B (zh) 2023-11-24

Family

ID=84705119

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211396624.0A Active CN115527610B (zh) 2022-11-09 2022-11-09 一种单细胞组学数据的聚类分析方法

Country Status (1)

Country Link
CN (1) CN115527610B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117854600B (zh) * 2024-03-07 2024-05-21 北京大学 基于多组学数据的细胞识别方法、装置、设备及存储介质

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105868791A (zh) * 2016-04-15 2016-08-17 上海交通大学 基于模糊聚类的多分辨率社区发现方法
CN107545133A (zh) * 2017-07-20 2018-01-05 陆维嘉 一种用于鉴别诊断慢性支气管炎的高斯模糊聚类计算方法
CN109979538A (zh) * 2019-03-28 2019-07-05 广州基迪奥生物科技有限公司 一种基于10x单细胞转录组测序数据的分析方法
CN110827921A (zh) * 2019-11-12 2020-02-21 玉林师范学院 一种单细胞聚类方法、装置、电子设备及存储介质
CN113808669A (zh) * 2021-09-28 2021-12-17 上海大学 一种宏基因组序列组装方法
CN113889192A (zh) * 2021-09-29 2022-01-04 西安热工研究院有限公司 一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法
WO2022058339A1 (en) * 2020-09-15 2022-03-24 Dhapola Parashar Method and system for subsampling of cells from single-cell genomics dataset
CN114864003A (zh) * 2022-03-17 2022-08-05 中国科学院深圳先进技术研究院 基于混合实验组和对照组单细胞样本的差异分析方法及系统
CN115293290A (zh) * 2022-08-29 2022-11-04 重庆理工大学 一种自动识别聚类数的层次聚类算法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105868791A (zh) * 2016-04-15 2016-08-17 上海交通大学 基于模糊聚类的多分辨率社区发现方法
CN107545133A (zh) * 2017-07-20 2018-01-05 陆维嘉 一种用于鉴别诊断慢性支气管炎的高斯模糊聚类计算方法
CN109979538A (zh) * 2019-03-28 2019-07-05 广州基迪奥生物科技有限公司 一种基于10x单细胞转录组测序数据的分析方法
CN110827921A (zh) * 2019-11-12 2020-02-21 玉林师范学院 一种单细胞聚类方法、装置、电子设备及存储介质
WO2022058339A1 (en) * 2020-09-15 2022-03-24 Dhapola Parashar Method and system for subsampling of cells from single-cell genomics dataset
CN113808669A (zh) * 2021-09-28 2021-12-17 上海大学 一种宏基因组序列组装方法
CN113889192A (zh) * 2021-09-29 2022-01-04 西安热工研究院有限公司 一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法
CN114864003A (zh) * 2022-03-17 2022-08-05 中国科学院深圳先进技术研究院 基于混合实验组和对照组单细胞样本的差异分析方法及系统
CN115293290A (zh) * 2022-08-29 2022-11-04 重庆理工大学 一种自动识别聚类数的层次聚类算法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Single cell transcriptomics of human epidermis identifies basal stem cell transition states;Shuxiong wang et al.;《Nature communications》;全文 *
基于WGCNA鉴定子宫颈鳞状细胞癌CD8+T细胞浸润相关的生物标志物;凡志祥 等;《中国现代医药杂志》;第24卷(第07期);全文 *
基于异构并行计算的单细胞测序数据聚类算法;谢林娟 等;《计算机工程与应用》;第58卷(第24期);全文 *

Also Published As

Publication number Publication date
CN115527610A (zh) 2022-12-27

Similar Documents

Publication Publication Date Title
Wirth et al. Mining SOM expression portraits: feature selection and integrating concepts of molecular function
US11954614B2 (en) Systems and methods for visualizing a pattern in a dataset
Kairov et al. Determining the optimal number of independent components for reproducible transcriptomic data analysis
Hanczar et al. Small-sample precision of ROC-related estimates
CN110827924B (zh) 基因表达数据的聚类方法、装置、计算机设备及存储介质
CN109891508A (zh) 单细胞类型检测方法、装置、设备和存储介质
Yun et al. Biclustering for the comprehensive search of correlated gene expression patterns using clustered seed expansion
US20130304783A1 (en) Computer-implemented method for analyzing multivariate data
CN112599199A (zh) 一种适用于10x单细胞转录组测序数据的分析方法
CN115527610B (zh) 一种单细胞组学数据的聚类分析方法
CN114864003A (zh) 基于混合实验组和对照组单细胞样本的差异分析方法及系统
Liu et al. Are dropout imputation methods for scRNA-seq effective for scATAC-seq data?
US20230352119A1 (en) Method and system for subsampling of cells from single-cell genomics dataset
Zou et al. An improved fast shapelet selection algorithm and its application to pervasive EEG
Pandey et al. Improved downstream functional analysis of single-cell RNA-sequence data using DGAN
Li et al. Simultaneous estimation of cluster number and feature sparsity in high-dimensional cluster analysis
Rossin et al. A framework for analytical characterization of monoclonal antibodies based on reactivity profiles in different tissues
CN103678709B (zh) 一种基于时序数据的推荐系统攻击检测方法
Alexander et al. Capturing discrete latent structures: choose LDs over PCs
CN109885685A (zh) 情报数据处理的方法、装置、设备及存储介质
Dash et al. Distance based feature selection for clustering microarray data
Giurcăneanu et al. Cluster structure inference based on clustering stability with applications to microarray data analysis
Kalinin et al. A versatile information retrieval framework for evaluating profile strength and similarity
López-Oriona et al. Analyzing categorical time series with the R package ctsfeatures
CN117992858B (zh) 一种免疫细胞亚群识别方法、系统、电子设备及存储介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant