Microbial trend analysis

作者: 3between7 | 来源:发表于2020-02-11 17:43 被阅读0次

介绍

2020年1月31日METHODOLOGY上线了一篇文章《Microbial trend analysis for common dynamic
trend, group comparison and classification in
longitudinal microbiome study》，该文章提供了一种MTA（Microbial trend analysis）分析方法，它可以完成3个目的的分析：

可以在菌群水平上捕捉常见的微生物动态趋势，并确定优势taxa;
可以分析这种微生物总体动态趋势在组间是否显著不同；
可以根据个体的纵向微生物图谱（longitudinal microbial profiling）对其进行分类。

安装

该方法对应的R包即MTA（manual），这两天经过我艰苦卓绝又冒着傻气的安装之后，我终于可以好好研究它一番了（安装过程详见安装R包MTA遇到的那些事儿）。

原理

微生物趋势分析

这一趴理解起来有些难，下面是我从文献的背景介绍中截取的几句介绍，大体上也能了解下这MTA到底是一种什么方法：

MTA算是在PTA的基础上扩展而来的一种方法。为了能够从一组样本中提取某种动态模式，MTA整合了基样条法（spline-based method）以及主成分分析法，其中，MTA使用矩阵分解以及lasso回归解决数据的降维问题；使用graph Laplacian penalty兼顾数据间的进化信息关系。以上几点都有助于MTA找到对趋势有贡献的优势taxa。

组间比较

MTA进行组间差异检验的算法有几个关键点，理解了他们基本上就明白是怎么回事了：

首先，无论是对照组还是处理组，其相对丰度数据均为一个 $N*P*T$ 的3维数组， $N$ 为样本数， $P$ 为taxa数， $T$ 为时间点数；

其次，该算法是围绕着差异数组 $\mathbf{G}^{*}$ 展开计算的：即以不放回的方式从对照组(X)与处理组(Z)中随机抽R个样，并构建差异数组 $\mathbf{G}^{*(r)} = \mathbf{Z}^{*(r)} - \mathbf{X}^{*(r)}$ ，其中无论是 $G_n^*$ 、 $Z_n^*$ 还是 $R_n^*$ 的维度都是 $P *T$

再有，零假设下， $G_n^*$ 应该是一个元素全部为0的 $P*T$ 的矩阵，换句话说，在备择假设下， $G_n^*$ 应当是与在零假设下的随机噪音信号不同。也就是说，组间比较的实质其实是比较从 $G^*$ 中获取的共同趋势与噪音趋势之间是否有显著差异（哈哈，感觉没描述清楚）

分类
关于这一部分，文献中有一句话是这么描述的：

we propose a distance-based classification algorithm to classify subjects based on their distances from the predicted microbial matrices of two different groups in the training set.

就我现在的水平，知道MTA是使用一种基于距离的分类算法对样本进行分类的即可，暂时不继续深入探索了。

MTA usage

在R中，可以做上述分析的函数即MTA()，其参数如下：

MTA(Y1,Y_new=NULL,phy.tree=NULL,timevec=NULL, transf="None",
M=NULL,proportion.explained=0.85,
k=5,lambda1.set=c(0.01,0.1,1,10),
lambda2.set=seq(0,2,0.1),lambda3.set=c(0),num.sam=10,alpha=0.05)

参数解释：

Y1：

Y1可以是一个 $N*P*T$ 的丰度数组，此时MTA()会为Y1分析微生物动态趋势并鉴定key taxa；
Y1也可以是一个列表，列表的元素与group对应，每个元素均是一个 $P*T$ 的矩阵，注意group数应该大于等于2；另外，分组名称可以作为列表的名字，默认是 1:length(Y1)。此时，MTA()会进行组间比较或者分类。

Y_new: 可选参数，用于MTA的classification分析中，Y_new是需要被分类的样本丰度数据，维度 $N*P*T$ ，默认为NULL。
phy.tree：P taxa的进化树，phylo-class类型。默认为NULL，即不对数据进行Laplacian penalty。
timevec：可选参数，记录连续时间结点的数字向量，与T等长的。
transf：控制选择哪种方法对原始丰度数据进行转化，“None”代表使用原始数据，“Log”、“ALR”和“CLR”分别代表对原始数据进行log变换、加性对数比变换和中心化对数比变换，默认为“None”。
M：想要得到前多少个共同趋势，默认为NULL。

