16分析之差异分析(limmar)
limmar package是一个集大成的包,对载入数据、数据前处理(背景矫正、组内归一化和组间归一化都有很多种选择方法)、差异基因分析都有很多的选择。但是这里我当然全部没有研究这这些功能,加油吧!
另外,本人唯一的实验数据结果又达不到预期,明天就是汇报,看来我只能做一个观众了!只能默默祭奠我的劳动吧,劳动光荣!
下面是整个分析过程,一次性贴出来吧,也不进行过多
#limma基于相对丰度计算差异
library(limma)
setwd("E:/Shared_Folder/HG_usearch_HG")
# 读入实验设计
design =read.table("map_lxdjhg.txt", header=T, row.names= 1,sep="\t")
# 读取OTU表
otu_table = read.delim("otu_table.txt",row.names= 1, sep="\t",header=T,check.names=F)
#本步骤采用otu。table基于抽平后的qiime文件,去掉第一行,和第二行首行的#符号,即可导入成功。
# 过滤数据并排序
idx =rownames(design) %in% colnames(otu_table)
idx
sub_design =design[idx,]
count =otu_table[, rownames(sub_design)]
head(count)
# 转换原始数据为百分比,
norm =t(t(count)/colSums(count,na=T))
#下面来筛选差异otu
design.mat =model.matrix(~ 0 + sub_design$SampleType)
colnames(design.mat)=levels(design$SampleType)
#可以同时设置好几组比较
contrast.matrix<- makeContrasts(GC1-G0, levels=design.mat)
#行线性模型拟合
fit <-lmFit(norm, design.mat)
#根据对比模型进行差值计算T-test对数据进行计算
fit2 <-contrasts.fit(fit, contrast.matrix)
#贝叶斯检验
fit2 <-eBayes(fit2)
#"fdr"(equivalent to "BH");lfc(log2-fold-change )
results<-decideTests(fit2,method="global", adjust.method="BH", p.value=0.05, lfc=0)
summary(results)
x<-topTable(fit2,coef=1, number=10000, adjust.method="BH", sort.by="p",resort.by=NULL)
head(x)
#提取各组相对丰度均值
Y=fit$coefficients
head(Y)
#提取我们需要的组
Y=data.frame(Y)
R=Y[c("GC1","G0")]
head(R)
#合并,将相对丰度放在前面
index =merge(R,x, by="row.names",all=F)
head(index)
#存为.xls文件
write.table(index,file="差异比较测试limma.xls", row.names=F, sep="\t")
#存为.txt文件
write.table(index,"差异比较测试limma.txt",quote= FALSE,row.names = F,
col.names = T,sep = "\t")
结果展示:
如果对差异进行可视化请选择:
学习永无止境,分享永不停歇!