[feature] Transition the phenotype associations to non-parametric tests

.github/workflows/pythonpackage.yml

@@ -16,7 +16,7 @@ jobs:
     strategy:
       matrix:
         os: ["ubuntu-latest", "windows-latest", "macos-latest"]
-        python-version: ["3.8", "3.11"]
+        python-version: ["3.9", "3.11"]
     runs-on: "${{matrix.os}}"
 
     steps:
@@ -30,12 +30,6 @@ jobs:
       run: |
         python -m pip install --upgrade pip cryptography
         pip install wheel numpy Cython biom-format
-    - name: Install highspy via dev wheel
-      if: matrix.os == 'macos-latest'
-      run: |
-        PYTHON_VERSION=${{ matrix.python-version }}
-        PYTHON_VERSION=${PYTHON_VERSION//./}
-        pip install .github/wheels/highspy-1.5.3-cp${PYTHON_VERSION}-cp${PYTHON_VERSION}-macosx_10_9_x86_64.whl
     - name: Install MICOM
       run: pip install -e .
     - name: install CPLEX
@@ -52,7 +46,7 @@ jobs:
         pip install pytest pytest-cov
         pytest --cov=micom --cov-report=xml
     - name: Upload to codecov
-      if: matrix.os == 'ubuntu-latest' && matrix.python-version == 3.9
+      if: matrix.os == 'ubuntu-latest' && matrix.python-version == '3.11'
       uses: codecov/[email protected]
 
   release:

docs/source/viz.ipynb

@@ -8,7 +8,7 @@
     "\n",
     "MICOM provides a a set of visualizations that can be used with the outputs from MICOM workflows. Those visualizations are the same as provided by the [MICOM Qiime 2 plugin](https://github.com/micom-dev/q2-micom) but are delivered as single HTML files that bundles interactive graphics and raw data.\n",
     "\n",
-    "To create some more interesting figures here we will use a realistic example data set which is the output of running the MICOM grow workflow on a data set of 8 healthy fecal samples and 8 fecal samples from individuals with colorectal cancer taken from https://doi.org/10.1158/1940-6207.CAPR-14-0129. To see the interactive evrsion of a visualization you can click on the provided previews. All visualization contain download buttons to download the raw data used to generate the plot. "
+    "To create some more interesting figures here we will use a realistic example data set which is the output of running the MICOM grow workflow on a data set of 10 healthy fecal samples from the [iHMP cohort](https://www.hmpdacc.org/ihmp/). To see the interactive version of a visualization you can click on the provided previews. All visualization contain download buttons to download the raw data used to generate the plot. "
    ]
   },
   {
@@ -19,8 +19,8 @@
    "source": [
     "import micom.data as mmd\n",
     "\n",
-    "crc_results = mmd.crc_results()\n",
-    "crc_tradeoff = mmd.crc_tradeoff()"
+    "results = mmd.test_results()\n",
+    "tradeoff = mmd.test_tradeoff()"
    ]
   },
   {
@@ -40,7 +40,7 @@
    "source": [
     "from micom.viz import plot_tradeoff\n",
     "\n",
-    "pl = plot_tradeoff(crc_tradeoff, filename=\"tradeoff.html\")"
+    "pl = plot_tradeoff(tradeoff, filename=\"tradeoff.html\")"
    ]
   },
   {
@@ -58,7 +58,7 @@
     {
      "data": {
       "text/plain": [
-       "<micom.viz.core.Visualization at 0x7f17f0774490>"
+       "<micom.viz.core.Visualization at 0x7ff72ce1d950>"
       ]
      },
      "execution_count": 3,
@@ -116,7 +116,7 @@
    "source": [
     "from micom.viz import plot_growth\n",
     "\n",
-    "pl = plot_growth(crc_results, filename=\"growth_rates.html\")"
+    "pl = plot_growth(results, filename=\"growth_rates.html\")"
    ]
   },
   {
@@ -145,7 +145,7 @@
    "source": [
     "from micom.viz import plot_exchanges_per_sample\n",
     "\n",
-    "pl = plot_exchanges_per_sample(crc_results, filename=\"consumption.html\")"
+    "pl = plot_exchanges_per_sample(results, filename=\"consumption.html\")"
    ]
   },
   {
@@ -170,11 +170,25 @@
    "cell_type": "code",
    "execution_count": 7,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">[13:02:03] </span><span style=\"color: #800000; text-decoration-color: #800000\">WARNING </span> Not enough samples. Adjusting T-SNE perplexity to <span style=\"color: #008080; text-decoration-color: #008080; font-weight: bold\">5</span>.                           <a href=\"file:///home/cdiener/code/micom/micom/viz/exchanges.py\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">exchanges.py</span></a><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">:</span><a href=\"file:///home/cdiener/code/micom/micom/viz/exchanges.py#127\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">127</span></a>\n",
+       "</pre>\n"
+      ],
+      "text/plain": [
+       "\u001b[2;36m[13:02:03]\u001b[0m\u001b[2;36m \u001b[0m\u001b[31mWARNING \u001b[0m Not enough samples. Adjusting T-SNE perplexity to \u001b[1;36m5\u001b[0m.                           \u001b]8;id=5441;file:///home/cdiener/code/micom/micom/viz/exchanges.py\u001b\\\u001b[2mexchanges.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=120284;file:///home/cdiener/code/micom/micom/viz/exchanges.py#127\u001b\\\u001b[2m127\u001b[0m\u001b]8;;\u001b\\\n"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
    "source": [
     "from micom.viz import plot_exchanges_per_taxon\n",
     "\n",
-    "pl = plot_exchanges_per_taxon(crc_results, filename=\"niche.html\")"
+    "pl = plot_exchanges_per_taxon(results, filename=\"niche.html\")"
    ]
   },
   {
@@ -190,29 +204,45 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Investigating the relationship with the phenotype\n",
+    "### Investigating associations with a phenotype\n",
     "\n",
-    "Finally we may want to which functions relate to an observed phenotype. This can be done with the `plot_fit` function which will:\n",
+    "Finally we may want to which fluxes relate to an observed phenotype. This can be done with the `plot_association` function which will:\n",
     "\n",
-    "1. calculate overall production fluxes (total metabolite amount produced by the microbiota)\n",
-    "2. log-transform and standardize production fluxes\n",
-    "3. train a LASSO regression (continuous response) or LASSO logistic regression (binary response)\n",
-    "4. present the coefficients and differences in production fluxes"
+    "1. calculate overall production or consumption fluxes (total metabolite amount produced by the microbiota)\n",
+    "2. run non-parametric tests for each metabolite against the phenotype\n",
+    "3. control the flase discovery rate and report significantly associated metabolite fluxes\n",
+    "4. log-transform and standardize production fluxes\n",
+    "5. train a LASSO regression (continuous response) or LASSO logistic regression (binary response)\n",
+    "6. present the overall performance of the fluxes inpredicting the phenotype\n",
+    "\n",
+    "So you will get data on local (metabolite) and global (all fluxes) associations. \n",
+    "\n",
+    "To illustrate this we will create a mock phenotype that is correlated with propionate production. We will allow very high q values here. In a real analysis the default of 0.05 is more appropriate,"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [],
    "source": [
-    "from micom.viz import plot_fit\n",
+    "from micom.viz import plot_association\n",
+    "from micom.measures import production_rates\n",
+    "import numpy as np\n",
     "import pandas as pd\n",
     "\n",
-    "samples = pd.Series(crc_results.growth_rates.sample_id.unique())\n",
-    "status = pd.Series(\"healthy\", index=samples)\n",
-    "status[status.index.str.startswith(\"Cancer\")] = \"colorectal cancer\"\n",
-    "pl = plot_fit(crc_results, phenotype=status, variable_type=\"binary\", filename=\"fit.html\", min_coef=1e-3)"
+    "prod = production_rates(results)\n",
+    "propionate = prod.loc[prod.metabolite == \"ppa[e]\", \"flux\"]\n",
+    "propionate.index = prod.sample_id.unique()\n",
+    "high_propionate = propionate > np.median(propionate)\n",
+    "\n",
+    "pl = plot_association(\n",
+    "    results,\n",
+    "    phenotype=high_propionate,\n",
+    "    variable_type=\"binary\",\n",
+    "    filename=\"association.html\",\n",
+    "    fdr_threshold=0.5,\n",
+    ")"
    ]
   },
   {
@@ -221,14 +251,14 @@
    "source": [
     "The output will look something like the following.\n",
     "\n",
-    "[![fit](_static/fit.png)](_static/fit.html)"
+    "[![fit](_static/association.png)](_static/association.html)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "So we see some higher sulfate production in cancer cases but also production of harmful components as TMAO. However, with the low sample numbers we don't get a very good prediction based on fluxes."
+    "So we see we recovered propionate in the analysis but we would need larger sample sizes here."
    ]
   },
   {
@@ -255,7 +285,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.2"
+   "version": "3.11.8"
   }
  },
  "nbformat": 4,

