POSTm_vignette.R
In POSTm: Phylogeny-Guided OTU-Specific Association Test for Microbiome Data

## ----setup, include=FALSE-----------------------------------------------------
library(POSTm)
knitr::opts_chunk$set(echo = TRUE)

## ----eval=FALSE---------------------------------------------------------------
#  post(y, OTU, tree = NULL, X = NULL, cValues = seq(from = 0, to = 0.05, by = 0.01))

## ----eval = FALSE-------------------------------------------------------------
#  p.adjust(p, ..., method = p.adjust.methods, n = length(x = p), alpha = 0.05)

## ----eval = FALSE-------------------------------------------------------------
#  plot(x, ..., siglevel = 0.05, method = "none", alpha = 0.05, subTree = TRUE)

## ----eval = FALSE-------------------------------------------------------------
#  print(x, ..., siglevel = 1.0)

## -----------------------------------------------------------------------------
library("POSTm")
data("POSTmData")

## -----------------------------------------------------------------------------
y <- as.integer(metadata[,"GC"] == "BV")

## -----------------------------------------------------------------------------
result1 <- post(y = y, OTU = otu[,1L:20L], tree = otutree)

## -----------------------------------------------------------------------------
class(x = result1)

## -----------------------------------------------------------------------------
print(x = result1)

## -----------------------------------------------------------------------------
print(x = result1, siglevel = 0.04)

## -----------------------------------------------------------------------------
p.adjust(p = result1)

## -----------------------------------------------------------------------------
p.adjust(p = result1, method = c("holm", "BH"))

## -----------------------------------------------------------------------------
plot(x = result1, siglevel = 0.2)

## -----------------------------------------------------------------------------
plot(x = result1, siglevel = 0.2, method = "BH")

## -----------------------------------------------------------------------------
plot(x = result1, siglevel = 0.2, subTree = FALSE)

## -----------------------------------------------------------------------------
X <- metadata$mRace

## -----------------------------------------------------------------------------
result2 <- post(y = y, OTU = otu[,1L:20L], X = X, tree = otutree)

## -----------------------------------------------------------------------------
Dmat <- ape::cophenetic.phylo(x = otutree)

## -----------------------------------------------------------------------------
# symmetric
all.equal(target = Dmat[lower.tri(x = Dmat)], 
          current = t(x = Dmat)[lower.tri(x = Dmat)])

# zero diagonal
sum(abs(x = diag(x = Dmat)))

# correct names
all(colnames(x = otu) %in% rownames(x = Dmat))
all(colnames(x = otu) %in% colnames(x = Dmat))

## -----------------------------------------------------------------------------
result3 <- post(y = y, OTU = otu[,1L:20L], X = X, tree = Dmat)

## -----------------------------------------------------------------------------
print(x = result3, siglevel = 0.04)

## -----------------------------------------------------------------------------
p.adjust(p = result3, method = "bonferroni")

## -----------------------------------------------------------------------------
try(expr = plot(x = result3))

## -----------------------------------------------------------------------------
result4 <- post(y = y, OTU = otu[,1L:20L])

## -----------------------------------------------------------------------------
print(x = result4)

## -----------------------------------------------------------------------------
try(expr = plot(x = result4))

## -----------------------------------------------------------------------------
result5 <- post(y = y, OTU = otu, X = X, tree = otutree, 
                cValues = c(0.03,0.045,0.06))

## -----------------------------------------------------------------------------
print(x = result5, siglevel = 0.01)

## -----------------------------------------------------------------------------
pv <- p.adjust(p = result5, method = "BH")
head(pv$adjPOST)
head(pv$adjSO)

## -----------------------------------------------------------------------------
plot(x = result5, method = "BY", alpha = 0.01)