其余的参数基本上保持默认即可，在这里就不再展开了（ps...实际上我也没办法展开介绍,啊哈哈哈！）

函数返回值

若Y1是数组，那么MTA则在菌群水平分析微生物动态趋势分析，函数返回值则将是一个含有两个元素的列表：

trend plot：从Y1中提取出的共同趋势图；
factor score：对所提取出的趋势有贡献的所有taxa的factor scores（贡献越大分值越高）；

当Y1是一个含有2组或2以上丰度数据的列表，但Y_new是NULL时，MTA则针对微生物总体动态趋势进行组间差异分析，其返回值将是记录有任一两组之间的比较结果的列表。且每一对比较的结果的返回内容包括以下5个元素：

P value：使用Wilcoxon rank-sum test进行组间比较返回的p值；
Trend plot：使用重新随机抽样得到的样本分析得到的趋势图；
Discover rate：随机抽样样本中所有taxa的Discover rate（只有当共同趋势组间差异显著时才会返回该值）；
Factor score：同上，但只有当共同趋势组间差异显著时才会返回；
Standard error：The standard error of the estimated factor scores for all taxa among R resamplings.（汉语解释会有点绕，直接上英文吧。同样的，共同趋势组间差异显著时才会返回）

最后当Y1是list，而Y_new是 $N*P*T$ 的数组时，MTA则会基于训练集Y1，对Y_new中的样本进行分类。返回结果是一个 $N*2$ 的矩阵，其中，第一列是样本ID，第二列是相应的分组结果；

Examples in manual

########### generation data
library(dirmult)
## sample size
N=20;
## time point
T=10;
##number of taxa
P=30;
beta0=abs(rnorm(P,10,20))
## key taxa
causal.ind=sample(1:20,3)
#### generating two groups
group1=array(NA, dim = c(N,P,T))
for(t in 1:T) group1[, ,t] = rdirichlet(N,(beta0))
group2=array(NA, dim = c(N,P,T))
for(t in 1:T){
betaT=rep(0,P)
effect.size1=c(-1,2, -1)*min(beta0[causal.ind])*sin(t)
betaT[causal.ind]=effect.size1
group2[, ,t] = rdirichlet(N,(beta0+betaT))}
########## MTA captures the common trends shared by all subjects in group1
### extracting the common trend from a group of subjects with explained variation >=85%
MTA(Y1=group1)
### extracting the common trend from a group of subjects with 2 common trends
MTA(Y1=group1,M=2)
######## comparing between group1 and group2
MTA(list(group1,group2))
######## classifying subjects based on the training set
MTA(list(group1,group2),group2[1:5,,])

实践一下

正好我手头上有一个项目的数据可以用来练习一下MTA分析，数据分为A，B两组，每组都有3个时间点，我的分析过程如下：

#载入分析基本的R包

library('MTA')
library('phyloseq')

#加载丰度数据
load('ps.species.rda')

#将每组丰度数据组织成N*P*T的数组
otu.tab <- as.data.frame(as.matrix(otu_table(ps.species)))
rownames(otu.tab) <- tax_table(ps.species)[,'Species']
meta <- sample_data(ps.species)

#过滤掉相对丰度全部过低的taxa
filter.index <- apply(otu.tab, 1, function(X){sum(X>0)>0.3*length(X)})
otu.filter <- otu.tab[filter.index, ]

a.meta <- subset(meta,group=="A")
arrayA <- array(NA,dim=c(5,dim(otu.filter)[1],3))
index=1
for (t in c(1,4,6)){
  tmp <- rownames(subset(a.meta,points==t))
  t.tab <- t(otu.filter[,match(tmp,colnames(otu.filter))])
  arrayA[,,index] <- t.tab
  index = index +1
}

b.meta <- subset(meta,group=="B")
arrayB <- array(NA,dim=c(3,dim(otu.filter)[1],3))
index=1
for (t in c(1,4,6)){
  tmp <- rownames(subset(b.meta,points==t))
  t.tab <- t(otu.filter[,match(tmp,colnames(otu.filter))])
  arrayB[,,index] <- t.tab
  index = index +1
}

