CN116486902A - 一种基于基因调控网络识别驱动调控因子的方法 - Google Patents

一种基于基因调控网络识别驱动调控因子的方法 Download PDF

Info

Publication number
CN116486902A
CN116486902A CN202310524599.8A CN202310524599A CN116486902A CN 116486902 A CN116486902 A CN 116486902A CN 202310524599 A CN202310524599 A CN 202310524599A CN 116486902 A CN116486902 A CN 116486902A
Authority
CN
China
Prior art keywords
gene
network
genes
node
driving
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.)
Pending
Application number
CN202310524599.8A
Other languages
English (en)
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.)
Tsinghua University
Original Assignee
Tsinghua 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 Tsinghua University filed Critical Tsinghua University
Priority to CN202310524599.8A priority Critical patent/CN116486902A/zh
Publication of CN116486902A publication Critical patent/CN116486902A/zh
Pending legal-status Critical Current

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
    • G16B20/00ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • 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
    • G16B40/00ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • 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

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Medical Informatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Bioethics (AREA)
  • Data Mining & Analysis (AREA)
  • Chemical & Material Sciences (AREA)
  • Molecular Biology (AREA)
  • Genetics & Genomics (AREA)
  • Artificial Intelligence (AREA)
  • Analytical Chemistry (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Databases & Information Systems (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Public Health (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及生物信息学领域,具体涉及一种基于基因调控网络识别驱动调控因子的方法。该方法针对单细胞转录组数据基于图注意力神经网络构建基因调控网络,利用该基因调控网络识别驱动基因,通过计算候选驱动基因的影响力得分,识别驱动调控因子以及由驱动调控因子控制的调控基因模块。本发明提供的识别方法准确性高,有利于进行与细胞命运决定相关的基因调控,对于细胞的分化、发育,疾病机制的研究等具有重要的理论意义和实用价值。

Description

一种基于基因调控网络识别驱动调控因子的方法
技术领域
本发明涉及生物信息学领域,具体涉及一种基于基因调控网络识别驱动调控因子的方法。
背景技术
细胞命运决定的机制关系到人体发育、维持稳态、癌症等疾病的发生发展。细胞通过分化产生功能各异的细胞类型,如骨髓造血干细胞受到激素刺激后可分化为淋巴干细胞,并进一步选择性分化为T淋巴细胞和B淋巴细胞。此外,细胞中基因变异等众多复杂因素,可导致正常细胞转变为癌细胞,使细胞命运发生根本性的变化。细胞的分化和重编程等过程是由复杂的基因调节控制的,每个细胞通过整合大量的信号并且执行复杂的基因调控变化来决定自己的命运。究竟什么因素决定了细胞命运,一直是生命科学领域一个基础且非常重要的问题。
从基因到细胞层面阐明细胞命运决定机制至关重要,随着单细胞测序技术的快速发展,越来越多的研究已经从细胞“群体”层面具体到单个细胞水平上,使得生命科学领域的研究更加精准化,然而通过计算方法从单细胞测序数据中挖掘基因调控关系,并检测决定细胞命运的关键调控因子,仍然是一个巨大的挑战。
大量研究发现,细胞命运选择是由一些关键的转录因子调控的,但是这些转录因子如何决定细胞的分化进程还并不完全清楚。很多计算研究通过重建转录因子驱动的基因调控网络来描述细胞命运的转变,这些最关键的转录因子通常被称为主要调控子(Masterregulators,MRs),主要调控子被认为负责控制相关表型下的细胞的整个调控程序。目前常用的主要调控子识别方法包括VIPER,ANANSE,SCENIC等方法,其中只有SCENIC是专门针对单细胞转录组数据设计的。然而,由于单细胞测序数据存在噪声高、缺失值较多等问题,当前基于单细胞转录组测序数据的基因调控网络的构建还不理想,并且目前缺乏专门针对单细胞数据的控制细胞命运决定的主要调控因子识别方法。一般来说,这些主要调控因子仅仅被局限为转录因子,且很难明确哪些调控因子是真正起到驱动作用的,因而迫切需要基于单细胞转录组测序数据的基因调控网络的构建以及对驱动调控因子准确高效地识别。
发明内容
本发明旨在至少在一定程度上解决相关技术中的技术问题之一。为此,本发明的一个目的在于提供一种基于基因调控网络识别驱动调控因子以及由驱动调控因子控制的调控基因模块的方法。本发明提供的方法充分利用单细胞转录组数据的特征,构建更加精确的基因调控网络,并且基于此网络获得起关键调控作用的驱动调控因子,有助于解析细胞命运决定的调控机制。
为此,本发明第一方面提供一种构建基因调控网络的方法,所述方法包括以下步骤;
S01:获取细胞内部基因相互作用网络作为背景基因网络;
S02:获取单细胞转录组数据,对所述数据进行预处理,获取单细胞发育轨迹信息;
S03:基于图注意力神经网络构建所述基因调控网络的编码器;
S04:输入所述背景基因网络和所述单细胞发育轨迹信息,基于图对比学习框架对所述编码器进行训练;
S05:获取基因相互关系,利用注意力系数为所述基因相互关系赋权;
S06:设定基因相互关系权重阈值,选取高于所述阈值的基因相互关系,得到所述基因调控网络。
其中背景基因网络作为先验网络,用来从中选择与特定细胞命运相关的基因相互作用。单细胞转录组数据进行预处理后获取细胞的不同发育轨迹信息和伪时间信息,将先验网络同上述信息输入基于图注意力神经网络构建的编码器中,并基于图对比学习框架对编码器进行训练,构建得到细胞谱系特异的基因调控网络。该基因调控网络还为后续驱动调控因子的识别以及调控基因模块的识别提供了有力的基础。
根据本发明的实施方案,步骤S01中所述细胞内部基因相互作用网络来自NicheNet(https://github.com/saeyslab/nichenetr)、OmniPath(https://omnipathdb.org/)、InbioMap(https://inbio-discover.com/)或PathwayCommons(https://www.pathwaycommons.org/)中的至少一个,其中所述NicheNet、OmniPath、InbioMap和PathwayCommons剔除掉细胞间的配体-受体相互作用关系,将无向边处理为双向边,由此获得有向的基因关系网络。
根据本发明的实施方案,所述步骤S02中所述预处理包括:
基于所述单细胞转录组数据中每个细胞表达的基因数量对细胞进行过滤,根据每个基因参与表达的细胞数目对基因进行过滤,去除低质量细胞和基因后,对基因表达量进行归一化处理,获取单细胞发育轨迹信息。
根据本发明的实施方案,所述单细胞发育轨迹信息包括采用Slingshot、DPT、Palantir中的至少一种方法获取。
根据本发明的实施方案,步骤S03中所述编码器包括:注意力机制函数、批量归一化、前馈神经网络和激活函数。其中,输入的先验基因网络被表示为一个有向图G=(V,E),其中节点V(|V|=N)表示基因,E表示边。图G的邻接矩阵表示为其中Aij=1表示有一条从vi到vj的边,否则Aij=0。
根据本发明的实施方案,所述编码器的组成方式包括:每一层图注意力神经网络之前进行所述批量归一化,每一层图注意力神经网络之后连接所述前馈神经网络。
根据本发明的实施方案,所述组成方式堆叠两层。
根据本发明的实施方案,所述注意力机制函数包括如下公式:
其中,att函数为注意力函数,度量节点vi和节点vj间的相关性;l表示所述图注意力神经网络中的层;和/>分别为与源节点和目标节点相关的权重矩阵;/>和/>)分别为节点vi和节点vj的基因特征表示,其中对于第一层/>为输入的基因表达谱数据,Dscore为基因差异表达得分;
eij通过softmax函数归一化得到注意力系数αij;τ为温度参数,且τ<1;
最终节点vi的输出特征是得到的注意力系数αij所对应的特征的线性组合,使用多头注意力机制同时捕捉多个表征子空间的基因特征表示得到σ是GELU激活函数;K是多头注意力的总头数;||表示连接操。
基于注意力机制的模型,通过对节点特征之间的相关性进行评分,来学习每个节点的邻居对其特征表示的重要性。节点vi和节点vj间的相关性可以用一个注意力函数att来度量,后通过softmax函数归一化得到一个概率值,即注意力系数αij,最终每个节点的输出特征是上述归一化注意力系数所对应的特征的线性组合,并使用多头注意力机制进行拼接。
在注意力函数att中,使用了基于余弦相似性的注意力得分函数,以更直接地考虑节点特征之间的相似性,即只考虑基因调控关系的强度,而忽略调控作用是否为激活或者是抑制的。同时此步骤考虑了基因在细胞命运变化过程中起始和结束状态之间的差异表达,如下式所示:
Dscore∈[0,1]为基因差异表达得分,用来放大细胞发育过程中显著变化的基因之间的关联性。
此外在softmax函数中使用了温度参数τ。当τ>1时,越大的值会导致更均匀的概率;而当τ<1时,越小的值会导致更有区分度的概率。因此此处限定τ<1使得概率分布更加尖锐,有助于注意力集中于少数更相关的邻居节点上。
根据本发明的实施方案,所述批量归一化包括采用如下公式对节点vi的基因表达谱进行处理:
其中,BN为Batch Normalization函数。
根据本发明的实施方案,所述前馈神经网络包括如下函数:
其中,w1和w2为系数矩阵,b1和b2为偏置项。
根据本发明的实施方案,所述激活函数如下式所示:
根据本发明的实施方案,所述步骤S03进一步包括:将背景基因网络划分为传入网络(incoming network)和传出网络(outgoing network),所述传入网络和传出网络上的计算同时进行,根据以下公式将节点vi的基因特征表示连接起来:
其中,和/>分别表示传入网络和传出网络上的基因特征表示,concat表示拼接操作。
背景基因网络的划分可分别考虑在图神经网络中不同方向的消息传递,有助于理解基因作为调控因子或被调控目标的重要性,由此可提高注意力的可解释性,使其更具有生物学意义。
根据本发明的实施方案,所述基因差异表达得分Dscore的计算方式如下:
其中,Lx表示基因x的表达量的log2FC,c和d是标量参数,所述c和d分别在每一层图注意力神经网络上具有相同的数值。
根据本发明的实施方案,步骤S04中所述图对比学习框架包括采用深度图信息最大化算法,所述深度图信息最大化算法的损失函数如下式所示:
其中,N是所述图注意力神经网络中的节点数目,E表示步骤S03中的编码器,(X,A)是一对真实的节点特征和网络结构,是通过随机打乱网络而得到的作为负样本的节点特征和网络结构,si代表全局的图级别的表示,/>表示被扰动后的嵌入向量,/>是一个权重矩阵,σ是sigmoid激活函数。
深度图信息最大化算法(Deep Graph Infomax,DGI)通过最大化局部图结构和整个图的全局表示之间的互信息来对比地获得网络中节点的表示。所述DGI的损失函数被定义为二元交叉熵损失。
根据本发明的实施方案,所述si的计算方式如下:
根据本发明的实施方案,步骤S05中为所述基因相互关系赋权包括利用如下函数进行赋权:
其中,Aij表示邻接矩阵,Aij=1表示有一条从节点vi到节点vj的边,否则Aij=0;和/>分别为由传入网络和传出网络上计算得到的注意力系数;。
在编码器模型收敛后,使用注意力系数来度量基因之间的相关性的强弱,本发明通过将注意力系数分别与中心节点在传入网络的入度相乘或在传出网络中的出度相乘来对其进行缩放,最终的相互作用权重βij是两个方向网络上缩放注意力系数的平均值。当模型采用两层图注意力神经网络时,最终通过将两层的缩放注意力系数结合起来。
根据本发明的实施方案,步骤S06中所述设定基因相互作用关系权重阈值的方法包括:选择阈值使得构建的基因调控网络中边的数目为基因数目的k倍,其中k为用户设定的参数,该参数k即为所构建的基因调控网络的平均度,一般k为选自5~12的任一整数,k越大网络越稠密。
本发明第二方面提供一种驱动调控因子的识别方法,所述识别方法包括以下步骤:
S07:基于基因调控网络识别驱动基因,得到候选驱动基因;
S08:计算所述候选驱动基因的影响力得分,并根据所述影响力得分排序;
S09:挑选出影响力得分排名前n的基因作为驱动调控因子;
其中,所述基因调控网络通过第一方面所述的方法构建;
n为选自1-200的任一整数。
基于本发明第一方面提供的基因调控网络构建方法构建基因调控网络,进而识别驱动基因,通过定义一个关于基因在调控网络中重要性的影响力得分,识别细胞命运决定的驱动调控因子。
根据本发明的实施方案,所述步骤S07进一步包括:采用最小反馈节点集和最小支配节点集分别识别驱动基因。
最小反馈节点集(MFVS)和最小支配节点集(MDS)属于网络控制理论方法。其中,反馈节点集是指图中的一个节点集合,使得去除该集合后的图没有反馈回路。支配节点集是指图中的一个节点集合,使得每个支配集以外的其它节点都是支配集成员的邻居。
根据本发明的实施方案,步骤S07中所述候选驱动基因为通过最小反馈节点集和最小支配节点集识别到的驱动基因的并集。
根据本发明的实施方案,步骤S08中所述影响力得分的计算方式如下:
其中,和/>分别为传入网络和传出网络上的候选驱动基因影响力得分,ISi为最终的候选基因影响力得分;/>和/>分别表示节点vi的前驱邻居集合和后继邻居集合;和/>分别为传出网络和传入网络上的缩放注意力系数。
根据本发明的实施方案,本发明第三方面提供一种由驱动调控因子控制的调控基因模块的识别方法,所述识别方法包括以下步骤:
S10.分别从传入网络和传出网络中识别驱动调控因子及其邻居基因构成的调控基因模块;
其中,所述驱动调控因子通过第二方面所述的识别方法识别得到。
根据本发明的实施方案,步骤S10中所述调控基因模块包括:i)在传入网络中,所述调控基因模块包括一个驱动调控因子和所述驱动调控因子的靶基因;ii)在传出网络中,所述调控基因模块包括一个驱动调控因子及调控所述驱动调控因子的其它调控基因。
根据本发明的实施方案,所述识别方法进一步包括:
S11.对所述调控基因模块在每个细胞中的相对活性进行度量,识别不同细胞类型或状态下具有活性的基因模块。
根据本发明的实施方案,步骤S11中所述相对活性的度量方法包括采用AUCell、GSEA或Pagoda2中的至少一种方法。
AUCell方法是SCENIC方法中提出的基于细胞全基因组排名中曲线下面积(AUC)来计算输入基因集的活性评分指标,可以用它来识别活性较高的细胞状态相关的基因模块。
GSEA方法是一种基因集合分析方法,通过比较一个预定义的基因集合在不同条件下的基因表达谱,来评估该基因集合在某一特定条件下的富集程度。
Pagoda2方法对每个细胞拟合一个误差模型,对细胞中每个基因的残差进行重新规范化,每个基因组的评分矩阵由其第一个加权主成分进行量化。
本发明相对于现有技术具有以下有益效果:
本发明提供的方法基于单细胞转录组测序数据构建了基因调控网络,通过该网络可高效准确识别控制细胞命运决定的主要调控因子和调控基因模块。本发明的方法准确可靠,可应用到多领域中,通过预先计算获得的基因调控关系及关键的调控基因等信息可以为实验人员提供参考,节省生物实验成本,提高环保效益;也可应用到疾病治疗中,筛选出的影响疾病发生发展的重要基因及其调控关系,为医生诊断提供帮助,并且能够提高药物治疗靶点筛选的效率等,具有广阔的前景。
本发明的附加方面和优点将在下面的描述中部分给出,部分将从下面的描述中变得明显,或通过本发明的实践了解到。
附图说明
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解,其中:
图1显示了本发明基于基因调控网络识别驱动调控因子和基因调控模块的流程示意图;
图2显示了本发明基于图注意力神经网络构建基因调控网络的编码器结构,其中最左侧为总的编码器结构,中间及右侧为编码器中图注意力神经网络的详细结构。图中abs表示取绝对值运算,cos表示余弦函数,Dscore为本发明所定义的基因差异表达得分;
图3显示了本发明实施例1中不同基因调控网络构建方法的AUPRC结果。其中,纵坐标为AUPRC值,横坐标为不同的数据集;
图4显示了本发明实施例1中不同基因调控网络构建方法的EPR结果。其中,纵坐标为EPR值,横坐标为不同的数据集;
图5显示了本发明实施例2中分别采用CEFCON、SCENIC、VIPER和ANANSE方法对驱动调控因子识别性能评估的精确度结果。其中,第一行图显示了驱动调控因子在mESC数据集上的精确度结果,第二行图为hESC数据集上的精确度结果;四列图分别对应四个验证数据集,即细胞命运定型、干细胞群体维持、内胚层发育和文献收集的关键调控因子;每一个子图的纵坐标为精确度指标,横坐标为基因排序;
图6显示了本发明实施例2中在人胚胎干细胞分化数据中识别的排名前20的驱动调控因子。其中,上方柱状图的纵坐标为影响力得分IS值,横坐标为20个驱动调控基因,柱状图下方对应四个基因的验证数据集;
图7显示了本发明实施例2中对人胚胎干细胞分化数据中基因嵌入的UMAP可视化及主要聚类的基因功能;
图8显示了本发明实施例3中对调控基因模块的GO和KEGG功能富集方面性能评估的精确度结果;
图9显示了本发明实施例4中对小鼠造血干细胞分化的三个发育谱系识别得到的前20个驱动调控因子;
图10显示了本发明实施例4中对小鼠造血干细胞分化的三个发育谱系上公共的驱动调控因子的表达趋势。其中纵坐标为规范化后的基因表达量,横坐标为伪时间;
图11显示了本发明实施例4中对小鼠造血干细胞分化的三个发育谱系上各自的一个最典型的驱动调控因子的表达趋势。其中纵坐标为规范化后的基因表达量,横坐标为伪时间;
图12显示了本发明实施例4中小鼠造血干细胞分化的三个发育谱系上的网络可控性指标。
具体实施方式
下面参考具体实施例,对本发明进行描述,需要说明的是,这些实施例仅仅是描述性的,而不以任何方式限制本发明。
需要说明的是,术语“第一”、“第二”仅用于描述目的,而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此,限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。进一步地,在本发明的描述中,除非另有说明,“多个”的含义是两个或两个以上。
在本文中所披露的范围的端点和任何值都不限于该精确的范围或值,这些范围或值应当理解为包含接近这些范围或值的值。对于数值范围来说,各个范围的端点值之间、各个范围的端点值和单独的点值之间,以及单独的点值之间可以彼此组合而得到一个或多个新的数值范围,这些数值范围应被视为在本文中具体公开。
为了更容易理解本发明,以下具体定义了某些技术和科学术语。除显而易见在本文件中的它处另有明确定义,否则本文中使用的所有其它技术和科学术语都具有本发明所属领域的一般技术人员通常理解的含义。
在本文中,术语“包含”或“包括”为开放式表达,即包括本发明所指明的内容,但并不排除其他方面的内容。
在本文中,术语“任选地”、“任选的”或“任选”通常是指随后所述的事件或状况可以但未必发生,并且该描述包括其中发生该事件或状况的情况,以及其中未发生该事件或状况的情况。
根据本发明的实施例,本发明第一方面提供一种构建基因调控网络的方法,所述方法包括以下步骤;
S01:获取细胞内部基因相互作用网络作为背景基因网络。
可使用来自于NicheNet、OmniPath、InbioMap或PathwayCommons的集成多种基因相互作用关系的网络作为先验基因网络,用来从中选择与特定细胞命运相关的基因相互作用。由于只关注细胞内部的基因相互作用关系,因此将细胞间的配体-受体相互作用剔除掉。通过将其中的无向边处理为双向的,最终获得有向的基因关系网络。
S02:获取单细胞转录组数据,对所述数据进行预处理,获取单细胞发育轨迹信息。
其中所述预处理包括:
基于所述单细胞转录组数据中每个细胞表达的基因数量对细胞进行过滤,根据每个基因参与表达的细胞数目对基因进行过滤,去除低质量细胞和基因后,对基因表达量进行归一化处理,采用Slingshot、DPT、Palantir中的至少一种方法获取单细胞发育轨迹信息。
提取所要研究的发育轨迹上细胞的数据,进行基因差异表达分析,具体为计算沿单细胞发育轨迹的基因差异表达的显著性p-value值和倍数变化log2FC值。
其中p-value值通过Benjamin-Hochberg方法对其进行多重检验调整,并考虑错误发现率(false discovery rate,FDR)阈值为0.01,即认为基因的FDR<0.01为显著差异表达。对于显著差异表达的基因,使用以2为底对数的倍数变化(log2 Fold Change,log2FC)的绝对值来衡量差异化表达水平。
S03:基于图注意力神经网络构建所述基因调控网络的编码器。
输入的先验基因网络被表示为一个有向图G=(V,E),其中节点V(|V|=N)表示基因,E表示边。图G的邻接矩阵表示为其中Aij=1表示有一条从节点vi到节点vj的边,否则Aij=0。在图神经网络的每一层l∈{1,2,…,L},节点特征表示为其中/>F是节点特征的数量,此处第一层神经网络的输入节点特征h(1)即基因的表达谱数据。
其中编码器包括以下组件:注意力机制函数、批量归一化、前馈神经网络和激活函数,每一层图注意力神经网络之前进行批量归一化,每一层图注意力神经网络之后连接前馈神经网络,按照此种组成方式堆叠两层。
其中,注意力机制函数具体包括:
给定一个节点vi,基于注意力机制的模型通过对节点特征之间的相关性进行评分来学习每个节点的邻居对其特征表示的重要性。节点vi和节点vj之间的相关性可以用一个注意力函数att来度量,即下式:
其中和/>分别为与源节点和目标节点相关的权重矩阵。同时,在本发明中注意力函数使用为基于余弦相似性的注意力得分函数,以更直接地考虑节点特征之间的相似性。具体来说,本发明使用了余弦相似性注意力得分的绝对值,即只考虑基因调控关系的强度,而忽略调控作用是否为激活或者是抑制的,并且考虑了基因在细胞命运变化过程中起始和结束状态之间的差异表达,如下式所示:
其中,Dscore∈[0,1]为基因差异表达得分,用来放大细胞发育过程中显著变化的基因之间的关联性。关于Dscore的计算,首先使用MSAT算法识别沿细胞发育轨迹的起始和末端状态之间有显著差异表达的基因,其中起始状态和末端状态分别被定义为根据细胞的伪时间排序的前四分之一和末端四分之一部分。具体计算如下:
其中,Lx表示基因x的表达量的log2 FC,c和d是可学习的标量参数,它们在所有图神经网络层上都共享相同的数值。
节点vi和它所有的邻居节点之间的注意力系数经softmax函数归一化处理,表示为一个概率,用αij来表示:
其中,在softmax函数中使用了一个温度参数τ(其中τ<1),以使得概率分布更加尖锐,有助于注意力集中于少数更相关的邻居节点上。当τ>1时,越大的值会导致更均匀的概率;而当τ<1时,越小的值会导致更有区分度的概率。
最终每个节点的输出特征是上述归一化注意力系数所对应的特征的线性组合,本发明使用多头注意力机制来同时捕捉多个表征子空间的基因特征表示,具体如下所示:
其中,σ是GELU非线性激活函数,K是多头注意力的总头数,||表示连接操。
为了提高注意力的可解释性,使其更具有生物学意义,本发明根据先验基因网络中边的方向性将网络划分为传入网络和传出网络,用于分别考虑在图注意力神经网络中不同方向的消息传递。这种操作有助于理解基因作为调控因子或被调控目标的重要性。两个方向网络上的计算是同时进行的,最终输出的基因表示向量被连接起来,如下:
其中和/>分别表示传入网络和传出网络上的基因特征表示,concat表示向量拼接操作。
批量归一化函数如下式所示:
其中,BN为Batch Normalization函数。
前馈神经网络函数如下式所示:
其中w1和w2为系数矩阵,b1和b2为偏置项。
激活函数函数采用GELU函数,如下式所示:
S04:输入所述背景基因网络和所述单细胞发育轨迹信息,基于图对比学习框架对所述编码器进行训练。
本发明使用一种无监督、对比学习的方法框架,即DGI来训练上述图注意力神经网络。该方法通过最大化局部图结构和整个图的全局表示之间的互信息来对比地获得网络中节点的表示。DGI的损失函数被定义为二元交叉熵损失,即下式所示:
其中,N是图注意力神经网络中的节点数目,E表示步骤S03中所提出的基于图注意力神经网络的编码器,(X,A)是一对真实的节点特征和网络结构,是通过随机打乱网络而得到的作为负样本的节点特征和网络结构。本发明中通过保持网络的拓扑结构不变,而将真实的基因表达谱随机分配给其他基因来进行随机打乱。/>代表全局的图级别的表示,/>表示被扰动后的嵌入向量,/>是一个可训练的权重矩阵。这里的非线性激活函数σ采用sigmoid函数。
S05:获取基因相互关系,利用注意力系数为所述基因相互关系赋权。
在模型收敛后,使用注意力系数来度量基因之间的相关性的强弱。本发明通过将注意力系数分别与中心节点在传入网络的入度相乘或在传出网络中的出度相乘来对其进行缩放,最终的相互作用权重βij是两个方向网络上缩放注意力系数的平均值,如下式所示:
其中,Aij表示邻接矩阵,Aij=1表示有一条从节点vi到节点vj的边,否则Aij=0;和/>分别为由传入网络和传出网络上计算得到的注意力系数。
由于本发明中模型采用两层图注意力神经网络,最终通过将两层的缩放注意力系数结合起来。
S06:设定基因相互关系权重阈值,选取高于所述阈值的基因相互关系,得到所述基因调控网络。
相互作用关系权重阈值的选取方式如下:选择阈值使得基因调控网络中边的数目为网络中基因数目的k倍,其中k为用户设定的参数,该参数k即为所构建的基因调控网络的平均度,一般k为选自5~12的任一整数,k越大网络越稠密。
根据本发明的实施例,本发明第二方面提供一种驱动调控因子的识别方法,包括以下步骤:
S07:基于基因调控网络识别驱动基因,得到候选驱动基因。
所述基因调控网络通过第一方面所述的方法构建。可采用网络控制理论中的MFVS和MDS方法识别驱动基因,得到的候选驱动基因为通过MFVS和MDS方法识别到的驱动基因的并集。
基于MFVS的方法需要同时控制网络中所有的源节点(即入度为0的节点)和最小的一个反馈节点集。其中,反馈节点集是指图中的一个节点集合,使得去除该集合后的图没有反馈回路。采用如下的0-1整数规划问题进行求解:
where 0≤wi≤n-1 and xi={0,1}
其中,E表示边集合,eij为节点vi和节点vj之间的一条边。
基于MDS的方法需要控制网络中最小的一个支配集。支配集是指图中的一个节点集合,使得每个支配集以外的其它节点都是支配集成员的邻居。采用如下0-1整数规划问题进行求解:
其中,Ni表示节点vi的后继邻居节点。
以上步骤的求解可以采用Gurobi、SCIP、CPLEX等优化求解器。
S08:计算所述候选驱动基因的影响力得分,并根据所述影响力得分排序。
取候选驱动调控基因集合,通过如下方式分别计算在入度网络和出度网络上基因的影响力得分,进一步合并两个影响力得分作为最终的基因的影响力得分,
其中,和/>分别为传入网络和传出网络上的基因影响力得分,ISi为最终的基因影响力得分;/>和/>分别表示节点vi的前驱邻居集合和后继邻居集合;/>和/>分别为传出网络和传入网络上的缩放注意力系数。
S09:挑选出影响力得分排名前n的基因作为驱动调控因子。
n为选自1-200的任一整数。基因的排名选择可根据实际需要进行选择,例如选择影响力得分排名前10的基因、影响力得分排名前20的基因、影响力得分排名前100的基因、影响力得分排名前200的基因等。
根据本发明的实施例,本发明第三方面提供一种由驱动调控因子控制的调控基因模块的识别方法,包括以下步骤:
S10:分别从传入网络和传出网络中识别驱动调控因子及其邻居基因构成的调控基因模块。
其中,驱动调控因子通过第二方面所述的识别方法识别得到。此外本发明分别从传入网络和传出网络中识别两种不同类型的基因模块:i)在传入网络中,识别包括一个驱动调控因子和它的靶基因构成的基因模块(即驱动调控基因及其出度邻居);ii)在传出网络中,识别驱动调控因子作为靶标被其它多个基因调控所构成的基因模块(即驱动调控基因及其入度邻居)。
S11:对所述调控基因模块在每个细胞中的相对活性进行度量,识别不同细胞类型或状态下具有活性的基因模块。
其中相对活性的度量方法包括但不限于采用AUCell、GSEA或Pagoda2方法,任何可以实现对调控基因模块在每个细胞中的相对活性的度量的方法均可用于此步骤中。
下面将结合实施例对本发明的方案进行解释。本领域技术人员将会理解,下面的实施例仅用于说明本发明,而不应视为限定本发明的范围。除非另有说明,否则本文使用的所有技术和科学术语具有本发明所述领域的常规技术人员通常理解的相同含义。虽然本发明仅描述了优选的方法和数据材料,但是在本发明的实施或测试中也可以使用与本文所述相似或等同的任何方法和数据。
以下实施例围绕三个方面展开:(1)针对单细胞转录组数据的基因调控网络的构建;(2)基于所构建的基因调控网络识别驱动调控因子;(3)由驱动调控因子控制的调控基因模块的识别及活性度量。总流程图如图1所示。
实施例1针对单细胞转录组数据的基因调控网络的构建
S01:获取人和小鼠的背景基因网络。
NicheNet提供的基因网络是用人类基因符号给出的,因此为了获得小鼠的基因网络,使用ENSEMBL中的对应的人与小鼠的同源关系将基因符号映射到小鼠上,并删除了具有多个对应信息的基因。最后,针对人类得到一个包含25,332个基因和5,290,993条边的先验基因相互作用网络,针对小鼠得到一个包含18,137个基因和5,018,393条边的先验基因相互作用网络;
S02:获取并处理包含细胞发育轨迹信息的单细胞转录组数据。
对单细胞转录组数据进行预处理,采用Scanpy对于原始数据为reads数目的数据用TPM或CPM方法进行规范化,并提取高变异的基因用于后续的计算和分析。采用Slingshot工具推断细胞的伪时间信息并构建细胞发育轨迹。采用MAST工具进行基因差异表达分析,对比细胞发育轨迹初始阶段和末端阶段两个状态,获取显著差异表达的基因并计算倍数变化,其中差异表达显著性选取多重校验后的p-value值小于0.05作为阈值;
S03:基于图注意力神经网络构建基因调控网络的编码器。
该编码器包括注意力机制函数、批量归一化、前馈神经网络和激活函数。该编码器的组成方式为每一层图注意力神经网络之前进行批量归一化,每一层图注意力神经网络之后连接前馈神经网络,按照此种方式堆叠两层。该编码器结构及图注意力神经网络构造如图2所示。
其中,注意力机制函数包括如下公式:
/>
在所述图注意力神经网络的每一层l中,节点特征表示为 F是节点特征的数量,此处第一层的/>为基因的表达谱数据,att函数为注意力函数,度量节点vi和节点vj间的相关性,/>和/>分别为与源节点和目标节点相关的权重矩阵,/>和/>)分别为节点vi和节点vj的基因表达谱数据,Dscore为基因差异表达得分;
eij通过softmax函数归一化得到注意力系数αij,τ为温度参数,且τ<1;
为节点vi的输出特征,σ是GELU激活函数,K是多头注意力的总头数,||表示连接操。
批量归一化包括采用如下公式对节点vi的基因表达谱进行处理:
BN为Batch Normalization函数。
前馈神经网络包括如下函数:
w1和w2为系数矩阵,b1和b2为偏置项。
激活函数如下式所示:
S04:输入步骤S01中的背景基因网络和步骤S02得到的单细胞发育轨迹信息,基于图对比学习框架对所述编码器进行训练。
选取二元交叉熵损失作为DGI的损失函数,如下式所示:
其中,N是图注意力神经网络中的节点数目,E表示步骤S03中所提出的基于图注意力神经网络的编码器,(X,A)是一对真实的节点特征和网络结构,是通过随机打乱网络而得到的作为负样本的节点特征和网络结构。本发明中通过保持网络的拓扑结构不变,而将真实的基因表达谱随机分配给其他基因来进行随机打乱。/>代表全局的图级别的表示,/>表示被扰动后的嵌入向量,/>是一个可训练的权重矩阵。这里的非线性激活函数σ采用sigmoid函数;
S05:获取基因相互关系,利用注意力系数为所述基因相互关系赋权。
通过将注意力系数分别与中心节点在传入网络的入度相乘或在传出网络中的出度相乘来对其进行缩放,最终的相互作用权重βij是两个方向网络上缩放注意力系数的平均值,如下式所示:
其中,Aij表示邻接矩阵,Aij=1表示有一条从节点vi到节点vj的边,否则Aij=0;和/>分别为由传入网络和传出网络上计算得到的注意力系数。最终通过/>将两层的缩放注意力系数结合起来;
S06:设定基因相互关系权重阈值,选取高于所述阈值的基因相互关系,得到所述基因调控网络。
选择阈值使得基因调控网络中边的数目为网络中基因数目的k倍,其中k为用户设定的参数,该参数k即为所构建的基因调控网络的平均度,k为选自5~12的任一整数。
实验结果验证:
为了验证所提供的构建基因调控网络方法(简称为CEFCON)的有效性,在不同类型的数据集上进行了测试,并且和经典的GRNBoost2、NetREX、SCINET方法以及在背景基因网络上随机选取与真实验证网络中相同数量边的方法(Random)进行了比较。
在本实施例中采用的单细胞转录组数据集包括了人类胚胎干细胞(hESC)、人类成熟肝细胞(hHep)、小鼠树突细胞(mDC)、小鼠胚胎干细胞(mESC)和三个小鼠造血干细胞谱系(mHSC),即红细胞(mHSC-E)、粒细胞-单核细胞(mHSC-GM)和淋巴细胞(mHSC-L)。对于用于验证的基准数网络数据集,采用对应细胞类型的特异性ChIP-seq数据。此外,还针对mESC数据集收集了另外两个真实的基准网络,以进行更全面的评估。第一个基准网络来自ESCAPE数据库的功能缺失/功能获得(lof/gof)的数据(mESC-lofgof),另一个基准网络来自于94个转录因子诱导表达的基因表达谱(mESC-TFind)。对于每个数据集,只考虑前1000个高变异基因(Highly VariableGene,HVG)进行评估。
采用AUPRC和EPR两个指标评价方法的性能:
AUPRC表示Precision-Recall曲线下的面积,该曲线的纵坐标为Precision(精确率),横坐标为Recall(召回率):
其中TP为真阳性样本数目,FP为假阳性样本数,FN为假阴性样本数;
EPR(Early Precision Rate)定义为早期精度值(Early Precision)与随机预测器的早期精度值的比率。早期精度值被定义为预测的前k条边中真阳性的部分,其中k为真实网络中的边的数量。
图3和图4分别为本发明提供的基因调控网络构建方法的AUPRC和EPR结果,可以看出本发明提供的网络构建方法具有优越性,在大多数数据集上的性能明显好于其他经典方法。
实施例2基于所构建的基因调控网络识别驱动调控因子
S07:基于所构建的基因调控网络识别驱动基因,得到候选驱动基因。
采用网络控制理论中的MFVS和MDS方法识别驱动基因。其中,基于MFVS的方法需要同时控制网络中所有的源节点(即入度为0的节点)和最小的一个反馈节点集。采用如下的0-1整数规划问题进行求解:
where 0≤wi≤n-1 and xi={0,1}
其中,E表示边集合,eij为节点i和节点j之间的一条边。
基于MDS的方法需要控制网络中最小的一个支配集。采用如下0-1整数规划问题进行求解:
其中,Ni表示节点i的后继邻居节点;
将通过最小反馈节点集和最小支配节点集识别到的驱动基因进行并集运算得到候选驱动基因;
S08:计算候选驱动基因的影响力得分,并根据该影响力得分排序。
取候选驱动调控基因集合,通过如下方式分别计算在入度网络和出度网络上基因的影响力得分,进一步合并两个影响力得分作为最终的基因的影响力得分,
其中,和/>分别为传入网络和传出网络上的基因影响力得分,ISi为最终的基因影响力得分,/>和/>分别表示节点vi的前驱邻居集合和后继邻居集合,/>和/>分别为传出网络和传入网络上的缩放注意力系数;
S09:挑选出影响力得分排名前20的基因作为最终的驱动调控因子。
实验结果验证:
为了评估本发明提供的驱动调控因子识别方法的性能,我们对于实施例1中实验结果验证中提到的mESC和hESC这两个数据集进行实验。对于验证数据集,从基因本体数据库(Gene Ontology,GO)中收集了三个与细胞命运和胚胎干血细胞发育相关的基因集合,以及从两个参考文献中收集的经过实验验证的基因集合,具体信息如下表1所示:
表1用于评估驱动调控因子识别方法的基准数据集信息
其中所用的GO基因集合是“细胞命运定型(GO_Cell Fate Commitiment)”、“干细胞群体维持(GO_Stem Cell Population Maintenace)”和“内胚层发育(GO_EndodermDevelopment)”相关的三个基因集合。“干细胞群体维持”是干细胞关于其自我更新和多能性的特征。使用“内胚层发育”是因为本发明使用的mESC和hESC数据集都是关于从胚胎干细胞分化为内胚层细胞的。两个参考文献中收集的经过实验验证的基因集合主要涉及细胞命运基因的发育表达和控制胚胎干细胞状态的转录因子。
为了评估本发明提供的驱动调控因子识别方法的准确性,与两个用于对传统群细胞RNA-seq数据上进行主要调控因子(Master regulators)预测的方法VIPER和ANANSE,以及一种基于单细胞数据用于发现关键转录因子及其调控子的方法SCENIC进行了比较,采用Top-k precision作为评估指标。Top-k precision被定义为前k个预测结果中属于基准基因集的比例,即:
其中,TPk是k个预测结果中正确的基因数,我们选择k的取值为1到20进行评估。
图5显示了各个方法识别得到的驱动调控基因在mESC和hESC数据集上对于四个真实基准基因集评估的性能结果。结果表明,本发明提供的驱动调控因子识别方法整体性能优越,在大多数基准基因集上得到最高的准确率,特别是排名靠前的基因基本都能被证明与细胞命运决定相关,例如,排名最高的基因在所有的基准基因集中都能被准确预测到(即Top-1 precision几乎都为1)。
而后进一步展示了在hESC数据中发现的排名前20的驱动调控因子,其中一半可以在真实的基准基因集中得到验证,如图6所示。其中三个转录因子,即NANOG、SOX2和POU5F1(OCT3/4),是人类和小鼠胚胎干细胞分化中著名的多能性因子,本发明提供的方法对这三个转录因子均给出了非常高的预测排名(分别为1/20,4/20和5/20)。此外,值得注意的是,另外两个预测排名较高的驱动调控因子GATA4(排名7/20)和GATA6(排名10/20),它们的过表达可以导致胚胎干细胞高度选择性地分化为原始内胚层细胞。事实上,导致细胞多能性或是分化的因素在胚胎干细胞中是精细平衡的,这种平衡的紊乱可以促进或抑制原始内胚层的分化。这些预测到的关键调控因子都是启动基因转录的先锋因子,能够驱动细胞命运的决定。
此外,本发明提供的方法可以识别转录因子以外的的驱动调控因子,如图6所示。在前20个驱动调控因子中,CYP26A1、ERBB4、EPCAM和CDH1并不是转录因子,但它们已被报道参与胚胎干细胞分化。其中CYP26A1通过调节视黄酸信号在脊椎动物的胚胎发育中起着关键作用,而视黄酸是胚胎干细胞分化的关键;CDH1编码E-钙黏蛋白(E-cadherin),它被证明能增强胚胎干细胞的多能性和维持自我更新;EPCAM被发现在未分化的胚胎干细胞中选择性地高表达,它直接调节几个重要的转录因子,包括POU5F1、NANOG和SOX2;ERBB4基因编码的蛋白是一种受体酪氨酸激酶,可调节细胞的增殖和分化。作为信号通路的源头,一些受体激酶可被视为在转录因子的上游发挥作用,也可能在细胞命运决定中发挥驱动作用。然而,目前大多数寻找主要调控因子的方法,包括ANANSE、SCENIC和VIPER,都仅限于找出重要的转录因子,无法获得上述这些基因。
尽管本发明提供的方法在识别调控关系和识别驱动调控因子时高度依赖于注意力系数,但作为一个基于图神经网络的方法,其模型输出部分的基因向量表示也十分重要,发明人发现模型学到的基因嵌入可以粗略地揭示出一定的生物意义。如图7所示,在对hESC数据上获得的基因嵌入使用UMAP可视化,并根据Leiden方法对其进行聚类。结果显示,由低分辨率的Leiden聚类方法得到的结果中,本发明提供的方法识别到的驱动调控因子主要分布在UMAP空间的一侧,其中前20的驱动调控因子中有13个都在聚类3中。为此,进一步对聚类3采用更高的分辨率,结果显示聚类3被分为了4个子类别,即图7中右侧的聚类7、8、9和10。值得注意的是,三个核心转录因子NANOG、SOX2和POU5F1都属于聚类7。图7给出了这些基因聚类显著富集到GO生物过程中最重要的功能条目,主要是关于内胚层的形成、发育过程和细胞表面受体信号通路等,表明CEFCON方法可以提供具有生物学意义的基因嵌入表示,并能揭示驱动调控因子的信息。
实施例3由驱动调控因子控制的调控基因模块的识别及活性度量
S10:对于通过实施例2识别到的驱动调控因子,分别从传入网络和传出网络中识别两种不同类型的基因模块:i)在传入网络中,识别一个驱动调控因子和它的靶基因构成的基因模块(即驱动调控基因及其出度邻居);ii)在传出网络中,识别驱动调控因子作为靶标被其它多个基因调控所构成的基因模块(即驱动调控基因及其入度邻居)。
只考虑与预测得到的驱动调控因子相关的基因模块,分别从传入网络和传出网络中获得影响力得分排名前50的驱动调控因子对应的基因模块,最后的结果只保留显著差异表达的基因的,并且去除规模小于10的基因模块作为最终的结果。
S11:对所识别到的调控基因模块在每个细胞中的相对活性进行度量,识别不同细胞类型或状态下具有活性的基因模块。
其中相对活性的度量方法包括采用AUCell方法,其是SCENIC中提出的基于细胞全基因组排名中曲线下面积(AUC)来计算输入基因集的活性评分指标,可以用它来识别活性较高的细胞状态相关的基因模块。
实验结果验证:
为了评估本发明提供的由驱动调控因子控制的调控基因模块的识别性能,我们从基因集合功能富集分析和细胞状态聚类两个层面做了分析,应用于具有已知细胞类型或细胞状态的scRNA-seq数据集,包括hESC、mESC和三个mHSC分化谱系,并且与SCENIC方法进行了比较。
首先,使用富集到基因本体(Gene Ontology,GO)中标注为“生物过程(biologicalprocess)”功能的基因集合以及KEGG(Kyoto Encyclopedia of Genes and Genomes)通路基因集合的精确度来评估各个算法得到的基因模块是否具有生物学意义。使用Fisher精确检验计算统计显着性p值,Fisher精确检验基于超几何分布,其概率公式是:
其中N是总的基因数目,M是已知在GO中注释为“生物过程”的基因数目,n是算法得到的该模块中基因的数目,k是M这个已知功能的基因和n个算法得到的基因交集的数目。最终选取经过Benjamin-Hochberg方法多重校验后p值小于0.05的那些基因模块作为显著的功能模块。可以使用gProfiler工具进行GO和KEGG通路富集分析。
结果如图8所示,本发明中在基因模块识别方面在所有测试数据上都明显优于SCENIC方法,表明本发明提供的方法可以获得更多生物学上合理且有意义的基因模块。
实施例4基于基因调控网络进行小鼠造血干细胞分化研究
造血过程是一个非常复杂且需要精细调节的过程,造血干细胞(Hematopoieticstem cells,HSCs)通过这一过程产生属于髓系或淋巴系的成熟血细胞。为了进一步说明本发明提出的方法对解析细胞命运决定调控机制的作用,因此应用于小鼠造血干细胞分化相关的单细胞RNA-seq数据集上,并希望深入探究造血干细胞向不同细胞命运分化的重要驱动调控因子,以解析细胞命运的调控机制。
使用Scanpy工具对原始的read counts的单细胞RNA-seq数据进行预处理,删去超过200个基因没有表达的细胞,删除在少于5个细胞上表达的基因。通过使用cell_ranger方法选择前3000个高变异基因,最后对表达值进行对数进行规范化及归一化处理。使用Slingshot方法计算伪时间信息并基于UMAP空间中的坐标计算轨迹推断。最后获得该造血干细胞分化过程中涉及的细胞类型、包含的主要三条分化轨迹以及伪时间信息,并使用MAST工具获得基因的差异表达信息。
从造血干细胞开始,细胞命运在三个方向上经过一些重要的细胞状态或类型,即红细胞样谱系(MEP)、粒细胞-单核细胞系谱系(GMP)、淋巴细胞样谱系(LMPP),为此分别应用本发明提供的方法计算了这三个主要分化谱系,以确定控制其不同细胞命运轨迹的驱动调控因子。图9为计算得到每个细胞谱系的排名前20个驱动调控因子,其中有6个基因为三个细胞谱系上公共的驱动调控因子。
对得到的驱动调控因子的在三个细胞命运谱系上的基因表达趋势进行分析,利用Palantir工具获得沿着伪时间的基因表达趋势。图10显示了三个细胞谱系上公共的驱动调控因子的基因表达趋势。发现这些基因沿着伪时间有相似的表达趋势,它们都在造血干细胞中高度表达,并且它们的表达随着分化水平的提高而降低。这些基因都被报道与维持造血干细胞自更新以及控制分化、调控细胞生长增殖相关。需要特别注意的是,本发明方法找到的Dusp1基因是6个共同驱动调控基因中唯一的一个非转录因子,它是一个著名的与细胞增殖相关的基因。该数据集中淋巴细胞样谱系的细胞相对较少,其与粒细胞-单核细胞的轨迹重合度比较高,它们共享更多的驱动调控基因,这与经典的造血分化模型一致。其次,各细胞谱系特异的驱动调控因子显示了它们在特定方向上引导细胞命运的能力。图11给出了三个细胞谱系上各自的一个最典型的驱动调控因子,如Klf1对于控制红细胞分化、Irf8控制粒细胞-单核细胞分化和Mef2c控制淋巴细胞分化,这些驱动调控因子明显表现出在发育中后期对于特定细胞命运方向上的表达增高的趋势。
进一步的,采用网络控制理论建模细胞命运动力学,进一步从网络可控性的角度分析造血干细胞的三个细胞谱系的分化情况。为此分别评估了基于MDS控制方法和基于MFVS控制方法的可控性得分(Controllability Score)、这两种控制方法得到的驱动调控因子的Jaccard系数(Jaccard Index),以及驱动调控因子在网络控制方法结果上的覆盖率(Driver Regulators Coverage),这些度量指标用于衡量控制网络状态的困难程度,它们的具体定义如下:
其中Dx和Dy分别表示通过方法x和方法y获得的驱动基因集合,G表示网络中所有基因的集合,CD表示通过本发明提供的方法最终识别的驱动调控因子的集合。这些指标值越大,表明控制网络的状态所需的基因越少,并且两种控制方法以及最终得到的关键驱动调控因子的一致性越高,进一步说明由基因调控网络状态代表的细胞状态发生变化的概率就越大。图12显示了基于本发明提供的方法得到的造血干细胞分化三个细胞谱系上获得驱动调控因子对应的可控性方面的度量指标,从该结果可以看出,控制造血干细胞向红细胞样分化需要较少比例的驱动基因,并且两种控制方法得到的驱动基因集合以及高影响力得分的基因重叠度比较高,说明整个调控网络的状态相对于其它分化谱系来说更容易控制。
在本说明书的描述中,参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、“一些实施方案”或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中,对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且,描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外,在不相互矛盾的情况下,本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管上面已经示出和描述了本发明的实施例,可以理解的是,上述实施例是示例性的,不能理解为对本发明的限制,本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (18)

