CN110808083B - 基于scRNA-seq及动态时间规整的基因调控网络构建方法 - Google Patents

基于scRNA-seq及动态时间规整的基因调控网络构建方法 Download PDF

Info

Publication number
CN110808083B
CN110808083B CN201911012024.8A CN201911012024A CN110808083B CN 110808083 B CN110808083 B CN 110808083B CN 201911012024 A CN201911012024 A CN 201911012024A CN 110808083 B CN110808083 B CN 110808083B
Authority
CN
China
Prior art keywords
gene
dtw
network
regulation
time
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
CN201911012024.8A
Other languages
English (en)
Other versions
CN110808083A (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.)
Nantong University
Original Assignee
Nantong 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 Nantong University filed Critical Nantong University
Priority to CN201911012024.8A priority Critical patent/CN110808083B/zh
Publication of CN110808083A publication Critical patent/CN110808083A/zh
Application granted granted Critical
Publication of CN110808083B publication Critical patent/CN110808083B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • Y—GENERAL 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
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information 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)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (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)
  • Molecular Biology (AREA)
  • Genetics & Genomics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Analytical Chemistry (AREA)
  • Physiology (AREA)
  • Artificial Intelligence (AREA)
  • Bioethics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Chemical & Material Sciences (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)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

本发明具体涉及基于scRNA‑seq及动态时间规整的基因调控网络构建方法，属于基因调控网络技术领域。该方法包括以下步骤：步骤1、基于单细胞数据的单细胞伪时间分析；步骤2、基于正负调控方式的双重DTW距离分析；步骤3、基于DTW距离的聚类分析以及GRN模块筛选；步骤4、基于DTW距离的基因间相似性矩阵构建；步骤5、通过随机置换检验，计算出基因与基因间网络连接度阈值C；步骤6、对GRN模块进行基因功能性分析；步骤7、pGRN网络的可视化得出GRN网络图。本发明方法的pGRN能够发现时间调控滞后的调控关系，同时处理时间滞后和时间伸缩的调控关系，为新的调控关系的发掘提供了途径。

Description

