【代谢组学】代谢组学与其他组学数据的整合

​主要内容：

1多组学数据整合概述

2主要分析软件汇总

3数据预处理的方法

1多组学数据整合概述

11 多组学数据整合的现实意义

在研究疾病和其他性状或扰动的分子机制时，在一个以上的组学平台上进行分子分析是一种常见的做法。 

意义：多组分平台数据的整合分析为生物系统提供了更全面的分子特征，有可能改进疾病分子机制的发现以及诊断和预后预测模型。 

目的：整合分析主要是为了更好地确定不同类型的生物分子是如何关联的（例如，代谢数量性状位点（mQTL）分析），为了确定与表型或疾病结局显著相关的分子途径，或者提高预测性能，例如，在生物标记物发现和诊断或预后应用的情况。 

如果数据集不是从同一个人那里收集的，那么跨组学数据集的整合分析的范围也是有限的。

2 分析使用的原材料

21 可以参与分析的“材料”

1软件：所有可以进行统计分析和数据可视化的软件都可以使用（例如，R）。 

2计算资源：原则上，所有分析都可以在功能相当强大的台式计算机上进行，尽管建议在某些分析中使用更强大的计算资源，例如计算服务器。 

3研究设计：主要考虑匹配样本研究设计，在多元分子表型（-omics）平台上分析来自相同个体的生物样本。例如，在匹配样本设计中，单个血液样本可以分成两份，其中一份在代谢平台上分析，另一份在转录平台上分析。然而，匹配样本设计可以基于同一生物样本小份体液，但也可以不一定要基于，例如在同一时间点从同一个人采集的血液样本和另一组织活检样本，即匹配。

4预处理。分析中使用的数据集需要适用于特定分子分析平台的标准方法先行进行预处理，包括质量控制和排除潜在异常值和/或非典型观测。

3 数据分析的方法

31 mQTL分析：代谢特性的遗传驱动因素

描述遗传变异（通常以单核苷酸多态性（SNPs）的形式）与代谢产物丰度之间的统计关联，为理解代谢的遗传驱动因素提供了手段。mQTL分析需要来自相同个体的代谢谱数据和基因分型数据（SNP阵列或DNA测序数据）。基于群体的mQTL研究采用类似于全基因组关联研究的统计方法，不同之处在于表型响应变量是代谢丰度。mQTL分析通常是通过全基因组的单变量关联分析进行的，通常假设一个附加的遗传模型，其遗传效应主要由等位基因的数量编码。对于每个代谢物或代谢特征，进行全基因组扫描，以测试与遗传变异的关联。mQTL分析的结果是关于影响代谢丰度的候选基因驱动因素的信息。
311 mQTL分析过程

1预处理和质控SNP分析数据，去除具有低质量位点和具有次等等位基因频（MAF）的变体；

2预处理和质控代谢组数据（取决于平台和样本类型）； 

3对于每一对代谢物—遗传变异位点： 

    a)拟合一个统计模型来检验代谢物-遗传变异的关联。通常使用线性模型，以代谢特征为响应变量，以遗传变异和相关协变量（如性别、年龄、批次）为预测变量； 

    b)对与遗传效应有关的模型参数进行统计零假设检验（如Wald检验或似然比检验），以确定p值

    c)保存与每个遗传变异-代谢物对儿的遗传效应相关的p值（和相关参数估计）。

4调整保存的p值向量以进行多次测试，例如基于错误发现率（FDR）的方法。 

5使用FDR调整的p值和可接受的FDR水平（通常<005）确定哪些mQTL模型具有统计显著性。

6对于重要的mQTL模型，对模型执行额外的质量控制（高杠杆点，检查残差分布），以确保没有异常值或其他非典型观测影响结果。 

7以表的形式展示结果，表中列出了重要的SNP-代谢物对儿，并展示了排名最高的重要模型。典型图表类型包括曼哈顿图，对于每个重要的代谢物-SNP对儿，绘制代谢物丰度与主要等位基因数（通常编码为0、1或2）的比较。 

8重要mQTL变异的功能解释可以通过确定变异体是否位于基因的编码区来进行，这将表明该基因与相关代谢物之间的功能关系。如果变异体不位于基因的编码区，则可以基于与mQTL变异的基因组距离来确定候选功能基因，在mQTL变异体附近定位的基因将被视为主要功能候选基因。
32 基于代谢途径的整合分析

代谢途径分析提供了一种方法来确定特定的分子路径或生物过程是否与特定的生物扰动相关，如疾病状态。 