1.一种构建基因调控网络的方法,其特征在于,所述方法包括以下步骤:
S01:获取细胞内部基因相互作用网络作为背景基因网络;
S02:获取单细胞转录组数据,对所述数据进行预处理,获取单细胞发育轨迹信息;
S03:基于图注意力神经网络构建所述基因调控网络的编码器;
S04:输入所述背景基因网络和所述单细胞发育轨迹信息,基于图对比学习框架对所述编码器进行训练;
S05:获取基因相互关系,利用注意力系数为所述基因相互关系赋权;
S06:设定基因相互关系权重阈值,选取高于所述阈值的基因相互关系,得到所述基因调控网络。
2.根据权利要求1所述的方法,其特征在于,步骤S01中所述细胞内部基因相互作用网络来自NicheNet、OmniPath、InbioMap、PathwayCommons中的至少一个,其中所述NicheNet、OmniPath、InbioMap和PathwayCommons剔除掉细胞间的配体-受体相互作用关系,将无向边处理为双向边。
3.根据权利要求1所述的方法,其特征在于,步骤S02中所述预处理包括:
基于所述单细胞转录组数据中每个细胞表达的基因数量对细胞进行过滤,根据每个基因参与表达的细胞数目对基因进行过滤,去除低质量细胞和基因后,对基因表达量进行归一化处理,获取单细胞发育轨迹信息;
任选地,所述单细胞发育轨迹信息包括采用Slingshot、DPT、Palantir中的至少一种方法获取。
4.根据权利要求1所述的方法,其特征在于,步骤S03中所述编码器包括:
注意力机制函数、批量归一化、前馈神经网络和激活函数。
5.根据权利要求4所述的方法,其特征在于,所述编码器的组成方式包括:每一层图注意力神经网络之前进行所述批量归一化,每一层图注意力神经网络之后连接所述前馈神经网络;
任选地,所述组成方式堆叠两层。
6.根据权利要求4所述的方法,其特征在于,所述注意力机制函数包括如下公式:
其中,att函数为注意力函数,度量节点vi和节点vj间的相关性;l表示所述图注意力神经网络中的层;和/>分别为与源节点和目标节点相关的权重矩阵;/>和/>分别为节点vi和节点vj的基因特征表示,其中对于第一层/>为输入的基因表达谱数据;Dscore为基因差异表达得分;
eij通过softmax函数归一化得到注意力系数αij;τ为温度参数,且τ<1;
最终节点vi的输出特征是得到的注意力系数αij所对应的特征的线性组合,使用多头注意力机制同时捕捉多个表征子空间的基因特征表示得到是GELU激活函数;K是多头注意力的总头数;||表示连接操;
任选地,所述批量归一化包括采用如下公式对节点vi的基因表达谱进行处理:
其中,BN为Batch Normalization函数;
任选地,所述前馈神经网络包括如下函数:
其中,w1和w2为系数矩阵,b1和b2为偏置项;
任选地,所述激活函数如下式所示:
7.根据权利要求6所述的方法,其特征在于,所述步骤S03进一步包括:将背景基因网络划分为传入网络和传出网络,所述传入网络和传出网络上的计算同时进行,根据以下公式将节点vi的基因特征表示连接起来:
其中,和/>分别表示传入网络和传出网络上的基因特征表示,concat表示拼接操作。
8.根据权利要求6所述的方法,其特征在于,所述基因差异表达得分Dscore的计算方式如下:
其中,Lx表示基因x的表达量的log2FC,c和d是标量参数,所述c和d分别在每一层图注意力神经网络上具有相同的数值。
9.根据权利要求1所述的方法,其特征在于,步骤S04中所述图对比学习框架包括采用深度图信息最大化算法,所述深度图信息最大化算法的损失函数如下式所示:
其中,N是所述图注意力神经网络中的节点数目;E表示步骤S03中的编码器;(X,A)是一对真实的节点特征和网络结构;是通过随机打乱网络而得到的作为负样本的节点特征和网络结构;si代表全局的图级别的表示;/>表示被扰动后的嵌入向量;/>是一个权重矩阵;σ是sigmoid激活函数;
任选地,所述si的计算方式如下:
10.根据权利要求7所述的方法,其特征在于,步骤S05中为所述基因相互关系赋权包括利用如下函数进行赋权:
其中,Aij表示邻接矩阵,Aij=1表示有一条从节点vi到节点vj的边,否则Aij=0;分别为由传入网络和传出网络上计算得到的注意力系数。
11.根据权利要求1所述的方法,其特征在于,步骤S06中所述设定基因相互作用关系权重阈值的方法包括:选择阈值使得构建的基因调控网络中边的数目为基因数目的k倍,其中k为选自5~12的任一整数。
12.一种驱动调控因子的识别方法,其特征在于,所述识别方法包括以下步骤:
S07:基于基因调控网络识别驱动基因,得到候选驱动基因;
S08:计算所述候选驱动基因的影响力得分,并根据所述影响力得分排序;
S09:挑选出影响力得分排名前n的基因作为驱动调控因子;
其中,所述基因调控网络通过权利要求1-11任一项所述的方法构建;
n为选自1-200的任一整数。
13.根据权利要求12所述的识别方法,其特征在于,所述步骤S07进一步包括:采用最小反馈节点集和最小支配节点集分别识别驱动基因;
任选地,步骤S07中所述候选驱动基因为通过最小反馈节点集和最小支配节点集识别到的驱动基因的并集。
14.根据权利要求12所述的识别方法,其特征在于,步骤S08中所述影响力得分的计算方式如下:
其中,和/>分别为传入网络和传出网络上的候选驱动基因的影响力得分,ISi为最终的候选驱动基因的影响力得分;/>和Ni s分别表示节点vi的前驱邻居集合和后继邻居集合;/>和/>分别为传出网络和传入网络上的缩放注意力系数。
15.一种由驱动调控因子控制的调控基因模块的识别方法,所述识别方法包括以下步骤:
S10.分别从传入网络和传出网络中识别驱动调控因子及其邻居基因构成的调控基因模块;
其中,所述驱动调控因子通过权利要求12-14任一项所述的识别方法识别得到。
16.根据权利要求15所述的识别方法,其特征在于,步骤S10中所述调控基因模块包括:i)在传入网络中,所述调控基因模块包括一个驱动调控因子和所述驱动调控因子的靶基因;ii)在传出网络中,所述调控基因模块包括一个驱动调控因子及调控所述驱动调控因子的其它调控基因。
17.根据权利要求15所述的识别方法,其特征在于,所述识别方法进一步包括:
S11.对所述调控基因模块在每个细胞中的相对活性进行度量,识别不同细胞类型或状态下具有活性的基因模块。
18.根据权利要求17所述的识别方法,其特征在于,步骤S11中所述相对活性的度量方法包括采用AUCell、GSEA、Pagoda2中的至少一种方法。
CN202310524599.8A 2023-05-10 2023-05-10 一种基于基因调控网络识别驱动调控因子的方法 Pending CN116486902A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310524599.8A CN116486902A (zh) 2023-05-10 2023-05-10 一种基于基因调控网络识别驱动调控因子的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310524599.8A CN116486902A (zh) 2023-05-10 2023-05-10 一种基于基因调控网络识别驱动调控因子的方法