基于scRNA-seq及动态时间规整的基因调控网络构建方法
技术领域
本发明具体涉及基于scRNA-seq及动态时间规整的基因调控网络构建方法，属于基因调控网络技术领域。
背景技术
基因调控网络(Gene Regulatory Network,GRN)对基因与基因之间调控关系进行系统性描述的重要方法,已被广泛应用于人类疾病、干细胞多能性及其改造、农作物育种等众多研究领域。在微观方面，GRN几乎参与了所有细胞的活动，包括细胞分裂、细胞信号传导、细胞代谢等活动。目前GRN的构建方法主要包括相关性模型、布尔模型、贝叶斯模型、常微分方程、信息论模型、神经网络模型等。上述方法主要针对基因表达芯片数据或者混合细胞转录组测序数据，尚未针对单细胞转录组测序数据(scRNA-seq)特征。随着单细胞转录组测序技术的发展及其相关数据的大量出现，开发针对单细胞转录组测序数据(scRNA-seq)的基因调控网络构建方法具有重要的应用价值。在过去几年，陆续有针对单细胞转录组测序数据(scRNA-seq)特征的基因调控网络方法出现，包括SingCellNet、SCNS、BTR、SCENIC、LEAP等。这些方法中，有的没有考虑单细胞转录组测序数据(scRNA-seq)伪时间信息，如SingCellNet、SCENIC、BTR，有的方法未能考虑时间滞后的调控关系，如LEAP之外的方法。考虑了时间滞后调控关系的LEAP方法，不能很好的处理有时间伸缩的调控关系，例如调控因子与其靶基因对应表达变化区段时间长度不同的调控关系。因此需要一种同时处理时间滞后和时间伸缩的调控关系的GRN新方法。
发明内容
本发明的目的是提供在一定程度上解决基因调控网络中时间滞后和时间伸缩问题的基于scRNA-seq及动态时间规整的基因调控网络构建方法。
为实现上述目的，本发明所用技术方案如下：
基于scRNA-seq及动态时间规整的基因调控网络构建方法，包括以下步骤：步骤1、基于单细胞数据的单细胞伪时间分析；步骤2、基于正负调控方式的双重DTW距离分析；步骤3、基于DTW距离的聚类分析以及GRN模块筛选；步骤4、基于DTW距离的基因间相似性矩阵构建；步骤5、通过随机置换检验，计算出基因与基因间网络连接度阈值C；步骤6、对GRN模块进行基因功能性分析；步骤7、pGRN网络的可视化得出GRN网络图。
进一步作为本发明的优选技术方案，步骤1的具体步骤如下：步骤1.1、分析单细胞转录组测序数据；步骤1.2、对于不同阶段的细胞赋予不同的伪时间值。
进一步作为本发明的优选技术方案，步骤2的具体步骤如下：步骤2.1、设有两组数组如下：
Xm＝(x1，x2，x3，…，xm)
Ym＝(y1，y2，y3，…，ym)
其中，数组Xm与Ym均为时间序列数组；
步骤2.2、定义DTW(m,n)为Xm与Yn之间最小的规整距离，通过如下动态规划算法：
DTW(i，j)＝Dist(i，j)+min[DTW(i-1，j)，DTW(i，j-1)，DTW(i-1，j-1)]
DTW(1，1)＝Dist(1，1)
DTW(i，1)＝Dist(i，1)+DTW(i-1，1)if(i>1)
DTW(1，j)＝Dist(1，j)+DTW(1，j-1)if(j>1)
其中，Dist(i，j)是Xi与Yj的DTW规整路径中所有对应于这两个数组中的点的欧几里德距离之和；
步骤2.3、设有基因A与基因B，随着伪时间表达量变化的数组为XA0与XB0，对原始数组做时间窗口平均取值得出XA与XB，XA与XB分别代表基因A与基因B叠加平移窗口内基因表达量的均值数组，其长度为窗口数量N；
步骤2.4、对表达数组XA与XB进行Z-score标准化，即为：
X′A＝(XA-Mean(XA))/(sd(XA))
X′B＝(XB-Mean(XB))/(sd(XB))
其中，X′A与X′B分别为XA与XB的Z-score标准化数组；
步骤2.5、对于基因A与基因B，计算双重DTW距离，即为：
D(A，B)＝dtw(X′A，X′B)/L
D(A，-B)＝dtw(X′A，-X′B)/L
其中，dtw()函数计算DTW距离，L为DTW匹配路径长度；
步骤2.6、将D(A，B)与D(A,-B)进行距离数值大小对比；当D(A，B)＞D(A,-B)，则调控关系为负值-1；当D(A，B)<D(A,-B)，则调控关系为正值+1；当D(A,B)＝D(A,-B)，则调控关系为零值0；步骤2.7、定义DTW(A，B)＝min(D(A,B)，D(A,-B))。
进一步作为本发明的优选技术方案，步骤3的具体步骤如下：步骤3.1、基于基因间DTW距离信息，对细胞的基因进行层次聚类分析；步骤3.2、设定聚类模块内最小基因数目阈值，使用dynamicTreeCut方法找出层次聚类树中的模块。
进一步作为本发明的优选技术方案，步骤4的具体步骤如下：步骤4.1、对于基因A和基因B；当DTW(A,B)距离数值高，则表示基因A和基因B随伪时间表达的变化趋势相关性弱；当DTW(A,B)距离数值低，则表示基因A和基因B随伪时间表达的变化趋势相关性强；步骤4.2、设定基因A和基因B之间的相似性数值为:
S(A，B)＝(1-DTW(A，B))k
其中，k取值为2到5之间的整数,S(A,B)代表了pGRN网络中基因A与基因B间的连接程度；
步骤4.3、将所有表达基因间的相似性数值构成基因表达相似性矩阵。
进一步作为本发明的优选技术方案，步骤5的具体步骤如下：步骤5.1、通过随机置换基因表达在伪时间中的顺序，采用同样的双重DTW分析方法，得到随机序列的DTW距离矩阵；
步骤5.2、基于随机序列的DTW距离矩阵，得到基因表达相似性矩阵即基因与基因间连接度矩阵；步骤5.3、设定一个基因与基因间连接度的阈值C，通过公式计算出基因调控网络的假阳性率FDR值，该公式为：
其中，M为随机置换的次数，N为基因相似表达矩阵中基因的个数，为基于真实数据在阈值C下的网络有效连接数，而/>为基于随机置换数据得到的在阈值C下的网络有效连接数；I为指示函数；当S(x,y)>C为真时，其数值为1；当S(x,y)>C为假时，其数值为0；步骤5.4、在一定范围内选取C值计算FDR随着C值变化的情况，进而确定在给定FDR阈值下基因与基因间连接度阈值C的取值。
进一步作为本发明的优选技术方案，步骤5.3与步骤5.4中基因调控网络的假阳性率FDR值为0.05。
进一步作为本发明的优选技术方案，步骤6的具体步骤如下：步骤6.1、筛选出GRN模块；步骤6.2、对于筛选出的GRN模块内部的基因进行GO富集性分析；步骤6.3、pGRN模块将整合topGO进行GO富集性分析。
进一步作为本发明的优选技术方案，步骤7的具体步骤如下：步骤7.1、基于基因相似表达性矩阵和模块信息以及RCy3，pGRN模块将通过R语句调用网络可视化工具Cytoscape；步骤7.2、得到各个模块的GRN网络图。
本发明的技术方案相对于现有技术的有益效果为：
本发明的pGRN能够发现时间调控滞后的调控关系，同时处理时间滞后和时间伸缩的调控关系，为新的调控关系的发掘提供了途径。pGRN通过应用DTW算法，能够处理时间伸缩的相似表达图谱查找，扩大了随时间的调控关系的查找范围，为进一步发现新的调控关系提供了可能。
附图说明
图1为本发明的流程示意图；
图2为本发明的单细胞测序数据聚类分析图；
图3为本发明的单细胞伪时间分析图；
图4为本发明的基因A与基因B基于DTW的矩阵比对图；
图5为本发明的基因A与基因B表达量随时间变化及其对于DTW比对点配对情况图；
图6为本发明的伪时间轴叠及平移窗口示意图；
图7为本发明的基于双重DTW距离的基因聚类分析图；
图8为本发明的pGRN基因调控网络模块图。
具体实施方式
下面结合附图对本发明做进一步的详细说明。
如图1所示，基于scRNA-seq及动态时间规整的基因调控网络构建方法，包括以下步骤：
步骤1、基于单细胞数据的单细胞伪时间分析；步骤2、基于正负调控方式的双重DTW距离分析；步骤3、基于DTW距离的聚类分析以及GRN模块筛选；步骤4、基于DTW距离的基因间相似性矩阵构建；步骤5、通过随机置换检验，计算出基因与基因间网络连接度阈值C；步骤6、对GRN模块进行基因功能性分析；步骤7、pGRN网络的可视化得出GRN网络图。
步骤1的具体步骤如下：步骤1.1、分析单细胞转录组测序数据；步骤1.2、对于不同阶段的细胞赋予不同的伪时间值。其中，单细胞测序数据来源于处在生物学动态过程的一群细胞，这些细胞处在了生物学动态过程中的不同阶段。通过单细胞测序数据的分析，如图2所示；能够对不同阶段的细胞赋以不同的伪时间值，如图3所示。伪时间值越小代表该细胞类型处在生物学过程的初始阶段，越大代表该细胞类型处在生物学过程的终端阶段。具体的，整合Monocle2、Monocle3、SLICER、DPT、Wishbone以及URD等伪时间分析方法，比较不同伪时间算法得到的伪时间构建的pGRN网络的差异性，并选用不同的伪时间方法构建pGRN基因调控网络。
步骤2的具体步骤如下：步骤2.1、设有两组数组如下：
Xm＝(x1，x2，x3，…，xm)
Ym＝(y1，y2，y3，…，ym)
其中，数组Xm与Ym均为时间序列数组；
步骤2.2、定义DTW(m,n)为Xm与Yn之间最小的规整距离，通过如下动态规划算法：
DTW(i，j)＝Dist(i，j)+min[DTW(i-1，j)，DTW(i，j-1)，DTW(i-1，j-1)]
DTW(1，1)＝Dist(1，1)
DTW(i，1)＝Dist(i，1)+DTW(i-1，1)if(i>1)
DTW(1，j)＝Dist(1，j)+DTW(1，j-1)if(j>1)
其中，Dist(i，j)是Xi与Yj的DTW规整路径中所有对应于这两个数组中的点的欧几里德距离之和；
步骤2.3、设有基因A与基因B，随着伪时间表达量变化的数组为XA0与XB0，对原始数组做时间窗口平均取值得出XA与XB，XA与XB分别代表基因A与基因B叠加平移窗口内基因表达量的均值数组，其长度为窗口数量N；
步骤2.4、对表达数组XA与XB进行z-score标准化，即为：
X′A＝(XA-Mean(XA))/(sd(XA))
X′B＝(XB-Mean(XB))/(sd(XB))
其中，X′A与X′B分别为XA与XB的z-score标准化数组；
步骤2.5、对于基因A与基因B，计算双重DTW距离，即为：
D(A，B)＝dtw(XA，XB)/L
D(A，-B)＝dtw(X′A，-X′B)/L
其中，dtw()函数计算DTW距离，L为DTW匹配路径长度；
步骤2.6、将D(A,B)与D(A,-B)进行距离数值大小对比；当D(A,B)>D(A,-B)，则调控关系为负值-1；当D(A,B)<D(A,-B)，则调控关系为正值+1；当D(A,B)＝D(A,-B)，则调控关系为零值0；步骤2.7、定义DTW(A,B)＝min(D(A,B),D(A,-B))。
其中，双重DTW距离分析是分别考虑了基因间正调控和负调控的DTW距离分析。DTW算法作为距离计算方法被广泛应用于时间序列数据的处理。步骤2.2中，两数组对应的DTW规整路径示例，如图4所示；通过规则路径，可以对两个时间序列数组对应的随时间变化的曲线做一一匹配，即使在两个相似时间变化曲线有时间滞后位移或者有伸缩的情况下，仍然能够找到两者之间的匹配关系，如图5所示。步骤2.3中，沿着伪时间轴从头到尾依次有叠加的设置时间窗口，然后再对窗口内的细胞表达量取均值，如图6所示。
步骤3的具体步骤如下：步骤3.1、基于基因间DTW距离信息，对细胞的基因进行层次聚类分析；步骤3.2、设定聚类模块内最小基因数目阈值，使用dynamicTreeCut方法找出层次聚类树中的模块。如图7所示，为基于双重DTW距离的基因聚类分析图。
步骤4的具体步骤如下：步骤4.1、对于基因A和基因B；当DTW(A,B)距离数值高，则表示基因A和基因B随伪时间表达的变化趋势相关性弱；当DTW(A,B)距离数值低，则表示基因A和基因B随伪时间表达的变化趋势相关性强；步骤4.2、设定基因A和基因B之间的相似性数值为:
S(A，B)＝(1-DTW(A，B))k
其中，k取值为2到5之间的整数,S(A,B)代表了pGRN网络中基因A与基因B间的连接程度；
步骤4.3、将所有表达基因间的相似性数值构成基因表达相似性矩阵。
步骤5的具体步骤如下：步骤5.1、通过随机置换基因表达在伪时间中的顺序，采用同样的双重DTW分析方法，得到随机序列的DTW距离矩阵；
步骤5.2、基于随机序列的DTW距离矩阵，得到基因表达相似性矩阵即基因与基因间连接度矩阵；步骤5.3、设定一个基因与基因间连接度的阈值C，通过公式计算出基因调控网络的假阳性率FDR值，该公式为：
其中，M为随机置换的次数，N为基因相似表达矩阵中基因的个数，为基于真实数据在阈值C下的网络有效连接数，而/>为基于随机置换数据得到的在阈值C下的网络有效连接数；I为指示函数；当S(x,y)>C为真时，其数值为1；当S(x,y)>C为假时，其数值为0；步骤5.4、在一定范围内选取C值计算FDR随着C值变化的情况，进而确定在给定FDR阈值下基因与基因间连接度阈值C的取值。步骤5.3与步骤5.4中基因调控网络的假阳性率FDR值为0.05。
步骤6的具体步骤如下：步骤6.1、筛选出GRN模块；步骤6.2、对于筛选出的GRN模块内部的基因进行GO富集性分析；步骤6.3、pGRN模块将整合topGO进行GO富集性分析。同时也可以根据实际需要，使用其它功能富集分析软件对pGRN模块基因进行分析，比如DAVID在线基因功能富集性分析网站以及Broad Institute的GSEA软件。
步骤7的具体步骤如下：步骤7.1、基于基因相似表达性矩阵和模块信息以及RCy3，pGRN模块将通过R语句调用网络可视化工具Cytoscape；步骤7.2、得到各个模块的GRN网络图。其中，RCy3即为R包，各个模块的GRN网络图，如图8所示。
具有生物学动态过程细胞的器官，如人类睾丸，在对人类睾丸组织的单细胞测序当中，能够获取处在精子发生过程不同阶段的细胞类型，有在精子发生过程时间上游的精原细胞、初级精母细胞，以及在精子发生下游的成熟精子细胞。pGRN能够对不同阶段的细胞赋以不同的伪时间值。结合单细胞转录组测序数据特征以及动态时间规整算法，更好的处理有时间滞后以及有调控区域时间伸缩的调控方式。将pGRN应用到异常精子发生小鼠模型中，期望发现与之相关的基因调控网络改变情况，并期望鉴别出关键的调控因子。将pGRN应用到正常人精子发生过程调控机理研究中，期望发现与之相关的新的基因调控模块。有望为男性不育提供基因调控网络层面的理论基础。
上述是为便于该技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对这些实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于这里的实施例，本领域技术人员根据本发明的揭示，不脱离本发明范畴所做出的改进和修改都应该在本发明的保护范围之内。

