We introduced adaANCOM for microbiome differential abundance analysis by incorporating phyloengy.

1. Installation and load the package

# devtools::install_github("ZRChao/adaANCOM")  
library(adaANCOM)

2. Usage

2.1 Data transfromation

  • Data generation from zero-inflated Dirichlet multinomial (ZIDM) distribution
N <- 100
D1 <- D2 <- round(runif(N, 1, 100))
Pi <- rzidirichlet(N=N, alpha =  c(3, 6), pi=0.1) #  ZID
which(Pi[,1]==0) # true zero
## [1]  7 27 42 47 54 63 66 68 72
X1 <- matrix(NA, N, 2)
for(i in 1:N) X1[i, ] <- rmultinom(1, D1[i], prob=Pi[i, ])
  • Parameter estimation
est0 <- mom.dm(X1) # moment methods for DM
est0
## [1] 1.140399 2.417518
est1 <- est.dm.NR(X1) # Newton-Raphson for DM
est1
## $alpha
## [1] 1.156717 2.611071
## 
## $loglik
## [1] -362.681
est2 <- est.zidm.EM(X1) # expectation-maximization (EM) for ZIDM
est2
## $alpha
## [1] 2.245509 4.234112
## 
## $pi
## [1] 0.08827392
## 
## $loglik
## [1] -351.9855
  • Model selection, return the p-value for testing \(MN \subset DM \subset ZIDM\), and the value of the log-likelihood
LikRTest(X1)
## $pvalue
## [1]  3.745636e-06 1.183943e-196
## 
## $loglik
## [1] -810.1934 -362.6810 -351.9855
  • Posterior mean transformation, return a matrix with the transformed relative data
prop <- Smooth_X(X1, type='MIX', rel=TRUE)
head(prop)
##           [,1]      [,2]
## [1,] 0.3280042 0.6719958
## [2,] 0.4820121 0.5179879
## [3,] 0.3776819 0.6223181
## [4,] 0.3000282 0.6999718
## [5,] 0.5235432 0.4764568
## [6,] 0.5000691 0.4999309
  • Outlier detection, return the index of the outliers and values of the log-ratio
res <- Outliers(X1, prop)
res
## $outliers
##  [1]  7 23 27 39 42 47 54 63 66 68 72
## 
## $ratios
##  [1] -0.7172258450 -0.0719827551 -0.4993991703 -0.8471634164  0.0942424353
##  [6]  0.0002763695 -0.9791278129 -0.0211519778 -0.5845057771 -0.1072084443
## [11]  0.4390277611 -1.0429226032 -0.4891798226 -0.9257007191 -1.8175248244
## [16] -0.1971520747 -0.1649203860 -0.3092861331 -0.2858987418 -2.0133557507
## [21] -1.3655117223 -1.2237462649 -1.1120660543 -1.2656675322 -2.5836617937
## [26] -1.2740635431 -0.9459402269  0.3387286151 -0.9190780371 -1.3335747102
## [31]  0.8265499988 -2.3847424613 -0.4237113170 -2.1096710885  0.0905030265
## [36]  0.1994493999  0.6751276723 -0.2434023866 -0.0856210638  0.4054651081
## [41] -1.1132537047 -0.1168530627  0.1575275415 -1.3258760251 -0.8472185835
## [46]  0.7647717148 -1.9773663933 -1.9326286807 -0.8067328426 -0.3045742815
## [51] -1.6615287267 -0.4893903835 -0.6714240512 -1.1081839931 -0.0207088779
## [56] -1.2150116138  0.0002519336 -1.2359178712 -0.5845057771  0.8028209818
## [61] -0.3999447058 -1.7713993555 -1.3195171089  0.5415849359 -0.5629767550
## [66] -0.9910045411 -1.1248226876 -1.1895183517 -0.5870921711  0.9561781232
## [71] -1.6823132393 -1.1485890451 -1.4195901595 -1.7388840212 -1.1520351954
## [76] -1.8533360276 -2.2714522682 -0.3646578600 -0.7036336801 -0.7509350499
## [81] -0.9685924573 -1.1308313608 -0.4959681068 -0.0198761336  0.6756435970
## [86]  0.3239736431 -1.0583663424 -1.4726076281 -0.9897511080
## 
## $abnormal
##  [1] -6.917767      -Inf -8.059501       Inf -6.215735 -5.405609 -8.024862
##  [8] -5.554645 -6.917767 -5.480930 -7.152760

2.2 Differential abundance test

  • Next, we illustrate the differential abundance testing based on the example data
library(phyloseq)
library(data.table)
data(example)
data <- example_data$data
fit <- adaANCOM(otu=data.frame(otu_table(data)), 
                Y=sample_data(data)$group, 
                tree=phy_tree(data), tfun = t.test, smooth=0.5)
data.table(fit$L)
##        otu      p.value   deg
##  1:  OTU_1 8.809676e-01 FALSE
##  2:  OTU_2 8.809676e-01 FALSE
##  3:  OTU_3 5.048874e-02 FALSE
##  4:  OTU_4 2.719568e-25 FALSE
##  5:  OTU_5 2.060812e-11 FALSE
##  6:  OTU_6 2.890715e-16 FALSE
##  7:  OTU_7 1.965896e-17 FALSE
##  8:  OTU_8 3.295955e-18 FALSE
##  9:  OTU_9 7.449108e-17 FALSE
## 10: OTU_10 1.833675e-17 FALSE
data.table(fit$V)
##    taxa      p.value        p.adj
## 1:   11 8.073647e-01 8.809676e-01
## 2:   12 5.048874e-02 2.271993e-01
## 3:   13 8.809676e-01 8.809676e-01
## 4:   14 2.911361e-32 2.620225e-31
## 5:   15 2.793830e-01 8.381489e-01
## 6:   16 5.430869e-01 8.767481e-01
## 7:   17 5.844988e-01 8.767481e-01
## 8:   18 3.801746e-01 8.553929e-01
## 9:   19 7.094116e-01 8.809676e-01

Any suggestions or problem, please contact Chao Zhou（) .