导读

NormalizeMets是一个R语言集成包，主要用于代谢组学研究中数据的归一化。这个包可以用于去除数据中的噪音，如大样本中存在的共性问题——质谱信号偏移。那么除此之外，这个包还可以进行图形的交互式可视化以及获得一些常规的统计结果，如生物标记物的发现，聚类及PCA分析，分类及相关性分析。

Pipeline

第一步 导入数据

rm(list = ls())
library(NormalizeMets)
data("alldata_eg")
featuredata_eg<-alldata_eg$featuredata

# 1.featuredata is a metabolomics data matrix taking the following format, with metabolites in columns and samples in rows. Unique sample names should be provided as row names.
dataview(featuredata_eg)

# 2.sampledata sampledata is a dataframe that contains sample specific information，行名是
# 样品名，列名是一些协变量信息，如性别、批次、年龄、BMI
sampledata_eg <- alldata_eg$sampledata
dataview(sampledata_eg)

# 3.metabolitedata 包含代谢物的特定信息，比如内标外标或者正负对照，行名是代谢物的名称
# ，其顺序要和featuredata一致
metabolitedata_eg<-alldata_eg$metabolitedata
dataview(metabolitedata_eg)

alldata_eg<-list(featuredata=featuredata_eg, sampledata=sampledata_eg,
                 metabolitedata=metabolitedata_eg)
dataview(alldata_eg$metabolitedata)

第二步 数据处理

1. log转换

代谢组学数据一般都呈现一个偏态分布（右偏），所以需要用一个合适的转换来使得数据的分布变得对称一些

logdata <- LogTransform(featuredata_eg,zerotona=TRUE) # zero=TRUE表示如果存在NA值则用数字0填充

2.缺失值的处理

代谢组学数据中一个常见的问题就是存在缺失值，那么尽可能多的减少缺失值是数据分析前一项非常有必要做的一件事，这里用的填充方法是"k次最近邻算法"，或者用矩阵中最小值的"一半"作为缺失值的填充值

imp <-  MissingValues(logdata$featuredata,sampledata_eg,metabolitedata_eg,
                      feature.cutof=0.8, sample.cutoff=0.8, method="knn")

3. 可视化

经过log转换的代谢物丰度数据可以通过诸多方式进行展示，这样可以直观的看出数据的变异情况聚类情况及离异值等

3.a 那么这里用的是根据个体不同批次或者整个代谢物的分布来看代谢物的一个relative log abundance(RLA)图来展示

RlaPlots(imp$featuredata, sampledata_eg[,1], cex.axis = 0.6,saveinteractiveplot = TRUE)
RlaPlots(t(imp$featuredata), groupdata=rep("group",dim(imp$featuredata)[2]),
         cex.axis = 0.6,saveinteractiveplot = TRUE,xlabel="Metabolites")

Fig.1 Rlaplot

3.b pca图，可以用于发现离异值

PcaPlots(imp$featuredata,sampledata_eg[,1],
         scale=FALSE, center=TRUE, multiplot = TRUE, varplot = TRUE)

3.c 热图展示（略）

4. 数据的归一化处理

这个包所采纳的数据归一化方法有4种：1)根据内标;2)根据QC样品;3)标度化方法;4)联合方法

4.a 如何根据QC样品来进行归一化，其是根据QC样品在进样是有规律的插入，然后基于LOESS(locally estimated scatterplot smoothing)信号校正方法，在statTarget包也有介绍。

这里用的是另外一个新的数据集，注意这里的参数lg，应该要在归一化后做log转换，所以lg参数应设置为lg=FALSE，示例方法

# NormQcsamples<- function(featuredata, sampledata, method = c("rlsc"), span = 0,
#                         deg = 2, lg = TRUE, saveoutput = FALSE,
#                         outputname = "qcsample_results", ...)
data(Didata)
dataview(Didata$sampledata)
Norm_rlsc<- NormQcsamples(sampledata=Didata$sampledata[order(Didata$sampledata$order),],
                    featuredata=Didata$featuredata[order(Didata$sampledata$order),],lg=FALSE)

4.b 评估以及选择最佳的归一化方法

通过鉴定生物标记物来评判归一化方法采集线性模型的数据归一化方法，并且能够鉴定与想要研究的目标条件相关的生物标记物

factormat<-model.matrix(~gender +Age +bmi, sampledata_eg)
ruv2Fit<-LinearModelFit(featuredata=imp$featuredata,
                        factormat=factormat,
                        ruv2=TRUE,k=2,
                        qcmets = which(metabolitedata_eg$IS ==1))
# Exploring metabolites associated with age
unadjustedFit<-LinearModelFit(featuredata=imp$featuredata,
                              factormat=factormat,
                              ruv2=FALSE)
Norm_is <-NormQcmets(imp$featuredata, method = "is", 
                     isvec = imp$featuredata[,which(metabolitedata_eg$IS ==1)[1]])
isFit<-LinearModelFit(featuredata=Norm_is$featuredata,
                      factormat=factormat,
                      ruv2=FALSE)
lcoef_age<-list(unadjusted=unadjustedFit$coefficients[,"Age"],
                is_age=isFit$coefficients[,"Age"],
                ruv2_age=ruv2Fit$coefficients[,"Age"])
lpvals_age<-list(unadjusted=unadjustedFit$p.value[,"Age"],
                 is=isFit$p.value[,"Age"],
                 ruv2=ruv2Fit$p.value[,"Age"])
negcontrols<-metabolitedata_eg$names[which(metabolitedata_eg$IS==1)]                   
CompareVolcanoPlots(lcoef=lcoef_age, 
                    lpvals_age, 
                    normmeth = c(":unadjusted", ":is", ":ruv2"),
                    xlab="Coef",
                    negcontrol=negcontrols)

Fig.2 火山图

# 线性模型拟合的残差RLA图
lresiddata<-list(unadjusted=unadjustedFit$residuals,
                 is=isFit$residuals,
                 ruv2=ruv2Fit$residuals)
CompareRlaPlots(lresiddata,groupdata=sampledata_eg$batch,
                yrange=c(-3,3),
                normmeth = c("unadjusted:","is:","ruv2:"))
# 不同方法之间与未校正的数据的比较，venn图
lnames<- list(names(ruv2Fit$coef[,"Age"])[which(ruv2Fit$p.value[,"Age"]<0.05)],
              names(unadjustedFit$coef[,"Age"])[which(unadjustedFit$p.value[,"Age"]<0.05)],
              names(isFit$coef[,"Age"])[which(isFit$p.value[,"Age"]<0.05)])

VennPlot(lnames, group.labels=c("ruv2","unadjusted","is"))

Fig.3 venn图

4.c 用于分类classification

svm<-SvmFit(featuredata=uv_ruvrandclust$featuredata, 
            groupdata=UVdata$sampledata$group,
            crossvalid=TRUE,
            k=5,
            rocplot = TRUE)

Fig.4 ROC图