介紹
2020年1月31日METHODOLOGY上線了一篇文章《Microbial trend analysis for common dynamic
trend, group comparison and classification in
longitudinal microbiome study》企垦，該文章提供了一種MTA（Microbial trend analysis）分析方法苇倡，它可以完成3個(gè)目的的分析：
- 可以在菌群水平上捕捉常見的微生物動(dòng)態(tài)趨勢(shì)万细，并確定優(yōu)勢(shì)taxa;
- 可以分析這種微生物總體動(dòng)態(tài)趨勢(shì)在組間是否顯著不同；
- 可以根據(jù)個(gè)體的縱向微生物圖譜（longitudinal microbial profiling）對(duì)其進(jìn)行分類籍滴。
安裝
該方法對(duì)應(yīng)的R包即MTA（manual），這兩天經(jīng)過我艱苦卓絕又冒著傻氣的安裝之后，我終于可以好好研究它一番了（安裝過程詳見安裝R包MTA遇到的那些事兒）涧至。
原理
- 微生物趨勢(shì)分析
這一趴理解起來有些難诅蝶，下面是我從文獻(xiàn)的背景介紹中截取的幾句介紹退个，大體上也能了解下這MTA到底是一種什么方法：
MTA算是在PTA的基礎(chǔ)上擴(kuò)展而來的一種方法。為了能夠從一組樣本中提取某種動(dòng)態(tài)模式调炬，MTA整合了基樣條法（spline-based method）以及主成分分析法语盈，其中，MTA使用矩陣分解以及l(fā)asso回歸解決數(shù)據(jù)的降維問題筐眷；使用graph Laplacian penalty兼顧數(shù)據(jù)間的進(jìn)化信息關(guān)系黎烈。以上幾點(diǎn)都有助于MTA找到對(duì)趨勢(shì)有貢獻(xiàn)的優(yōu)勢(shì)taxa。
- 組間比較
MTA進(jìn)行組間差異檢驗(yàn)的算法有幾個(gè)關(guān)鍵點(diǎn)匀谣，理解了他們基本上就明白是怎么回事了：
首先照棋，無論是對(duì)照組還是處理組，其相對(duì)豐度數(shù)據(jù)均為一個(gè)的3維數(shù)組武翎，
為樣本數(shù)烈炭，
為taxa數(shù)，
為時(shí)間點(diǎn)數(shù)宝恶；
其次符隙，該算法是圍繞著差異數(shù)組展開計(jì)算的：即以不放回的方式從對(duì)照組(X)與處理組(Z)中隨機(jī)抽R個(gè)樣趴捅，并構(gòu)建差異數(shù)組
，其中無論是
霹疫、
還是
的維度都是
再有拱绑，零假設(shè)下，應(yīng)該是一個(gè)元素全部為0的
的矩陣丽蝎，換句話說猎拨，在備擇假設(shè)下，
應(yīng)當(dāng)是與在零假設(shè)下的隨機(jī)噪音信號(hào)不同屠阻。也就是說红省，組間比較的實(shí)質(zhì)其實(shí)是比較從
中獲取的共同趨勢(shì)與噪音趨勢(shì)之間是否有顯著差異（
哈哈，感覺沒描述清楚）
- 分類
關(guān)于這一部分国觉，文獻(xiàn)中有一句話是這么描述的：
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.
就我現(xiàn)在的水平吧恃，知道MTA是使用一種基于距離的分類算法對(duì)樣本進(jìn)行分類的即可，暫時(shí)不繼續(xù)深入探索了麻诀。
MTA usage
在R中痕寓，可以做上述分析的函數(shù)即MTA()，其參數(shù)如下：
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)
參數(shù)解釋：
- Y1：
- Y1可以是一個(gè)
的豐度數(shù)組 针饥，此時(shí)
MTA()會(huì)為Y1分析微生物動(dòng)態(tài)趨勢(shì)并鑒定key taxa厂抽； - Y1也可以是一個(gè)列表，列表的元素與group對(duì)應(yīng)丁眼，每個(gè)元素均是一個(gè)
的矩陣筷凤，注意group數(shù)應(yīng)該大于等于2；另外苞七，分組名稱可以作為列表的名字藐守，默認(rèn)是
1:length(Y1)。此時(shí)蹂风，MTA()會(huì)進(jìn)行組間比較或者分類卢厂。
- Y_new: 可選參數(shù)，用于MTA的classification分析中惠啄，Y_new是需要被分類的樣本豐度數(shù)據(jù)慎恒，維度
- phy.tree：P taxa的進(jìn)化樹融柬，
phylo-class類型。默認(rèn)為NULL趋距，即不對(duì)數(shù)據(jù)進(jìn)行Laplacian penalty粒氧。 - timevec：可選參數(shù)，記錄連續(xù)時(shí)間結(jié)點(diǎn)的數(shù)字向量节腐，與T等長的外盯。
- transf：控制選擇哪種方法對(duì)原始豐度數(shù)據(jù)進(jìn)行轉(zhuǎn)化摘盆，“None”代表使用原始數(shù)據(jù)，“Log”饱苟、“ALR”和“CLR”分別代表對(duì)原始數(shù)據(jù)進(jìn)行l(wèi)og變換孩擂、加性對(duì)數(shù)比變換和中心化對(duì)數(shù)比變換，默認(rèn)為“None”掷空。
- M：想要得到前多少個(gè)共同趨勢(shì)肋殴，默認(rèn)為NULL。
其余的參數(shù)基本上保持默認(rèn)即可坦弟，在這里就不再展開了（ps...實(shí)際上我也沒辦法展開介紹,啊哈哈哈！）
函數(shù)返回值
若Y1是數(shù)組官地，那么MTA則在菌群水平分析微生物動(dòng)態(tài)趨勢(shì)分析酿傍，函數(shù)返回值則將是一個(gè)含有兩個(gè)元素的列表：
- trend plot：從Y1中提取出的共同趨勢(shì)圖；
- factor score：對(duì)所提取出的趨勢(shì)有貢獻(xiàn)的所有taxa的factor scores（貢獻(xiàn)越大分值越高）驱入；
當(dāng)Y1是一個(gè)含有2組或2以上豐度數(shù)據(jù)的列表赤炒，但Y_new是NULL時(shí)，MTA則針對(duì)微生物總體動(dòng)態(tài)趨勢(shì)進(jìn)行組間差異分析亏较，其返回值將是記錄有任一兩組之間的比較結(jié)果的列表莺褒。且每一對(duì)比較的結(jié)果的返回內(nèi)容包括以下5個(gè)元素：
- P value：使用Wilcoxon rank-sum test進(jìn)行組間比較返回的p值；
- Trend plot：使用重新隨機(jī)抽樣得到的樣本分析得到的趨勢(shì)圖雪情；
- Discover rate：隨機(jī)抽樣樣本中所有taxa的Discover rate（只有當(dāng)共同趨勢(shì)組間差異顯著時(shí)才會(huì)返回該值）遵岩；
- Factor score：同上，但只有當(dāng)共同趨勢(shì)組間差異顯著時(shí)才會(huì)返回巡通；
- Standard error：The standard error of the estimated factor scores for all taxa among R resamplings.（漢語解釋會(huì)有點(diǎn)繞尘执，直接上英文吧。同樣的宴凉，共同趨勢(shì)組間差異顯著時(shí)才會(huì)返回）
最后當(dāng)Y1是list誊锭，而Y_new是的數(shù)組時(shí)，MTA則會(huì)基于訓(xùn)練集Y1弥锄，對(duì)Y_new中的樣本進(jìn)行分類丧靡。返回結(jié)果是一個(gè)
的矩陣，其中籽暇，第一列是樣本ID温治，第二列是相應(yīng)的分組結(jié)果；
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,,])
實(shí)踐一下
正好我手頭上有一個(gè)項(xiàng)目的數(shù)據(jù)可以用來練習(xí)一下MTA分析图仓，數(shù)據(jù)分為A罐盔，B兩組，每組都有3個(gè)時(shí)間點(diǎn)救崔，我的分析過程如下：
#載入分析基本的R包
library('MTA')
library('phyloseq')
#加載豐度數(shù)據(jù)
load('ps.species.rda')
#將每組豐度數(shù)據(jù)組織成N*P*T的數(shù)組
otu.tab <- as.data.frame(as.matrix(otu_table(ps.species)))
rownames(otu.tab) <- tax_table(ps.species)[,'Species']
meta <- sample_data(ps.species)
#過濾掉相對(duì)豐度全部過低的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) ##這一步的運(yùn)算還挺慢......
上面分析過程最重要的一點(diǎn)是要對(duì)相對(duì)豐度過低的菌種進(jìn)行過濾惶看，否則捏顺，在計(jì)算隨機(jī)選擇的兩組樣本的相對(duì)豐度差時(shí)，很有可能就會(huì)得到某列全部為0的數(shù)據(jù)纬黎，這樣在進(jìn)行后面的某步奇異值分解時(shí)就會(huì)遇見報(bào)錯(cuò)７尽！
唉本今，原以為不會(huì)再出問題了拆座，可：
我看了半天MTA package的原始代碼，猜測(cè)問題出現(xiàn)在我的樣本量過少從而導(dǎo)致交叉驗(yàn)證步驟出了問題:
MTA默認(rèn)k值為5冠息，而我的樣本數(shù)一組為5挪凑，一組為3，樣本量太少逛艰，從而導(dǎo)致交叉驗(yàn)證報(bào)錯(cuò)躏碳，而且這一報(bào)錯(cuò)在本項(xiàng)目中無解（實(shí)在太少......），我按照原始代碼中抽樣個(gè)數(shù)的計(jì)算規(guī)則進(jìn)行計(jì)算后散怖，發(fā)現(xiàn)若按k=5來算菇绵，樣本數(shù)最小應(yīng)為7才能保證在交叉驗(yàn)證時(shí)不會(huì)報(bào)錯(cuò)！
為了驗(yàn)證我的想法镇眷，我翻箱倒柜的找到了另外一個(gè)項(xiàng)目的數(shù)據(jù)咬最，那個(gè)項(xiàng)目的樣本數(shù)目多一些，現(xiàn)在剛把腳本跑上欠动，等等看會(huì)不會(huì)有啥報(bào)錯(cuò)～～
--------------------------------------手動(dòng)分割線（2.18）-----------------------------------------
可以罵人么......又報(bào)錯(cuò)了：
這個(gè)報(bào)錯(cuò)的的意思是永乌，我的 Allresults長度為1，但后面的向量里卻有2個(gè)翁垂，長度不等所以報(bào)錯(cuò)铆遭。但是為啥我的 Allresults長度為1呢？待我一探究竟沿猜！
--------------------------------------手動(dòng)分割線（2.19）-----------------------------------------
我找到原因了枚荣！我能說......原來是MTA原始腳本有錯(cuò)么，話不多說啼肩，咱們來一起看一看MTA_comparison腳本的具體內(nèi)容：
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) }
看過原始代碼后就能發(fā)現(xiàn)橄妆，出錯(cuò)的地方屬于判斷是否有差異顯著的trend的部分，他的邏輯大概是這樣的：如果有p小于等于0.05的trend祈坠，那么輸出的結(jié)果中將返回"P value","Trend plot","Discover rate", "Factor score", "Standard error"這5個(gè)值害碾，若沒有則僅返回"P value","Trend plot"。
可是上面的腳本的else后卻是這么寫的：
......
else {Allresults=list(pvalue.tpc);
names(Allresults)=c("P value","Trend plot")}
Allresults=list(pvalue.tpc)僅有1個(gè)元素赦拘，但下面賦名時(shí)卻提供了兩個(gè)慌随，所以肯定會(huì)報(bào)錯(cuò)，而查看MTA manual：
也是能確定如果沒有common trend是差異顯著的，MTA確實(shí)想要返回P value和Trend Plot兩個(gè)結(jié)果阁猜，所以無疑是腳步有問題沒錯(cuò)了丸逸。
按照這個(gè)邏輯，我把繪制Trend Plot的代碼行稍加修改移到了if (length(index.sig)!=0) ...... else ......條件判斷語句的外面了剃袍，這樣無論是否顯著黄刚，均會(huì)返回Trend Plot。下圖則是我使用新的數(shù)據(jù)返回的Trend Plot：
只是可惜民效，這兩個(gè)commend trend的p value都異常的大：
到此憔维，MTA包的學(xué)習(xí)算是告一段落啦！開心～
