library(GUniFrac)
library(MASS)
data(throat.tree)
data(throat.otu.tab)
data(throat.meta)
## Simulate covariate data
set.seed(123)
n = nrow(throat.otu.tab)
Sex <- throat.meta$Sex
Smoker <- throat.meta$SmokingStatus
anti <- throat.meta$AntibioticUsePast3Months_TimeFromAntibioticUsage
Male = (Sex == "Male")**2
Smoker = (Smoker == "Smoker") **2
Anti = (anti != "None")^2
cova = cbind(1, Male, Smoker, Anti)
## Simulate microbiome data
otu.tab.rff <- Rarefy(throat.otu.tab)$otu.tab.rff
unifracs <- GUniFrac(otu.tab.rff, throat.tree, alpha=c(0, 0.5, 1))$unifracs
# Distance matrices
D.weighted = unifracs[,,"d_1"]
D.unweighted = unifracs[,,"d_UW"]
# Kernel matrices
K.weighted = D2K(D.weighted)
K.unweighted = D2K(D.unweighted)
if (requireNamespace("vegan")) {
library(vegan)
D.BC = as.matrix(vegdist(otu.tab.rff, method="bray"))
K.BC = D2K(D.BC)
}
# Simulate phenotype data
rho = 0.2
Va = matrix(rep(rho, (2*n)^2), 2*n, 2*n)+diag(1-rho, 2*n)
Phe = mvrnorm(n, rep(0, 2*n), Va)
K.y = Phe %*% t(Phe) # phenotype kernel
# Simulate genotype data
G = matrix(rbinom(n*10, 2, 0.1), n, 10)
K.g = G %*% t(G) # genotype kernel
## Unadjusted analysis (microbiome and phenotype)
KRV(y = Phe, kernels.otu = K.weighted, kernel.y = "Gaussian") # numeric y
KRV(kernels.otu = K.weighted, kernel.y = K.y) # kernel y
## Adjusted analysis (phenotype only)
KRV(kernels.otu = K.weighted, y = Phe, kernel.y = "linear", X = cova, adjust.type = "phenotype")
if (requireNamespace("vegan")) {
## Adjusted analysis (adjust both kernels; microbiome and phenotype)
KRV(kernels.otu = K.BC, kernel.y = K.y, X = cova, adjust.type='both')
## Adjusted analysis (adjust both kernels; microbiome and genotype)
KRV(kernels.otu = K.BC, kernel.y = K.g, X = cova, adjust.type='both')
}
Run the code above in your browser using DataLab