vignettes/sexContrastMicrobiomeAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sex-related differences in the human microbiome using curatedMetagenomicData 3 and the python 3 language"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Paolo Manghi, (paolo.manghi@unitn.it)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains the instructions to run a meta-analysis of sex-related contrasts in the human gut microbiome, using `curatedMetagenomicDataTerminal` and a set of freely-available python programs. If not running in docker, see the [installation instructions](installation.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described here, we are now going to\n",
    "\n",
    "1. create a folder called `species_abundances_from_cMD3Terminal`\n",
    "2. go in that directory\n",
    "3. download all the taxonomic profiles using `curatedMetagenomicDataTerminal`\n",
    "\n",
    "Step 3 may take some time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "mkdir $HOME/species_abundances_from_cMD3Terminal\n",
    "cd $HOME/species_abundances_from_cMD3Terminal\n",
    "curatedMetagenomicData -m \"*relative_abundance\"\n",
    "cd $HOME"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The flag \"-m\" will attach the per-sample metadata available in curatedMetagenomicData 3 to their taxonomic profiles.\n",
    "We now switch to a python 3 set of instructions that can be used to perform the main analysis of Figure 2, panel a, of the paper \"***\". \n",
    "\n",
    "The following instructions are meant as an example and also to clarify some aspect of the inner-code: the whole analysis is avaiable, as explained, via command-line programs at [https://github.com/waldronlab/curatedMetagenomicAnalyses/tree/main/python_tools] and can be all completed from the shell.\n",
    "\n",
    "## Imports section\n",
    "\n",
    "So we must first import the basic modules to work with cMDTerminal in a meta-analytical way. Note: If your `python_modules` and `python_tools` directories are in a different location, correct the paths below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "sys.path.append(os.path.join(os.environ[\"HOME\"], \"curatedMetagenomicAnalyses/python_modules\"))\n",
    "sys.path.append(os.path.join(os.environ[\"HOME\"], \"curatedMetagenomicAnalyses/python_tools\"))\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "## The following two modules handle either \n",
    "## a meta-analysis based on a std. linear model\n",
    "## OR a std. differential abundance analysis.\n",
    "## This is useful as some biological questions are\n",
    "## commonly explorable just in one dataset.\n",
    "from meta_analysis_runner import meta_analysis_with_linear_model\n",
    "from meta_analysis_runner import analysis_with_linear_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construction of the dataset\n",
    "\n",
    "Run the following code to see its help page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "python3 $HOME/curatedMetagenomicAnalyses/python_tools/meta_analysis_data.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "## The function `select` builds a dataset according to the requirements/constraints passed by the user. \n",
    "## Now, the `select` function can be imported as a module and used:\n",
    "## To to this, it is enough to pass it a dictionary of the parameters,\n",
    "## as follows:\n",
    "\n",
    "from meta_analysis_data import select\n",
    "\n",
    "## Following from the \"Help\" page, we next try to build \n",
    "## a dataset to study the Sex-contrast in the healthy, human, gut microbiome:\n",
    "\n",
    "params = {\n",
    "    'input_folder': os.path.join(os.environ[\"HOME\"], \"species_abundances_from_cMD3Terminal\"),\n",
    "    \"output_dataset\": os.path.join(os.environ[\"HOME\"], \"a_dataset_for_the_sex_contrast_in_gut_species.tsv\"),\n",
    "    \"min\": [\"age:16\"],\n",
    "    \"max\": [],\n",
    "    \"cat\": [\"study_condition:control\", \"body_site:stool\"], \n",
    "    \"multiple\": -1,\n",
    "    \"min_perc\": [\"gender:25\"],\n",
    "    \"cfd\":[\"BMI\"], \n",
    "    \"iqr\": [],\n",
    "    \"minmin\": \"gender:40\",\n",
    "    \"study_identifier\": \"study_name\", \n",
    "    \"verbose\": False, \n",
    "    \"debug\": False,\n",
    "    \"binary\": [],\n",
    "    \"search\": [],\n",
    "    \"exclude\": []\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the parameters that are needed to narrow this dataset. Normally, these are inserted via command-line, \n",
    "but python is flexible on this.\n",
    "\n",
    "Now we write a dataset table: a dataset-table is meant here as a table with Samples as columns-ID, and metadata + features as index (row-names). Features are then distingushed from metadata based on the \"feat-id\" keyword."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "select(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, you should have saved a file named `a_dataset_for_the_sex_contrast_in_gut_species.tsv`. \n",
    "\n",
    "At June, 2021, this dataset should store 4007 sample-IDs (columns) + 1 column index.\n",
    "The current dataset is not yet usable to perform a meta-analysis as the compositional \n",
    "data returned by MetaPhlAn3 software are not recommended for work with most statistical methods.\n",
    "\n",
    "Among the transformation more often applied, there are the *Centered-Log-Ratio* (CLR) and the *arcsin-square-root* of \n",
    "the relative-abundance-corresponding proportions. This last method is widely applied in meta-analyses on \n",
    "microbiome-related questions as it considers equally all the zeros independent from the dataset.\n",
    "\n",
    "We thus will apply the arcsin-square root, but a utility script is currently available to apply both. Run the following code to read its help page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "python3 $HOME/curatedMetagenomicAnalyses/python_tools/apply_a_transform.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## The general procedure to transform only the **features**\n",
    "## in a dataset-merged-table is quite simple.\n",
    "## We must first read the table we have created:\n",
    "sex_contrast_dataset = pd.read_csv(os.path.join(os.environ[\"HOME\"], \"a_dataset_for_the_sex_contrast_in_gut_species.tsv\"), sep=\"\\t\", \\\n",
    "        header=0, index_col=0, low_memory=False, engine=\"c\").fillna(\"NA\")\n",
    "\n",
    "## We then identify it's features\n",
    "feats = [j for j in sex_contrast_dataset.index.tolist() if (\"s__\" in j)]\n",
    "\n",
    "## We transform the data\n",
    "for sample in sex_contrast_dataset.columns.tolist():\n",
    "    sex_contrast_dataset.loc[feats, sample] = sex_contrast_dataset.loc[feats, sample].values.astype(float)\n",
    "    sex_contrast_dataset.loc[feats, sample] /= np.sum(sex_contrast_dataset.loc[feats, sample].values)\n",
    "    sex_contrast_dataset.loc[feats, sample] = np.arcsin(np.sqrt(\\\n",
    "        sex_contrast_dataset.loc[feats, sample].values.astype(float)))\n",
    "    \n",
    "## Now that we have transformed the data we can save a table which is suitable for the analyses:\n",
    "sex_contrast_dataset.index.name = \"sample_id\"\n",
    "sex_contrast_dataset.to_csv(os.path.join(os.environ[\"HOME\"], \"a_dataset_for_the_sex_contrast_in_gut_species_arcsin.tsv\"), sep=\"\\t\", \\\n",
    "        header=True, index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's now time to use the dataset we have created to perform a meta-analysis.\n",
    "The meta-analysis we'll perform is a standard meta-analysis (e.g. a weighted-average of several \n",
    "effect-sizes, computed over several, independent populations).\n",
    "\n",
    "The effect-size will be in the class of the differences of means (Cohen's d) and will characterize \n",
    "the difference between males and females (men and women).\n",
    "\n",
    "The only difference, with respect to a canonical paradigm of meta-analysis, is that we will extract, \n",
    "each time, **the Cohen'd effect size from the sex-relative coefficient of an ordinary least squares (OLS)\n",
    "model**. Beside being averaged over 15 populations, these effect-sizes will be so adjusted by age and by BMI\n",
    "in addition.\n",
    "\n",
    "Now, the program allowing you to run directly the meta-analysis is called `metaanalyze.py`, and \n",
    "asks for a short set of parameters in order to understand how to perform the desired meta-analysis.\n",
    "\n",
    "We will now run, **manually**, the basic steps of this program, in order to show them explicitly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "## First, we reproduce part of the program's import section\n",
    "from meta_analyses import paule_mandel_tau\n",
    "from meta_analyses import RE_meta_binary\n",
    "from meta_analyses import RE_meta\n",
    "\n",
    "## we specify the important parameters\n",
    "outfile = \"sex_4007_individuals_meta_analysis\"\n",
    "heterogeneity = \"PM\" ## this corresponds to the default of the program\n",
    "study_id = \"study_name\" ## this is needed in order to segregate the populations\n",
    "type_of_analysis = \"CLS\" ## CLASSIFY (SEE AFTER)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to specifiy a `formula` for the **model** that will be optimized before to compute the effect-size.\n",
    "The model will be applied to all the features in `feats`. Given this, we don't have to specify the Y, but only the **X(s)**.\n",
    "\n",
    "The first X will be the main predictor, and the effect-sizes will be computed on that.\n",
    "The following X will be used to adjuste the model.\n",
    "\n",
    "NOTE THAT:\n",
    "1) Though this is not needed for the main one, **categorical** covariates must be explicitly indicated via **C(name of c.)**.\n",
    "\n",
    "2) Though the first variable is automatically detected as categorical (if `type_of_analysis == \"CLS\"`), we must\n",
    "    specify the **positive** and the **negative** direction of the analysis we are going to perform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "formula = \"gender + age + BMI\"\n",
    "control_side, case_side = \"female\", \"male\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now can run the meta-analysis main module, `metaanalyze.py`.\n",
    "The module can be run from command-line with a minimum parameters, provided that a suitable dataset has \n",
    "been built.\n",
    "\n",
    "Run the following to view the help page for `metaanalyze.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "python3 $HOME/curatedMetagenomicAnalyses/python_tools/metaanalyze.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "ma = meta_analysis_with_linear_model(\\\n",
    "    sex_contrast_dataset, \n",
    "    formula, \\\n",
    "    study_id, \\\n",
    "    feats, \\\n",
    "    outfile + \".tsv\", \\\n",
    "    type_of_analysis, \\\n",
    "    heterogeneity, \\\n",
    "    pos=case_side, \\\n",
    "    neg=control_side\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ma.random_effect_regression_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have computed a meta-analysis, we want to plot it. Meta-analyses and Metagenomics are both dealing\n",
    "each with a problem of multidimensionality of the results, so coupling them requires to do our best to fit the \n",
    "different layers of results in one figure.\n",
    "\n",
    "In particular, we can plot many forest-plots in a single figure with the script: `draw_figure_with_ma.py`.\n",
    "Run the following to see its help page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "python3 $HOME/curatedMetagenomicAnalyses/python_tools/draw_figure_with_ma.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to call the program from python, we have to prepare a dictionary, faking its command-line arguments.\n",
    "Also, the program is able to plot several meta-analysis in the same set of axes, though we will for now \n",
    "limit ourselves to the **basic usage**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure_params = {\n",
    "    \"metaanalysis\": [os.path.join(os.environ[\"HOME\"], \"sex_4007_individuals_meta_analysis.tsv\")],\n",
    "    \"narrowed\": True,\n",
    "    \"boxes\": True,\n",
    "    \"relative_abundances\": [os.path.join(os.environ[\"HOME\"], \"a_dataset_for_the_sex_contrast_in_gut_species.tsv\")],\n",
    "    \"imp\": 30, \n",
    "    \"how\": \"first\",\n",
    "    \"positive_direction\": \"sex: male\",\n",
    "    \"negative_direction\": \"sex: female\",\n",
    "    \"x_axis\": \"standardized mean difference\",\n",
    "    \"y_axis\": \"\",\n",
    "    \"title\": \"meta-analysis of sex-related microbial species\",\n",
    "    \"e_suff\": [\"_Effect\"],\n",
    "    \"q_suff\": [\"_Qvalue\"],\n",
    "    \"prevalence\": [\"\"],\n",
    "    \"min_prevalence\": [0.01],\n",
    "    \"min_ab\": [0.000],\n",
    "    \"min_studies\": [4],\n",
    "    \"markers\": False,\n",
    "    \"outfile\": None,\n",
    "    \"random_effect\": [\"RE_Effect\"],\n",
    "    \"confint\": [\"RE_conf_int\"],\n",
    "    \"random_effect_q\": [\"RE_Effect_Qvalue\"],\n",
    "    \"color_red\": [\"goldenrod\"],\n",
    "    \"color_blue\": [\"dodgerblue\"],\n",
    "    \"color_black\": [\"black\"],\n",
    "    \"diam_marker\": [\"D\"],\n",
    "    \"important_lines\": [0.20],\n",
    "    \"a_single\": 0.2,\n",
    "    \"a_random\": 0.05, \n",
    "    \"dotsize\": 9,\n",
    "    \"diamsize\": 17,\n",
    "    \"neg_max_rho\": 0.8,\n",
    "    \"pos_max_rho\": 0.8,\n",
    "    \"legloc\": \"best\"\n",
    "}\n",
    "\n",
    "## We then import the main function\n",
    "from draw_figure_with_ma import draw_figure\n",
    "\n",
    "## We run it\n",
    "draw_figure(figure_params, show=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}