代谢途径分析提供了一个将更广泛的生物学功能分配给分子层面的机会，并且可以帮助对研究结果进行生物学解释。 

进行代谢途径分析的两个最常见的框架是基于过度代表（OR）（或）或基于秩的假设检验。 

这里我们采用基于秩的检验方法，其中基因集富集分析（GSEA）是在转录组数据背景下进行代谢路径分析的一个常见例子，它也可以应用于代谢组数据。使用GSEA或OR分析的途径分析也可以通过结合代谢组学和转录组学数据来进行，好处是从这两种数据类型中收集信息以确定代谢途径富集程度。 

代谢途径分析的结果是基于代谢组学和转录组学数据中的证据，给出关于哪些分子途径与所研究的表型相关的信息。

321 代谢途径分析过程

1使用平台特定的方法对转录组和代谢组数据进行预处理和质控。 

2将代谢物标识符分配给轮廓代谢物，然后可以将其映射到生物途径。我们假设转录组数据已经有了带有变量标准标识符的注释信息（例如，Ensembl、Entrez或基因ID），它们也可以映射到路径。 

3对每个带注释的代谢物和每个转录组变量进行单变量关联分析，以确定它们与感兴趣的表型或结果的关联，例如病例-对照状态。在基于秩的（GSEA）分析中，秩是由每个变量的估计效应大小（例如，固定效应模型中的系数）决定。在OR分析中，重要的代谢物和/或基因集合由每个变量的统计零假设检验确定。

4分别使用转录组和代谢组数据进行途径富集分析，例如GSEA，并存储与每个路径相关的p值。 

5结合来自转录组和代谢组数据代谢途径丰度证据，确定两个数据集的组合路径显著性。可以使用基于排列的测试来确定显著性。 

6以表的形式显示结果，表中列出了与p值和FDR调整p值相关的重要途径。
33 结合代谢组学和其他组学数据进行预测建模

预测建模，例如分类或回归，是生物医学研究中的一个共同目标，可以针对疾病诊断、亚型或预后等的预测。 

有时，这种模型的预测性能可以通过包含一种以上的分子表型（-omics）数据来提高。

如果额外的数据（即额外的分子表型）有助于提供补充与预测感兴趣的结果相关的信息，则预期会出现改进的预测性能。如果预测性能没有得到改善，这意味着添加的数据要么根本不是预测性的，要么只是在第一个数据集中捕获的内容上提供冗余信息。 

这项分析的结果是哪个分子表型平台提供了最好的预测信息，如果差异在统计学上是显著的，并且如果两种分子层面数据的联合（组合）建模提供了一个改进的预测模型。

331 预测模型过程

1预处理和质控代谢组学数据和从同一个体收集的额外分子数据，例如转录组数据。

2选择适合高维数据的多元预测模型（如PLS、OPLS-DA或lasso）。在后续步骤中使用此模型。 

3应用（嵌套）交叉验证优化模型参数并评估预测性能，或使用外部测试集评估预测性能。在分类的情况下，受试者工作特征曲线（ROC）和ROC曲线下面积（AUC）通常用于确定分类性能。分别基于各分子表型平台的数据评价模型的预测性能。 

4基于来自两个分子分析平台的数据优化和评估模型的预测性能（参见步骤3），在这两个平台上，数据通过变量的直接连接进行组合。
5比较单个数据集和组合数据集之间ROC曲线和ROC-AUC的差异。如果需要，统计零假设检验可用于确定ROC-AUC估计值是否存在显著差异。 

6可用于未来数据预测的最终预测模型使用优化模型参数（步骤3和4）使用数据集中的所有观测值进行拟合，而不是基于交叉验证训练集拟合的模型。 

7以表格的形式呈现结果，表中分别包含每个数据和组合模型的交叉验证（或外部测试集）的AUC估计值。另外还包括相应ROC曲线的可视化和与AUC值比较相关的假设检验p值。

近年来，由于细胞的异质性及发育分化等相关的问题越来越被研究者们所关注，单细胞转录组分析为研究异质细胞群的复杂生物学过程提供了方法和工具。每一个细胞进行转录组测序时就是细胞发育过程中的快照，单细胞拟时间分析软件Monocle2是基于R语言的安装包，其功能基于单细胞转录组的表达矩阵，通过无监督学习（Reversed Graph Embedding算法）的方式将细胞置于发育轨迹的不同分支上，从而模拟细胞群体生物学过程。也就是我们经常说的拟时序（pseudotime）分析，又称细胞轨迹（cell trajectory）分析。通过拟时分析可以推断出发育过程细胞的分化轨迹或细胞亚型的演化过程，在发育相关研究中使用频率较高。
模拟细胞的分化轨迹的软件，最常用的软件为Monocle2。该算法不仅能模拟细胞的发育轨迹，同时也能对细胞进行聚类（t-SNE）。通过聚类获得不同状态下的差异基因，分析影响分支形成的关键基因及其功能，对研究相关生物学问题有指导性的作用。
Monocle2主要基于关键基因的表达模式，通过学习每个细胞必须经历的基因表达变化的序列，根据拟时间值中对单个细胞进行排序，模拟出时间发育过程的动态变化。而这个排序技术表现是一种在低维空间排布高维数据的降维技术。（具体描叙请参考： >

