代谢组连载-MetaboDiff包-上
写在前面
代谢组目前在人体上的发展比较快也比较早，很多数据库当然都比较完善，总是在发展，环境代谢组机会还是蛮多的。
这次介绍的MetaboDiff包拥有一套对代谢组懂数据预处理到结果分析和可视化完整解决方案，今天我先带大家学习数据预处理部分，这部分做的好，后面你才能所向披靡，因为后面无论是什么组学，都是类似的多元统计和机器学习部分，前面才是一个组学拥有自己特性的部分。
开始
安装R包，默认数据
MetaboDiff包自带的示例数据来自于这篇文献AKT1 and MYC Induce Distinctive Metaboli c Fingerprints in Human Prostate Cancer。代谢组数据来自于61个前列腺癌病人和25个正常人的前列腺组织。先查看一下这个三个数据。
# library("devtools")
# install_github("andreasmock/MetaboDiff")
library(MetaboDiff)
使用本包自带数据进行代谢组的分析；一般而言代谢组前期工作较少，但是比较麻烦，到我们处理这一个步骤，得到的最重要的一个表格就是···代谢组峰面积矩阵··表格。
这里MetaboDiff自带的数据就是如此，met包提供了一个封装的函数，用于MultiAssayExperiment对象的生成。就这样用也挺方便。只要记住，核心都是SummarizedExperiment，只要我们懂SummarizedExperiment即可自己构造这类S4对象。毕竟在多组学整合情况下有一个统一的框架还是不错的。
# rowData # 分泌物注释表格
# colData # 样本信息注释表格
met <- create_mae(assay,rowData,colData)
如果HMDB、KEGG或ChEBI id是rowData数据集的一部分，则可以从小分子通路数据库(SMPDB)检索进行代谢产物注释。
介绍代谢组数据库 SMPDB - http://smpdb.ca/
不过对于环境来讲，无论是人体数据库还是植物数据库，都无法对环境代谢组进行一个好的通路注释。所以未来就环境微生物代谢组数据库还需要做太多的工作。
SMPDB为“小分子通路数据库（The Small Molecule Pathway Database）”的缩写，它是一个互动性较强且具有数据可视化的数据资源库，含有超过618种发现于人体中的小分子通路。其中70%的通路（超过433条）为该数据独有，在其他数据库中无法找到。这些通路包括代谢、药物以及疾病通路。
SMPDB设计创建的初衷是为了满足并支持代谢组学、转录组学、蛋白质组学和系统生物学中的通路鉴定与通路发现。还有部分数据提供了有关人类代谢通路、代谢疾病通路、代谢物信号通路和药物作用通路的超级精美的超链接图表。
所有的SMPDB通路包括了有关器官、亚细胞、蛋白质复合物的辅助因子、蛋白复合物的位置、代谢物的位置、化学结构和蛋白质复合物的四级结构信息。每一个小分子有相应的链接，点击后可以弹出关于该小分子的详细内容。
所有的SMPDB通路都伴有详细的描述和参考资料，每张图表上都有通路预览图、条件或过程概述。本数据库浏览起来非常方便，支持全文本格式、序列和化学结构搜索。
用户可以使用代谢物名称、药物名称、基因/蛋白质复合物名称、SwissProt ID以及GenBank ID列表进行查询。查询结果可以形成一个与之相匹配的通路列表，并重点给出每一个通路图上的匹配分子。绘图界面还可以得到基因、代谢物和蛋白质复合物数据的可视化图表。SMPDB上所有的图片和表格都可以下载。
注意具体注释的过程中由于会注释到多个通路，所以会有好多通路和描述信息，这里都会添加到rowdata数据后面，后续分析还需要进一步处理。其次这里的注释是调用本地数据库完成的，所以R包的版本会影响注释结果，
met <- get_SMPDBanno(met,column_kegg_id=6,column_hmdb_id=7,column_chebi_id=NA)
#--提取各部分的数据
# 查看有那些子数据集
experiments(met)
#--这里只有一个raw数据集，都是S4对象，这里通过下面方式选择
microbe <- met[["raw"]]
GC_table <- as.data.frame(SummarizedExperiment::assays(microbe))
GC_tax <- as.data.frame(SummarizedExperiment::rowData(microbe))
head(GC_table)
head(GC_tax)
# #--提取样本的分组文件
# colData(met)
就这个热图来看，这部分用于可视化NA值，因为某些样本峰面积缺失，所以NA值的处理是第一个重要问题。这里通过将NA值转换为0,1 矩阵然后通过热图可视化展示，我们可以清楚的看到有哪些代谢物空缺比较多，空缺多的就是要去除的，哪些样本有较多的空缺值，也是要去除的。
通过基础包的热图很容易并且很快的就能统计出来这些信息，右侧的样本或者底部的代谢物都是要去除的候选目标。
# 这部分用于可视化NA值，因为某些样本峰面积缺失，所以NA值的处理是第一个重要问题。这里通过将NA值转换为0,1 矩阵然后通过热图可视化展示
# 使用基础的heatmap进行数据缺失的统计。
na_heatmap(met,group_factor="tumor_normal",label_colors=c("darkseagreen","dodgerblue"))
#--转换代码如下
# met_na = is.na(assay(met)) * 1
空缺代谢组的处置：这个步骤去除了在超过40%的样本中存在空缺值的代谢组，然后对于空缺值在少于40%样本中出现的代谢组使用k-临近法进行空缺值的填充。
注意，这里并没有对原始数据进行修改，只是将修改的数据作为补充的S4对象添加到整体数据中。
# 这个步骤去除了在超过40%的样本中存在空缺值的代谢组，然后对于空缺值在少于40%样本中出现的代谢组使用k-临近法进行空缺值的填充
met = knn_impute(met,cutoff=0.4)
#--如何使用k-临近法填充空缺值，详见下面代码impute.knn
cutoff = 0.4
met_temp = met[["raw"]]
met_temp = met_temp[!rowSums(is.na(assay(met_temp))) > (ncol(assay(met_temp)) *
cutoff), ]
assay(met_temp) = impute::impute.knn(assay(met_temp))$data
异常值热图
# 下面展示样本空缺值情况，和聚类情况，用于后续不好的样本的去除
outlier_heatmap(met,group_factor="tumor_normal",label_colors=c("darkseagreen","dodgerblue"), k=2)
根据上述热图，设置了k=2, 热图形成了cluster1和cluster2，cluster1相对cluster2便是异常值，我们将剔除cluster1。
这是什么原理呢，为什么cluster1就是异常值呢，原始在进行聚类过程中，将代谢物丰度矩阵进行了如下处理：每个代谢物除以中位数，得到样本中代谢组在总体数据中的分布百分数，如果百分比过小或者过大，都会是离散过大的情况， 在后续聚类过程中就会被聚为另外一类。
下面去除cluster为1的样本，从而完成样本的过滤。
# #统计每个样本空缺值的百分比
# col_na = apply(met_na, 2, sum)/nrow(met_na)
#
# #查看每个代谢组部分在数据的位置，百分比
# mat = log2(assays(met)[["imputed"]]/apply(assays(met)[["imputed"]],
# 1, median))
met <- remove_cluster(met,cluster=1)
数据标准化
调用justvsn进行标准化
met <- normalize_met(met)
数据标准化质控
将前面四个部分的数据进行可视化箱线图，看下整体的归一性，发现这一套流程后归一性变得非常好。
quality_plot(met, group_factor="tumor_normal",label_colors=c("darkseagreen","dodgerblue"))