是否能设置相对丰度、绝对丰度？

[DuXinChina]

在trans_diff函数的lefsef分析中的物种丰富度是相对丰度，似乎不能调节为绝对丰度。应用相对的确是斯特根原文中的用法，但如果能设置一个选择是否将丰度数据转化为相对丰度的参数（relative_abun=T/F）可能会使方法更加灵活。

[ChiLiubio]

您好，lefse虽然默认用的是相对丰度，但是没有强制去检测输入数据是否是相对丰度。所以要实现的话最好的方式是通过改变数据来实现，而且这种操作很简单。下面是示例数据的一个用法。

library(microeco)
data(dataset)
t1 <- trans_diff$new(dataset = dataset, method = "lefse", group = "Group", alpha = 0.05)
# 改变cal_abund函数的rel参数，使用count数据
dataset$cal_abund(rel = FALSE)
t2 <- trans_diff$new(dataset = dataset, method = "lefse", group = "Group", alpha = 0.05)
t1$plot_diff_bar()
t2$plot_diff_bar()
# 结果组别平均丰度也能从res_abund中看出
head(t1$res_abund)
head(t2$res_abund)

[DuXinChina]

谢谢您，感谢您的解惑。还有一个问题，如果将丰度设置为绝对丰度，是否也需要将*10^6参数进行更改？lefse_norm = 1e+06是否可以更改为10^6/（各局域群落中最大总丰度）？这样处理应能使LDA得分的极大值为6。

[ChiLiubio]

理论上是这样的，可以试下。这样确实会更好约束最终获取的LDA的范围。

[DuXinChina]