micom/data/__init__.py

@@ -1,5 +1,6 @@
 """Submodule including some common data sets."""
 
+from ..workflows import load_results
 from os.path import split, join
 from numpy.random import randint
 import pandas as pd
@@ -68,7 +69,7 @@ def test_data(n_samples=4, uses_db=True):
     return data
 
 
-def crc_results():
+def test_results():
     """Return growth results for a realistic example data sets.
 
     This will give you access to the results as produced by the community tutorial
@@ -77,12 +78,12 @@ def crc_results():
     Returns
     -------
     GrowthResults
-        The results for running the `grow` workflow with a tradeoff of 0.3.
+        The results for running the `grow` workflow with a tradeoff of 0.8.
     """
-    return pickle.load(open(join(this_dir, "artifacts", "crc_results.pickle"), "rb"))
+    return load_results(join(this_dir, "artifacts", "results.zip"))
 
 
-def crc_tradeoff():
+def test_tradeoff():
     """Return tradeoff results for a realistic example data set.
 
     This will give you access to the tradeoffs as produced by the community tutorial

micom/data/templates/growth.html

@@ -6,11 +6,9 @@
     <title>Growth rates</title>
     <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/[email protected]/css/bulma.min.css">
     <link href="https://fonts.googleapis.com/css?family=Roboto:300,400,700&display=swap" rel="stylesheet">
-    <script defer src="https://use.fontawesome.com/releases/v5.3.1/js/all.js"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]/dist/d3.min.js"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
+    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
+    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
+    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
 
     <style>
     #vis {
@@ -44,33 +42,42 @@ <h1 class="title">
     <script>
       // Assign the specification to a local variable vlSpec.
       var vlSpec = {
-        $schema: "https://vega.github.io/schema/vega-lite/v4.json",
-        width: {{width}},
-        height: {{height}},
+        $schema: "https://vega.github.io/schema/vega-lite/v5.json",
+        width: 300,
+        height: {{18 * n_taxa}},
         data: {
             values: {{data | safe}}
         },
         selection: {
           pts: {type: "single", fields: ["sample_id"]}
         },
+        transform: [{
+          calculate: "random()-0.5",
+          as: "jitter"
+        }],
         mark: {
             type: "point",
             filled: true,
             size: 40,
             tooltip: {content: "data"}
         },
         encoding: {
-          y: {
+          x: {
             field: "growth_rate",
             type: "quantitative",
             scale: {type: "log"},
             axis: {title: "growth rate [1/h]", grid: false}
           },
-          x: {
+          y: {
             field: "taxon",
             type: "nominal",
-            sort: {op: "median", field: "growth_rate"},
-            axis: {title: null}
+            sort: {op: "median", field: "growth_rate", order: "descending"},
+            axis: {title: null, grid: true}
+          },
+          yOffset: {
+            field: "jitter",
+            type: "quantitative",
+            scale: {domain: [-1, 1]}
           },
           color: {
             condition: {
@@ -93,8 +100,9 @@ <h1 class="title">
 
       // Embed the visualization in the container with id `vis`
       vegaEmbed("#vis", vlSpec, {renderer: "svg"}).then(function(r) {
-        svg = d3.select("#vis svg");
-        svg.attr("width", "100%").attr("height", "auto");
+        svg = document.querySelector("#vis svg");
+        svg.setAttribute("width", "100%");
+        svg.setAttribute("height", "auto");
       });
 
       {% include "blobify.js" %}

micom/data/templates/reduced.html

@@ -6,11 +6,9 @@
     <title>Taxa niche/production overlap</title>
     <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/[email protected]/css/bulma.min.css">
     <link href="https://fonts.googleapis.com/css?family=Roboto:300,400,700&display=swap" rel="stylesheet">
-    <script defer src="https://use.fontawesome.com/releases/v5.3.1/js/all.js"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]/dist/d3.min.js"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
+    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
+    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
+    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
 
     <style>
     #vis {
@@ -27,7 +25,7 @@ <h1 class="title">
         Taxa niche/production overlap
       </h1>
       <p class="subtitle">
-        Each point denotes the reduced import/export fluxes for a single taxon
+        Each point denotes the reduced {{direction}} fluxes for a single taxon
         in a single sample. If two points are close to each other they consume
         (for direction="import") or produce (for direction="export") the same
         metabolites.<br>
@@ -47,7 +45,7 @@ <h1 class="title">
     <script>
       // Assign the specification to a local variable vlSpec.
       var vlSpec = {
-        $schema: "https://vega.github.io/schema/vega-lite/v4.json",
+        $schema: "https://vega.github.io/schema/vega-lite/v5.json",
         width: {{width}},
         height: {{height}},
         data: {
@@ -91,8 +89,9 @@ <h1 class="title">
 
       // Embed the visualization in the container with id `vis`
       vegaEmbed("#vis", vlSpec, {renderer: "svg"}).then(function(r) {
-        svg = d3.select("#vis svg");
-        svg.attr("width", "100%").attr("height", "auto");
+        svg = document.querySelector("#vis svg");
+        svg.setAttribute("width", "100%");
+        svg.setAttribute("height", "auto");
       });
 
       {% include "blobify.js" %}

micom/data/templates/sample_heatmap.html

@@ -6,11 +6,9 @@
     <title>Sample metabolite utilization/production</title>
     <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/[email protected]/css/bulma.min.css">
     <link href="https://fonts.googleapis.com/css?family=Roboto:300,400,700&display=swap" rel="stylesheet">
-    <script defer src="https://use.fontawesome.com/releases/v5.3.1/js/all.js"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]/dist/d3.min.js"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
-    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
+    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
+    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
+    <script src="https://cdn.jsdelivr.net/npm/[email protected]"></script>
 
     <style>
     #vis {
@@ -44,7 +42,7 @@ <h1 class="title">
     <script>
       // Assign the specification to a local variable vlSpec.
       var vlSpec = {
-        $schema: "https://vega.github.io/schema/vega-lite/v4.json",
+        $schema: "https://vega.github.io/schema/vega-lite/v5.json",
         width: {{width}},
         height: {{height}},
         data: {
@@ -76,8 +74,9 @@ <h1 class="title">
 
       // Embed the visualization in the container with id `vis`
       vegaEmbed("#vis", vlSpec, {renderer: "svg"}).then(function(r) {
-        svg = d3.select("#vis svg");
-        svg.attr("width", "100%").attr("height", "auto");
+        svg = document.querySelector("#vis svg");
+        svg.setAttribute("width", "100%");
+        svg.setAttribute("height", "auto");
       });
 
       {% include "blobify.js" %}