R语言：unifrac的计算

Python010

R语言：unifrac的计算,第1张

Unifrac是一个十分常用的衡量不同群落之间谱系结构差异的指标。在R语言中，计算unifrac的函数不只一种，不同函数之间有什么差别呢？本文目的就是对几个常用的计算unifrac的函数的使用方法做个记录。

比较对象

首先，每个包的安装方法：

其次，每个函数的使用方法：

为了更容易区分函数是来自哪个包，每个函数前面都添加了包的名字。

最后，每个函数的运行效率：

虽然phyloseq的安装稍微有些麻烦，并且在计算unfirac之前还需要先转换一下数据类型，但其计算unfirac的效率最高。

对于每个函数的计算时间是否是随OTU数目和样品数目线性增加，还需要进一步探究。

需要注意的是：

picante和PhyloMeasures只能计算unweighted unfirac。

GUniFrac不仅同时计算weighted和unweighted unfirac，同时还能通过参数调节丰度加权的程度。

phyloseq通过控制参数weighted=T或F，可以计算weighted和unweighted unifrac。

虽然一般的16S或者宏基因组等分析流程当中都会包含PCoA分析，但如果自己想要更改分组的形状，或者挑选特定的OTU进行分析，那么自己进行操作会高效很多。

PCoA的作图主要分为三个步骤：

选择特定的相似性距离并计算距离矩阵。距离的选择可以有Bray-curits、Unifrac等，不同的距离有不同的作用和意义（具体可以参考 微生物β多样性常用计算方法比较）。相似性距离可以利用R的GUniFrac和vegan等包计算，也可以利用QIIME计算。

进行PCoA分析，也就是利用表征分析选择最能表示样本距离的坐标轴。这个可以利用R的ape包的pcoa（）命令完成。

PCoA图形展示。图形可以用ordiplot（）命令展示，但如果需要比较美观的图形，建议用ggplot来画。

下面我们以R为基础，展示如何根据Unweighted Unifrac距离来画PCoA图：

----------------------代码开始了-----------------------

###导入需要的R包

library(GUniFrac) #用于计算Unifrac距离

library(ape) # 用于pcoa分析

library(ggplot2) #用于画图

##读文件

Otu_tab <- read.table("otu_table",row.names=1,header=T,sep="t",check.names=FALSE)

Tree <- read.tree("Muscle.align.tree.nhx") #输入OTU表格

Otu_tab <- as.data.frame(t(Otu_tab))

Otu_tab_rff <- Rarefy(Otu_tab)$otu.tab.rff

unifracs <- GUniFrac(Otu_tab_rff,Tree,alpha=c(0, 0.5, 1))

du <- unifracs$unifracs[, , "d_UW"] # 计算Unweighted UniFrac距离

Group <- c('A','B','C') #按照目的输入样本

shape <- c("A" =16,"B" =17,"C" =16) #定义点形状

color <- c("A" ='#CCFF33',"B" ='#CCFF33',"C" ='#CCFF33') #定义点颜色

PCOA <- pcoa(du, correction="none", rn=NULL) #利用PCOA()指令做pcoa分析

result <-PCOA$values[,"Relative_eig"]

pro1 = as.numeric(sprintf("%.3f",result[1]))*100

pro2 = as.numeric(sprintf("%.3f",result[2]))*100

x = PCOA$vectors

sample_names = rownames(x)

pc = as.data.frame(PCOA$vectors)

pc$names = sample_names

legend_title = ""

group = Group

pc$group = group

xlab=paste("PCOA1(",pro1,"%)",sep="")

ylab=paste("PCOA2(",pro2,"%)",sep="")

pca=ggplot(pc,aes(Axis.1,Axis.2)) + #用ggplot作图

geom_point(size=3,aes(color=group,shape=group)) +

# geom_text(aes(label=names),size=4,vjust=-1) +

labs(x=xlab,y=ylab,title="PCOA",color=legend_title,shape=legend_title) +

geom_hline(yintercept=0,linetype=4,color="grey") +

geom_vline(xintercept=0,linetype=4,color="grey") +

scale_shape_manual(values=shape) +

scale_color_manual(values=color) +

theme_bw()

-------------------------代码结束------------------------

结果图展示：

输入OTU表格之后，运行上面代码，就可以出来图形（当然结果的数据输入是经过一定修改，自己根据需求定义样本数目，点的形状和颜色就可以）。