非靶向代谢组学分析连载（第一篇：缺失数据处理、归一化、标准化）
点击上方「微生信生物」关注我们！专业干货推送！
引子
今天开始我们一起来学习非靶向代谢组学（GC-TOF-MS）的数据处理过程。我是从注释到化学式带有相对丰度的表格开始处理。相信这也是目前大家最为迫切的过程。前面的峰提取，等操作，后面我们会结合MetaboAnalyst给大家补充，但是不会首先做这部分工作。
《非靶向代谢组学》系列推送有事要说
本完整的教程涉及R，MetaboAnalyst（网站工具，也是R包）基本完整的教程也可以仅仅需要R即可；
本推送中所需的代谢组数据、R语言分析等中间流程，均可以在百度网盘下载。链接：https://pan.baidu.com/s/1uwMaOz1vezcpWT8_XSowRg 密码：lh4x
第一篇. 缺失数据处理、归一化、标准化
本节课程，需要完成《非靶向代谢组学》系列之前的阅读
1. 非靶向代谢组学数据分析连载（第零篇引子）
缺失值在我们使用excel打开数据框中的表现就是空的格子，在R语言中缺失值通常以NA表示，判断是否缺失值的函数是is.na；
GC-MS测得数据进过峰对齐，相同峰即为认定相同的物质，但是有时候会有少量样品的峰值对齐的并不好，这个时候我们要注意，如果一个处理的样品大多数的峰都在一个位置上，只有个别峰没有检测到，这个时候没有检测到的峰的样品中很可能也是拥有这个峰的。按照这个思考，我们使用一下方式来对数据进行缺失值的处理：
即在同一个处理大部分样本测得该物质峰，在同一个处理中却又少数样品峰面积数值缺失，此时我们使用大部分样品此物质峰面积均值来替代少数样品的缺失值；
同理，如果一个处理中只有少数样品测得某种物质的峰面积，则表明这种物质很可能不存在，是峰没有对齐产生的错误，此时将该处理中的所有的这种物质峰面积值均去掉；
基于以上思路我们在R语言中简单编程，得到这样的函数：
for(i in 1:nrow(a)){
if(sum(a[i,1:3][1,]==0)==1){
c[i,1:3]= a[i,1:3]
c[i,1:3][,c[i,1:3][1,]==0]=sum(a[i,1:3]/2)
}else if(sum(b=a[i,1:3][1,]==0)==2){
c[i,1:3]= a[i,1:3]
c[i,1:3][,!c[i,1:3][1,]==0]=0
}else{
c[i,1:3]=a[i,1:3]
}
}
实战
我们将文件导入，并且做一个简单的整理
######################################################
将得到的代谢组数据表格原始表格导入
########
#清空内存
rm(list=ls())
#########下面来分类统计分泌物的信息
setwd("E:/微生信生物推送内容编辑/非靶向代谢组数据分析全套")
# 读取OTU表，这里我选择的是整个otu表格，但是一般没有必要全部做差异的啊，相对丰度高的做做就可以了
RE_table = read.delim("RET_a1原始数据.txt",sep="\t",header=T,check.names=F)
RE_table1=RE_table
RE_table$id=NULL
head(RE_table)
a=RE_table
a=as.data.frame(a)
#####将缺失值替换为数值0，方便后续我们进行数值判断
a[is.na(a)] <- 0.0
str(a)
head(a)
使用我们编辑好的脚本进行缺失值处理
c=data.frame(a=c(rep(0,nrow(a))),b=c(rep(0,nrow(a))),c=c(rep(0,nrow(a))),d=c(rep(0,nrow(a))),e=c(rep(0,nrow(a))),f=c(rep(0,nrow(a)))
)
colnames(c)=colnames(a)
row.names(c)=row.names(a)
c
str(c)
head(c)
for(i in 1:nrow(a)){
if(sum(a[i,1:3][1,]==0)==1){
c[i,1:3]= a[i,1:3]
c[i,1:3][,c[i,1:3][1,]==0]=sum(a[i,1:3]/2)
}else if(sum(b=a[i,1:3][1,]==0)==2){
c[i,1:3]= a[i,1:3]
c[i,1:3][,!c[i,1:3][1,]==0]=0
}else{
c[i,1:3]=a[i,1:3]
}
}
for(i in 1:nrow(a)){
if(sum(a[i,4:6][1,]==0)==1){
c[i,4:6]=a[i,4:6]
c[i,4:6][,c[i,4:6][1,]==0]=sum(a[i,4:6]/2)
}else if(sum(b=a[i,4:6][1,]==0)==2){
c[i,4:6]=a[i,4:6]
c[i,4:6][,!c[i,4:6][1,]==0]=0
}else{
c[i,4:6]=a[i,4:6]
}
}
head(c)
c$id=RE_table1$id
dim(c)
head(c,n = 50L)
##最后又将0替换为NA，方便excel使用
#c[c==0] = NA
write.table(c,"RET_a2缺失值处理后分泌物原始数据.txt",row.names= T,
col.names = T,sep = "\t")
笔记：
1， 我们这里是有三个重复，所以代码编写的会简单一些，四个重复，五个对应我们可以修改代码；这里的代码会不能顺利使用，切记！
2. 对缺失值处理之后可能会有很多物质全部为缺失，所以我们需要将这些物质去除掉
保存文件：相信大家都会遇到这样的问题，就是R语言保存的文件首行总是错位一个，这里大家在Excel中打开的时候手动移动即可；（解决办法之后我们会用到）
数据标准化和归一化
缺失值处理结束，我们其实就可以对自己的数据进行多元方差分析了，这个时候我们其实也不怎么必须要将数据标准化，这部分内容可能在之后我们多元统计结果不好的时候会用到，这里我仅仅只介绍几种数据标准化（归一化）的方法
这里我们明确两个概念：
数据标准化：将每个数据减去均值并将结果除以标准差，对应的参数：scale=T
数据中心化：即是将数据减去均值，对应的参数：center=T
scale(data,center=T,scale=F)
会还有一些别的归一化方法，比如使用中位数归一化，使用某一个参考数据进行归一化，关于归一化和标准化的方法如果大家有想法可以留言，我们讨论；
其次还有数据转化
常见的我们可以进行对数和指数变换
R语言函数为log10，sqrt
#进行log转化
c[1:6]<--log10(c[1:6]+0.000001)
##进行平方根转化
c[1:6]<-sqrt(c[1:6])
head(c)
将全部为缺失的物质去掉代码：
##########去掉全部为NA的代谢物
d=c[1:6]
n=ncol(d)
d[is.na(d)] <- 0.0
head(d)
#增加一行，为整列的均值，计算每一列的均值，2就是表示列
d[n+1]=apply(d[c(1:nrow(d)),],1,sum)
#
count=d[d[n+1] > 0,1:n]
head(count)
dim(count)
#整合输出
index = merge(count,c[7],by="row.names",all=F)
head(index)
dim(index)
index$Row.names=NULL
write.table(index,"RET_a21缺失值处理后并去除全部为NA的物质.txt",row.names= T,
col.names = T,sep ="\t")
到此我们的数据预处理便完成了，下面我们就可以开始对数据进行多元统计分析了敬请期待下一篇：多元统计分析