inst/doc/IFAA.R

In IFAA: Robust Inference for Absolute Abundance in Microbiome Analysis

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(IFAA)

## ----eval=FALSE---------------------------------------------------------------
#  install.packages("IFAA", repos = "http://cran.us.r-project.org")

## ----eval=FALSE---------------------------------------------------------------
#  require(devtools)
#  devtools::install_github("gitlzg/IFAA")

## -----------------------------------------------------------------------------
library(IFAA)

set.seed(1)

## create an ID variable for the example data
ID <- seq_len(20)

## generate three covariates x1, x2, and x3, with x2 binary
x1 <- rnorm(20)
x2 <- rbinom(20, 1, 0.5)
x3 <- rnorm(20)
dataC <- data.frame(cbind(ID, x1, x2, x3))

## Coefficients for x1, x2, and x3 among 10 taxa.
beta_1 <- c(0.1, rep(0, 9))
beta_2 <- c(0, 0.2, rep(0, 8))
beta_3 <- rnorm(10)
beta_mat <- cbind(beta_1, beta_2, beta_3)

## Generate absolute abundance for 10 taxa in ecosystem.
dataM_eco <- floor(exp(10 + as.matrix(dataC[, -1]) %*% t(beta_mat) + rnorm(200, sd = 0.05)))

## Generate sequence depth and generate observed abundance
Ci <- runif(20, 0.01, 0.05)
dataM <- floor(apply(dataM_eco, 2, function(x) x * Ci))
colnames(dataM) <- paste0("rawCount", 1:10)

## Randomly introduce 0 to make 25% sparsity level.
dataM[sample(seq_len(length(dataM)), length(dataM) / 4)] <- 0

dataM <- data.frame(cbind(ID, dataM))

## The following steps are to create a SummarizedExperiment object.
## If you already have a SummarizedExperiment format data, you can
## ignore the following steps and directly feed it to the IFAA function.

## Merge two dataset by ID variable
data_merged <- merge(dataM, dataC, by = "ID", all = FALSE)

## Seperate microbiome data and covariate data, drop ID variable from microbiome data
dataM_sub <- data_merged[, colnames(dataM)[!colnames(dataM) %in% c("ID")]]
dataC_sub <- data_merged[, colnames(dataC)]

## Create SummarizedExperiment object
test_dat <- SummarizedExperiment::SummarizedExperiment(
  assays = list(MicrobData = t(dataM_sub)), colData = dataC_sub
)

## ---- eval=TRUE---------------------------------------------------------------
set.seed(123) # For full reproducibility
results <- IFAA(
  experiment_dat = test_dat,
  testCov = c("x1", "x2"),
  ctrlCov = c("x3"),
  sampleIDname = "ID",
  fdrRate = 0.05,
  nRef = 2,
  paraJobs = 2
)

## ----eval=TRUE----------------------------------------------------------------
summary_res <- results$full_results
sig_results <- subset(summary_res, sig_ind == TRUE)
sig_results

## ---- eval=FALSE--------------------------------------------------------------
#  set.seed(123) # For full reproducibility
#  
#  results <- IFAA(
#    experiment_dat = test_dat,
#    microbVar = c("rawCount1", "rawCount2", "rawCount3"),
#    testCov = c("x1", "x2"),
#    ctrlCov = c("x3"),
#    sampleIDname = "ID",
#    fdrRate = 0.05,
#    nRef = 2,
#    paraJobs = 2
#  )

## ---- eval=TRUE---------------------------------------------------------------
set.seed(123) # For full reproducibility

resultsRatio <- MZILN(
  experiment_dat = test_dat,
  microbVar = c("rawCount5","rawCount8"),
  refTaxa = c("rawCount10"),
  allCov = c("x1", "x2", "x3"),
  sampleIDname = "ID",
  fdrRate = 0.05,
  paraJobs = 2
)

## -----------------------------------------------------------------------------
summary_res_ratio <- resultsRatio$full_results
summary_res_ratio

## ---- eval=T------------------------------------------------------------------
set.seed(123) # For full reproducibility
resultsAllRatio <- MZILN(
  experiment_dat = test_dat,
  refTaxa = c("rawCount10"),
  allCov = c("x1", "x2", "x3"),
  sampleIDname = "ID",
  fdrRate = 0.05,
  paraJobs = 2
)

## -----------------------------------------------------------------------------
summary_res_ratio <- resultsAllRatio$full_results
summary_res_ratio[summary_res_ratio$taxon == "rawCount5", , drop = FALSE]

## -----------------------------------------------------------------------------
sessionInfo()