Example data sets illustration
The CAFT package includes filtered data from two
studies: URT (upper respiratory tract), originally from [1]
and obtained from the R package LOCOM [2]; and
Colon from the R package
curatedMetagenomicData [3] (which can be installed
using Bioconductor). These processed microbiome data
sets can be loaded directly from CAFT package.
The URT data are a subset of data from a study of the microbiome of the upper respiratory tract [1], also used in the LOCOM package [2]. Following [2], we used only left oropharyngeal observations, and followed the same filtering steps as in the LOCOM paper [2]. We excluded three samples from individuals who had used antibiotics within three months prior to sample collection, and subset both the OTU count and taxonomy tables accordingly. The filtered OTU table have 57 rows (samples) and 856 columns (taxa), and the average censoring rate (the mean proportion of zero counts across taxa) is 89.35%.
The Colon data in CAFT are described in
the R package curatedMetagenomicData [3]. For our
analyses here, we retain only adult participants, stool specimens, and
samples with disease status labeled as colorectal cancer (CRC), adenoma,
or healthy.
Binary outcome phenotypes example
We use the URT microbiome data set to illustrate the simplest use of CAFT for a single binary outcome. We first load the data and apply the filters described previously.
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)
#> [1] 57 856
cens.prop = colMeans(throat.otu.table.filter == 0, na.rm = T)
mean(cens.prop)
#> [1] 0.8935276Setting filter.thresh = 0.10 and
fdr.nominal = 0.10, we conduct the taxa differential
abundant analysis by entering \(x_{test}\) and \(x_{adj}\) as separate arguments, using
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)With the default multiple comparison adjustment method of “BH”, we can determine the significantly differential abundant taxa at the nominal FDR level as shown below.
If we want to see the p- or q-values for these taxa, we can use the code (only first 5 shown)
res.CAFT.throat$p.otu[res.CAFT.throat$q.detected.otu][1:5]
#> 4440 2831 2434 4363 1490
#> 1.975870e-03 2.776166e-03 1.604953e-03 2.861630e-03 2.939009e-05Note these are in the order that the taxa are given in the otu table.
Note also that discoveries at different FDR thresholds can be obtained by simply selecting values of res.CAFT.throat$q.otu < FDR threshold.
Setting the grouping label of each sample, we can obtain a
ggplot2-based boxplot showing the relative abundance of a specified OTU
using the function otuboxplot() in our
CAFT package.
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),]
#> [1] "Root;Bacteria;Bacteroidetes;Bacteroidetes;Bacteroidales;Prevotellaceae;Prevotella"
otuboxplot(plot.otu=boxplot.otu, count.data=throat.otu.table.filter,
plot.title = "<i>Prevotellaceae Prevotella</i>",groups=groups)
Multicategory outcome phenotypes example
We use the Colon data to illustrate multivariate covariates such as multigroup categorial data. Here we also use the phyloseq package to extract the data: note that phyloseq is not required to run CAFT.
age_category == "adult"body_site == "stool"disease %in% c("CRC", "healthy", "adenoma")
The filtered OTU table have 856 rows (samples) and 935 columns (taxa). For convenience of display in this vignette, we also set the names of each taxon to its column number.
data(Colon)
count.tab = Colon$otu
sample.tab = Colon$meta
tax.tab = Colon$tax
dim(count.tab)
#> [1] 856 935
colnames(count.tab)=1:ncol(count.tab)We also removed samples with missing age values (no samples were missing gender) and then filtered the dataset to retain only taxa present in more than one sample (i.e., taxa with positive counts in at least two samples), subsetting both the count and taxonomy tables accordingly. After this additional filtering, the OTU table has 849 rows (samples) and 708 columns (taxa), and the average censoring rate (the mean proportion of zero counts across taxa) is 87.20%.
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)
#> [1] 849 708
cens.prop = colMeans(count.tab == 0, na.rm = T)
mean(cens.prop)
#> [1] 0.8720479Then we put the covariates we wish to test into the columns of
x.test, setting “healthy” as the reference level and
including two indicators for “CRC” and “adenoma”. We also put the
covariates we wish to adjust for (the variables “Age” and “Gender”) into
the columns of x.adj.
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)We can now run CAFT exactly as in the univariate example:
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
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0 0 0
#> [2,] 0 1 0 0
res.CAFT.Colon$q.detected.otu
#> [1] "1" "2" "3" "5" "7" "8" "9" "11" "12" "13" "14" "16"
#> [13] "17" "18" "19" "21" "22" "24" "25" "27" "28" "29" "31" "32"
#> [25] "34" "36" "38" "48" "49" "54" "60" "64" "74" "79" "80" "82"
#> [37] "83" "87" "94" "99" "112" "113" "123" "140" "153" "163" "172" "174"
#> [49] "195" "196" "216" "222" "228" "243" "252" "306" "312" "331" "381"We see that 59 taxa were detected at the FDR=0.1 nominal level. If we want to see the p- or q-values for these taxa, we can use the code (only first 10 shown)
res.CAFT.Colon$p.otu[res.CAFT.Colon$q.detected.otu][1:10]
#> 1 2 3 5 7 8
#> 1.098651e-21 5.461108e-08 6.405993e-07 5.423462e-03 2.634203e-03 4.515287e-08
#> 9 11 12 13
#> 5.357031e-06 1.454428e-09 1.506928e-06 1.835382e-02Again, note these are in the order the taxa appear in the otu table, not the size of the p-value.
Separate estimation and testing workflow
The CAFT analysis can also be carried out in two separate steps. The advantage of this approach is computational, and is explained after we describe the two steps.
For the first step, we use caft_estimate to obtain
the \(\widehat{\beta}_j\)s, the
unrestricted estimates of parameters \(\beta_j\). Because no hypothesis is
specified, data should be entered using the full data matrix x. If data
are stored as x_test' andx_adj’, then simply combine them
using `x=cbind(x_test,x_adj)’.
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)
#> CAFT unrestricted estimation
#> Samples: 849
#> Taxa: 708
#> Covariates: 4
#> Taxa used: 217
#> Taxa skipped by rarity filter: 491
#> Taxa failed in unrestricted fit: 0
head(est.CAFT$beta.est)
#> b1.est b2.est b3.est b4.est
#> 1 0.12129121 0.10417073 0.0009675841 0.01772128
#> 2 0.14476170 -0.08198201 -0.0127392812 0.12949872
#> 3 0.08282945 0.05023630 -0.0039049851 0.08690792
#> 4 -0.31959266 0.17216115 0.0032853348 -0.54353098
#> 5 -0.15146946 -0.08585543 -0.0005585079 -0.09654635
#> 6 -0.08010923 0.18089992 0.0128099367 0.22558895For the second step, we specify a hypothesis we want to test, either
by specifying values of x.test and x.adj or by
specifying Gamma and b. Using these inputs,
caft_test will calculate the restricted parameter
estimates \(\widehat{\beta}_j^{(r)}\)
that satisfy the null hypothesis \(\Gamma
\widehat{\beta}_j^{(r)}=b\), and conduct the hypothesis test to
determine whether to accept or reject the hypothesis \(\Gamma \beta=b\). Since the unrestricted
estimates are the same regardless of the null hypothesis we wish to test
(as long as the variables in x are not changed), the
multi-step analysis can be more efficient. In the example below, we test
the hypothesis \(\Gamma \beta = b\).
Because we do not specify \(b\), the
default choice \(b=\text{median}(\Gamma
\widehat{\beta})\) is used, where \(\widehat{\beta}\) are the unrestricted
estimators obtained from caft_estimate.
For example, we can test the significance of the two coefficients corresponding to disease category, controlling for age and gender using
Gamma=cbind(diag(2), matrix(0,nrow=2, ncol=2) )
Gamma
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0 0 0
#> [2,] 0 1 0 0which we can see is the same as
from the previous section. Then, we run caft_test to finish the analysis:
test.CAFT = caft_test(est.CAFT, Gamma=Gamma, fdr.nominal=0.10, adjust.method="BH")
print(test.CAFT)
#> CAFT restricted score test
#> Samples: 849
#> Taxa: 708
#> Covariates: 4
#> Tested contrasts: 2
#> FDR level: 0.1
#> Adjustment method: BH
#> Taxa with p < 0.1 : 89
#> Taxa with q < 0.1 : 59
#> Taxa failed in restricted test: 0
test.CAFT$betahat.median
#> [1] -0.33370892 0.06554493
test.CAFT$b.null
#> [1] -0.33370892 0.06554493
test.CAFT$q.detected.otu
#> [1] "1" "2" "3" "5" "7" "8" "9" "11" "12" "13" "14" "16"
#> [13] "17" "18" "19" "21" "22" "24" "25" "27" "28" "29" "31" "32"
#> [25] "34" "36" "38" "48" "49" "54" "60" "64" "74" "79" "80" "82"
#> [37] "83" "87" "94" "99" "112" "113" "123" "140" "153" "163" "172" "174"
#> [49] "195" "196" "216" "222" "228" "243" "252" "306" "312" "331" "381"
test.CAFT$p.otu[test.CAFT$q.detected.otu][1:10]
#> 1 2 3 5 7 8
#> 1.098651e-21 5.461108e-08 6.405993e-07 5.423462e-03 2.634203e-03 4.515287e-08
#> 9 11 12 13
#> 5.357031e-06 1.454428e-09 1.506928e-06 1.835382e-02These results should be identical to those obtained in the previous section.
The advantage of the separate analysis path is that we can now test for a difference in the effect of the disease categories CRC and adenoma without having to run caft_estimate again (i.e., using the same unrestricted null estimates), by executing the following steps:
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)
#> CAFT restricted score test
#> Samples: 849
#> Taxa: 708
#> Covariates: 4
#> Tested contrasts: 1
#> FDR level: 0.2
#> Adjustment method: BH
#> Taxa with p < 0.2 : 69
#> Taxa with q < 0.2 : 39
#> Taxa failed in restricted test: 0
test.CAFT.2$betahat.median
#> [1] -0.4017934
test.CAFT.2$b.null
#> [1] -0.4017934
test.CAFT.2$q.detected.otu
#> [1] "1" "2" "3" "5" "7" "9" "11" "14" "16" "17" "18" "19"
#> [13] "24" "27" "28" "29" "60" "79" "80" "87" "94" "99" "112" "116"
#> [25] "123" "129" "140" "172" "195" "196" "200" "203" "216" "222" "252" "312"
#> [37] "331" "372" "381"
length(test.CAFT.2$q.detected.otu)
#> [1] 39
test.CAFT.2$p.otu[test.CAFT.2$q.detected.otu][1:10]
#> 1 2 3 5 7 9
#> 7.506836e-08 3.769040e-05 5.299649e-03 8.241583e-03 1.473945e-02 6.265867e-03
#> 11 14 16 17
#> 3.485724e-04 2.784984e-03 2.540297e-03 5.011159e-04We see that there are 39 taxa at which there is evidence that \(\beta_1=\beta_{\text{CRC}}\) is different from \(\beta_2=\beta_{\text{adenoma}}\).
Note that, if the x.test and x.adj input
format is used in caft_test, the order of the variables
must be the same as was used when running
caft_estimate.
The object test.CAFT$betahat.median gives the value of
\(\text{median}(\Gamma
\widehat{\beta})\), default value of \(b\), the median-based compositional
reference value. The user may specify a different value of \(b\) (e.g., for sensitivity analyses or to
choose a different reference, such as the mode of \(\Gamma \widehat{\beta}\) rather than the
median), which is returned in test.CAFT$b.null.
Parallel computation
We can also implement multi-threads parallel computing to speed up
the computation by setting the values of n.cores based on
available computation resources.
Bootstrap implementation of CAFT
CAFT also provides bootstrap version to calculate taxon-level
p-values when either the number of samples or the number of taxa is
small. The bootstrap calculation can be requested by setting
boot.B >= 2. In practice, a larger number of bootstrap
replicates, such as boot.B = 1000, is recommended, as the
smallest non-zero \(p\)-value that can
be obtained is 1/(boot.B+1). Note that the bootstrap version of CAFT can
be very slow.
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.otuBootstrap calibration can also be parallelized over bootstrap replicates. The following code is not evaluated in this vignette because the optimal number of cores depends on the user’s computing environment.