16s分析之差异OTU 挑选（edgeR）
16s分析到多样性后我们就开始对差异OTU 进行挑选
基于原始count数目，我们参考两个包来做差异分析：
library(DESeq2)
library(edgeR)
这里当时是首先使用edgeR包做的，下面我们来使用这个包做一次差异分析
现在我们读入进入工作目录并读入文件
#使用edgr统计差异otu
setwd("E:/Shared_Folder/HG_usearch_HG")
# 读入实验设计
design =read.table("map_lxdjhg.txt", header=T, row.names= 1,sep="\t")
# 读取OTU表,全部otu表没有抽平
otu_table =read.delim("otu_table.txt", row.names= 1,sep="\t",header=T,check.names=F)
为什么需要做这一步，edgeR并不是很智能，需要变量一样，并且连排列顺序都一样
# 过滤数据并排序
idx =rownames(design) %in% colnames(otu_table)
#提取分组文件
sub_design =design[idx,]
#提取OTUtable文件
count =otu_table[, rownames(sub_design)]
head(count)
这次我们载入这个包
library(edgeR)
两个重点，第一个是数据整理到DGEList中，这是一个列表，注意列表提取数据方式
# create DGElist
d =DGEList(counts=count, group=sub_design$SampleType)
d =calcNormFactors(d)#默认为TMM标准化
另一个就是完成实验分组矩阵的构建
# 生成实验设计矩阵
#生成矩阵
design.mat =model.matrix(~ 0 + d$samples$group)
#矩阵命名
colnames(design.mat)=levels(design$SampleType)
d2 =estimateGLMCommonDisp(d, design.mat)
d2 =estimateGLMTagwiseDisp(d2, design.mat)
fit = glmFit(d2,design.mat)
# 设置比较组
BvsA <-makeContrasts(contrasts = "GF5-GC5",levels=design.mat)#注意是以GF1为对照做的比较
# 组间比较,统计Foldchange, Pvalue
lrt =glmLRT(fit,contrast=BvsA)
# FDR检验，控制假阳性率小于5%
de_lrt =decideTestsDGE(lrt, adjust.method="fdr", p.value=0.05,lfc=0)#lfc=0这个是默认值
summary(de_lrt)
# 导出计算结果
x=lrt$table
x$sig=de_lrt
#这里加一列调整后的p值
x <- cbind(x,padj = p.adjust(x$PValue, method = "fdr"))
enriched =row.names(subset(x,sig==1))
depleted =row.names(subset(x,sig==-1))
x$level =as.factor(ifelse(x$sig==1, "enriched",ifelse(x$sig==-1,"depleted","nosig")))
做一张完整的差异分析表格一定要有OTU分类信息
setwd("E:/Shared_Folder/HG_usearch_HG/uclust_assigned_taxonomy")
#在加入这个文件taxonomy时，去除后面两列不相干的列
# 读取taxonomy，并添加各列名称
taxonomy =read.delim("rep_seqs_tax.txt", row.names= 1,header=F,sep="\t")
head(taxonomy)
#发现这个注释文件并不适用于直接作图。
#采用excel将其分列处理，并且删去最后一列，才可以运行
colnames(taxonomy)=c("kingdom","phylum","class","order","family","genus","species")
# 标记差异OTU类型
x$otu =rownames(x)
#使用dplyr包，很好用的数据整理包
library(dplyr)
taxonomy$id=rownames(taxonomy)
rownames(taxonomy)= taxonomy$id
idx =rownames(taxonomy) %in% x$otu
tax =taxonomy[idx, ]
# 手动筛选显著的组
x =x[rownames(tax), ]
x$phylum =gsub("p__","",tax$phylum,perl=TRUE)
x$class =gsub("c__","",tax$class,perl=TRUE)
x$order =gsub("o__","",tax$order,perl=TRUE)
x$family =gsub("f__","",tax$family,perl=TRUE)
x$genus =gsub("g__","",tax$genus,perl=TRUE)
x$species =gsub("s__","",tax$species,perl=TRUE)
head(x)
x$otu <- NULL
输出一张完整的OTU差异表格
write.table(x,"otu差异统计表格GC5-GF5.txt",quote= FALSE,row.names = T,
col.names = T,sep = "\t")
这个当然不算完成，加上我们的相对丰度
#加上相对丰度值，同样没有计算均值，通过excel解决
WT =read.table("otu差异统计表格GF5-GC5.txt", header=T, row.names= 1, sep="\t")
head(WT)
setwd("E:/Shared_Folder/HG_usearch_HG")
otu_table =read.delim("otu_table.txt", row.names= 1,sep="\t",header=T,check.names=F)
design =read.table("map_lxdjhg.txt", header=T, row.names= 1,sep="\t")
# 过滤数据并排序
idx =rownames(design) %in% colnames(otu_table)
idx
sub_design =design[idx,]
count =otu_table[, rownames(sub_design)]
# 转换原始数据为百分比
norm =t(t(count)/colSums(count,na=T))# * 100 # normalization to total 100
index =merge(norm,WT, by="row.names",all=F)
head(index)
write.table(index,"otu差异统计表格GF5-GC5.txt",quote= FALSE,row.names = T,
col.names = T,sep = "\t")
去掉不显著的otu
# 去除不显著的OTU
edgr =read.table("otu差异统计表格GF5-GC5.txt", header=T, row.names= 1, sep="\t")
head(edgr)
#匹配第一行中含有nosig的行
edgr1<-edgr[grep(pattern="nosig",edgr[,25]),]
#做一个判别，将没有NA的行输出为TRUE
idx = !edgr[,25]%in% edgr1[,25]#前一个数据集的行名和都一个数据集的列名做匹配，输出逻辑变量
#选择T的行
edgr2=edgr[idx,]
head(edgr2)
write.table(edgr2,"otu差异统计表格GF5-GC5.txt",quote= FALSE,row.names = T,
col.names = T,sep = "\t")
学习永无止境，分享永不停歇！