CN110211634A - 一种多组学数据联合分析的方法 - Google Patents

一种多组学数据联合分析的方法 Download PDF

Info

Publication number
CN110211634A
CN110211634A CN201810111865.3A CN201810111865A CN110211634A CN 110211634 A CN110211634 A CN 110211634A CN 201810111865 A CN201810111865 A CN 201810111865A CN 110211634 A CN110211634 A CN 110211634A
Authority
CN
China
Prior art keywords
value
data
module
group
multiple groups
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.)
Granted
Application number
CN201810111865.3A
Other languages
English (en)
Other versions
CN110211634B (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.)
BGI Technology Solutions Co Ltd
Original Assignee
BGI Technology Solutions Co Ltd
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 BGI Technology Solutions Co Ltd filed Critical BGI Technology Solutions Co Ltd
Priority to CN201810111865.3A priority Critical patent/CN110211634B/zh
Publication of CN110211634A publication Critical patent/CN110211634A/zh
Application granted granted Critical
Publication of CN110211634B publication Critical patent/CN110211634B/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
    • 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
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics

Landscapes

  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Bioethics (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Public Health (AREA)
  • Software Systems (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种多组学数据联合分析的方法。本发明所提供的多组学数据联合分析的方法包括如下步骤：(A)对待分析的多组学数据中的每个单一组学的指标数据进行共表达网络分析，找到各自的表达模块；(B)根据不同组学数据各自的表达模块之间的重叠关系，筛选出所述待分析的多组学数据之间显著相关的互作模块。本发明所提供的多组学数据联合分析方法不受组学数据的数量限制，理论上可以是任意多个组；同时也不依赖输入数据的来源，只要是能衡量相应组学的指标数据(如基因表达量值、表观的甲基化程度，SNP突变率等)，都可以作为输入数据。

Description

一种多组学数据联合分析的方法
技术领域
本发明涉及生物信息学领域，具体涉及一种多组学数据联合分析的方法。
背景技术
随着科学技术的不断进步，高通量的组学数据开始变得容易获取，他们提供了细胞中几乎所有的成员和相互作用的综合描述。Joyce等将这些数据分成3类：成员、相互作用和功能状态数据。成员数据描述细胞分子的属性；相互作用数据记录分子成员之间的作用关系；功能状态数据指的是整体的细胞表型，揭示所有组学数据作用的整体表现。已有的组数数据描述了从基因组到代谢组的生物信号流。首先，DNA(基因组)被转录为mRNA(转录组)，然后mRNA翻译为蛋白质(蛋白质组)，蛋白质催化反应生成代谢物、糖蛋白和寡糖，以及不同的脂类(脂类组)。其中大部分成员可以在细胞中标记和定位(定位组)。产生和改变这些细胞成分的过程通常取决于分子相互作用(相互作用组)，如转录过程中的蛋白质-DNA相互作用、翻译后的蛋白质相互作用以及酶相互作用等。最后，由代谢通路组成整合的网络或流量图(代谢组)，决定细胞动物或表型(表型组)。
对于复杂的生物过程，单一的组学研究已经很难对其进行深入的解读。多组学数据的分析，为挖掘科学研究热点，提供了新的思路。然而目前，还没有一套成熟的分析方法对多组学进行联合分析，尤其是超过2个组学的数据的分析。
目前，常用的多组学数据的联合分析方法，一般是基于位置关系，像lncRNA与mRNA之间位置关系、miRNA与基因之间的靶向关系、基因与组蛋白结合位置等，或者基于基因表达量的相关性，计算两个组学数据之间的相互关系。
现有的分析方法，通用性较差，不同的项目，可能都需要调整方法、参数来分析；方法局限性比较高，对于超过两个组学的数据，很难满足其科研需求。
发明内容
为了有效的解决以上问题，本发明开发了一套适用于任意多组的组学数据的联合分析方法。
本发明所提供的一种多组学数据联合分析的方法，具体可包括如下步骤：
(A)对待分析的多组学数据中的每个单一组学的指标数据进行共表达网络分析，找到各自的表达模块；
(B)根据步骤(A)得到的不同组学数据各自的表达模块之间的重叠关系，筛选出所述待分析的多组学数据中的互作模块。
其中，这里的所述“共表达网络”是基于某单一组学数据中不同指标的相似性而构建的网络图，图中的节点代表指标，具有共性的指标被连接起来形成网络。所述“共表达网络”中被连接起来的一组具有共性的指标即为一个所述“表达模块”。如“基因共表达网络”，是基于基因间表达数据的相似性而构建的网络图，图中的节点代表基因(此处的基因即对应所述指标)，具有相似表达谱的基因被连接起来形成网络。对于所述“基因共表达网络”而言，所述“表达模块”则为所述“基因共表达网络”中被连接起来的一组具有类似表达趋势的基因，如果某些基因在一个生理过程或不同组织中总是具有相类似的表达变化，那么则有理由认为这些基因在功能上是相关的，把它们定义为一个表达模块(module)。
这里所述的“重叠关系”有两层含义，其一是说不同组学数据的表达模块之间共有的指标数，其二是说这个共有的情况是否显著(基于显著模型判断)。
在所述方法中，步骤(A)之前还包括获得所述待分析多组学数据的步骤。
进一步地，步骤(A)中，是基于聚类方法对所述单一组学的指标数据进行共表达网络分析，从而找到所述单一组学的指标数据的所述表达模块的。如可采用加权基因共表达网络构建(Weighted Gene Co-Expression Network Analysis,WGCNA)的方法对所述单一组学的指标数据进行共表达网络分析，从而找到所述单一组学的指标数据的所述表达模块的。
进一步地，步骤(B)中，是利用超几何分布检验的方法根据不同组学数据各自的所述表达模块之间的重叠关系，筛选出所述待分析的多组学数据之间显著相关的互作模块。
更进一步地，在本发明中，步骤(A)具体是按照包括如下步骤的方法对所述单一组学的指标数据进行共表达网络分析，从而找到所述单一组学的指标数据的所述表达模块的：
(a1)计算所述单一组学的指标数据中任意两个指标之间的Pearson相关性系数(皮尔逊相关系数)，得到相关性系数矩阵；
(a2)按照无尺度网络的标准选择邻接矩阵的权重参数β值；
(a3)根据步骤(a2)中得到的β值，计算步骤(a1)中相关性系数矩阵的邻接矩阵；
(a4)用1减去步骤(a3)中的邻接矩阵所得数值作为距离，构建系统聚类树；然后根据混合动态剪切树(dynamicTreeCut)确定分类模块，即得到所述单一组学的指标数据的所述表达模块。
更加具体的，步骤(a1)中，用于计算任意两个指标之间的Pearson相关性系数的所述单一组学的指标数据为经过预处理后得到预处理后数据。所述预处理包括删除缺失率高于设定阈值的数据(删除缺失率高于设定阈值的整行数据)。在本发明的一个实施例中，所述设定阈值具体为0.2(表示所述设定阈值为“缺失率为20％”)，即删除缺失率高于20％的整行数据。
其中，所述“缺失率”指某个指标数据在所统计的样本中，被漏记的比例；譬如10个样本，只有8个样品有该指标数据的数值，则缺失率为20％。
更加具体的，步骤(a1)中，是按照包括如下步骤的方法得到所述相关性系数矩阵的：当所述“共表达网络”的网路类型为unsigned时，指标i和指标j之间的Pearson相关性系数Sij＝|cor(i,j)|；当所述“共表达网络”的网路类型为signed时，指标i和指标j之间的Pearson相关性系数Sij＝|(1+cor(i,j))/2|；由此，得到所述相关性系数矩阵S＝[Sij]。
更加具体的，步骤(a2)中，是按照包括如下步骤的方法确定所述β值的：
(i)β值分别取1到30的正整数，分别计算各β值对应的R2值；所述R2值是共表达网络中某节点连接度(k)的对数(log(k))和该节点出现的概率的对数(log(p(k)))之间的相关性系数的平方值。所述R2值越大，网络越逼近无网络尺度的分布。
(ii)按照如下确定所述β值：
第一种情况：如果存在所述R2值大于等于0.9，则取第一个出现的所述R2值大于等于0.9时对应的β值(即所述R2值大于等于0.9时对应的最小的β值)；
第二种情况：如果所有的所述R2值都小于0.9，则将所述R2值和所述β值进行局部多项式回归分析，取第一个饱和点的所述R2值对应的β值；
第三种情况：如果不存在满足所述第一种情况和所述第二种情况的β值，则取β值为30。
更加具体的，步骤(a3)中，是按照包括如下步骤的方法计算得到所述邻接矩阵的：根据步骤(a1)中计算得到的指标i和指标j之间的Pearson相关性系数Sij以及步骤(a2)中得到的β值，计算得到指标i和指标j之间的邻接系数aij＝|Sij|β；由此，得到所述邻接矩阵A＝[aij]。
在本发明的一个具体实施例中，步骤(a4)中，根据混合动态剪切树确定所述分类模块时定义每个模块中指标的最少数目为50；deepSplit为0。
更进一步地，在本发明中，步骤(B)中，具体是按照如下方法筛选出所述待分析的多组学数据之间显著相关的互作模块的：
(b1)假定所述待分析的多组学数据为X个组学数据，来自于所述X个组学数据的X个表达模块组合成模块组(不失一般性，X为大于等于2的正整数；所述模块组中包含的表达模块分别记为M1，M2，…，Mx)；利用超几何分布检验的方法，按照如下所示计算所述模块组之间相互关联的显著性P值；
假定所述模块组中的表达模块Mi(i属于1到X的正整数)包含的指标数量为M，所述模块组中所有表达模块之间的交集的指标数量为m，所述模块组中除了表达模块Mi之外的所有表达模块(即M1，…，Mi-1，Mi+1，…，Mx)的并集的指标数量为n，分析时输入的所有的指标数量为N，则对应的显著性P值为：
(b2)分别以所述模块组中的表达模块M1，M2，…，Mx替代步骤(b1)中的表达模块Mi，进行步骤(b1)，共得到X个P值，从中选择最小的P值；
(b3)如果所述m大于等于设定阈值，且步骤(b2)中所得的所述最小的P值小于等于设定阈值，则视所述模块组之间存在显著相关关系(同时输出m个指标数据列表，用于后续互作网络图等分析)，所述模块组即为所述待分析的多组学数据之间显著相关的互作模块。
进一步地，在本发明的一个实施例中，所述m的设定阈值具体为10；所述P值的设定阈值具体为0.05。
其中，所述m和所述P值的阈值设定的一般原则是根据物种的基因组注释情况和组学数据的类别来设置，譬如物种注释很全面，可以适当严格些，保留相对更可信的结果；组学数据中有miRNA时，目前来说miRNA的靶点很多，那么建立的链接可能就会比较多，相应的指标也可以严格些。当然，也可以设置宽松些，然后根据结果自行过滤。再有，如果基因组注释比较差，有可能出现没有显著差异的结果；此时则需要降低交集的判读阈值(所述m和所述P值)。这些是本领域技术人员可以根据实际情况进行把控的。
根据需要，所述方法在步骤(B)之后还可包括如下步骤(C1)或(C2)：
(C1)对构成所述互作模块的各表达模块之间的交集部分的指标数据进行功能富集分析(包括Geno Ontology、Pathway富集)，得到显著富集的功能条目。
其中，GO是Gene ontology的缩写，GO数据库分别从功能、参与的生物途径及细胞中的定位对基因产物进行了标准化描述，即对基因产物进行简单注释，通过GO富集分析可以粗略了解目标基因富集在哪些生物学功能、途径或者细胞定位。Pathway指代谢通路，对目标基因进行pathway分析，可以了解实验条件下显著改变的代谢通路，在机制研究中显得尤为重要。
(C2)根据所述互作模块的关联信息，构建有权重的互作网络图。这样可以直观的展示组学数据之间相互作用的关系，也可以直观挑选一些影响比较大的基因以便后续做进一步深入研究。
进一步地，步骤(C2)中，是利用visNetwork构建的所述有权重的互作网络图，其中权重信息为各指标数据的平均值。
本发明还提供了一种能够用于多组学数据联合分析的系统。
本发明所提供的能够用于多组学数据联合分析的系统，包括基本组成单元A和基本组成单元B；所述基本组成单元A具有能够实现前文所述方法中的步骤(A)的功能；所述基本组成单元B具有能够实现前文所述方法中的步骤(B)的功能。
另外，本发明还要求保护如下任一应用：
(1)所述系统在多组学数据联合分析中的应用；
(2)前文所述方法在农业育种、组织发育和药物靶点筛选等各项研究中的应用。
本发明所提供的多组学数据联合分析方法不受组学数据的数量限制，理论上可以是任意多个组；同时也不依赖输入数据的来源，只要是能衡量相应组学的指标数据(如基因表达量值、表观的甲基化程度，SNP突变率等)，都可以作为输入数据。
在本发明的一个实施例中，所述多组学数据具体为mRNA组学数据、small RNA组学数据。在本发明的另一个实施例中，所述多组学数据具体为mRNA组学数据、small RNA组学数据和DNA甲基化组学数据。可以理解的，根据分析目的的不同，有针对性的选择生物样本和组学数据，例如，如果需要筛选药物靶点，可以获取病人和健康人的样本，通过分析转录组和/或蛋白组和/或甲基化组等数据，找到差异表达位点，可获得候选的药物靶点，最终通过验证实验确定药物靶点。
本发明的优势在于：(1)本发明可以适用于不同组学的数值(可转换为数值的)数据，进行联合分析；(2)本发明不限制组学数据的数量，理论上可以是任意多组；(3)本发明不仅从统计分类的角度给出了潜在的关联关系，同时还从生物功能、代谢过程的角度(通过功能富集分析)给出了潜在的关联因素。
附图说明
图1为实施例1中mRNAMblack_smRNAMbrown互作网络图。
图2为实施例2中mRNAMturquoise_smRNAMturquoise_methyMturquoise互作网络图。
具体实施方式
下述实施例中所使用的实验方法如无特殊说明，均为常规方法。
下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。
实施例1、大鼠mRNA和small RNA数据联合分析
一、实验样本
大鼠的6个不同的发育时间阶段的正常组织样本。
二、数据采集
采用常规的RNA-Seq和Small RNA标准信息分析，从大鼠的6个不同的发育时间阶段的正常组织样本中获得mRNA和small RNA的输入数据。
最终获得：mRNA的数据包括25289个基因的表达量，small RNA的数据是631个miRNA的表达量。miRNA与mRNA之间的靶向关系，是根据TargetScan和miRanda靶基因预测软件的预测结果，包含了76069对靶向关系。
三、数据预处理
删除缺失率高于设定阈值(为默认阈值0.2，即缺失率为20％)的数据(删除相应的行)。
四、共表达网络分析
按照如下步骤对步骤三预处理后的mRNA数据和small RNA数据分别进行共表达网络分析：
1、计算预处理后数据中任意两个指标之间的Pearson相关性系数(皮尔逊相关系数)，得到相关性系数矩阵：
设定“共表达网络”的网路类型为unsigned，指标i和指标j之间的Pearson相关性系数Sij＝|cor(i,j)|；由此，得到所述相关性系数矩阵S＝[Sij]。
2、按照无尺度网络的标准选择邻接矩阵的权重参数β值：
(i)β值分别取1到30的正整数，分别计算各β值对应的R2值；所述R2值是共表达网络中某节点连接度(k)的对数(log(k))和该节点出现的概率的对数(log(p(k)))之间的相关性系数的平方值。所述R2值越大，网络越逼近无网络尺度的分布。
(ii)按照如下确定所述β值：
第一种情况：如果存在所述R2值大于等于0.9，则取第一个出现的所述R2值大于等于0.9时对应的β值(即所述R2值大于等于0.9时对应的最小的β值)；
第二种情况：如果所有的所述R2值都小于0.9，则将所述R2值和所述β值进行局部多项式回归分析，取第一个饱和点的所述R2值对应的β值；
第三种情况：如果不存在满足所述第一种情况和所述第二种情况的β值，则取β值为30。
3、根据步骤2中得到的β值，计算步骤1中相关性系数矩阵的邻接矩阵：
根据步骤1中计算得到的指标i和指标j之间的Pearson相关性系数Sij以及步骤2中得到的β值，计算得到指标i和指标j之间的邻接系数aij＝|Sij|β；由此，得到所述邻接矩阵A＝[aij]。
4、用1减去步骤3中的邻接矩阵所得数值作为距离，构建系统聚类树；然后根据混合动态剪切树确定分类模块(定义每个模块中指标的最少数目为50；deepSplit为0)，即得到所述单一组学的指标数据的所述表达模块。
分析完成后，mRNA部分得到48个模块，small RNA部分得到6个模块。
五、筛选mRNA和small RNA之间的显著相关互作模块
1、利用超几何分布检验的方法，按照如下方法计算表达模块A(来自于mRNA数据的模块)和表达模块B(来自于small RNA数据的模块)的之间相互关联的显著性P值：
假定所述表达模块A包含的指标数量为M，所述表达模块A与所述表达模块B交集的指标数量为m，所述表达模块B的指标数量为n，分析时输入的所有的指标数量为N(N是指背景的总数，如无特殊说明，则N为该物种全部的基因数。注：不管涉及的是什么组学数据，最终都会按照基因为参照去转换，所以如无特殊说明，则N可以理解为是整个基因组所涉及的基因集)，则对应的所述显著性P值为：
2、如果所述m大于等于设定阈值(为默认阈值10)，且所述P值小于等于设定阈值(为默认阈值0.05)，则视所述表达模块A与所述表达模块B为显著相关的互作模块(输出所述表达模块A与所述表达模块B的组合信息)；反之，则视所述表达模块A与所述表达模块B不为显著相关的互作模块。
通过筛选，得到116个mRNA和small RNA显著关联的模块对，部分结果见表1：
表1显著关联的模块
注：表中第一列为模块名称，不同组学的模块用下划线分隔；第二列为显著性P值；第三列为模块之间的交集的指标数。
六、Geno Ontology富集
Geno Ontology富集是将目标转录本向GO数据库(http://www.geneontology.org/)的各条目(term)映射，并计算每个term的转录本数，从而得到具有某个GO功能的转录本列表及转录本数目统计。然后应用超几何检验，找出与整个转录本组背景相比，在目标转录本中显著富集的GO条目，该假设检验的p-value计算公式为：
式中，N为所有转录本中具有GO注释的转录本数目；n为N中目标转录本的数目；M为所有转录本中注释为某特定GO term的转录本数目；m为注释为某特定GOterm的目标转录本数目。计算得到的pvalue通过Bonferroni校正之后，以corrected-pvalue≤0.05为阈值，满足此条件的GO term定义为在目标转录本中显著富集的GO term。
对表1中第一个模块组合mRNAMblack_smRNAMbrown的交集基因(交集基因是基于miRNA的靶基因与mRNAMblack进行判断的)进行功能富集(GO数据库版本是201710(下载日期))，找到216条显著富集的条目，部分结果见表2。
表2功能富集条目
七、构建有权重的互作网络图
根据模块组合mRNAMblack_smRNAMbrown中的关联关系，用visNetwork构建有权重的互作网络图，其中权重信息为各指标数据的平均值。结果如图1所示。
通过本方案的分析，在mRNAMblack_smRNAMblue显著模块对中发现了很多与器官发育相关的基因，包括细胞发育(涉及25个基因)、神经系统发育(35个)、神经元投影的发展(17个)等，这为后续系统的研究大鼠的神经发育提供了精确的探索方向。
实施例2、苹果甲基化、mRNA和small RNA数据联合分析
一、实验样本
苹果的4个不同处理的样本。
二、数据采集
采用常规的RNA-Seq、全基因组Bisulfite甲基化和Small RNA标准信息分析，从苹果的4个不同样本中获得mRNA、甲基化和small RNA的输入数据。
最终获得：mRNA的数据包括31964个基因的表达量，small RNA的数据是167个miRNA的表达量，甲基化是34889个甲基化区域的甲基化率结果。miRNA与mRNA之间的靶向关系，是根据psRobot和targetfinder靶基因预测软件的预测结果，包含了9033对靶向关系。
三、数据预处理
删除缺失率高于设定阈值(为默认阈值0.2，即缺失率为20％)的数据(删除相应的行)。
四、共表达网络分析
按照如下步骤对步骤三预处理后的甲基化数据、mRNA数据和small RNA数据分别进行共表达网络分析：
1、计算预处理后数据中任意两个指标之间的Pearson相关性系数(皮尔逊相关系数)，得到相关性系数矩阵：
设定“共表达网络”的网路类型为unsigned，指标i和指标j之间的Pearson相关性系数Sij＝|cor(i,j)|；由此，得到所述相关性系数矩阵S＝[Sij]。
2、按照无尺度网络的标准选择邻接矩阵的权重参数β值：
(i)β值分别取1到30的正整数，分别计算各β值对应的R2值；所述R2值是共表达网络中某节点连接度(k)的对数(log(k))和该节点出现的概率的对数(log(p(k)))之间的相关性系数的平方值。所述R2值越大，网络越逼近无网络尺度的分布。
(ii)按照如下确定所述β值：
第一种情况：如果存在所述R2值大于等于0.9，则取第一个出现的所述R2值大于等于0.9时对应的β值(即所述R2值大于等于0.9时对应的最小的β值)；
第二种情况：如果所有的所述R2值都小于0.9，则将所述R2值和所述β值进行局部多项式回归分析，取第一个饱和点的所述R2值对应的β值；
第三种情况：如果不存在满足所述第一种情况和所述第二种情况的β值，则取β值为30。
3、根据步骤2中得到的β值，计算步骤1中相关性系数矩阵的邻接矩阵：
根据步骤1中计算得到的指标i和指标j之间的Pearson相关性系数Sij以及步骤2中得到的β值，计算得到指标i和指标j之间的邻接系数aij＝|Sij|β；由此，得到所述邻接矩阵A＝[aij]。
4、用1减去步骤3中的邻接矩阵所得数值作为距离，构建系统聚类树；然后根据混合动态剪切树确定分类模块(定义每个模块中指标的最少数目为50；deepSplit为0)，即得到所述单一组学的指标数据的所述表达模块。
分析完成后，mRNA部分得到67个模块，small RNA部分得到3个模块，甲基化部分得到12个模块。
五、筛选mRNA、甲基化和small RNA之间的显著相关互作模块
按照如下步骤筛选筛选mRNA、甲基化和small RNA之间的显著相关互作模块：
1、利用超几何分布检验的方法，按照如下方法计算表达模块A(来自于mRNA数据的模块)、表达模块B(来自于small RNA数据的模块)和表达模块C(来自于甲基化数据的模块)之间相互关联的显著性P值；
2、假定所述表达模块A包含的指标数量为M，所述表达模块A、表达模块B与表达模块C的交集的指标数量为m，所述表达模块B和表达模块C的并集指标数量为n，分析时输入的所有的指标数量为N(N是指背景的总数，如无特殊说明，则N为该物种全部的基因数。注：不管涉及的是什么组学数据，最终都会按照基因为参照去转换，所以如无特殊说明，则N可以理解为是整个基因组所涉及的基因集)，则对应的所述显著性P值为：
3、重复步骤2，将A、B和C的P值都计算出来，取最小的P值；
4、如果所述m大于等于设定阈值(为默认阈值10)，且步骤3所得的所述最小的P值小于等于设定阈值(为默认阈值0.05)，则视所述表达模块A、所述表达模块B和所述表达模块C之间存在显著相关关系；
5、输出交集m中所有模块信息用于后续分析
通过筛选，得到7个mRNA、甲基化和small RNA显著关联的模块对，结果见表3：
表3苹果显著关联的模块
注：表中第一列为模块名称，不同组学的模块用下划线分隔；第二列为显著性P值；第三列为模块之间的交集的指标数。
六、Geno Ontology富集
Geno Ontology富集是将目标转录本向GO数据库(http://www.geneontology.org/)的各条目(term)映射，并计算每个term的转录本数，从而得到具有某个GO功能的转录本列表及转录本数目统计。然后应用超几何检验，找出与整个转录本组背景相比，在目标转录本中显著富集的GO条目，该假设检验的p-value计算公式为：
式中，N为所有转录本中具有GO注释的转录本数目；n为N中目标转录本的数目；M为所有转录本中注释为某特定GO term的转录本数目；m为注释为某特定GOterm的目标转录本数目。计算得到的pvalue通过Bonferroni校正之后，以corrected-pvalue≤0.05为阈值，满足此条件的GO term定义为在目标转录本中显著富集的GO term。
对模块组合mRNAMturquoise_smRNAMturquoise_methyMturquoise的交集基因(交集基因是基于miRNA的靶基因与mRNAMblack进行判断的。注：甲基化部分，是按照基因区来统计的甲基化率；所以甲基化可以直接对应到基因。)进行功能富集(GO数据库版本是201710(下载日期))，找到5条显著富集的条目，结果见表4。
表4苹果功能富集条目
七、构建有权重的互作网络图
根据模块组合mRNAMturquoise_smRNAMturquoise_methyMturquoise中的关联关系，用visNetwork构建有权重的互作网络图，其中权重信息为各指标数据的平均值。结果如图2所示。

Claims (10)

1.一种多组学数据联合分析的方法，包括如下步骤：
(A)对待分析的多组学数据中的每个单一组学的指标数据进行共表达网络分析，找到各自的表达模块；
(B)根据步骤(A)得到的不同组学数据各自的表达模块之间的重叠关系，筛选出所述待分析的多组学数据之间显著相关的互作模块。
2.根据权利要求1所述的方法，其特征在于：步骤(A)中，是基于聚类方法对所述单一组学的指标数据进行共表达网络分析，从而找到所述单一组学的指标数据的所述表达模块的。
3.根据权利要求1或2所述的方法，其特征在于：步骤(B)中，是利用超几何分布检验的方法根据不同组学数据各自的所述表达模块之间的重叠关系，筛选出所述待分析的多组学数据中的互作模块的。
4.根据权利要求1-3中任一所述的方法，其特征在于：步骤(A)中，是按照包括如下步骤的方法对所述单一组学的指标数据进行共表达网络分析，从而找到所述单一组学的指标数据的所述表达模块的：
(a1)计算所述单一组学的指标数据中任意两个指标之间的Pearson相关性系数，得到相关性系数矩阵；
(a2)按照无尺度网络的标准选择邻接矩阵的权重参数β值；
(a3)根据步骤(a2)中得到的β值，计算步骤(a1)中相关性系数矩阵的邻接矩阵；
(a4)用1减去步骤(a3)中的邻接矩阵所得数值作为距离，构建系统聚类树；然后根据混合动态剪切树确定分类模块，即得到所述单一组学的指标数据的所述表达模块。
5.根据权利要求4所述的方法，其特征在于：步骤(a1)中，用于计算任意两个指标之间的Pearson相关性系数的所述单一组学的指标数据为经过预处理后得到预处理后数据；
进一步地，所述预处理包括删除缺失率高于设定阈值的数据；
更进一步地，所述设定阈值为0.2。
6.根据权利要求4或5所述的方法，其特征在于：步骤(a2)中，是按照包括如下步骤的方法确定所述β值的：
(i)β值分别取1到30的正整数，分别计算各β值对应的R2值；所述R2值是共表达网络中某节点连接度的对数和该节点出现的概率的对数之间的相关性系数的平方值；
(ii)按照如下确定所述β值：
第一种情况：如果存在所述R2值大于等于0.9，则取第一个出现的所述R2值大于等于0.9时对应的β值；
第二种情况：如果所有的所述R2值都小于0.9，则将所述R2值和所述β值进行局部多项式回归分析，取第一个饱和点的所述R2值对应的β值；
第三种情况：如果不存在满足所述第一种情况和所述第二种情况的β值，则取β值为30。
7.根据权利要求1-6中任一所述的方法，其特征在于：步骤(B)中，是按照如下方法筛选出所述待分析的多组学数据之间显著相关的互作模块的：
(b1)所述待分析的多组学数据为X个组学数据，X为大于等于2的正整数；来自于所述X个组学数据的X个表达模块组合成模块组；将所述模块组中所包含的X个表达模块分别记为M1，M2，…，Mx；然后利用超几何分布检验的方法，按照如下所示计算所述模块组之间相互关联的显著性P值；
所述模块组中的表达模块Mi包含的指标数量为M，所述模块组中所有表达模块之间的交集的指标数量为m，所述模块组中除了表达模块Mi之外的所有表达模块的并集的指标数量为n，分析时输入的所有的指标数量为N，则对应的显著性P值为：
(b2)分别以所述模块组中的表达模块M1，M2，…，Mx替代步骤(b1)中的所述表达模块Mi，进行步骤(b1)，共得到X个P值，从中选择最小的P值；
(b3)如果所述m大于等于设定阈值，且步骤(b2)中得到的所述最小的P值小于等于设定阈值，则所述模块组为所述待分析的多组学数据之间显著相关的互作模块。
8.根据权利要求1-7中任一所述的方法，其特征在于：所述方法在步骤(B)之后还包括如下步骤(C1)或(C2)：
(C1)对构成所述显著相关的互作模块的各表达模块之间的交集部分的指标数据进行功能富集分析，得到显著富集的功能条目；
(C2)根据所述互作模块的关联信息，构建有权重的互作网络图；
进一步地，步骤(C2)中，是利用visNetwork构建的所述有权重的互作网络图，其中权重信息为各指标数据的平均值。
9.一种能够用于多组学数据联合分析的系统，包括基本组成单元A和基本组成单元B；所述基本组成单元A具有能够实现权利要求1-8任一所述方法中的步骤(A)的功能；所述基本组成单元B具有能够实现权利要求1-8任一所述方法中的步骤(B)的功能。
10.如下任一应用：
(1)权利要求9所述的系统在多组学数据联合分析中的应用；
(2)权利要求1-8中任一所述方法在农业育种、组织发育或药物靶点筛选中的应用。
CN201810111865.3A 2018-02-05 2018-02-05 一种多组学数据联合分析的方法 Active CN110211634B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810111865.3A CN110211634B (zh) 2018-02-05 2018-02-05 一种多组学数据联合分析的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810111865.3A CN110211634B (zh) 2018-02-05 2018-02-05 一种多组学数据联合分析的方法

Publications (2)

Publication Number Publication Date
CN110211634A true CN110211634A (zh) 2019-09-06
CN110211634B CN110211634B (zh) 2022-04-05

Family

ID=67778649

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810111865.3A Active CN110211634B (zh) 2018-02-05 2018-02-05 一种多组学数据联合分析的方法

Country Status (1)

Country Link
CN (1) CN110211634B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111709219A (zh) * 2020-04-28 2020-09-25 上海欧易生物医学科技有限公司 单组学及多组学KEGG PATHWAY map表达热图个性化展示的方法及应用
CN113469816A (zh) * 2021-09-03 2021-10-01 浙江中科华知科技股份有限公司 基于多组学技术的数字货币识别方法、系统和存储介质

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1777686A (zh) * 2003-03-28 2006-05-24 科根泰克股份有限公司 差别表达的基因的调节因子结合位点的统计分析
CN103473416A (zh) * 2013-09-13 2013-12-25 中国人民解放军国防科学技术大学 蛋白质相互作用的模型建立方法和装置
WO2015200823A1 (en) * 2014-06-26 2015-12-30 Institute For Systems Biology Markers and therapeutic indicators for glioblastoma multiforme (gbm)
CN105209918A (zh) * 2013-05-09 2015-12-30 宝洁公司 生物标记鉴定方法和系统
WO2016138488A2 (en) * 2015-02-26 2016-09-01 The Broad Institute Inc. T cell balance gene expression, compositions of matters and methods of use thereof
US20160341731A1 (en) * 2015-05-21 2016-11-24 Ge Healthcare Bio-Sciences Corp. Method and system for classification and quantitative analysis of cell types in microscopy images
CN106709231A (zh) * 2016-10-19 2017-05-24 王�忠 评价生物分子网络中药物对模块间关系的影响的方法
CN107066835A (zh) * 2017-01-19 2017-08-18 东南大学 一种利用公共数据资源发现并整合直肠癌相关基因及其功能分析的方法及系统和应用
US20170277826A1 (en) * 2016-03-27 2017-09-28 Insilico Medicine, Inc. System, method and software for robust transcriptomic data analysis
CN107480467A (zh) * 2016-06-07 2017-12-15 王�忠 一种判别或比较药物作用模块的方法

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1777686A (zh) * 2003-03-28 2006-05-24 科根泰克股份有限公司 差别表达的基因的调节因子结合位点的统计分析
CN105209918A (zh) * 2013-05-09 2015-12-30 宝洁公司 生物标记鉴定方法和系统
CN103473416A (zh) * 2013-09-13 2013-12-25 中国人民解放军国防科学技术大学 蛋白质相互作用的模型建立方法和装置
WO2015200823A1 (en) * 2014-06-26 2015-12-30 Institute For Systems Biology Markers and therapeutic indicators for glioblastoma multiforme (gbm)
WO2016138488A2 (en) * 2015-02-26 2016-09-01 The Broad Institute Inc. T cell balance gene expression, compositions of matters and methods of use thereof
US20160341731A1 (en) * 2015-05-21 2016-11-24 Ge Healthcare Bio-Sciences Corp. Method and system for classification and quantitative analysis of cell types in microscopy images
US20170277826A1 (en) * 2016-03-27 2017-09-28 Insilico Medicine, Inc. System, method and software for robust transcriptomic data analysis
CN107480467A (zh) * 2016-06-07 2017-12-15 王�忠 一种判别或比较药物作用模块的方法
CN106709231A (zh) * 2016-10-19 2017-05-24 王�忠 评价生物分子网络中药物对模块间关系的影响的方法
CN107066835A (zh) * 2017-01-19 2017-08-18 东南大学 一种利用公共数据资源发现并整合直肠癌相关基因及其功能分析的方法及系统和应用

Non-Patent Citations (10)

* Cited by examiner, † Cited by third party
Title
CHAO XU 等: "A Systemic Analysis of Transcriptomic and Epigenomic Data To Reveal Regulation Patterns for Complex Disease", 《GENES|GENOMES|GENETICS》 *
EDWARD LAU 等: "Integrated omics dissection of proteome dynamics during cardiac remodeling", 《NATURE COMMUNICATIONS》 *
KARTHICKEYAN CHELLA KRISHNAN 等: "Integration of Multi-omics Data from Mouse Diversity PanelHighlights Mitochondrial Dysfunction in Non-Alcoholic Fatty Liver Disease", 《CELL SYST》 *
PETER LANGFELDER 等: "WGCNA: an R package for weighted correlation network analysis", 《BMC BIOINFORMATICS》 *
RAQUEL L. COSTA 等: "Distinct co-expression networks using multi-omic data reveal novel interventional targets in HPV-positive and negative head-and-neck squamous cell cancer", 《BIORXIV 》 *
VASILEIOS STATHIAS 等: "Identifying Glioblastoma Gene Networks Based on Hypergeometric Test Analysis", 《PLOS ONE》 *
宋长新 等: "基于 WGCNA 算法的基因共表达网络构建理论及其 R 软件实现", 《基因组学与应用生物学》 *
张淑梅 等: "癌症中 DNA 甲基化基因模块筛选", 《生物信息学》 *
李云飞: "基于加权网络的精神分裂症基因脑区共表达研究", 《万方数据》 *
杨旭 等: "基于新一代高通量技术的人类疾病组学研究策略", 《遗传》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111709219A (zh) * 2020-04-28 2020-09-25 上海欧易生物医学科技有限公司 单组学及多组学KEGG PATHWAY map表达热图个性化展示的方法及应用
CN111709219B (zh) * 2020-04-28 2021-06-18 上海欧易生物医学科技有限公司 单组学及多组学KEGG PATHWAY map表达热图个性化展示的方法及应用
CN113469816A (zh) * 2021-09-03 2021-10-01 浙江中科华知科技股份有限公司 基于多组学技术的数字货币识别方法、系统和存储介质

Also Published As

Publication number Publication date
CN110211634B (zh) 2022-04-05

Similar Documents

Publication Publication Date Title
Lin et al. lncLocator 2.0: a cell-line-specific subcellular localization predictor for long non-coding RNAs with interpretable deep learning
Caudai et al. AI applications in functional genomics
Rockman Reverse engineering the genotype–phenotype map with natural genetic variation
Saint et al. Single-cell imaging and RNA sequencing reveal patterns of gene expression heterogeneity during fission yeast growth and adaptation
Alhamdoosh et al. Easy and efficient ensemble gene set testing with EGSEA
Oellrich et al. An ontology approach to comparative phenomics in plants
CN103778349B (zh) 一种基于功能模块的生物分子网络分析的方法
Sazzini et al. Complex interplay between neutral and adaptive evolution shaped differential genomic background and disease susceptibility along the Italian peninsula
Wei et al. Multiple testing in genome-wide association studies via hidden Markov models
Alexander et al. Quantifying age-dependent extinction from species phylogenies
CN113597645A (zh) 用于重建药物应答和疾病网络的方法和系统以及其用途
CN107368707A (zh) 基于us‑elm的基因芯片表达数据分析系统及方法
Mendoza-Revilla et al. Disentangling signatures of selection before and after European colonization in Latin Americans
CN110211634A (zh) 一种多组学数据联合分析的方法
Ioannidis et al. Knowledge integration in cancer: current landscape and future prospects
Gao et al. Comprehensive network modeling from single cell RNA sequencing of human and mouse reveals well conserved transcription regulation of hematopoiesis
Ma et al. Combined Analysis of BSA-Seq Based Mapping, RNA-Seq, and metabolomic unraveled candidate genes associated with panicle grain number in rice (Oryza sativa L.)
Aigner et al. Functional genomics of osteoarthritis: on the way to evaluate disease hypotheses.
Ai et al. Key genes in the liver fibrosis process are mined based on single-cell transcriptomics
Zhou et al. EVlncRNA-Dpred: Improved prediction of experimentally validated lncRNAs by deep learning
Jia et al. iMPT-FRAKEL: a simple multi-label web-server that only uses fingerprints to identify which metabolic pathway types compounds can participate in
CN113921084B (zh) 疾病相关非编码rna调控轴多维靶向预测方法及系统
Jin et al. CellDrift: inferring perturbation responses in temporally sampled single-cell data
Zamani Dadaneh et al. Covariate-dependent negative binomial factor analysis of RNA sequencing data
Witteveen et al. Publicly available privacy-preserving benchmarks for polygenic prediction

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
REG Reference to a national code

Ref country code: HK

Ref legal event code: DE

Ref document number: 40010423

Country of ref document: HK

GR01 Patent grant
GR01 Patent grant