Claims (8)

1.基于scRNA-seq及动态时间规整的基因调控网络构建方法，其特征在于，包括以下步骤：
步骤1、基于单细胞数据的单细胞伪时间分析；
步骤2、基于正负调控方式的双重DTW距离分析；
步骤3、基于DTW距离的聚类分析以及GRN模块筛选；
步骤4、基于DTW距离的基因间相似性矩阵构建；
步骤5、通过随机置换检验，计算出基因与基因间网络连接度阈值C；
步骤6、对GRN模块进行基因功能性分析；
步骤7、pGRN网络的可视化得出GRN网络图；
步骤2的具体步骤如下：
步骤2.1、设有两组数组如下：
Xm＝(x1，x2，x3，…，xm)
Ym＝(y1，y2，y3，…，ym)
其中，数组Xm与Ym均为时间序列数组；
步骤2.2、定义DTW(m,n)为Xm与Yn之间最小的规整距离，通过如下动态规划算法:
DTW(i，j)＝Dist(i，j)+min[DTW(i-1，j)，DTW(i,j-1)，DTW(i-1,j-1)]
DTW(1，1)＝Dist(1，1)
DTW(i，1)＝Dist(i，1)+DTW(i-1，1)if(i＞1)
DTW(1，j)＝Dist(1，j)+DTW(1，j-1)if(j＞1)
其中，Dist(i,j)是Xi与Yj的DTW规整路径中所有对应于这两个数组中的点的欧几里德距离之和；
步骤2.3、设有基因A与基因B,随着伪时间表达量变化的数组为XA0与XB0，对原始数组做时间窗口平均取值得出XA与XB，XA与XB分别代表基因A与基因B叠加平移窗口内基因表达量的均值数组，其长度为窗口数量N；
步骤2.4、对表达数组XA与XB进行z-score标准化，即为：
X′A＝(XA-Mean(XA))/(sd(XA))
X′B＝(XB-Mean(XB))/(sd(XB))
其中，X'A与X'B分别为XA与XB的z-score标准化数组；
步骤2.5、对于基因A与基因B,计算双重DTW距离，即为：
D(A，B)＝dtw(X′A，X′B)/L
D(A，-B)＝dtw(X′A，-X′B)/L
其中，dtw()函数计算DTW距离，L为DTW匹配路径长度；
步骤2.6、将D(A,B)与D(A,-B)进行距离数值大小对比；当D(A,B)>D(A,-B)，则调控关系为负值-1；当D(A,B)<D(A,-B)，则调控关系为正值+1；当D(A,B)＝D(A,-B)，则调控关系为零值0；
步骤2.7、定义DTW(A,B)＝min(D(A,B),D(A,-B))。
2.根据权利要求1所述的基于scRNA-seq及动态时间规整的基因调控网络构建方法，其特征在于，步骤1的具体步骤如下：
步骤1.1、分析单细胞转录组测序数据；步骤1.2、对于不同阶段的细胞赋予不同的伪时间值。
3.根据权利要求1所述的基于scRNA-seq及动态时间规整的基因调控网络构建方法，其特征在于，步骤3的具体步骤如下：
步骤3.1、基于基因间DTW距离信息，对细胞的基因进行层次聚类分析；
步骤3.2、设定聚类模块内最小基因数目阈值，使用dynamicTreeCut方法找出层次聚类树中的模块。
4.根据权利要求1所述的基于scRNA-seq及动态时间规整的基因调控网络构建方法，其特征在于，步骤4的具体步骤如下：
步骤4.1、对于基因A和基因B，当DTW(A,B)距离数值高，则表示基因A和基因B随伪时间表达的变化趋势相关性弱；当DTW(A,B)距离数值低，则表示基因A和基因B随伪时间表达的变化趋势相关性强；
步骤4.2、设定基因A和基因B之间的相似性数值为:
S(A，B)＝(1-DTW(A，B))k
其中，k取值为2到5之间的整数,S(A,B)代表了pGRN网络中基因A与基因B间的连接程度；
步骤4.3、将所有表达基因间的相似性数值构成基因表达相似性矩阵。
5.根据权利要求1所述的基于scRNA-seq及动态时间规整的基因调控网络构建方法，其特征在于，步骤5的具体步骤如下：
步骤5.1、通过随机置换基因表达在伪时间中的顺序，采用同样的双重DTW分析方法，得到随机序列的DTW距离矩阵；
步骤5.2、基于随机序列的DTW距离矩阵，得到基因表达相似性矩阵即基因与基因间连接度矩阵；
步骤5.3、设定一个基因与基因间连接度的阈值C，通过公式计算出基因调控网络的假阳性率FDR值，该公式为：
其中，M为随机置换的次数，N为基因相似表达矩阵中基因的个数，为基于真实数据在阈值C下的网络有效连接数，而/>为基于随机置换数据得到的在阈值C下的网络有效连接数；I为指示函数；当S(x,y)>C为真时，其数值为1；当S(x,y)>C为假时，其数值为0；
步骤5.4、在一定范围内选取C值计算FDR随着C值变化的情况，进而确定在给定FDR阈值下基因与基因间连接度阈值C的取值。
6.根据权利要求5所述的基于scRNA-seq及动态时间规整的基因调控网络构建方法，其特征在于，步骤5.3与步骤5.4中基因调控网络的假阳性率FDR值为0.05。
7.根据权利要求1所述的基于scRNA-seq及动态时间规整的基因调控网络构建方法，其特征在于，步骤6的具体步骤如下：
步骤6.1、筛选出GRN模块；
步骤6.2、对于筛选出的GRN模块内部的基因进行GO富集性分析；
步骤6.3、pGRN模块将整合topGO进行GO富集性分析。
8.根据权利要求1所述的基于scRNA-seq及动态时间规整的基因调控网络构建方法，其特征在于，步骤7的具体步骤如下：
步骤7.1、基于基因相似表达性矩阵和模块信息以及RCy3，pGRN模块将通过R语句调用网络可视化工具Cytoscape；
步骤7.2、得到各个模块的GRN网络图。
CN201911012024.8A 2019-10-23 2019-10-23 基于scRNA-seq及动态时间规整的基因调控网络构建方法 Active CN110808083B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911012024.8A CN110808083B (zh) 2019-10-23 2019-10-23 基于scRNA-seq及动态时间规整的基因调控网络构建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911012024.8A CN110808083B (zh) 2019-10-23 2019-10-23 基于scRNA-seq及动态时间规整的基因调控网络构建方法

