技术领域
[0001] 本申请涉及人工智能技术领域,具体涉及一种基因共表达网络建立方法以及相关装置。
相关背景技术
[0002] 在生物领域,基因共表达分析对于理解复杂的生物过程至关重要,其中,单细胞测序技术能够在单细胞分辨率下分析基因表达,揭示细胞之间的异质性。然而由于单细胞测序数据的高维度、高噪声以及高稀疏性,增加了网络构建的复杂性,使得研究者们从单细胞数据推断基因表达面临着巨大挑战。
具体实施方式
[0030] 下面将结合本申请实施例中的附图,对本申请实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅是本申请一部分实施例,而不是全部的实施例。基于本申请中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本申请保护的范围。
[0031] 如前述,基因共表达分析对于理解复杂生物过程至关重要。比如,在微阵列和常规转录组测序(也即RNA‑seq)数据的研究中,基因共表达分析在基因模块检测、疾病发现和生物相关基因的检测方面表现卓越性能。尤其是,单细胞RNA测序(也即scRNA‑seq)技术的发展,使得能够在单个细胞分辨率下分析基因表达,揭示细胞之间的异质性。目前,多种方法或指标已经被提出用于从scRNA‑seq数据中推断基因共表达关系,包括但不限于皮尔逊(pearson)相关系数,斯皮尔曼(spearman)相关系数,GENIE3,PIDC,minet,sclink和glasso。当从scRNA‑seq数据推断基因共表达时,前述这些指标或方法考虑了不同的关联度量。比如,现在流行的加权基因共表达网络分析(WGCNA)就是采用pearson作为共表达相关性指标,该方法主要应用于常规转录组基因共表达网络的构建。PIDC和Minet利用互信息(MI)作为共表达测量方法。GENIE3基于树的集成方法来估计基因之间的相关性。sclink通过在单细胞数据中明确建模,修正了高斯图形模型,从而进行基因共表达的估计。然而,由于scRNA‑seq数据的高维度和高稀疏性,前述的指标或方法在进行基因共表达网络的构建中有较低的准确性。此外,前述的指标或方法在处理多来源、异质性高的大规模转录组整合数据时,效果较差,无法识别到相关性较高的基因模块。并且,对于从大型数据集中推断细胞类型特异性共表达网络的方法,如sclink,CS‑CORE,直接筛选细胞类型高表达的基因用于构建网络,但是这些方法会遗漏掉一些细胞类型低表达的基因,使得最终构建的网络不准确。因此,迫切需要开发新的计算方法,从这些大型scRNA‑seq数据中有效提取和整合基因共表达信息,以提高基因共表达预测的准确性。并改进目前基于大样本中构建特定细胞类型的基因共表达网络方式,推断更准确的细胞类型特异性网络,揭示其潜在的基因共表达模式。
[0032] 为此,本申请实施例提供一种共表达网络建立方法,通过获取目标基因表达量矩阵,并且对该目标基因表达量矩阵进行特征提取,获得第一基因特征矩阵,然后,利用该第一基因特征矩阵构建全局的基因共表达网络,和/或,利用该第一基因特征矩阵和该目标基因表达量矩阵构建对应的细胞类型特异性共表达网络。以此,使用基因特征来构建全局的基因共表达网络和/或细胞类型特异性共表达网络,可以避免直接基于稀疏的单细胞测序数据预测单细胞基因共表达导致假阳性偏高的问题,不仅提高了单细胞基因共表达预测的准确性,而且能够预测更显著相关的细胞类型特异性基因模块。
[0033] 下面结合附图详细阐述本申请描述的实施例。
[0034] 如图1所示,其示出本申请实施例提供的一种基因共表达网络建立方法的流程示意图。在图1中,该基因共表达网络建立方法,可以包括:步骤101:获取目标基因表达量矩阵;所述目标基因表达量矩阵的元素用于衡量某一基因在某一单细胞内的表达水平。
[0035] 需要说明的是,所谓基因表达(gene expression)是指将来自基因的遗传信息合成功能性基因产物的过程,其中,基因表达产物通常是蛋白质,所有已知的生命均是利用基因表达来合成生命的大分子。基因表达过程可以包括以下三个主要步骤:转录、RNA加工以及翻译。所谓的基因表达量矩阵可以是生物信息学和转录学研究中的核心数据结构,它可以通过高通量测序技术如可以是指scRNA‑seq数据集生成。这种基因表达量矩阵的行代表基因,列代表细胞,其中,该基因表达量矩阵中每一个元素(和/或单元格)包含了响应基因在对应的细胞中的表达水平。一个基因表达量矩阵实例如图2所示。
[0036] 在一些实施例中,对于步骤101,可以包括:获取至少一个生物样本的scRNA‑seq数据集,并且根据每一个scRNA‑seq数据集生成对应生物样本包含基因的子表达量矩阵;对至少一个所述子表达量矩阵进行预处理,获得所述目标基因表达量矩阵。
[0037] 需要说明的是,本申请所说的目标基因表达量矩阵可以是从至少一个生物样本中提取的多个基因表达组成的矩阵。换句话说,该目标基因表达量矩阵中的基因或细胞可以来自至少一个生物样本。所说的生物样本可以是任何具有遗传基因的生物,比如,小鼠、人的某些组织(比如,皮层组织、外周血单核细胞等)。其中,每一样生物样本经scRNA‑seq之后,会获得一个scRNA‑seq数据集;每一个scRNA‑seq数据集会生成一个对应子表达量矩阵。
[0038] 在获得每一个子表达量矩阵后,对该至少一个子表达量矩阵进行预处理,获得该目标基因表达量矩阵。其中,所述预处理可以包括以下至少之一:质量控制、筛选高变异基因、数据归一化。
[0039] 需要说明的是,所谓质量控制可以是指由于经scRNA‑seq之后得到的scRNA‑seq数据集的测序质量会严重影响后续的分析,因此,可以对得到的scRNA‑seq数据集进行质量控制,以过滤不满足条件的基因或细胞。例如,在一些实施例中,该质量控制的目的主要是去除低质量的细胞和基因,因此,可以采用python包中的scanpy对得到的每一scRNA‑seq数据集进行质量控制。其中,scanpy是一个专门用于处理单细胞数据的python包。也即,scanpy包的filter_cells函数和filter_genes函数去除每一scRNA‑seq数据集中表达量低的细胞和基因,得到过滤后的表达量数据。
[0040] 这里,所谓的筛选高可变基因可以是指从经前述的质量控制后的scRNA‑seq数据集中选择在不同细胞间表达差异较大的基因进行研究。例如,在本申请中,可以经过质量控制的scRNA‑seq数据集中筛选出1000个高可变基因进行分析。
[0041] 这里,所谓的归一化可以是指为了避免神经网络训练出现梯度消失或爆炸,可以对筛选高可变基因的数据进行归一化操作,以将筛选高可变基因的数据限制在一个合适的范围。例如,可以采用离差标准化(min‑max normalization)方法,将筛选高可变基因的数据映射到0‑1之间,这样,可以加速基因聚类模型的训练速度,提高训练效率。
[0042] 经上述预处理后,获得后续基因表达量矩阵。如前述图2所示为一个示例性获得的基因表达量矩阵,其中,每一列表示一个细胞对应的所有基因,每一行表示一个基因对应的所有细胞,例如,图2中的基因表达量数据中总共包括n个细胞,每个细胞包括m个基因,n和m均为正整数。每一个元素表征的是一个基因在一个细胞中表达程度,每一个元素的数值在0‑1之间。
[0043] 步骤102:对所述目标基因表达量矩阵进行特征提取,获得第一基因特征矩阵;所述第一基因特征矩阵中的元素用于衡量所述目标基因表达量矩阵中的基因在不同维度的特征值。
[0044] 在一些实施例中,对于步骤102,可以包括:利用训练好的基因聚类模型包含的特征提取单元对所述目标基因表达量矩阵进行特征提取,获得所述第一基因特征矩阵;其中,所述特征提取单元包括输入层、卷积层、激活函数层、池化层及全连接层,所述池化层为空间金字塔池化层,所述池化层的输出结果作为全连接层的输入;所述卷积层包括N个,并且,第i个卷积层的输出结果为第i+1个卷积层的输入,第N个卷积层的输出结果经激活函数层处理后为池化层的输入,其中,N和i为正整数,i小于N。
[0045] 也就是说,该特征提取单元可以由深度神经网络实现的。一种示例性实施方案是,该特征提取单元可以是卷积神经网络(CNN,Convolutional Neural Network),比如,该特征提取单元可以是AlexNet网络。
[0046] 示例性地,如图3所示,该特征提取单元可以包括输入层(input)、卷积层(Cnov)、激活函数层(Relu)、池化层(Pool)及全连接层(FC)。其中,示例性地,该特征提取单元(比如n特征提取单元中的卷积层)可以满足以下一项或多项:卷积核的个数2 、卷积核的大小为3*
1、步长(也即卷积的步长)设置为1、填充设置为1;其中,n为正整数。这里,卷积核(也可以称之为滤波器或特征检测器),其是一个小型矩阵,它在输入数据(如图像)上滑动(或卷积),并且在每个位置计算卷积核与输入数据的元素乘积之和,这个过程会产生一个新的二维数组,通常称之为特征(feature map)或激活图(activation map),它代表了输入数据中某些特征的强度和位置。其中,该步长可以是指卷积核在输入数据上移动的间隔,也就是卷积核矩阵本身。填充可以是指在输入数据的边缘上添加额外的数值,使得输入和输出大小保持一致。需要说明的是,在本申请实施例中,也可以根据实际需要更改特征提取单元中的卷积核个数、卷积核的大小、卷积的步长、填充、卷积层数等。
[0047] 在一些实施例中,该特征提取单元中的池化层可以为空间金字塔池化层(SPPNet,Spatial Pyramid Pooling Net)。示例性地,池化层的输出结果可作为全连接层的输入。这样,通过将不同尺寸的基因表达量矩阵转换成固定大小的特征向量,可以在不改变基因聚类模型的情况下处理不同尺寸的数据,从而可以提高基因聚类模型的适应性与泛化性。也就是说,该已训练好的基因聚类模型可以适用于不同大小的基因表达量矩阵。
[0048] 示例性地,可以将SPPNet的输出结果的维度设置为44维,使得不同尺寸的基因表达量矩阵经SPPNet处理后的输出结果的维度均为44维,然后将SPPNet的输出结果输入全连接层中进行处理,得到基因特征矩阵,这样,可以提高基因聚类模型在不同样本量下特征提取的能力,从而提高基因聚类模型的适应性与泛化性。
[0049] 在一些实施例中,该特征提取单元中除最后一个卷积层外,其他卷积层的后面可以均没有池化层。例如,该特征提取单元可以包括N个卷积层,第i个卷积层的输出结果可以为第i+1个卷积层的输入,第N个卷积层的输出结果经激活函数层处理后可以为池化层的输入,其中,N和i为正整数,i小于N。这样,可以有效地降低基因聚类模型的计算量。
[0050] 经过上述的已训练好的基因聚类模型中的特征提取单元对该目标基因表达量矩阵进行特征提取,获得第一基因特征矩阵,该第一基因特征矩阵中的元素表征着目标基因表达量矩阵中基因的特定维度的表达强度,也就是,该第一基因特征矩阵中的元素用于衡量该目标基因表达量矩阵中某一基因在不同维度的特征值。
[0051] 在经步骤102得到的第一基因特征矩阵,可以如图4所示。在图4中,每一行是同一基因在不同维度的特征值,每一列是同一维度不同基因的特征值。其中,例如,图4中的第一基因特征矩阵中总共包括p维,也即每一个基因包含p个维度上的特征值,每一维度上有m个基因,每一个基因对应一个特征值。p和m均为正整数。在一些实施例中,p小于n(n为细胞数量)。
[0052] 步骤103:根据所述第一基因特征矩阵和所述目标基因表达量矩阵构建对应的细胞类型特异性共表达网络;和/或,根据所述第一基因特征矩阵构建全局的基因共表达网络。
[0053] 需要说明的是,在得到该第一基因特征矩阵后,可以根据需求获得全局的基因共表达网络和/或获得对应的细胞类型特异性共表达网络。
[0054] 在一些实施例中,对于步骤103中的所述根据所述第一基因特征矩阵和所述目标基因表达量矩阵构建对应的细胞类型特异性的共表达网络,如图5所示,可以包括:步骤501:根据所述目标基因表达量矩阵和所述第一基因特征矩阵确定每一细胞类型的第二基因特征矩阵;所述第二基因特征矩阵中的每一行是所述第一基因特征矩阵中的基因在不同维度的特征值,每一列是不同基因在同一维度的特征值;
步骤502:根据每一个细胞类型的所述第二基因特征矩阵获得对应的细胞类型特异性模块;
步骤503:对每一个细胞类型特异性模块中的基因进行层级聚类,获得对应的细胞类型特异性共表达网络。
[0055] 在一些实施例中,对于步骤501,可以包括:根据所述目标基因表达量矩阵获得每一个细胞类型包含的差异表达基因的基因表达量矩阵;计算每一个细胞类型包含的差异表达基因的基因表达量矩阵与所述第一基因特征矩阵之间的相关性,形成对应细胞类型的细胞‑嵌入基因相关性矩阵;其中,所述细胞‑嵌入基因相关性矩阵中的每一行表征不同特征在所述对应细胞中的重要程度,每一列表征同一特征在所述对应的细胞中的重要程度;对于每一细胞类型,根据对应细胞类型的细胞‑嵌入基因相关性矩阵和所述第一基因特征矩阵确定对应细胞类型的所述第二基因特征矩阵。
[0056] 这里,所说的细胞类型可以是指基于特定的形态、功能、结构和分子标记物来分类的细胞群。在生物学中,细胞类型的多样性是生物体复杂性的基础。其中,按照不用分类方式有不同细胞类型,比如,按照形态结构分类,可以包括上皮细胞群、肌肉细胞群、神经细胞群等等。再比如,按照功能分类,可以包括:光感受细胞群、嗅觉细胞群、味觉细胞群等等。所说的差异表达基因(Differential Gene Expression)可以是指在不同条件下(如不同组织类型、不同发育阶段、不同环境压力或不同病理状态等),基因表达水平发生显著变化的基因,也就是表达差异比较大的基因。在一些实施例中,可以设置阈值来筛选差异表达基因,比如,基因表达相似度小于一定阈值时,定义为差异表达基因;再比如,若基因表达相似度大于一定阈值定义为表达相似基因,其余的为差异表达基因等等。
[0057] 也就是说,对于步骤501,先从该目标基因表达量矩阵中获得每一个细胞类型对应的差异表达基因的基因表达量矩阵;然后,计算每一个细胞类型对应的差异表达基因的基因表达量矩阵与前述计算的第一基因特征矩阵之间的相关性,获得每一个细胞类型对应的细胞‑嵌入基因(cell‑embedding)相关性矩阵。所述细胞‑嵌入基因相关性矩阵中的每一行表征不同特征在所述对应细胞中的重要程度,每一列表征同一特征在所述对应的细胞中的重要程度;之后,对于每一细胞类型,根据对应细胞类型的细胞‑嵌入基因相关性矩阵和所述第一基因特征矩阵确定对应细胞类型的所述第二基因特征矩阵。
[0058] 其中,在一些实施例中,所说的根据所述目标基因表达量矩阵获得每一个细胞类型包含的差异表达基因的基因表达量矩阵可以是指从目标基因表达量矩阵中筛选出每一个细胞类型包含的差异性表达基因的基因表达量矩阵。
[0059] 其中,在一些实施例中,所说的计算每一个细胞类型包含的差异表达基因的基因表达量矩阵与所述第一基因特征矩阵之间的相关性,形成对应细胞类型的细胞‑嵌入基因相关性矩阵可以是指通过细胞类型包含的差异表达基因的基因表达量矩阵与所述第一基因特征矩阵计算对应细胞类型包含的基因对之间的相关系数(如皮尔逊相关系数)来评估它们之间的相关性,然后根据相关性,将各细胞类型包含的差异表达基因的基因表达量矩阵和共表达网络的相关性矩阵结合起来,就可以得到细胞‑嵌入基因相关性矩阵。
[0060] 其中,在一些实施例中,对于每一细胞类型,所述根据对应细胞类型的细胞‑嵌入基因相关性矩阵和所述第一基因特征矩阵确定对应细胞类型的所述第二基因特征矩阵,可以包括:根据所述细胞‑嵌入基因相关性矩阵从所述第一基因特征矩阵中提取设定维数的特征作为对应细胞类型的特征;根据提取的对应细胞类型的特征和所述第一基因特征矩阵中包含的各个基因形成所述第二基因特征矩阵。
[0061] 也就是,在得到细胞‑嵌入基因相关性矩阵后,结合第一基因特征矩阵,来确定第二基因特征矩阵。具体地,就是根据细胞‑嵌入基因相关性矩阵从第一基因特征矩阵中提取设定维数的特征作为对应细胞类型的特征,比如,提取100维作为对应的细胞类型的特征;然后,根据提取的对应细胞类型的特征和第一基因特征矩阵中包含的各个基因形成所述第二基因特征矩阵,比如,第一基因特征矩阵包含的基因是1000个,那么得到第二基因特征矩阵的维数也就是(1000,100),其中,1000为基因个数;100为特征维数。
[0062] 之后,在一些实施例中,对于步骤502,如图6所示,可以包括:步骤601:计算所述第二基因特征矩阵中的各个基因之间的余弦相似性;
步骤602:根据所述各个基因之间的余弦相似性对所述第二基因特征矩阵进行聚类处理,形成第一子细胞类型特异性模块和第二子细胞类型特异性模块;
步骤603:计算所述第一子细胞类型特异性模块的第一得分和所述第二子细胞类型特异性模块的第二得分;
步骤604:比较所述第一得分和第二得分;其中,响应所述第一得分大于所述第二得分,确定所述第一子细胞类型特异性模块为所述对应的细胞类型特异性模块;响应于所述第一得分小于所述第二得分,确定所述第二子细胞类型特异性模块为所述对应的细胞类型特异性模块;响应于所述第一得分等于所述第二得分,确定所述第一子细胞类型特异性模块或第二子细胞类型特异性模块为所述对应的细胞类型特异性模块。
[0063] 需要说明的是,对于步骤502中的前两步:获得第一子细胞类型特异性模块和第二子细胞类型特异性模块的方式与后述的步骤701和702的步骤相似,在此不再赘述详细细节。仅有一点需说明,这里的聚类数设置为2,也即聚类形成第一子细胞类型特异性模块和第二子细胞类型特异性模块。
[0064] 在一些实施例中,对于步骤502中的所述计算所述第一子细胞类型特异性模块的第一得分和所述第二子细胞类型特异性模块的第二得分,可以包括:计算所述第一子细胞类型特异性模块中每一个差异表达基因在不同细胞中的第一相对表达水平;并根据每一个所述第一相对表达水平和所述第一子细胞类型特异性模块中包含的基因数量计算所述第一得分;以及计算所述第二子细胞类型特异性模块中每一个差异表达基因在不同细胞中的第二相对表达水平;并根据每一个所述第二相对表达水平和所述第二子细胞类型特异性模块中包含的基因数量计算所述第二得分。
[0065] 需要说明的是,所谓的细胞类型特异性是指不同细胞类型在基因表达水平上的差异,这些差异反映了不同细胞类型的功能和特性。细胞类型特异性的基因表达模式对于理解细胞发育、组织功能以及疾病机制具有重要意义。
[0066] 本申请实施例定义了一个新的细胞类型特异性得分来筛选细胞类型特异性模块,具体步骤如下:第一步,首先计算基因在不同细胞中的相对表达水平。例如,基因g的相对表达水平由基因g在细胞c的表达水平除以基因g在该细胞类型中的所有细胞中的表达水平总和。
计算公式如下:
(1)
其中, 表示基因g在细胞c中的表达水平,C表示某一细胞类型中包含的细胞总
数。 表示基因g在该细胞c中的相对表达水平。
[0067] 第二步,计算不同子细胞类型特异性模块中的细胞类型特异性得分,也即第一得分或第二得分。具体地,将子细胞类型特异性模块中所有基因的细胞类型特异性得分总和除以子细胞类型特异性模块中基因的数量得到该子细胞类型特异性模块的细胞类型特异性得分。计算公式如下: (2)
其中, 表示基因g在细胞类型t中的表达水平, 也即前面计算的基因在不同
细胞中的相对表达水平之和。∣m∣表示模块m中基因的总数量。 表示子细胞类型特异性模块的细胞类型特异性得分,也即第一得分或第二得分。也就是,该第一得分和第二得分均按照上述两个步骤计算得到。
[0068] 在得到第一子细胞类型特异性模块的第一得分和第二子细胞类型特异性模块的第二得分后,比较该第一得分和该第二得分的大小,在该第一得分大于该第二得分时,确定该第一子细胞类型特异性模块为所述对应的细胞类型特异性模块;若该第一得分小于该第二得分时,确定该第二子细胞类型特异性模块为所述对应的细胞类型特异性模块。若该第一得分等于该第二得分时,确定该第一子细胞类型特异性模块或第二子细胞类型特异性模块为所述对应的细胞类型特异性模块。也即,比较不同子细胞类型特异性模块的细胞类型特异性得分,筛选较高得分的子细胞类型特异性模块作为该种细胞类型的特异性模块。
[0069] 之后,对筛选细胞类型特异性模块进行层级聚类,得到对应细胞类型特异性的共表达网络。其中,该细胞类型特异性的共表达网络是指在特定细胞类型中,基因表达模式相似的基因集形成的网络。这种网络可以帮助理解不同细胞类型在基因表达层面的独特性和功能。
[0070] 在另一些实施例中,所谓的全局的基因共表达网络(Global Gene Co‑expression Network)可以是一种生物信息学分析方法,用于研究基因之间的相互作用和协同表达模式。这种方法可以帮助我们理解基因在不同生物学过程中的功能和调控网络。在全局的基因共表达网络中,每个基因被表示为网络中的一个节点,而节点之间的边(连接线)表示基因之间的共表达关系。这些关系通常通过计算基因表达数据的相关系数来确定,如皮尔逊相关系数或斯皮尔曼相关系数。网络中的边的权重可以反映共表达关系的强度,权重越高,表示两个基因的表达模式越相似。也就是,根据第一基因特征矩阵获得该基因表达量矩阵中各个基因之间的共表达网络。
[0071] 在一些实施例中,对于步骤103中的所述根据所述第一基因特征矩阵构建全局基因共表达网络,如图7所示,可以包括:步骤701:根据所述第一基因特征矩阵计算各个基因之间的共表达关系,得到基因共表达相关矩阵;
步骤702:对所述基因共表达相关矩阵中的基因进行层级聚类构建所述全局基因共表达网络;所述全局基因共表达网络呈树状结构,其中,该树状结构的每一个分支包括一个基因共表达模块;所述基因共表达模块包含表达模式相似的一组基因。
[0072] 也就是,基于该第一基因特征构建该全局基因共表达网络可以包括以下两步:根据前面获得的第一基因特征矩阵计算前述目标基因共表达矩阵中的各个基因之间的共表达关系,得到基因共表达相关矩阵,也即步骤701;之后,对该基因共表达相关矩阵中的基因进行层级聚类得到该全局基因共表达网络,也即步骤702。其中,所谓层级聚类(Hierarchical Clustering)是一种常用的聚类分析方法,它不需要事先指定簇的数量,而是通过逐步合并或分裂的方式构建一个由层次组成的聚类树(称为树状图或dendrogram),也就是说,该全局基因共表达网络呈树状结构,并且,该树状结构的每一个分支包括一个基因共表达模块;该基因共表达模块包括表达模式相似比较高的一组基因。在实现时,可以通过阈值来确定表达模式相似比较高的一组基因。
[0073] 这里,所述第一基因特征矩阵的每一行表征同一基因在不同维度的特征值,每一列表征同一维度下不同基因的特征值;对于步骤701可以包括:计算所述第一基因特征矩阵中基因间的余弦相似性;所述余弦相似性用于衡量基因之间的共表达关系;根据获得的各个基因之间的余弦相似性获得所述基因共表达相关矩阵。
[0074] 也就是,利用余弦相似性计算基因之间的共表达关系,以获得该基因共表达相关矩阵。所谓余弦相似性(Cosine Similarity)是一种衡量两个向量之间角度的相似性度量方法。它通过计算两个向量的点积和它们模的乘积的比值来确定。余弦相似性的范围从‑1到1,其中1表示向量完全相同的方向,‑1表示完全相反的方向,而0表示向量之间是正交的,即无相关性。在本申请中两个向量可以是指该第一基因特征矩阵中的任意两行。
[0075] 其中,余弦相似性的计算公式如下:余弦相似性(Cosine Similarity)=(A⋅B)/(||A||*||B||)。
[0076] 其中:A⋅B 是向量 A 和B 的点积;||A||是向量 A 的模(或范数);||B||是向量 B 的模(或范数)。
[0077] 在一些实施例中,对于步骤602,可以包括:设定最小基因模块中包含的基因数目,按照动态剪切树策略逐步对所述基因共表达相关矩阵中的数据进行处理,形成至少一个基因共表达模块;记录每一个得到的基因共表达模块,以构建所述全局基因共表达网络。
[0078] 这里,所说的设定最小基因模块中包含的基因数目可以是指定基因模块的最小单元,比如说,该设定最小基因模块中包含的基因数目为10,那么,聚类出来的基因模块中至少包含10个基因。该动态剪切树策略是一种在加权基因共表达网络分析(WGCNA)中用于识别基因模块的方式。因此,对于步骤602可以根据动态剪切树策略构建全局的基因共表达网络。
[0079] 本申请实施例提供的基因共表达网络建立方法,通过将一个高维且稀疏的单细胞基因表达量矩阵转换成低维的基因特征矩阵,来构建全局基因共表达网络或细胞类型特异性的共表达网络,实现了矩阵解耦,以避免直接基于稀疏的单细胞数据进行网络构建,提升共表达计算准确性。也即,本申请实施例提供一种解耦特征学习框架(DeepCSCN)来构建细胞类型特异性的基因共表达网络(也即前述的细胞类型特异性的共表达网络)。其中,DeepCSCN通过解耦细胞类型特征(也即分离不同细胞类型相关的特征维度),可以从大规模的样本中推断细胞类型特异性的共表达网络。DeepCSCN具有以下创新和优势:(1)有效的特征提取:DeepCSCN利用深度聚类模型(也即GeneCluster,也即前述的已训练好的基因聚类模型)进行基因特征提取,解决了scRNA‑seq数据的稀疏性问题,提高了基因相关性计算的准确性。(2)更高的准确性:与现有方法相比,DeepCSCN的共表达预测准确性更高,能够识别更多的细胞类型特异性基因模块。(3)全局到局部的网络构建:DeepCSCN先利用所有基因的嵌入特征构建全局的基因共表达网络;然后解耦特定细胞类型的嵌入进行特定细胞类型共表达网络的构建。
[0080] 为了理解本申请,以下提供实例进行说明。
[0081] 其中,基因共表达网络建立方法流程可以如下:(1)数据收集与预处理
在本申请中,我们从Single Cell Protal一共获取了8个单细胞数据集,分别是4个小鼠皮层数据和4个人外周血单核数据(PBMC)。在下面的描述中,我们主要以人的PBMC1 drop‑seq数据集为例,该数据集包括9种细胞类型,共计10802个细胞和22792个基因。将该数据集首先经过python包scanpy来初步过滤基因和细胞,然后筛选1000个高可变基因,并对原始的count矩阵进行归一化,得到基因表达量矩阵。然后将基因表达量矩阵中的数据输入GeneCluster进行训练,获得第一基因特征矩阵。
[0082] (2)模型训练,获取第一基因特征矩阵将预处理后的基因表达量矩阵中的数据输入GeneCluster模型进行特征提取,以获的第一基因特征矩阵,其中,该第一基因特征矩阵的维度为(1000,2048)。
[0083] (3)全局基因共表达网络构建将GeneCluster模型提取的PBMC1 drop‑seq数据集的基因特征矩阵基于余弦相似性计算基因间的相关性矩阵,并对相关性矩阵运行层级聚类,得到PBMC1 drop‑seq数据集全局的基因共表达模块。最终经过实验,该数据集一共识别到25个全局的基因共表达模块。
然后我们将不同的基因模块与细胞类型相关联,发现模块1与CD16+ monocyte、Dendritic cell显著相关,模块2与Megakaryocyte显著相关,模块3与Cytotoxic T cell、Dendritic cell和Natural killer cell显著相关,模块5与Plasmacytoid dendritic cell显著相关,模块8与B cell显著相关。
[0084] (4)细胞类型特异性共表达网络的构建尽管全局的基因共表达网络能够将基因模块与细胞类型相关联,但是无法准确地分配每一个模块与特定细胞类型相关联。PBMC1 drop‑seq数据集的细胞类型特异性网络的构建主要分为三步。第一步,我们首先生成一个相关性矩阵来量化基因表征在识别不同细胞类型的重要性。其通过首先提取该数据集在不同细胞类型的差异表达基因,计算这些差异表达基因的基因表达矩阵与基因特征矩阵的相关性,得到cell‑embedding的相关性矩阵,维度为(10802,2048);第二步,提取特定细胞类型(如B细胞,细胞数量为477)top100维的特征作为B细胞的特征,维度为(477,100)。基于所有基因在特定细胞类型B细胞的基因特征(1000,100)计算余弦相似性,并聚类,利用细胞类型特异性得分筛选B细胞特异性模块,最后我们筛选的B细胞特异性模块共包含182个基因。然后将该模块的基因进行再次聚类(182,100),最终一共识别到4个特异性的B细胞共表达模块。将每个模块运行GO富集分析,发现模块1主要富集在“MHC protein complex binding”,模块2主要富集在“B cell activation”,模块3主要富集在“immunoglobulin receptor binding functions”。将相同的数据集应用于sclink和CS‑CORE以进行比较,在B细胞基因共表达模块的识别上,CS‑CORE一共识别到4个B细胞功能模块,但是只有一个模块与B细胞相关;sclink没有识别到B细胞相关的功能模块。因此,与其它细胞类型特异性网络构建方法相比,DeepCSCN能够识别到更多而且更相关的功能通路。
[0085] (5)不同基因共表达方法的比较我们对不同的基因共表达方法,包括pearson、spearman、GENIE3、minet、PIDC、sclink、glasso和CS‑CORE进行比较。同样以PBMC1 drop‑seq数据集作为示例,针对不同共表达方法,我们利用AUPRC和AUROC作为评估指标。AUROC是假阳性率(FDR)与真阳性率(TPR)的曲线图,AUPRC是召回率和精准率围成的曲线下面积。我们首先筛选400个高变异基因来进行比较,计算不同共表达方法的基因相关性矩阵,基于此矩阵计算不同方法的AUPRC和AUROC得分,发现在PBMC1 drop‑seq数据集中,pearson、spearman、GENIE3、minet、PIDC、sclink、glasso,CS‑CORE和DeepCSCN的AUPRC的得分分别是0.50,0.48,0.40,0.07,0.48,
0.51,0.18,0.20,0.60。它们的AUROC的得分分别是0.89,0.89,0.74,0.30,0.88,0.89,
0.54,059,0.91。因此,当DeepCSCN与其它方法相比时,DeepCSCN都有最高的AUPRC和AUROC得分,有效提升了共表达预测的准确性。
[0086] 如图8所示,本申请实施例提供一种基因共表达网络建立装置的结构示意图。在图8中,该共表达网络建立装置800,可以包括:
获取模块801,用于获取目标基因表达量矩阵;所述目标基因表达量矩阵中的元素用于衡量某一基因在某一细胞内的表达水平;
提取模块802,用于对所述目标基因表达量矩阵进行特征提取,获得第一基因特征矩阵;所述第一基因特征矩阵中的元素用于衡量所述目标基因表达量矩阵中的基因在不同维度的特征值;
构建模块803,用于根据所述第一基因特征矩阵和所述目标基因表达量矩阵构建对应的细胞类型特异性共表达网络;和/或,根据所述第一基因特征矩阵构建全局的基因共表达网络。
[0087] 在一些实施例中,所述构建模块803,还用于:根据所述目标基因表达量矩阵和所述第一基因特征矩阵确定每一细胞类型的第二基因特征矩阵;所述第二基因特征矩阵中的每一行是所述第一基因特征矩阵中的基因在不同维度的特征值,每一列是不同基因在同一维度的特征值;根据每一个细胞类型的所述第二基因特征矩阵获得对应的细胞类型特异性模块;对每一个细胞类型特异性模块中的基因进行层级聚类,获得对应的细胞类型特异性共表达网络。
[0088] 在一些实施例中,所述构建模块803,还用于:根据所述目标基因表达量矩阵获得每一个细胞类型包含的差异表达基因的基因表达量矩阵;计算每一个细胞类型包含的差异表达基因的基因表达量矩阵与所述第一基因特征矩阵之间的相关性,形成对应细胞类型的细胞‑嵌入基因相关性矩阵;其中,所述细胞‑嵌入基因相关性矩阵中的每一行表征不同特征在所述对应细胞中的重要程度,每一列表征同一特征在所述对应的细胞中的重要程度;对于每一细胞类型,根据对应细胞类型的细胞‑嵌入基因相关性矩阵和所述第一基因特征矩阵确定对应细胞类型的所述第二基因特征矩阵。
[0089] 在一些实施例中,所述构建模块803,还用于:根据所述细胞‑嵌入基因相关性矩阵从所述第一基因特征矩阵中提取设定维数的特征作为对应细胞类型的特征;根据提取的对应细胞类型的特征和所述第一基因特征矩阵中包含的各个基因形成所述第二基因特征矩阵。
[0090] 在一些实施例中,所述构建模块803,还用于:计算所述第二基因特征矩阵中各个基因之间的余弦相似性;根据所述各个基因之间的余弦相似性对所述第二基因特征矩阵进行聚类处理,形成第一子细胞类型特异性模块和第二子细胞类型特异性模块;计算所述第一子细胞类型特异性模块的第一得分和所述第二子细胞类型特异性模块的第二得分;比较所述第一得分和第二得分;其中,响应所述第一得分大于所述第二得分,确定所述第一子细胞类型特异性模块为所述对应的细胞类型特异性模块;响应于所述第一得分小于所述第二得分,确定所述第二子细胞类型特异性模块为所述对应的细胞类型特异性模块;响应于所述第一得分等于所述第二得分,确定所述第一子细胞类型特异性模块或第二子细胞类型特异性模块为所述对应的细胞类型特异性模块。
[0091] 在一些实施例中,所述构建模块803,还用于:计算所述第一子细胞类型特异性模块中每一个基因在不同细胞中的第一相对表达水平;并根据每一个所述第一相对表达水平和所述第一子细胞类型特异性模块中包含的基因数量计算所述第一得分;以及计算所述第二子细胞类型特异性模块中每一个基因在不同细胞中的第二相对表达水平;并根据每一个所述第二相对表达水平和所述第二子细胞类型特异性模块中包含的基因数量计算所述第二得分。
[0092] 在一些实施例中,所述构建模块803,还用于:根据所述第一基因特征矩阵计算各个基因之间的共表达关系,得到基因共表达相关矩阵;对所述基因共表达相关矩阵中的基因进行层级聚类构建所述全局基因共表达网络;所述全局基因共表达网络呈树状结构,其中,该树状结构的每一个分支包括一个基因共表达模块;所述基因共表达模块包含表达模式相似的一组基因。
[0093] 在一些实施例中,所述第一基因特征矩阵的每一行表征同一基因在不同维度的特征值,每一列表征同一维度下不同基因的特征值;所述构建模块803,还用于:计算所述第一基因特征矩阵中基因间之间的余弦相似性;所述余弦相似性用于衡量基因之间的共表达关系;根据获得的各个基因之间的余弦相似性获得所述基因共表达相关矩阵。
[0094] 在一些实施例中,所述构建模块803,还用于:设定最小基因模块中包含的基因数目,按照动态剪切树策略逐步对所述基因共表达相关矩阵中的数据进行处理,形成至少一个基因共表达模块;记录每一个得到的基因共表达模块,以构建所述全局基因共表达网络。
[0095] 在一些实施例中,所述获取模块801,具体用于:获取至少一个生物样本的单细胞转录组测序(scRNA‑seq)数据集,并且根据每一个scRNA‑seq数据集生成对应生物样本包含基因的子表达量矩阵;对至少一个所述子表达量矩阵进行预处理,获得所述目标基因表达量矩阵。
[0096] 在一些实施例中,所述提取模块802,具体用于:利用训练好的基因聚类模型包含的特征提取单元对所述目标基因表达量矩阵进行特征提取,获得所述第一基因特征矩阵;其中,所述特征提取单元包括输入层、卷积层、激活函数层、池化层及全连接层,所述池化层为空间金字塔池化层,所述池化层的输出结果作为全连接层的输入;所述卷积层包括N个,并且,第i个卷积层的输出结果为第i+1个卷积层的输入,第N个卷积层的输出结果经激活函数层处理后为池化层的输入,其中,N和i为正整数,i小于N。
[0097] 需要说明的是,本申请实施例提供的基因共表达网络建立装置用于实现前述的共表达网络建立方法的装置,因此,这里出现的技术特征在上述描述共表达网络建立方法时已经详细说明,可参考前面描述进行理解,为了节省篇幅,在此不再意义赘述。
[0098] 本申请一方面提供了一种基因共表达网络建立设备,包括处理器,用于存储器中调用程序,以使所述设备执行如前述任一项所述的方法。
[0099] 本申请一方面提供了一种芯片,包括处理器,用于从存储器调用程序,使得安装有所述芯片的设备执行如前述任一项所述的方法。
[0100] 本申请一方面提供了一种计算机可读存储介质,其上存储有程序,所述程序是的计算积极执行如前述任一项所述的方法。
[0101] 图9是本申请一个实施例的设备900的示意性框图。图9所示的设备900包括存储器901、处理器902、通信接口903以及总线904。其中,存储器901、处理器902、通信接口903通过总线904实现彼此之间的通信连接。
[0102] 存储器901可以是图形处理单元(GPU,Graphics Processing Unit)存储体系,只读存储器(ROM,Read Only Memory),静态存储设备,动态存储设备和/或随机存取存储器(RAM,Random Access Memory)。存储器901可以存储程序,当存储器901中存储的程序被处理器902执行时,处理器902用于执行本申请实施例的方法的各个步骤,例如,可以执行前述图1、图5、图6、图7等所示实施例的各个步骤。
[0103] 处理器902可以采用图形处理单元(GPU,Graphics Processing Unit) ,神经网络处理单元(NPU,Neural Network Processing Unit),微处理器,应 用专用集成电路(ASIC,Application Specific Integrated Circuit),和/或一个或多个集成电路,用于执行相关程序,以实现本申请方法实施例的方法。
[0104] 处理器902还可以是一种集成电路芯片,具有信号的处理能力。在实现过程中,本申请实施例的方法的各个步骤可以通过处理器902中的硬件的集成逻辑电路和/或软件形式的指令完成。
[0105] 上述处理器902还可以是通用处理器、数字信号处理器(DSP,Digital Signal Processing)、专用集成电路、现场可编程门阵列(FPGA,Field Programmable Gate array)和/或其他可编程逻辑器件、分立门和/或晶体管逻辑器件、分立硬件组件。可以实现和/或执行本申请实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器和/或该处理器也可以是任何常规的处理器等。
[0106] 结合本申请实施例所公开的方法的步骤可以直接体现为硬件译码处理器执行完成,和/或用译码处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器,闪存、只读存储器,可编程只读存储器和/或电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器901,处理器902读取存储器901中的信息,结合其硬件完成本申请实施例中的各装置包括的单元所需执行的功能 ,和/或,执行本申请中的各方法实施例的方法,例如,可以执行图1、图5、图6及图7等所示实施例的各个步骤/功能。
[0107] 通信接口903可以使用但不限于收发器一类的收发装置,来实现设备900与其他设备或通信网络之间的通信。
[0108] 总线904可以包括在设备900各个部件(例如,存储器901、处理器902、通信接口903)之间传送信息的通路。
[0109] 应理解,本申请实施例所示的设备900可以是处理器或芯片,以用于执行本申请实施例所述的方法。
[0110] 应理解,本申请实施例中,该处理器可以为图形处理单元(GPU,Graphics Processing Unit),该处理器还可以是、数字信号处理器(DSP)、专用集成电路(ASIC)、现场可编程门阵列(FPGA)和/或其他可编程逻辑器件、分立元器件门电路和/或晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器和/或该处理器也可以是任何常规的处理器等。
[0111] 应理解,在本申请实施例中,“与A相应的B”表示B与A相关联,根据A可以确定B。但还应理解,根据A确定B并不意味着仅仅根据A确定B,还可以根据A和/或其它信息确定B。
[0112] 应理解,本文中术语“和/或”,仅仅是一种描述关联对象的关联关系,表示可以存在三种关系,例如,A和/或B,可以表示:单独存在A,同时存在A和B,单独存在B这三种情况。另外,本文中字符“/”,一般表示前后关联对象是一种“或”的关系。
[0113] 应理解,在本申请的各种实施例中,上述各过程的序号的大小并不意味着执行顺序的先后,各过程的执行顺序应以其功能和内在逻辑确定,而不应对本申请实施例的实施过程构成任何限定。
[0114] 在本申请所提供的几个实施例中,应该理解到,所揭露的系统、装置和方法,可以通过其它的方式实现。例如,以上所描述的装置实施例仅仅是示意性的,例如,所述单元的划分,仅仅为一种逻辑功能划分,实际实现时可以有另外的划分方式,例如多个单元或组件可以结合和/或可以集成到另一个系统,或一些特征可以忽略,或不执行。另一点,所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口,装置或单元的间接耦合或通信连接,可以是电性,机械或其它的形式。
[0115] 所述作为分离部件说明的单元可以是和/或也可以不是物理上分开的,作为单元显示的部件可以是和/或也可以不是物理单元,即可以位于一个地方,和/或也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分和/或全部单元来实现本实施例方案的目的。
[0116] 另外,在本申请各个实施例中的各功能单元可以集成在一个处理单元中,也可以是各个单元单独物理存在,也可以两个或两个以上单元集成在一个单元中。
[0117] 在上述实施例中,可以全部或部分地通过软件、硬件、固件和/或其任意组合来实现。当使用软件实现时,可以全部或部分地以计算机程序产品的形式实现。所述计算机程序产品包括一个或多个计算机指令。在计算机上加载和执行所述计算机程序指令时,全部或部分地产生按照本申请实施例所述的流程或功能。所述计算机可以是通用计算机、专用计算机、计算机网络、和/或其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中,和/或从一个计算机可读存储介质向另一个计算机可读存储介质传输,例如,所述计算机指令可以从一个网站站点、计算机、服务器或数据中心通过有线(例如同轴电缆、光纤、数字用户线(DSL,Digital Subscriber Line)或无线(例如红外、无线、微波等)方式向另一个网站站点、计算机、服务器或数据中心进行传输。所述计算机可读存储介质可以是计算机能够读取的任何可用介质和/或是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介质,(例如,软盘、硬盘、磁带)、光介质(例如,数字通用光盘(DVD,Digital Video Disk)和/或半导体介质(例如,固态硬盘(SSD,Solid State Disk)等。
[0118] 以上所述仅为本申请的较佳实施例而已,并不用以限制本申请,凡在本申请的精神和原则之内,所作的任何修改、等同替换等,均应包含在本申请的保护范围之内。