Publications (1)

Publication Number Publication Date
CN116486902A true CN116486902A (zh) 2023-07-25

Family

ID=87219459

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310524599.8A Pending CN116486902A (zh) 2023-05-10 2023-05-10 一种基于基因调控网络识别驱动调控因子的方法

Country Status (1)

Country Link
CN (1) CN116486902A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116825204A (zh) * 2023-08-30 2023-09-29 鲁东大学 一种基于深度学习的单细胞rna序列基因调控推断方法
CN116844645A (zh) * 2023-08-31 2023-10-03 云南师范大学 一种基于多视角分层超图的基因调控网络推断方法
CN117497062A (zh) * 2023-11-15 2024-02-02 广州瑞能精准医学科技有限公司 一种特发性肺纤维化浆细胞特征基因预后模型构建方法
CN117746995A (zh) * 2024-02-21 2024-03-22 厦门大学 基于单细胞rna测序数据的细胞类型识别方法、装置及设备

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116825204A (zh) * 2023-08-30 2023-09-29 鲁东大学 一种基于深度学习的单细胞rna序列基因调控推断方法
CN116825204B (zh) * 2023-08-30 2023-11-07 鲁东大学 一种基于深度学习的单细胞rna序列基因调控推断方法
CN116844645A (zh) * 2023-08-31 2023-10-03 云南师范大学 一种基于多视角分层超图的基因调控网络推断方法
CN116844645B (zh) * 2023-08-31 2023-11-17 云南师范大学 一种基于多视角分层超图的基因调控网络推断方法
CN117497062A (zh) * 2023-11-15 2024-02-02 广州瑞能精准医学科技有限公司 一种特发性肺纤维化浆细胞特征基因预后模型构建方法
CN117746995A (zh) * 2024-02-21 2024-03-22 厦门大学 基于单细胞rna测序数据的细胞类型识别方法、装置及设备
CN117746995B (zh) * 2024-02-21 2024-05-28 厦门大学 基于单细胞rna测序数据的细胞类型识别方法、装置及设备

