From c0cc17104523963a60e8d63be76dca03c2ef804a Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Mon, 1 Sep 2025 12:44:17 -0700 Subject: [PATCH 1/3] draft of vignette --- .gitignore | 1 + DESCRIPTION | 8 ++- R/example_data.R | 18 ++++++ data/example_data.rda | Bin 0 -> 1598 bytes man/example_data.Rd | 30 +++++++++ vignettes/.gitignore | 2 + vignettes/intro_enviromtx.Rmd | 113 ++++++++++++++++++++++++++++++++++ 7 files changed, 171 insertions(+), 1 deletion(-) create mode 100644 R/example_data.R create mode 100644 data/example_data.rda create mode 100644 man/example_data.Rd create mode 100644 vignettes/.gitignore create mode 100644 vignettes/intro_enviromtx.Rmd diff --git a/.gitignore b/.gitignore index 1968f99..1dab18d 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ .httr-oauth .DS_Store enviromtx.Rproj +inst/doc diff --git a/DESCRIPTION b/DESCRIPTION index fa8e8c7..d338f5c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: enviromtx Title: Predicting Species-Species Interactions from Environmental Metatranscriptomics Data -Version: 1.0.0 +Version: 1.1.0 Authors@R: c(person("Amy D", "Willis", email = "adwillis@uw.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-2802-4317")), person("Sarah", "Teichman", role = "aut"), @@ -8,6 +8,8 @@ Authors@R: Description: Tools to answer the question: Does the abundance of a given species impact another species' expression of a gene? License: MIT + file LICENSE Suggests: + knitr, + rmarkdown, testthat (>= 3.0.0) Config/testthat/edition: 3 Encoding: UTF-8 @@ -21,3 +23,7 @@ Imports: raoBust (>= 1.1.1), tibble Remotes: statdivlab/raoBust +Depends: + R (>= 3.5) +LazyData: true +VignetteBuilder: knitr diff --git a/R/example_data.R b/R/example_data.R new file mode 100644 index 0000000..d72d8b9 --- /dev/null +++ b/R/example_data.R @@ -0,0 +1,18 @@ +#' Dataset for use in vignettes: example_data +#' +#' A dataset containing metatranscriptome data for responder taxon *Pseudoalteromonas sp.* and companion taxon *Chaetoceors dichaeta*. +#' +#' @format A data frame with 22 rows and 9 variables: +#' \describe{ +#' \item{sample_name_r}{Sample id with replicate information (character)} +#' \item{xx}{Observed abundance of responder taxon (integer)} +#' \item{xstar}{Observed abundance of companion taxon (integer)} +#' \item{sample_name}{Sample id (character)} +#' \item{temp}{Temperature in Celsius of ocean from which sample was taken (integer)} +#' \item{K02520_counts_sum}{Observed abundance of gene K02520 expressed by the responder taxon (integer)} +#' \item{K01006_counts_sum}{Observed abundance of gene K01006 expressed by the responder taxon (integer)} +#' \item{K01497_counts_sum}{Observed abundance of gene K01497 expressed by the responder taxon (integer)} +#' \item{K03106_counts_sum}{Observed abundance of gene K03106 expressed by the responder taxon (integer)} +#' } +#' @source Bartolek et al. "Functional patterns of microbial interaction in the North Pacific revelaed through statistical modeling of environmental metatranscriptomes." (2025+) +"example_data" diff --git a/data/example_data.rda b/data/example_data.rda new file mode 100644 index 0000000000000000000000000000000000000000..24e098f8c1c0345fe14a3b0e296125f4245ab1ca GIT binary patch literal 1598 zcmV-E2EqA4T4*^jL0KkKSwOl-&;SEa|NsC0|Nnpg|NsBL{{8>||Ns8~|NsAgfB*l# z|G)qK_y5oX{(7#5z1;P8n-6t}g%D~sk%2K9Xwyt5sKml#H1!Q0iK7vd(r6lDG{6M) zGz^+y4IZN)^qL>2hK8P@rjJN^fb|bZ$?82dH1v#205u1xrk+sr)FM+70F4X;44PnQ z0F4@GV3{<;z#0GmA)`zHBM8$I)Wpe}CMEzA00h7QX*AO#OhzDNFd9r6q{T5k38ZQ1 zq#lUV)byU9$?9!QsisVxrbcQSXgx-N!UjeVXkk4QCPN9PjDfU6)C1Hs88i(s9*}5c z#-JGip`#Fb1R4}202*oNnrWw~8X9ENKxESCU@00SV%+7m#?7=R57AOO(Q zAjB9=0Ac}B5JAi^*J00yRIF)eK^Z{_egA%TQQ#>*mLoGGaeO3&jdl@`Vb)||_y#N(Xd@ypYCz<{ zowUsn5fKd05iJ5F?BHHZf$@M0!V@6~*WuhISP55w!mgp*S%1C2KtNox|H)t}ZsfDG zl3vjVQGJAg5pCj&1?zb>3zzy*$P5+0+db&T5a5#C+1w6PH{d!UIc;qTNR#~S4@)?l zj@mdtLFkzs@y){eRtf&2?iE_61d{*;Vps#f2}OQ1#Pe(Q;2i4p9$F;2MWe)1zAj{N zr_TK&^I3mC=Z-Gh?4G|U`Z@n5c)8Nnz4bb6@ed31kgB2O5D5~Xk!b{r-p>Oe28p2t z%zZ$mE&5?7qY(Nfq@EIb9l=dz-%tg9)W3izVP@9Xy}#tVUax?W5fCoCJAjGDACwL3 zqu^9mM`>`flw|;6fWMFpF%(}6{J8vCDa)Io;Dv#SAOCo)J`S5R#+buS=A$ObOvku+ zj2bSO{k;N8R7Di2!a<^ci;s=s1enau=--T=Yf^K#vZA;R4yu`4#bH+!l|A0!BF&kO zYe=OF`z}w+A!0`SCegQA-`!+E3&E9KN+U;lc#jKgCX1y(-k7AzOS^dh5QF@~oOq`= zv_{G=xY8`S`bhlb)xxf$6Bssz8gWPf3^semO4E+wXVxaY&CFEciGu_)GcdypGGxh< z3^HWNmL^J0m0Euv7UXpFBS9jj2^B|yA`M2hgzECJi#r&8vb|~3|Dt8$RO{@-BO(Rm zFc^}Gs;a82YN&{+fkH55W({B*HjN7~X#|cx`=@>h64DpB8}!k*y+wN(;TUeGBXN(U z>oY|&!|2g?K4O;1m_3>pFo=dy_anOI3fObtbUJ^Bu{8SYNCf4GE?e3y^6wS3_uPO5 zgLnJ;K}nGjlwl3agGw+AAQ%VGXeF*Y@G&2^lR`fF$SXtzl)V0v0c$>QWI$B_w1A}1 z2|^^xat8p_YikAF>JT0hsBBNz-I78yJ-!&=uq6qI_cY8MnWTL+k^z+QnRDPtln$@K z0f+Y%20R5G9E(EGU=#b@P`u>ZzrUJ8D$H&6}U~%2u(ZEXdS4 zlrI6qSK#-?Jm(37tc-wQ%j1U33v&m+?R-IG*j@|j#fqwp>bSXtTXfIDUU@pRtEJa> zXAdFnv zF!zhvjU@-GQW7~#){)vg;s;jr<5*TYg?pC266zc8%H8yNRCVxjGdk7A+ztq=lXX%V z{)bwJGqFTYP$=DBl}^SOUaAbCx5gjla4nXjaAtg*M|5Ji0Q+e%gvmwufEC73qw)JikgMN(tik0HO}_iMl9Eb&gN`b=H|ih$0E`*#l}v%cYK3+SJ*FY6PE% wz``${Tu6jz65KHgIIT<(+!vExJ~E`?N|Az~=t8S?0^jj>BvXY61Pi2%05nJArvLx| literal 0 HcmV?d00001 diff --git a/man/example_data.Rd b/man/example_data.Rd new file mode 100644 index 0000000..7c3e511 --- /dev/null +++ b/man/example_data.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/example_data.R +\docType{data} +\name{example_data} +\alias{example_data} +\title{Dataset for use in vignettes: example_data} +\format{ +A data frame with 22 rows and 9 variables: +\describe{ +\item{sample_name_r}{Sample id with replicate information (character)} +\item{xx}{Observed abundance of responder taxon (integer)} +\item{xstar}{Observed abundance of companion taxon (integer)} +\item{sample_name}{Sample id (character)} +\item{temp}{Temperature in Celsius of ocean from which sample was taken (integer)} +\item{K02520_counts_sum}{Observed abundance of gene K02520 expressed by the responder taxon (integer)} +\item{K01006_counts_sum}{Observed abundance of gene K01006 expressed by the responder taxon (integer)} +\item{K01497_counts_sum}{Observed abundance of gene K01497 expressed by the responder taxon (integer)} +\item{K03106_counts_sum}{Observed abundance of gene K03106 expressed by the responder taxon (integer)} +} +} +\source{ +Bartolek et al. "Functional patterns of microbial interaction in the North Pacific revelaed through statistical modeling of environmental metatranscriptomes." (2025+) +} +\usage{ +example_data +} +\description{ +A dataset containing metatranscriptome data for responder taxon \emph{Pseudoalteromonas sp.} and companion taxon \emph{Chaetoceors dichaeta}. +} +\keyword{datasets} diff --git a/vignettes/.gitignore b/vignettes/.gitignore new file mode 100644 index 0000000..097b241 --- /dev/null +++ b/vignettes/.gitignore @@ -0,0 +1,2 @@ +*.html +*.R diff --git a/vignettes/intro_enviromtx.Rmd b/vignettes/intro_enviromtx.Rmd new file mode 100644 index 0000000..1b1cfdb --- /dev/null +++ b/vignettes/intro_enviromtx.Rmd @@ -0,0 +1,113 @@ +--- +title: "Introduction to enviromtx" +output: rmarkdown::html_vignette +author: "Sarah Teichman and Amy Willis" +date: "`r Sys.Date()`" +vignette: > + %\VignetteIndexEntry{Introduction to enviromtx} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +First we will install `enviromtx` if haven't already. + +```{r, eval = FALSE} +# if (!require("remotes", quietly = TRUE)) +# install.packages("remotes") +# +# remotes::install_github("statdivlab/enviromtx") +``` + +Next we will load `enviromtx`, as well as relevant packages in the `tidyverse` suite. + +```{r setup} +library(enviromtx) +``` + +## Introduction + +This vignette serves as an introduction to `enviromtx`, a statistical method to model microbial interactions and their transcriptional signatures from metatranscriptome data. We will demonstrate this method on data from Bartolek et al., collected from the North Pacific Ocean. + +The main function in `enviromtx` is `fit_mgx_model()`. This function considers a pair of taxa, one designated the "responder" taxon and the other the "companion" taxon, and tests whether the responder taxon shows a difference in gene expression associated with the abundance of the companion taxon. To use this function, you will need a data frame with the following variables: + +- `yy`: expression data for a gene expressed by the responder taxon. +- `xx`: observed abundance of the responder taxon. +- `xstar`: observed abundance of the companion taxon. + +Optionally, the data frame can also include the following variables: + +- `replicates`: technical replicate information for samples, if applicable. +- `wts`: optional non-negative weights for samples. These could be sequencing depths, to put more emphasis on deeply sequenced samples. +- enviromental covariates: any relevant environmental covariates for the samples. + +```{r} +data("example_data", package = "enviromtx") +head(example_data) +``` + +This `example_data` involves *Pseudoalteromonas sp.* as the responder taxon and *Chaetoceors dichaeta* as the companion taxon. We can see that this `example_data` contains each of these required variables including expression data for four KOfams (representing genes), as well as `sample_name_r` which gives replicate information and the environmental covariate `temp`. + +## Estimation and testing + +Now that we've looked at the data, lets fit a model with `fit_mgx_model()` for KOfam K02520. + +```{r} +model_K02520 <- fit_mgx_model(enviro_df = example_data, + yy = "K02520_counts_sum", + formula = ~ temp, + replicates = "sample_name") +model_K02520 +``` + +Here we can see the output of `fit_mgx_model()`. We have a row for `predictor` and a row for `temp`, as well as a row for `correlation:alpha` which we will not focus on in this vignette. We have columns for `Estimate`, `Robust Std Error`, upper and lower boundaries of a 95\% confidence interval, a robust Wald p-value, and a robust score p-value. + +You may notice that while `temp` was included as a covariate in the model, `predictor` was not. `predictor` is included in all `enviromtx` models. It is a transformation of `xx` and `xstar`, specifically `predictor` = log(`xx` / `xstar`). This is the log ratio of abundance of the responder taxon to the companion taxon. + +Let's interpret our estimates. For `predictor` we have an estimate of $-0.18$. This means that the average expression of function K02520 per unit abundance of the responder taxon has a fold-difference of $\exp(-0.18)$ = $0.84$ associated with a one unit increase in the relative abundance of the companion taxon to the responder taxon, holding temperature constant. + +For `temp` we have an estimate of $-0.27$. This means that the average expression of function K02520 per unit abundance of the responder taxon has a fold-difference of $\exp(-0.27) = 0.76$ associated with a one degree increase in temperature, holding the relative abundance of companion taxon to responder taxon constant. + +Next, let's consider the p-values. By default `fit_mgx_model()` runs both a robust Wald test and a robust score test. Both test the same hypothesis: that the average expression of the gene per unit abundance of the responder taxon does not change with a one unit increase in the relative abundance of the companion taxon to the responder taxon (for the test of the `predictor` covariate). We recommend doing inference using the robust score tests, due to their better error rate control with small sample sizes and sparse gene expression data. + +Let's now fit models for the other three KOfams. + +```{r} +model_K01006 <- fit_mgx_model(enviro_df = example_data, + yy = "K01006_counts_sum", + formula = ~ temp, + replicates = "sample_name") +model_K01006 +model_K01497 <- fit_mgx_model(enviro_df = example_data, + yy = "K01497_counts_sum", + formula = ~ temp, + replicates = "sample_name") +model_K01497 +model_K03106 <- fit_mgx_model(enviro_df = example_data, + yy = "K03106_counts_sum", + formula = ~ temp, + replicates = "sample_name") +model_K03106 +``` + +We can see that K01006 has similar estimates for `predictor` and `temp` as K02520, while KOfams K01497 and K03106 have estimates with smaller magnitudes. Although robust score test p-values are smaller for KOfams K01006 and K02520, none are significant at an alpha threshold of $0.05$ (even before accounting for multiple testing). + +Now you know how to use the `enviromtx` package for your metatranscriptome analyses! If you have additional questions, feel free to post an [issue](https://github.com/statdivlab/enviromtx/issues). If you want to know more of the details of how `enviromtx` fits these models, continue to the technical appendix below. + +## Technical Appendix + +The package `enviromtx` calls the package `raoBust`, which fits gee or glm models (for data with and without replicates) and runs robust score tests. When there are not replicates, `raoBust` calls the `glm()` function from the `stats` package. However, when there are replicates, `raoBust` calls `geelm()` from the `geeasy` package, and if this fails then it calls `glm()`. + +When `geelm()` fails and `glm()` is called, robust score tests cannot be fit by default. When `geelm()` fails, the output for `fit_mgx_model()` will contain `NA` values for the entire `Robust Score p` column of the output. However, the robust score tests can be fit if the user inputs the argument `cluster_corr_coef` to `fit_mgx_model()`. Our suggested workflow is running `fit_mgx_model()` for all taxon pairs of interest, and then for pairs for which `geelm()` successfully runs, finding the average `correlation:alpha` value. Then, `fit_mgx_model()` can be rerun for pairs for which `geelm()` failed, inputting the average `correlation:alpha` value as `cluster_corr_coef`. + +## Citation + +If you use `enviromtx`, please cite the preprint: + +Bartolek et al. "Functional patterns of microbial interaction in the North Pacific revelaed through statistical modeling of environmental metatranscriptomes." (2025+) From 79370324c7667dc4c41930ff9c677198052f22ee Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Mon, 1 Sep 2025 12:47:46 -0700 Subject: [PATCH 2/3] fix small error in citation --- R/example_data.R | 2 +- man/example_data.Rd | 2 +- vignettes/intro_enviromtx.Rmd | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/example_data.R b/R/example_data.R index d72d8b9..16119c3 100644 --- a/R/example_data.R +++ b/R/example_data.R @@ -14,5 +14,5 @@ #' \item{K01497_counts_sum}{Observed abundance of gene K01497 expressed by the responder taxon (integer)} #' \item{K03106_counts_sum}{Observed abundance of gene K03106 expressed by the responder taxon (integer)} #' } -#' @source Bartolek et al. "Functional patterns of microbial interaction in the North Pacific revelaed through statistical modeling of environmental metatranscriptomes." (2025+) +#' @source Bartolek et al. "Functional patterns of microbial interaction in the North Pacific revealed through statistical modeling of environmental metatranscriptomes." (2025+) "example_data" diff --git a/man/example_data.Rd b/man/example_data.Rd index 7c3e511..5d16dd3 100644 --- a/man/example_data.Rd +++ b/man/example_data.Rd @@ -19,7 +19,7 @@ A data frame with 22 rows and 9 variables: } } \source{ -Bartolek et al. "Functional patterns of microbial interaction in the North Pacific revelaed through statistical modeling of environmental metatranscriptomes." (2025+) +Bartolek et al. "Functional patterns of microbial interaction in the North Pacific revealed through statistical modeling of environmental metatranscriptomes." (2025+) } \usage{ example_data diff --git a/vignettes/intro_enviromtx.Rmd b/vignettes/intro_enviromtx.Rmd index 1b1cfdb..b32e37a 100644 --- a/vignettes/intro_enviromtx.Rmd +++ b/vignettes/intro_enviromtx.Rmd @@ -110,4 +110,4 @@ When `geelm()` fails and `glm()` is called, robust score tests cannot be fit by If you use `enviromtx`, please cite the preprint: -Bartolek et al. "Functional patterns of microbial interaction in the North Pacific revelaed through statistical modeling of environmental metatranscriptomes." (2025+) +Bartolek et al. "Functional patterns of microbial interaction in the North Pacific revealed through statistical modeling of environmental metatranscriptomes." (2025+) From 475147cf2c17624ab1453575d98cb8a9e52d05c6 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Fri, 5 Sep 2025 09:42:58 -0700 Subject: [PATCH 3/3] increment version for minor version update --- DESCRIPTION | 2 +- NEWS.md | 8 ++++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index d338f5c..ef78d3b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: enviromtx Title: Predicting Species-Species Interactions from Environmental Metatranscriptomics Data -Version: 1.1.0 +Version: 1.2.0 Authors@R: c(person("Amy D", "Willis", email = "adwillis@uw.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-2802-4317")), person("Sarah", "Teichman", role = "aut"), diff --git a/NEWS.md b/NEWS.md index da2d63a..0dd178c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,11 @@ +# enviromtx 1.2.0 + +This is a minor release that provides a vignette and an example dataset. + +## Minor changes + +* `example_data` is added and documented, "intro_enviromtx.Rmd" vignette is added. + # enviromtx 1.1.0 This is a minor release that provides the option of centering covariates.