#组间比较
tab <- list(arrayA,arrayB)
names(tab) <- c('A','B')
compairions <- MTA(tab)  ##这一步的运算还挺慢......

上面分析过程最重要的一点是要对相对丰度过低的菌种进行过滤，否则，在计算随机选择的两组样本的相对丰度差时，很有可能就会得到某列全部为0的数据，这样在进行后面的某步奇异值分解时就会遇见报错！！

唉，原以为不会再出问题了，可：

2020-02-17 17-32-07屏幕截图.png

我看了半天MTA package的原始代码，猜测问题出现在我的样本量过少从而导致交叉验证步骤出了问题:

MTA默认k值为5，而我的样本数一组为5，一组为3，样本量太少，从而导致交叉验证报错，而且这一报错在本项目中无解（实在太少......），我按照原始代码中抽样个数的计算规则进行计算后，发现若按k=5来算，样本数最小应为7才能保证在交叉验证时不会报错！

为了验证我的想法，我翻箱倒柜的找到了另外一个项目的数据，那个项目的样本数目多一些，现在刚把脚本跑上，等等看会不会有啥报错～～

--------------------------------------手动分割线（2.18）-----------------------------------------

可以骂人么......又报错了：

2020-02-18 17-52-50屏幕截图.png

这个报错的的意思是，我的 Allresults长度为1，但后面的向量里却有2个，长度不等所以报错。但是为啥我的 Allresults长度为1呢？待我一探究竟！

--------------------------------------手动分割线（2.19）-----------------------------------------

我找到原因了！我能说......原来是MTA原始脚本有错么，话不多说，咱们来一起看一看MTA_comparison脚本的具体内容：

MTA_comparison=function(Controls, Cases,k,Laplacian.matrix,timevec,lambda1.set,
                        lambda2.set,lambda3.set,num.sam,alpha)
{

"......此处有省略......"

index.sig=which(pvalue.tpc<=alpha)
if (length(index.sig)!=0) {
    Alltrend=NULL
    for(j in index.sig) {

        "......此处有省略......"

        qq=ggplot(plot.data, aes(x=time, y=value, group = variable))+
        geom_line() +geom_point()+theme_bw()+scale_x_continuous(breaks=unique(plot.data$time))+
        facet_wrap(~trend,scales="free_y")

"......此处有省略......"

Allresults=list(pvalue.tpc,qq,AllDiscover,Effect,SE)

names(Allresults)=c("P value","Trend plot","Discover rate", "Factor score", "Standard error")

} else {Allresults=list(pvalue.tpc);names(Allresults)=c("P value","Trend plot")}
  return(Allresults) }

看过原始代码后就能发现，出错的地方属于判断是否有差异显著的trend的部分，他的逻辑大概是这样的：如果有p小于等于0.05的trend，那么输出的结果中将返回"P value","Trend plot","Discover rate", "Factor score", "Standard error"这5个值，若没有则仅返回"P value","Trend plot"。

可是上面的脚本的else后却是这么写的：

......
 else {Allresults=list(pvalue.tpc);
     names(Allresults)=c("P value","Trend plot")}

Allresults=list(pvalue.tpc)仅有1个元素，但下面赋名时却提供了两个，所以肯定会报错，而查看MTA manual：

2020-02-19 14-26-49屏幕截图.png

也是能确定如果没有common trend是差异显著的，MTA确实想要返回P value和Trend Plot两个结果，所以无疑是脚步有问题没错了。

按照这个逻辑，我把绘制Trend Plot的代码行稍加修改移到了if (length(index.sig)!=0) ...... else ......条件判断语句的外面了，这样无论是否显著，均会返回Trend Plot。下图则是我使用新的数据返回的Trend Plot：

2020-02-19 14-29-43屏幕截图.png

只是可惜，这两个commend trend的p value都异常的大：

2020-02-19 14-31-24屏幕截图.png

到此，MTA包的学习算是告一段落啦！开心～

Microbial trend analysis

介绍

安装

原理

MTA usage

参数解释：

函数返回值

Examples in manual

实践一下

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读