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

## ----include=FALSE------------------------------------------------------------
library(CAFT)

## ----eval=FALSE---------------------------------------------------------------
# pak::pak("CAFT_0.1.0.tar.gz") # local directory

## ----eval=FALSE---------------------------------------------------------------
# remotes::install_github("mli171/CAFT", dependencies = TRUE)

## ----eval=FALSE---------------------------------------------------------------
# browseVignettes("CAFT")

## ----eval=FALSE---------------------------------------------------------------
# vignette("vignette", package = "CAFT")

## ----eval=FALSE---------------------------------------------------------------
# caft(
#   otu.table,
#   x.test = NULL,
#   x.adj = NULL,
#   x = NULL,
#   Gamma = NULL,
#   b = NULL,
#   filter.thresh = 0.05,
#   fdr.nominal = 0.20,
#   adjust.method = "BH",
#   regularize = TRUE,
#   test.method = "rank",
#   n.cores = 1L,
#   boot.B = 0L,
#   boot.parallel = FALSE,
#   boot.n.cores = 1L,
#   boot.seed = NULL,
#   boot.return.dist = FALSE,
#   verbose = FALSE,
#   return.mr.resid = FALSE
# )

## -----------------------------------------------------------------------------
data("URT")

throat.otu.table    <- URT$otu
throat.meta         <- URT$meta
throat.otu.taxonomy <- URT$tax

filter.out.sam = which(throat.meta$AntibioticUsePast3Months_TimeFromAntibioticUsage != "None")
throat.otu.table.filter = throat.otu.table[-filter.out.sam,]
throat.meta.filter = throat.meta[-filter.out.sam,]

dim(throat.otu.table.filter)

cens.prop = colMeans(throat.otu.table.filter == 0, na.rm = T)
mean(cens.prop)

## -----------------------------------------------------------------------------
x.test = ifelse(throat.meta.filter$SmokingStatus == "NonSmoker", 0, 1)
x.adj  = ifelse(throat.meta.filter$Sex == "Male", 0, 1)
res.CAFT.throat = caft(otu.table=throat.otu.table.filter, x.test=x.test, x.adj=x.adj, 
                filter.thresh = 0.10, fdr.nominal = 0.10)

## -----------------------------------------------------------------------------
res.CAFT.throat$q.detected.otu

## -----------------------------------------------------------------------------
res.CAFT.throat$p.otu[res.CAFT.throat$q.detected.otu][1:5]

## ----fig.width=7, fig.height=5, out.width="80%", fig.align='center'-----------
groups = interaction(x.test, x.adj, sep = "_", drop = TRUE)
groups = factor(groups,
                labels = c("Non-Smoker\nMale", "Smoker\nMale",
                           "Non-Smoker\nFemale", "Smoker\nFemale"))

boxplot.otu = "2831"
throat.otu.taxonomy[which(colnames(throat.otu.table.filter) == boxplot.otu),]
otuboxplot(plot.otu=boxplot.otu, count.data=throat.otu.table.filter, 
           plot.title = "<i>Prevotellaceae Prevotella</i>",groups=groups)

## -----------------------------------------------------------------------------
data(Colon)

count.tab  = Colon$otu
sample.tab = Colon$meta
tax.tab    = Colon$tax

dim(count.tab)
colnames(count.tab)=1:ncol(count.tab)

## -----------------------------------------------------------------------------
pNA = which(is.na(sample.tab$age))
if(length(pNA) > 0){
  count.tab = count.tab[-pNA, ]
  sample.tab = sample.tab[-pNA,]
}
# No samples have missing values for gender

## otu presence filtering
p_otu = which(rowSums(t(count.tab) > 0) > 1)
count.tab = count.tab[,p_otu]
tax.tab = tax.tab[p_otu,]

dim(count.tab)

cens.prop = colMeans(count.tab == 0, na.rm = T)
mean(cens.prop)

## -----------------------------------------------------------------------------
Disease1 = Disease2 = rep(0, NROW(sample.tab)) # healthy
Disease1[sample.tab$disease == "CRC"] = 1
Disease2[sample.tab$disease == "adenoma"] = 1

Age = as.numeric(sample.tab$age)
Gender = as.numeric(factor(sample.tab$gender)) - 1

x.test = cbind(CRC = Disease1, adenoma = Disease2)
x.adj  = cbind(Age = Age, Gender = Gender)

## -----------------------------------------------------------------------------
res.CAFT.Colon = caft(otu.table = count.tab,   x.test=x.test, x.adj=x.adj, filter.thresh = 0.10, fdr.nominal = 0.10)

res.CAFT.Colon$Gamma
res.CAFT.Colon$q.detected.otu

## -----------------------------------------------------------------------------
res.CAFT.Colon$p.otu[res.CAFT.Colon$q.detected.otu][1:10]

## -----------------------------------------------------------------------------
est.CAFT = caft_estimate(otu.table=count.tab, x=cbind(x.test, x.adj), filter.thresh=0.10, regularize=TRUE, n.cores=1L) 
print(est.CAFT) 
head(est.CAFT$beta.est) 

## -----------------------------------------------------------------------------
Gamma=cbind(diag(2), matrix(0,nrow=2, ncol=2) )
Gamma

## -----------------------------------------------------------------------------
res.CAFT.Colon$Gamma

## -----------------------------------------------------------------------------
test.CAFT = caft_test(est.CAFT, Gamma=Gamma, fdr.nominal=0.10, adjust.method="BH") 
print(test.CAFT) 
test.CAFT$betahat.median 
test.CAFT$b.null 
test.CAFT$q.detected.otu
test.CAFT$p.otu[test.CAFT$q.detected.otu][1:10]

## -----------------------------------------------------------------------------
Gamma=matrix( c(1,-1,0,0), nrow=1)
test.CAFT.2 = caft_test(est.CAFT, Gamma=Gamma, fdr.nominal=0.20, adjust.method="BH") 
print(test.CAFT.2) 
test.CAFT.2$betahat.median 
test.CAFT.2$b.null 
test.CAFT.2$q.detected.otu
length(test.CAFT.2$q.detected.otu)
test.CAFT.2$p.otu[test.CAFT.2$q.detected.otu][1:10]

## ----eval=FALSE---------------------------------------------------------------
# res.CAFT = caft(otu.table=count.tab, x.test=x.test, x.adj=x.adj, n.cores=2)

## ----eval=FALSE---------------------------------------------------------------
# res.CAFT.boot = caft(otu.table=count.tab, x.test=x.test, x.adj=x.adj, boot.B=1000, boot.seed=1)
# res.CAFT.boot$q.boot.detected.otu
# res.CAFT.boot$q.boot.chi.detected.otu

## ----eval=FALSE---------------------------------------------------------------
# res.CAFT.boot = caft(otu.table=count.tab, x.test=x.test, x.adj=x.adj,
#   boot.B=1000, boot.parallel=TRUE, boot.n.cores=2, boot.seed=1)

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