Similar Documents

Publication Publication Date Title
CN116486902A (zh) 一种基于基因调控网络识别驱动调控因子的方法
Talukder et al. Interpretation of deep learning in genomics and epigenomics
Krishnasamy et al. A hybrid approach for data clustering based on modified cohort intelligence and K-means
Hvidsten et al. Predicting gene function from gene expressions and ontologies
Li et al. Grouped gene selection of cancer via adaptive sparse group lasso based on conditional mutual information
Liu et al. Robust PCA based method for discovering differentially expressed genes
Wang et al. Unsupervised cluster analysis and gene marker extraction of scRNA-seq data based on non-negative matrix factorization
Dhyaram et al. RANDOM SUBSET FEATURE SELECTION FOR CLASSIFICATION.
Wang et al. Poisson-based self-organizing feature maps and hierarchical clustering for serial analysis of gene expression data
Chen et al. DGEPN-GCEN2V: a new framework for mining GGI and its application in biomarker detection
Yan et al. Identification of cell-type marker genes from plant single-cell RNA-seq data using machine learning
TW202121223A (zh) 訓練類神經網路以預測個體基因表現特徵的方法及系統
Shah et al. An experiment on ab initio discovery of biological knowledge from scRNA-Seq data using machine learning
Li et al. Multiclass nonnegative matrix factorization for comprehensive feature pattern discovery
CN115206423A (zh) 基于标签指导的蛋白质作用关系预测方法
Ma et al. EnsembleKQC: an unsupervised ensemble learning method for quality control of single cell RNA-seq sequencing data
Gong et al. Interpretable single-cell transcription factor prediction based on deep learning with attention mechanism
Tripathy et al. A Healthcare Data Analysis Approach for Breast Cancer Gene expression
Qiao et al. A personalized low-rank subspace clustering method based on locality and similarity constraints for scRNA-seq data analysis
Budiarto et al. Explainable supervised method for genetics ancestry estimation
Karaaslanli et al. Multiview Graph Learning for single-cell RNA sequencing data
Bouazza et al. Classifying Leukemia through DNA Expression Data Mining Techniques
Yu et al. CNLLRR: a novel low-rank representation method for single-cell RNA-seq data analysis
Kawaguchi et al. Exploiting marker genes for robust classification and characterization of single-cell chromatin accessibility
Bell et al. Development of novel methodology for gene identification-based classification of leukaemia disorder

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