感谢您，还想请教您一个关于LDA的分的分组的问题。
情况是这样的，我计算了许多个局域群落中的Rao Q指数，并根据两各个群落中Rao Q的均值对群落进行了二分类（Big Rao_Q group，Small Rao_Q group）。同时我计算了每个局域群落中两两物种对间相对丰度的乘积，这类似于微生物共现网络中的边。我在一个函数中套用了您R包中的lefse分析代码，以进行绘图，绘图结果用到了您的t1$plot_diff_abund与t1$plot_diff_bar函数。绘图结果中，在Small Rao_Q group种均值较高的种对的LDA的分被分入了Big Rao_Q group组，我对LDA得分的分组机制并不十分了解，希望您能对一方面给出指教。感谢您。以下是集成了您软件包函数的源代码、输出图片及原始数据。
plot_lefse_RaoQ_sp_pair=function(mc_sp,trait,p_adjust_method,alpha,lefse_norm,add_sig)
{
library(FD)
library(vegan)
library(ggcor)
library(tidygraph)
library(ggraph)
library(netET)
library(microeco)
library(magrittr)

#windowsFonts(Times_New_Roman=windowsFont("Times New Roman"))

Rao_decomposition=function(mc_sp)
{
Rao_decomposition_mat=function(mc_sp)
{
row_divided=function(a){a/rowSums(mc_sp)}
mc_sp=apply(mc_sp,2,row_divided)
MAT=list()

  for (i in 1:nrow(mc_sp))
  {
    mat=as.matrix(mc_sp[i,],,1) %*% t(as.matrix(mc_sp[i,],,1))
    mat[upper.tri(mat)] = NA
    diag(mat)=NA
    #assign(paste0("mc_",i),mat)
    MAT[[i]]=mat
  }  
  MAT
}

a=Rao_decomposition_mat(mc_sp)

n_col=ncol(mc_sp)
n_row=nrow(mc_sp)

vaule=a[[1]][1:n_col^2]
for (i in 2:n_row)
{
  vaule_i=a[[i]][1:n_col^2]
  vaule=rbind(vaule,vaule_i)
}

row_name=c()
for (i in 1:n_row)
{
  row_name[i]=paste0("com_",i)
}

col_name=c()
for (i in 1:n_col)
{
  for(j in 1:n_col)
  {
    a=paste0("sp",i,"_sp",j)
    col_name[(i-1)*n_col+j]=a
  }
}

colnames(vaule)=col_name
rownames(vaule)=row_name
vaule=t(vaule)
vaule=vaule[complete.cases(vaule),]
vaule=t(vaule)
vaule

}

spe_pair=Rao_decomposition(mc_sp)#计算种对间相对丰度的乘积。

RaoQ=function(mc_sp,trait)
{
RaoQ_decomposition_mat_sp.dist=function(mc_sp,trait)
{
row_divided=function(a){a/rowSums(mc_sp)}
trait=gowdis(trait)
mc_sp=apply(mc_sp,2,row_divided)
MAT=list()
for (i in 1:nrow(mc_sp))
{
mat=as.matrix(mc_sp[i,],,1) %% t(as.matrix(mc_sp[i,],,1))
mat[upper.tri(mat)] = NA
diag(mat)=NA
#assign(paste0("mc_",i),mat)
mat=matas.matrix(trait)
MAT[[i]]=mat
}
MAT

}
a=RaoQ_decomposition_mat_sp.dist(mc_sp,trait)
n_col=ncol(mc_sp)
n_row=nrow(mc_sp)
vaule=a[[1]][1:n_col^2]
for (i in 2:n_row)
{
  vaule_i=a[[i]][1:n_col^2]
  vaule=rbind(vaule,vaule_i)
}

row_name=c()
for (i in 1:n_row)
{
  row_name[i]=paste0("com_",i)
}

col_name=c()
for (i in 1:n_col)
{
  for(j in 1:n_col)
  {
    a=paste0("sp",i,"_sp",j)
    col_name[(i-1)*n_col+j]=a
  }
}
colnames(vaule)=col_name
rownames(vaule)=row_name
vaule=t(vaule)
vaule=vaule[complete.cases(vaule),]
vaule=t(vaule)
vaule
vaule=rowSums(vaule)*2

}

Rao_Q=RaoQ(mc_sp,trait)#计算Rao指数

spe_pair_Rao_Q=cbind(spe_pair,Rao_Q)%>%as.data.frame()
spe_pair_Rao_Q$Group=NA
big_row=which(spe_pair_Rao_Q$Rao_Q>=mean(spe_pair_Rao_Q$Rao_Q))
small_row=which(spe_pair_Rao_Q$Rao_Q<mean(spe_pair_Rao_Q$Rao_Q))
spe_pair_Rao_Q[big_row,]$Group="Big Rao_Q group"
spe_pair_Rao_Q[small_row,]$Group="Small Rao_Q group"
#准备数据
sample_table=data.frame(SampleID=row.names(spe_pair_Rao_Q),Group=spe_pair_Rao_Q$Group)
row.names(sample_table)=row.names(spe_pair_Rao_Q)
otu_table=spe_pair_Rao_Q[,!colnames(spe_pair_Rao_Q)%in%c("Rao_Q","Group")]%>%t()%>%as.data.frame()
colnames(otu_table)=row.names(sample_table)
tax_table=data.frame(class=row.names(otu_table))
row.names(tax_table)=row.names(otu_table)
m1=microtable$new(sample_table=sample_table,otu_table=otu_table,tax_table=tax_table)
m1$cal_abund(rel = FALSE)#该步为以绝对值计算，不转换为相对丰度
if (missing(p_adjust_method)){p_adjust_method = "none"}
if (missing(alpha)){alpha = 0.05}
if (missing(lefse_norm)){lefse_norm=1e+06/0.25}#0.25是sp*sp的极大值
if (missing(add_sig)){add_sig=T}
t1 = trans_diff$new(dataset = m1, method = "lefse", group = "Group",p_adjust_method=p_adjust_method ,alpha=alpha,lefse_norm=lefse_norm)
n=length(t1$res_diff$P.unadj)
a=t1$plot_diff_bar(color_values =c("lightblue","pink"), width = 0.4,use_number = 1:n)+
theme_bw()+
ggtitle("LDA Effect Size of Com_dominance degree \nbetween different Rao_Q group")+
theme(legend.position="none")+
theme(axis.title=element_text(size=12))+
xlab("Species pair")
b=t1$plot_diff_abund(color_values =c("lightblue","pink"),barwidth=0.9,add_sig=add_sig,group_order=c("Big Rao_Q group","Small Rao_Q group"),add_sig_label= 'P.adj')+
ggtitle("Bar chart of Com_dominance degree \nbetween different Rao_Q group")+
theme_bw()+
ylab("Common dominance degree")+
theme(legend.position="bottom")+
theme(axis.title=element_text(size=12))+
xlab("Species pair")

require(grid)
grid.newpage() ###新建图表版面
pushViewport(viewport(layout = grid.layout(1,2)))
vplayout <- function(x,y){viewport(layout.pos.row = x, layout.pos.col = y)}
print(a, vp = vplayout(1,1))
print(b, vp = vplayout(1,2))

result=list(a=a,b=b)
result

}
CAN_SP_trait_mc5_seed=read.csv("路径
mc_sp=CAN_SP_trait_mc5_seed[,1:11]
traits=read.csv("C:\大论文\数据\第6章 阔叶红松林下小群落构建机制\traitseed.csv",row.names=1)
trait=scale(traits,center=TRUE,scale=TRUE)
trait=as.data.frame(trait)
trait=trait[1:11,]
plot_RaoQ_sp_pair(mc_sp,trait)
Uploading Rplot04.pdf…
CAN_SP_trait_mc5_seed.csv

[DuXinChina]

traitseed.csv
Rplot04.pdf

[ChiLiubio]

[DuXinChina]

感谢您，祝您生活事业一切顺利！