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")


学习永无止境,分享永不停歇!

(0)

相关推荐