Publications (2)

Publication Number Publication Date
CN110808083A CN110808083A (zh) 2020-02-18
CN110808083B true CN110808083B (zh) 2023-08-29

Family

ID=69488954

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911012024.8A Active CN110808083B (zh) 2019-10-23 2019-10-23 基于scRNA-seq及动态时间规整的基因调控网络构建方法

Country Status (1)

Country Link
CN (1) CN110808083B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112768001A (zh) * 2021-01-27 2021-05-07 湖南大学 一种基于流形学习和主曲线的单细胞轨迹推断方法
CN113611368B (zh) * 2021-07-26 2022-04-01 哈尔滨工业大学（深圳） 基于2d嵌入的半监督单细胞聚类方法、装置、计算机设备
CN116129992A (zh) * 2023-04-17 2023-05-16 之江实验室 基于图神经网络的基因调控网络构建方法及系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109979538A (zh) * 2019-03-28 2019-07-05 广州基迪奥生物科技有限公司 一种基于10x单细胞转录组测序数据的分析方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2007110707A2 (en) * 2005-10-12 2007-10-04 Gni, Ltd. Methods and systems for modeling gene networks

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109979538A (zh) * 2019-03-28 2019-07-05 广州基迪奥生物科技有限公司 一种基于10x单细胞转录组测序数据的分析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
薛劼等.一种动态时间弯曲距离的时延调控基因相似度量聚类方法.福州大学学报.2013,第41卷(第2期),158-163. *

