R语言完整的Lefse分析

写在前面

怎么说呢，我在之前的推送中，在R语言中实现了LEFse分析过程。但是我没有将数据前处理和结果的物种分类图形做出来。这是一定要在R语言中操作的，要不然，仅仅是实现了lefse算法，也是没有什么作用。

所以这篇贴子主要的工作在于：

使用扩增子注物种注释文件构造Lefse分析输入文件
使用R语言实现Lefse核型计算过程
使用R语言实现物种分类树额绘制，比graphlan更加强大。
全部分析提供原始数据，可重现

注意：

物种注释不同数据库格式不同，如果你的注释文件同参考不同，转化为和示例文件相同的格式，方可进行后续操作。

所需R包

library("tidyverse") library(microbiomeViz) library(ggtree) library(tidyverse) library(phyloseq)


##--功能函数#---------
# 功能函数

### 提取OTU表格

vegan_otu <- function(physeq){

OTU <- otu_table(physeq)

if(taxa_are_rows(OTU)){

OTU <- t(OTU)

}

return(as(OTU,"matrix"))

}

### 添加OTU注释信息

vegan_tax <- function(physeq){

tax <- tax_table(physeq)

return(as(tax,"matrix")) }

整理LEFSE分析的数据

这里如果大家熟悉phyloseq对象，直接用这个phyloseq对象就是了

#导入数据 ps = readRDS("../../data/ps_liu.rds")

如果不熟悉phyloseq可以使用普通数据导入

# 导入otu表格


otu = read.delim("../../data/otutab.txt",row.names = 1)

head(otu)

#导入注释文件

tax = read.delim("../..//data/taxonomy.txt",row.names = 1)

head(tax)

#导入分组文件

map = read.delim("../../data/metadata.tsv",row.names = 1)

head(map)
# head(otu)

otu = as.matrix(otu)

str(otu)
tax = as.matrix(tax)

# taxa_names(tax)