B站上有大量生信干货, 学生信看B站足以 , 别问我B站是什么哈 ;
生信技能树论坛( >1）ROI

了解过三代测序数据分析的人，对于CCS 环形一致性序列的概念肯定不会陌生，在iso-seq中，提出了比CCS 更加灵活的一个概念：ROI

ROI , 全称 reads of insert，可以理解为插入片段，首先看下三代测序文库构建阶段的reads示意图：

对于上述的文库片段，测序产生的reads 示意图如下：

由于是一个环状分子， 随着测序反应的进行，会循环测序；如果把插入片段的正负链都测了一次，就做1个full pass;

对于CCS 而言，要求至少有2个full pass , 才能去生成CCS reads; 三代测序的特点就是读长很长，可以达到十几kb, 对于短的插入片段而言，CCS这样定义当然没有问题，但是对于全长转录本

而言，转录本长度很长，比如转录本长度1kb, 读长3kb, 此时在一个零模波导孔（ZMW）中测序的reads 就不可能达到2个full pass , 也就产生不了CCS reads, 为了解决这个问题，提高reads的利用率，提出了ROI 的概念，ROI 指的就是插入片段，上图测序reads 产生的ROI 如下：

ROI 不要求满足2个full pass， 相对CCS 而言，更加适合全长转录本的分析；

2）artifacts

这种序列是由于文库制备阶段，adapter 序列错误的将两条转录本的序列链接构成了一个环状分子，这个和adapter 浓度有关，通常这种reads 产生的比例很少，小于05%， 在后续的分析中，这部分reads 需要去除

在PCR 反应中，由于不完全延伸的产物作为了下次扩增反应的引物，导致出现嵌合体序列，直观上看，就是PCR产物来源于两条或者多条reads；

PCR 产生的嵌合体序列，在PCR 反应体系中，这种序列是不可避免的，大约有3%的比例，在后续的分析过程中，可以借助软件去除这部分reads；

FL , Full-length reads, 全长转录本

从raw data 到 ROI , 在从ROI 去除 artifacts reads 之后，我们就得到了用于后续分析的clean reads;

clean reads 就已经是转录本的序列了，我们首先看一下clean reads 当中，哪些是全长转录本；哪些不是全长转录本，这个 *** 作就是分类，classify

全长转录本的示意图如下：

对于全长转录本而言，其ROI reads 中包含5‘ primer 和 3‘ primer; 而且会出现polyA 为结构；（polyA 针对mRNA和部分lncRNA）

对于不同大小的文库，其全长转录本的比例也不同：

可以看到，文库片段越长，全长转录本的比例越低；

4） consensus transcript isoforms

一致性转录本序列，一个ZMW 产生一个转录本的reads, 肯定会有冗余的reads 出现，这是通过聚类（cluster）的方式，就全长转录本序列进行聚类，可以得到一致性的转录本序列；

数据分析流程： >初学RNA-seq，用于有参原核转录组的分析，主要参照DESeq2说明书：(Analyzing RNA-seq data with DESeq2)和(RNA-seq workflow: gene-level exploratory
analysis and differential expression)。reads的count矩阵来源于featureCounts的结果，为原始mapping上的reads数，其格式如下：

接下来构建DESeq2分析所需的分组信息，分组信息包括了实验的分组情况和平行样的情况。比如在我使用的数据中R0_1和R0_2是同一个处理的两个平行样，而R0，R16，R24和R32是不同的处理（就是不同培养时间的样本）。那么分组信息可以按照如下格式构建为dataframe

格式如下：

在构建DESeq数据集时，使用design参数告诉DESeq分组信息：

至此，完成了从featureCounts原始数据到R中DESeq2分析所需数据集的建立。可以使用colData命令查看分组是否正确：

也可以直接运行dds显示数据集的信息：

dim：5846 8 #数据集共5846行，8列
assays(1): counts # 分析数据为readscount
colData names(2): batch condition #分组信息的名字

---