We introduce a Statistical Approach via Pseudo-value Information and Estimation (SOHPIE; pronounced as “Sofie”). This is a regression modeling method for differential network (DN) analysis that can include covariate information in analyzing microbiome data. SOHPIE is a software extension of our methodology work [1].
Please install these R packages prior to use SOHPIE.
# library(robustbase) # To fit a robust regression.
# library(parallel) # To use mclapply() when reestimating the association matrix.
# library(dplyr) # For the convenience of tabulating p-values, coefficients, and q-values.
# library(fdrtool) # For false discovery rate control.
# library(gtools) # To estimate an association matrix via SparCC.Two sample datasets are available in this package. One
(combinedamgut) is from the American Gut Project [2] and
contains 138 taxa and 268 subjects. In this vignette, the first 30 out
of 138 taxa will be used for the simple demonstration purpose. The other
(combineddietswap) is from the geographical epidemiology
study of diet swap intervention [3] that includes 112 taxa with 37
subjects (20 African Americans from Pittsburgh and 17 rural South
Africans). The full data of each study are available in the SpiecEasi
and microbiome R packages, respectively.
The main grouping variable will be the indicator variable for the status of living with a dog. After the data processing, the indices of subjects will be available for each ‘Not living with a dog (Group A)’ vs. ‘Living with a dog (Group B).’ We need these indices for the estimation of group-specific p × p association matrices (and re-estimation of association matrices for pseudo-value calculations later).
# Note: Again, we will use a toy example with the first 30 out of 138 taxa.
OTUtab = combinedamgut[ , 8:37]
# Clinical/demographic covariates (phenotypic data):
# Note: All of these covariates in phenodat below will be included in the regression
# when you use SOHPIE_DNA function later. Please make sure
# phenodat below include variables that will be analyzed only.
phenodat = combinedamgut[, 1:7] # first column is ID, so not using it.# Obtain indices of each grouping factor.
# In this example, a variable indicating the status of living with a dog was chosen (i.e. bin_dog).
# Accordingly, Groups A and B imply living without and with a dog, respectively.
newindex_grpA = which(combinedamgut$bin_dog == 0)
newindex_grpB = which(combinedamgut$bin_dog == 1)Upon our data processing step above is complete, we can then fit a
pseudo-value regression using SOHPIE_DNA function. An
important note! Please provide the object name of each OTU table and
clinical/demographic data (i.e. metadata) separately in the function. In
addition, you must indicate the object names of the indices for each
group of a binary indicator variable that is used as a main predictor
variable (e.g. living with a dog vs. without a dog). Lastly, you must
enter a trimming proportion c, which ranges from 0.5 to
1.
Now, I would like to show you that SOHPIE has some convenient
tools/functions after fitting a pseudo-value regression. There are
functions that you can quickly extract names of taxa that are
significantly differentially connected (DC; DCtaxa_tab), as
well as p-values (pval and pval_specific_var),
adjusted p-values (q-values; qval and
qval_specific_var), coefficient estimates
(coeff and coeff_specific_var), and standard
errors (stderrs and stderrs_specific_var) of
all variables that are considered in the regression or a specific
variable.
# qval() function will get you a table with q-values.
qval(SOHPIEres)
#> bin_dog age sex bin_floss bin_exercise cat_alcohol1
#> 326792 0.73680874 2.180102e-02 0.6887811 0.0006955236 0.6559827 0.9506848
#> 348374 0.52163971 1.051708e-05 0.7805108 0.0626509883 0.4383095 0.1077605
#> 181016 0.57984394 9.028940e-02 0.3425181 0.1100483595 0.1740821 0.8426311
#> 191687 0.09757318 8.106589e-02 0.4677735 0.0719793426 0.1740821 0.9181658
#> 305760 0.18227112 2.936511e-02 0.4631969 0.0924203987 0.5657920 0.9528519
#> 326977 0.77732873 1.538560e-02 0.5404417 0.0540482806 0.6797655 0.9574773
#> 194648 0.75750536 2.617974e-01 0.7993647 0.0506919457 0.6849604 0.9481556
#> 28186 0.09757318 1.216955e-01 0.7376420 0.0495120780 0.6004570 0.9624822
#> 541301 0.42078021 7.415186e-02 0.7409960 0.0800154670 0.1740821 0.9423535
#> 198941 0.21665777 2.831475e-01 0.7758507 0.0531993063 0.2788675 0.9264058
#> 353985 0.65826787 1.983552e-01 0.3425181 0.0951237175 0.4892058 0.9620659
#> 187524 0.56751793 1.548248e-02 0.8102886 0.0489607405 0.1740821 0.1443733
#> 182054 0.09757318 8.023735e-03 0.3425181 0.0915912303 0.6067564 0.9581964
#> 175537 0.58774123 5.701687e-02 0.7585698 0.1067463972 0.1983390 0.9683967
#> 9753 0.75811340 3.091259e-01 0.6809264 0.0917822880 0.1740821 0.9233998
#> 194211 0.53938145 6.364261e-02 0.3425181 0.1095740486 0.1942578 0.9625979
#> 188518 0.77525433 1.321364e-02 0.7693914 0.0543396473 0.7045958 0.9633357
#> 189396 0.09757318 7.749439e-02 0.3425181 0.0758115149 0.6915717 0.9548474
#> 90487 0.62075598 7.583828e-02 0.3425181 0.0718001845 0.1983335 0.6202351
#> 203708 0.63472910 7.837297e-02 0.3630674 0.0668562166 0.5132607 0.9302210
#> 173965 0.79119486 1.835456e-01 0.5529293 0.1001370390 0.4846947 0.9671396
#> 194661 0.42236986 3.109899e-01 0.5297498 0.0725314269 0.1740821 0.9244561
#> 512309 0.77950254 3.636737e-02 0.3425181 0.0878097258 0.2893775 0.9596493
#> 170124 0.63627241 3.874625e-03 0.6320331 0.0331773059 0.1740821 0.9666128
#> 216862 0.09757318 2.982981e-01 0.7150792 0.0530007735 0.2745939 0.9386284
#> 352304 0.78498326 7.202617e-02 0.3425181 0.0950248398 0.4606316 0.9704272
#> 191306 0.68825922 1.051708e-05 0.6082806 0.0747296741 0.2793558 0.9684887
#> 191541 0.22301675 8.855131e-02 0.3425181 0.1086789858 0.5557820 0.9545399
#> 191547 0.16732644 1.376309e-01 0.7669113 0.0898095077 0.4748078 0.9313576
#> 195493 0.77994877 1.507536e-01 0.7597142 0.1183160423 0.6396434 0.9334123
#> cat_alcohol2 bin_migraine
#> 326792 1 0.93515350
#> 348374 1 0.85830796
#> 181016 1 0.71535938
#> 191687 1 0.93415070
#> 305760 1 0.81966766
#> 326977 1 0.91030297
#> 194648 1 0.90424419
#> 28186 1 0.19858613
#> 541301 1 0.92754595
#> 198941 1 0.66239978
#> 353985 1 0.92784606
#> 187524 1 0.58108418
#> 182054 1 0.50152223
#> 175537 1 0.90490747
#> 9753 1 0.92661893
#> 194211 1 0.88107913
#> 188518 1 0.93515855
#> 189396 1 0.07214349
#> 90487 1 0.91362981
#> 203708 1 0.91955312
#> 173965 1 0.90582383
#> 194661 1 0.79456765
#> 512309 1 0.67914530
#> 170124 1 0.83947467
#> 216862 1 0.92887565
#> 352304 1 0.87415708
#> 191306 1 0.89040123
#> 191541 1 0.90980992
#> 191547 1 0.89663536
#> 195493 1 0.93106150qval_specific_var function will be useful to retrieve
the q-values of a specific variable, bin_dog in this
example.
# Create an object to keep the table with q-values.
qvaltab <- qval(SOHPIEres)
# Retrieve a vector of q-values for a single variable of interest.
qval_specific_var(qvaltab = qvaltab, varname = "bin_dog")
#> bin_dog
#> 326792 0.73680874
#> 348374 0.52163971
#> 181016 0.57984394
#> 191687 0.09757318
#> 305760 0.18227112
#> 326977 0.77732873
#> 194648 0.75750536
#> 28186 0.09757318
#> 541301 0.42078021
#> 198941 0.21665777
#> 353985 0.65826787
#> 187524 0.56751793
#> 182054 0.09757318
#> 175537 0.58774123
#> 9753 0.75811340
#> 194211 0.53938145
#> 188518 0.77525433
#> 189396 0.09757318
#> 90487 0.62075598
#> 203708 0.63472910
#> 173965 0.79119486
#> 194661 0.42236986
#> 512309 0.77950254
#> 170124 0.63627241
#> 216862 0.09757318
#> 352304 0.78498326
#> 191306 0.68825922
#> 191541 0.22301675
#> 191547 0.16732644
#> 195493 0.77994877DCtaxa_tab will return a list containing of (1) names
and q-values of taxa that are significantly DC between two biological
conditions and (2) names of DC taxa only.
# Please do NOT forget to provide the name of variable in DCtaxa_tab(groupvar = )
# and the level of significance (0.3 in this example).
DCtaxa_tab <- DCtaxa_tab(qvaltab = qvaltab, groupvar = "bin_dog", alpha = 0.3)
#> Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
#> ℹ Please use one dimensional logical vectors instead.
#> ℹ The deprecated feature was likely used in the SOHPIE package.
#> Please report the issue to the authors.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
DCtaxa_tab
#> $DCtaxa_complete_tab
#> bin_dog
#> 191687 0.09757318
#> 305760 0.18227112
#> 28186 0.09757318
#> 198941 0.21665777
#> 182054 0.09757318
#> 189396 0.09757318
#> 216862 0.09757318
#> 191541 0.22301675
#> 191547 0.16732644
#>
#> $DCtaxa_names_only
#> [1] "191687" "305760" "28186" "198941" "182054" "189396" "216862" "191541"
#> [9] "191547"data(combineddietswap)
OTUtab = combineddietswap[ , 5:ncol(combineddietswap)]
phenodat = combineddietswap[ , 1:4] # first column is ID, so not using it.# Obtain indices for each groups
# (i.e. African-Americans from Pittsburgh (AAM) vs. Africans from rural South Africa (AFR))
# at baseline (time = 1) and at 29-days (time = 6)
# Group A1 for AAM at baseline.
newindex_A1 = which(combineddietswap$timepoint == 1 & combineddietswap$nationality == "AAM")
# Group A6 for AAM at 29-days.
newindex_A6 = which(combineddietswap$timepoint == 6 & combineddietswap$nationality == "AAM")
# Group B1 for AFR at baseline.
newindex_B1 = which(combineddietswap$timepoint == 1 & combineddietswap$nationality == "AFR")
# Group A6 for AFR at 29-days.
newindex_B6 = which(combineddietswap$timepoint == 6 & combineddietswap$nationality == "AFR")We are done with loading the data and obtaining indices for each group for each time point. We then move onto the analysis step!
The first step is to estimate and re-estimate association matrices
for each of groups A1, A6, B1, and A6 shown above. Of note, a list
output comprises data.frame objects as assomat
and reest.assomat.
est_asso_matA1 = asso_mat(OTUdat=OTUtab, group=newindex_A1)
est_asso_matA6 = asso_mat(OTUdat=OTUtab, group=newindex_A6)
est_asso_matB1 = asso_mat(OTUdat=OTUtab, group=newindex_B1)
est_asso_matB6 = asso_mat(OTUdat=OTUtab, group=newindex_B6)
## For each group, we take the difference of estimated
## association matrices between time points (29-days minus baseline).
asso_mat_diffA61 = est_asso_matA6$assomat - est_asso_matA1$assomat
asso_mat_diffB61 = est_asso_matB6$assomat - est_asso_matB1$assomat
## Similarly, for each group, we take the difference of re-estimated
## association matrices between time points (29-days minus baseline).
asso_mat_drop_diffA61 = mapply('-', est_asso_matA6$reest.assomat,
est_asso_matA1$reest.assomat, SIMPLIFY=FALSE)
asso_mat_drop_diffB61 = mapply('-', est_asso_matB6$reest.assomat,
est_asso_matB1$reest.assomat, SIMPLIFY=FALSE)Then, for each group, we calculate changes in network centrality of each taxa between the association matrices estimated from the whole data and between the association matrices re-estimated from the leave-one-out sample.
## For changes in network centrality between association matrices estimated from the whole data.
thetahat_grpA = thetahats(asso_mat_diffA61)
thetahat_grpB = thetahats(asso_mat_diffB61)
## For changes in network centrality between association matrices
## re-estimated from the leave-one-out sample.
thetahat_drop_grpA = sapply(asso_mat_drop_diffA61, thetahats)
thetahat_drop_grpB = sapply(asso_mat_drop_diffB61, thetahats)Next, we calculate jackknife pseudo-values.
# Sample sizes for each group.
n_A <- length(newindex_A1)
n_B <- length(newindex_B1)
# Jackknife pseudo-values.
thetatilde_grpA = thetatildefun(thetahat_grpA, thetahat_drop_grpA, n_A)
thetatilde_grpB = thetatildefun(thetahat_grpB, thetahat_drop_grpB, n_B)
thetatilde = rbind(thetatilde_grpA, thetatilde_grpB)Fit a robust regression regressing covariates on pseudo-values.
Please note that the metadata/phenotype data should contain a set of
predictors to be fitted only, so choose them wisely. Additionally, a
trimming proportion c of the least trimmed squares
estimator should be entered between 0.5 and 1.
Lastly, we can extract summary results from the fitted model from
fitmod object above. A list of data.frame
objects for coefficient estimates, p-values, and q-values will be
available.
summary.result = pseudoreg.summary(pseudo.reg.res=fitmod,
taxanames=colnames(OTUtab))
## In this study, the main grouping variable was nationality.
# We can use qval_specific_var to see the q-values only.
# of nationality. Alternatively, we further use DCtaxa_tab to see
# the DC taxa with their p-values.
qval_specific_var(summary.result$q_values, varname = "nationalityAFR")
#> nationalityAFR
#> Akkermansia 0.79686693
#> Alcaligenes faecalis et rel. 0.78674054
#> Allistipes et rel. 0.69836946
#> Anaerostipes caccae et rel. 0.14770650
#> Anaerotruncus colihominis et rel. 0.78537649
#> Anaerovorax odorimutans et rel. 0.14442740
#> Aquabacterium 0.77443040
#> Atopobium 0.71429048
#> Bacillus 0.73201348
#> Bacteroides fragilis et rel. 0.27790855
#> Bacteroides intestinalis et rel. 0.61071818
#> Bacteroides ovatus et rel. 0.64464014
#> Bacteroides plebeius et rel. 0.07859087
#> Bacteroides splachnicus et rel. 0.46914444
#> Bacteroides stercoris et rel. 0.20345191
#> Bacteroides uniformis et rel. 0.69551114
#> Bacteroides vulgatus et rel. 0.78929317
#> Bifidobacterium 0.10172409
#> Bilophila et rel. 0.65720446
#> Brachyspira 0.09363041
#> Bryantella formatexigens et rel. 0.70026602
#> Bulleidia moorei et rel. 0.62221009
#> Burkholderia 0.47840776
#> Butyrivibrio crossotus et rel. 0.73837070
#> Campylobacter 0.62643267
#> Catenibacterium mitsuokai et rel. 0.14789014
#> Clostridium (sensu stricto) 0.75159988
#> Clostridium cellulosi et rel. 0.10899752
#> Clostridium colinum et rel. 0.34550748
#> Clostridium difficile et rel. 0.79133223
#> Clostridium leptum et rel. 0.69717012
#> Clostridium nexile et rel. 0.12349428
#> Clostridium orbiscindens et rel. 0.11254808
#> Clostridium ramosum et rel. 0.33322716
#> Clostridium sphenoides et rel. 0.79114988
#> Clostridium stercorarium et rel. 0.33616140
#> Clostridium symbiosum et rel. 0.79877018
#> Collinsella 0.49226761
#> Coprobacillus catenaformis et rel. 0.46200718
#> Coprococcus eutactus et rel. 0.11076336
#> Corynebacterium 0.28908833
#> Desulfovibrio et rel. 0.10706825
#> Dialister 0.13381918
#> Dorea formicigenerans et rel. 0.11269705
#> Eggerthella lenta et rel. 0.73760977
#> Enterobacter aerogenes et rel. 0.62545132
#> Enterococcus 0.55591673
#> Escherichia coli et rel. 0.15027605
#> Eubacterium biforme et rel. 0.76046164
#> Eubacterium cylindroides et rel. 0.76680001
#> Eubacterium hallii et rel. 0.51621105
#> Eubacterium limosum et rel. 0.68990948
#> Eubacterium rectale et rel. 0.13169676
#> Eubacterium siraeum et rel. 0.50574749
#> Eubacterium ventriosum et rel. 0.78189507
#> Faecalibacterium prausnitzii et rel. 0.10713588
#> Fusobacteria 0.07845342
#> Haemophilus 0.69275215
#> Helicobacter 0.59934168
#> Klebisiella pneumoniae et rel. 0.73587124
#> Lachnobacillus bovis et rel. 0.39268182
#> Lachnospira pectinoschiza et rel. 0.30986631
#> Lactobacillus catenaformis et rel. 0.20969332
#> Lactobacillus gasseri et rel. 0.45110003
#> Lactobacillus plantarum et rel. 0.79859532
#> Lactobacillus salivarius et rel. 0.58507766
#> Lactococcus 0.73449794
#> Leminorella 0.76922522
#> Megamonas hypermegale et rel. 0.77376977
#> Megasphaera elsdenii et rel. 0.70981786
#> Mitsuokella multiacida et rel. 0.60600242
#> Moraxellaceae 0.80030969
#> Oceanospirillum 0.76075212
#> Oscillospira guillermondii et rel. 0.70775199
#> Outgrouping clostridium cluster XIVa 0.77089468
#> Oxalobacter formigenes et rel. 0.09185572
#> Papillibacter cinnamivorans et rel. 0.07889832
#> Parabacteroides distasonis et rel. 0.65980220
#> Peptococcus niger et rel. 0.76132551
#> Peptostreptococcus micros et rel. 0.07074831
#> Phascolarctobacterium faecium et rel. 0.66087262
#> Prevotella melaninogenica et rel. 0.56635133
#> Prevotella oralis et rel. 0.46333464
#> Prevotella ruminicola et rel. 0.52120651
#> Prevotella tannerae et rel. 0.50786635
#> Propionibacterium 0.69965920
#> Proteus et rel. 0.69413630
#> Roseburia intestinalis et rel. 0.79403689
#> Ruminococcus bromii et rel. 0.79215704
#> Ruminococcus callidus et rel. 0.10546188
#> Ruminococcus gnavus et rel. 0.56863660
#> Ruminococcus lactaris et rel. 0.39711515
#> Ruminococcus obeum et rel. 0.05178277
#> Serratia 0.75896883
#> Sporobacter termitidis et rel. 0.80126709
#> Staphylococcus 0.07671802
#> Streptococcus bovis et rel. 0.78008032
#> Streptococcus intermedius et rel. 0.28745234
#> Streptococcus mitis et rel. 0.75125468
#> Subdoligranulum variable at rel. 0.79803774
#> Sutterella wadsworthia et rel. 0.76291552
#> Tannerella et rel. 0.78041183
#> Uncultured Bacteroidetes 0.77365671
#> Uncultured Clostridiales I 0.77466109
#> Uncultured Clostridiales II 0.11381679
#> Uncultured Mollicutes 0.76629779
#> Uncultured Selenomonadaceae 0.54903006
#> Veillonella 0.59795559
#> Vibrio 0.22748351
#> Weissella et rel. 0.64872353
#> Xanthomonadaceae 0.66913221
#> Yersinia et rel. 0.51093007
DCtaxa_tab(summary.result$q_values, groupvar = "nationalityAFR", alpha=0.05)
#> $DCtaxa_complete_tab
#> [1] nationalityAFR
#> <0 rows> (or 0-length row.names)
#>
#> $DCtaxa_names_only
#> character(0)[1] Ahn S, Datta S. (2023). Differential Co-Abundance Network Analyses for Microbiome Data Adjusted for Clinical Covariates Using Jackknife Pseudo-Values. Under Review at BMC Bioinformatics.
[2] McDonald D. et al. (2018). American Gut: an Open Platform for Citizen Science Microbiome Research. mSystems. 3(3), e00031–18
[3] O’Keefe SJ. et al. (2015). Fat, fibre and cancer risk in African Americans and rural Africans. Nat Commun. 6, 6342