導(dǎo)讀

NormalizeMets是一個(gè)R語(yǔ)言集成包，主要用于代謝組學(xué)研究中數(shù)據(jù)的歸一化。這個(gè)包可以用于去除數(shù)據(jù)中的噪音，如大樣本中存在的共性問題——質(zhì)譜信號(hào)偏移缰儿。那么除此之外，這個(gè)包還可以進(jìn)行圖形的交互式可視化以及獲得一些常規(guī)的統(tǒng)計(jì)結(jié)果散址，如生物標(biāo)記物的發(fā)現(xiàn)乖阵，聚類及PCA分析宣赔，分類及相關(guān)性分析。

Pipeline

第一步 導(dǎo)入數(shù)據(jù)

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瞪浸，行名是
# 樣品名儒将，列名是一些協(xié)變量信息，如性別对蒲、批次钩蚊、年齡、BMI
sampledata_eg <- alldata_eg$sampledata
dataview(sampledata_eg)

# 3.metabolitedata 包含代謝物的特定信息蹈矮，比如內(nèi)標(biāo)外標(biāo)或者正負(fù)對(duì)照砰逻，行名是代謝物的名稱
# ，其順序要和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)

第二步 數(shù)據(jù)處理

1. log轉(zhuǎn)換

代謝組學(xué)數(shù)據(jù)一般都呈現(xiàn)一個(gè)偏態(tài)分布（右偏）泛鸟，所以需要用一個(gè)合適的轉(zhuǎn)換來使得數(shù)據(jù)的分布變得對(duì)稱一些

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

2.缺失值的處理

代謝組學(xué)數(shù)據(jù)中一個(gè)常見的問題就是存在缺失值蝠咆，那么盡可能多的減少缺失值是數(shù)據(jù)分析前一項(xiàng)非常有必要做的一件事，這里用的填充方法是"k次最近鄰算法"北滥，或者用矩陣中最小值的"一半"作為缺失值的填充值

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

3. 可視化

經(jīng)過log轉(zhuǎn)換的代謝物豐度數(shù)據(jù)可以通過諸多方式進(jìn)行展示刚操，這樣可以直觀的看出數(shù)據(jù)的變異情況聚類情況及離異值等

3.a 那么這里用的是根據(jù)個(gè)體不同批次或者整個(gè)代謝物的分布來看代謝物的一個(gè)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圖，可以用于發(fā)現(xiàn)離異值

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

3.c 熱圖展示（略）

4. 數(shù)據(jù)的歸一化處理

這個(gè)包所采納的數(shù)據(jù)歸一化方法有4種：1)根據(jù)內(nèi)標(biāo);2)根據(jù)QC樣品;3)標(biāo)度化方法;4)聯(lián)合方法

4.a 如何根據(jù)QC樣品來進(jìn)行歸一化再芋，其是根據(jù)QC樣品在進(jìn)樣是有規(guī)律的插入赡茸，然后基于LOESS(locally estimated scatterplot smoothing)信號(hào)校正方法，在statTarget包也有介紹祝闻。

這里用的是另外一個(gè)新的數(shù)據(jù)集，注意這里的參數(shù)lg遗菠，應(yīng)該要在歸一化后做log轉(zhuǎn)換联喘，所以lg參數(shù)應(yīng)設(shè)置為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 評(píng)估以及選擇最佳的歸一化方法

通過鑒定生物標(biāo)記物來評(píng)判歸一化方法采集線性模型的數(shù)據(jù)歸一化方法辙纬，并且能夠鑒定與想要研究的目標(biāo)條件相關(guān)的生物標(biāo)記物

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:"))
# 不同方法之間與未校正的數(shù)據(jù)的比較豁遭，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圖