ps <- phyloseq(otu_table(otu, taxa_are_rows=TRUE), sample_data(map) , tax_table(tax) # phy_tree(tree) )

没有不要用全部OTU，我按照序列数过滤，当然可以不过滤

# otu数量很多，所以选择一部分展示，一般树展示200个作用最为合适 ps = filter_taxa(ps, function(x) sum(x ) > 400 , TRUE) ps

提取注释文件

这一步是我我们准备Lefse物种注释表格整理的开始。

#-提取otu表和tax表格#-------- otu_table = as.data.frame(t(vegan_otu(ps ))) tax_table = (vegan_tax(ps ))

整理物种注释文件

首先去除未注释出来的物种，这里为：Unassigned。

其次将物种不同级别之间的分隔符号修改为|。

再有将每一级的名称都加上前面级别的全部注释名称，形成结构。

最后按照不同级别数据合并otu表格，然后将注释文件和otu文件合并。

初步整理，去除Unassigned

tax_table[tax_table == "Unassigned"] = NA

tax_table = as.data.frame( tax_table )

物种数据再整理

其次将物种不同级别之间的分隔符号修改为|。

# design = as.data.frame(sample_data(ps)) #合并otu表格和tax表格#--------- otu_tax = merge(otu_table,tax_table,by = "row.names",all = F) # head(otu_tax)


#--------添加门类标记，未注释的结果去掉#----

otu_tax$Kingdom = paste("k",otu_tax$Kingdom,sep = "__")

otu_tax$Kingdom[otu_tax$Kingdom == "k__NA"] =""
otu_tax$Phylum = paste("|p",otu_tax$Phylum,sep = "__")

otu_tax$Phylum[otu_tax$Phylum == "|p__NA"] = ""
otu_tax$Class = paste("|c",otu_tax$Class,sep = "__")

otu_tax$Class[otu_tax$Class == "|c__NA"] = ""
otu_tax$Order = paste("|o",otu_tax$Order,sep = "__")

otu_tax$Order[otu_tax$Order == "|o__NA"] = ""
otu_tax$Family = paste("|f",otu_tax$Family,sep = "__")

otu_tax$Family [otu_tax$Family == "|f__NA"] = ""
otu_tax$Genus = paste("|g",otu_tax$Genus,sep = "__")

otu_tax$Genus[ otu_tax$Genus == "|g__NA"] = ""
otu_tax$Species = paste("|s",otu_tax$Species,sep = "__")

otu_tax$Species[otu_tax$Species == "|s__NA"] = ""

otu_tax$Row.names = paste("|t",otu_tax$Row.names,sep = "__") #合并得到结合全部门类的OTU名称#---- OTU_name = paste(otu_tax$Kingdom,otu_tax$Phylum,otu_tax$Class,otu_tax$Order,otu_tax$Family, otu_tax$Genus,sep = "")

再整理

再有将每一级的名称都加上前面级别的全部注释名称，形成结构。这一步，我们去除了一些特殊的符号。主要是怕影响美观

#----不同等级注释修改#---- otu_tax$Phylum=paste(otu_tax$Kingdom,otu_tax$Phylum,sep = "") otu_tax$Class=paste(otu_tax$Phylum,otu_tax$Class,sep = "") otu_tax$Order=paste(otu_tax$Class,otu_tax$Order,sep = "") otu_tax$Family=paste(otu_tax$Order,otu_tax$Family,sep = "") otu_tax$Genus=paste(otu_tax$Family,otu_tax$Genus,sep = "") otu_tax$Species=paste(otu_tax$Genus,otu_tax$Species,sep = "")

#替换两个括号等特殊符号#-------- library("tidyverse") OTU_name = str_replace(OTU_name,"[(]","") OTU_name = str_replace(OTU_name,"[)]","") # as.character(OTU_name ) # OTU_name = gsub("(","",OTU_name[311]) # row.names(otu_tax) = OTU_name # #将otu表格和tax文件分离#----- # otu_table = otu_tax[2:(ncol(otu_table)+1)] # tax_table = otu_tax[(ncol(otu_table)+2):(ncol(otu_table)+2+6)] #

再整理

最后按照不同级别数据合并otu表格，然后将注释文件和otu文件合并。

#-------将不同分类登记的也添加上去 HA = otu_table # 按Kingdom合并 grp <- otu_tax[rownames(otu_tax), "Kingdom", drop=F] merge=cbind(HA, grp) HA_Kingdom = merge %>% group_by(Kingdom) %>% summarise_all(sum) colnames(HA_Kingdom)[1]="Class"


# 按Phylum合并

grp <- otu_tax[rownames(otu_tax), "Phylum", drop=F]

merge=cbind(HA, grp)

HA_Phylum = merge %>% group_by(Phylum) %>% summarise_all(sum)

colnames(HA_Phylum)[1]="Class"
# 按Class合并

grp <- otu_tax[rownames(otu_tax), "Class", drop=F]

merge=cbind(HA, grp)

HA_Class = merge %>% group_by(Class) %>% summarise_all(sum)

colnames(HA_Class)[1]="Class"
# 按Order合并

grp <- otu_tax[rownames(otu_tax), "Order", drop=F]

merge=cbind(HA, grp)

HA_Order = merge %>% group_by(Order) %>% summarise_all(sum)

colnames(HA_Order)[1]="Class"
# 按Family合并

grp <- otu_tax[rownames(otu_tax), "Family", drop=F]

merge=cbind(HA, grp)

HA_Family = merge %>% group_by(Family) %>% summarise_all(sum)

colnames(HA_Family)[1]="Class"
# 按Genus合并

grp <- otu_tax[rownames(otu_tax), "Genus", drop=F]

merge=cbind(HA, grp)

HA_Genus = merge %>% group_by(Genus) %>% summarise_all(sum)

colnames(HA_Genus)[1]="Class"

# colnames(otu_tax)

# # 按Species合并

# grp <- otu_tax[rownames(otu_tax), "Species", drop=F]

# merge=cbind(HA, grp)

# HA_Species = merge %>% group_by(Species) %>% summarise_all(sum)

# colnames(HA_Species)[1]="Class"

## OTU水平

# merge=cbind(HA, OTU_name)

# head(otu_table)

# HA_OTU = otu_table
# HA_OTU = data.frame(Class = row.names(otu_table),otu_table)

# colnames(HA_OTU )

# head(HA_OTU)

# 合并6个分类级

all = rbind(HA_Kingdom, HA_Phylum, HA_Class, HA_Order, HA_Family, HA_Genus)

dim(all)

# head(OTU_name)

因为我们去除了未注释出来的，所以会有重复

这里只需要简单去除就可以了。我们还构建了phyloseq对象，也就是说基于这个群落我们做其他分析也是是十分方便。

# 去除重复 all = distinct(all, Class, .keep_all = TRUE) dim(all)


#-------------lefse分析构建

all1 = all

row.names(all1) = all1$Class

all1$Class = NULL

all1 = as.matrix(all1)

# head(all1)
#-构建phylose对象

ps_G_graphlan = phyloseq(otu_table(as.matrix(all1),taxa_are_rows = TRUE),

sample_data(ps))

ps_G_graphlan
#----提取OTU表格

otu_table = as.data.frame((vegan_otu(ps_G_graphlan)))

otu_table[otu_table==0] <- 1
# row.names(otu_table)
# head(design)

design = as.data.frame(sample_data(ps_G_graphlan))

如果我们要使用功能python那个脚本运行，保存就可以使用以下代码

这里我就不做运行了，之前的帖子相信大家都看过了，我是做了对比的，发现算法都是一致。

# filename = paste(lefsepath,"/LEFSE_to_run_G_level.txt",sep = "") # write.table(otu_table, filename,append = F, quote = F,col.names= F,sep = "\t") #文件预处理 # format_input.py LEFSE_to_run_G_level.txt pri_lefse.in -c 1 -u 2 -o 1000000 # ~/src/nsegata-lefse/run_lefse.py pri_lefse.in pri_lefse_2.res -l 4 # plot_res.py pri_lefse_2.res lefse_barplot.pdf --format pdf # plot_cladogram.py pri_lefse_2.res lefse_tree.pdf --format pdf

# # 注意这里 –c用来指定分组信息-u 1指定样品信息 # #文件分析,这里-l设置LDA阈值，默认为2，我们使用4 会更加严格 # ~/src/nsegata-lefse/run_lefse.py pri_lefse.in pri_lefse_2.res -l 4 # #柱状图绘制 # plot_res.py pri_lefse_2.res lefse_barplot.pdf --format pdf # #树状图绘制 # plot_cladogram.py pri_lefse_2.res lefse_tree.pdf --format pdf # #做每个差异的柱状图 # mkdir biomarkers_raw_images ## plot_features.py pri_lefse.in pri_lefse_2.res biomarkers_raw_images/

R语言lefse分析

这部分如果大家不懂得原理建议查看之前的推送，基于R语言实现Lefse分析。

#-------------R语言进行lefse分析 # ps_sub <- subset_samples(ps_G_graphlan,Group %in% c("July-Orchard","July-cropland"));ps_sub ps_sub = ps_G_graphlan


otu_table = as.data.frame((vegan_otu(ps_sub)))

otu_table[otu_table==0] <- 1
# head(otu_table)
map = as.data.frame(sample_data(ps_sub))

otu = (otu_table)

claslbl= map$Group

set.seed(56290);

#KW rank sum test

dat3t = otu_table

# head(otu)

# dat3t = as.matrix(dat3t)

# str(dat3t)

rawpvalues <- apply(dat3t, 2, function(x) kruskal.test(x, claslbl)$p.value);

#--得到计算后得到的p值

ord.inx <- order(rawpvalues);

rawpvalues <- rawpvalues[ord.inx];

clapvalues <- p.adjust(rawpvalues, method ="fdr");

# p.adjust

dat3t <- dat3t[,ord.inx];

dim(dat3t)
wil_datadf <- as.data.frame(dat3t);

# head(wil_datadf)

#if no subclass within classes then no wilcoxon rank sum test

#Linear Discriminant analysis (LDA)

library( MASS)

ldares <- lda(claslbl~ .,data = wil_datadf);

# ldares

ldamean <- as.data.frame(t(ldares$means));

# ldamean

class_no <<- length(unique(claslbl));

ldamean$max <- apply(ldamean[,1:class_no],1,max);

ldamean$min <- apply(ldamean[,1:class_no],1,min);
#---计算LDA

ldamean$LDAscore <- signif(log10(1+abs(ldamean$max-ldamean$min)/2),digits=3);
# head(ldamean)
a = rep("A",length(ldamean$max))

i = 1

for (i in 1:length(ldamean$max)) {

name =colnames(ldamean[,1:class_no])

a[i] = name[ldamean[,1:class_no][i,] %in% ldamean$max[i]]

}

ldamean$class = a

# head(ldamean)
ldamean$Pvalues <- signif(rawpvalues,digits=5);

ldamean$FDR <- signif(clapvalues,digits=5);

resTable <- ldamean;

# head(ldamean)

# it seems lda add ` around names containing dash "-", need to strip this off

rawNms <- rownames(resTable);

rownames(resTable) <- gsub("`", '', rawNms);

# head(resTable)

# if(pvalOpt == "raw"){

# de.Num <- sum(rawpvalues<=p.lvl & ldamean$LDAscore>=lda.lvl)

# }else{

p.lvl =0.05

lda.lvl = 2

de.Num <- sum(clapvalues<=p.lvl & ldamean$LDAscore>=lda.lvl)

# }
if(de.Num == 0){

current.msg <<- "No significant features were identified with given criteria.";

}else{

current.msg <<- paste("A total of", de.Num, "significant features with given criteria.")

}

current.msg

# sort by p value

ord.inx <- order(resTable$Pvalues, resTable$LDAscore);

resTable <- resTable[ord.inx, ,drop=FALSE];

#p-values column to appear first; then FDR and then others

resTable <- resTable[,c(ncol(resTable),1:(ncol(resTable)-1))];

resTable <- resTable[,c(ncol(resTable),1:(ncol(resTable)-1))];

write.csv(resTable,"cs.csv")

R语言绘制物种分类树

这部分我首先进行了颜色的映射数据整理，然后使用microbiomeviz做了圈图的输出。

taxtree = resTable[clapvalues<=p.lvl & ldamean$LDAscore>=lda.lvl,]


delimiter = "\\|"

tax_split <- strsplit(row.names(taxtree), delimiter)

row.names(taxtree)<- vapply(tax_split, tail, n = 1, "")
# head(taxtree)
#-提取所需要的颜色

colour = c('darkgreen','red',"blue")
selececol = colour[1:length(levels(as.factor(taxtree$class)))]
names(selececol) = levels(as.factor(taxtree$class))

A = rep("a",length(row.names(taxtree)))
i = 1

for (i in 1:length(row.names(taxtree))) {

A[i] = selececol [taxtree$class[i]]

}
# A
lefse_lists = data.frame(node=row.names(taxtree),

color=A,

stringsAsFactors = FALSE

)
# str(all)
## 计算均值用于呈现结点大小

dat <- data.frame(V1=all[,1], V2=rowMeans(all[,-1]), stringsAsFactors = FALSE)
# head(dat)

# dim(dat)

# write.csv(dat,"./tree_tax.csv")

# dat$V2 = NULL

# 用物种和丰度生成树骨架

tr <- parseMetaphlanTSV(dat, node.size.offset=2, node.size.scale=0.8)
# tr
# 构造树#------

p =tree.backbone(tr, size=1,layout= 'circular') p # ?tree.backbone

# 注释树 p <- clade.anno(p, lefse_lists, alpha=0.3) p

ggsave("./cs.pdf",p,width = 10,height = 10)

R语言实现LEfse分析从数据整理到树形图绘制