Also Published As

Publication number Publication date
CN110808083A (zh) 2020-02-18

Similar Documents

Publication Publication Date Title
CN110808083B (zh) 基于scRNA-seq及动态时间规整的基因调控网络构建方法
Huang et al. Large-scale regulatory network analysis from microarray data: modified Bayesian network learning and association rule mining
CN105117618B (zh) 一种基于牵连犯罪原则和网络拓扑结构特征的药物‑靶标相互作用识别方法
CN107885971B (zh) 采用改进花授粉算法识别关键蛋白质的方法
CN108682457A (zh) 患者长期预后定量预测和干预系统及方法
CN112599187B (zh) 一种基于双流神经网络预测药物与靶标蛋白结合分数的方法
CN111951886A (zh) 一种基于贝叶斯归纳型矩阵补全的药物重定位预测方法
CN108197431B (zh) 染色质相互作用差异的分析方法和系统
CN105590039B (zh) 一种基于bso优化的蛋白质复合物识别方法
Seo et al. Cancer-related gene signature selection based on boosted regression for multilayer perceptron
CN116580848A (zh) 一种基于多头注意力机制的分析癌症多组学数据方法
Peres da Silva et al. TUGDA: task uncertainty guided domain adaptation for robust generalization of cancer drug response prediction from in vitro to in vivo settings
Stegmayer et al. A biologically inspired validity measure for comparison of clustering methods over metabolic data sets
Couto et al. Effects of threshold on the topology of gene co-expression networks
CN109671468B (zh) 一种特征基因选择及癌症分类方法
Bilen et al. A new hybrid and ensemble gene selection approach with an enhanced genetic algorithm for classification of microarray gene expression values on leukemia cancer
Rezaeenour et al. Prediction of protein thermostability by an efficient neural network approach
Jubair et al. Gxenet: Novel fully connected neural network based approaches to incorporate gxe for predicting wheat yield
CN107526937A (zh) 一种基于协同滤波的MiRNA‑疾病关联预测方法
Mendoza-Revilla et al. A Foundational Large Language Model for Edible Plant Genomes
Yu et al. Predicting phenotypes from high-dimensional genomes using gradient boosting decision trees
Fraiman et al. Statistics of dynamic random networks: a depth function approach
Mondal et al. Clustering Single-Cell RNA Sequence Data Using Information Maximized and Noise-Invariant Representations
Kojima et al. Identifying regulational alterations in gene regulatory networks by state space representation of vector autoregressive models and variational annealing
Shamaiah et al. Graphical models and inference on graphs in genomics: challenges of high-throughput data analysis

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