generic_metagenomic_demo.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/dcolinmorgan/grph/blob/main/generic_metagenomic_demo.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ozrca88hza85"
      },
      "source": [
        "# Tutorial: not-so-quick metagenomic analysis using GPU UMAP analysis & visualization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tLRHg2VEzoYy"
      },
      "source": [
        "UMAP analysis & visualization, metagenomic samples can be compared and explored.\n",
        "\n",
        "*   Task: Analyze metagenomic samples for similarity\n",
        "*   Data: time series samples\n",
        "**   563 samples collected from 84 donors, producing 4 dense long-term time series (up to 1 sample every other day during 18 months)\n",
        "*   [data](https://www.ebi.ac.uk/ena/browser/view/PRJNA544527)\n",
        "*   [metadata](https://static-content.springer.com/esm/art%3A10.1038%2Fs41591-019-0559-3/MediaObjects/41591_2019_559_MOESM3_ESM.xlsx)\n",
        "*   [paper](https://www-nature-com.eproxy.lib.hku.hk/articles/s41591-019-0559-3)\n",
        "\n",
        "**Insight/ Result:**\n",
        "\n",
        "2342s to umap and dbscan, which [the analysis reference paper](https://journals.asm.org/doi/full/10.1128/msystems.00118-23) does 12x per dataset\n",
        "\n",
        "(See also: [GPU acceleration via graphistry](https://github.com/dcolinmorgan/grph/blob/main/accelearting_metagenomic_demo.ipynb))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PjnS_PCWaClg"
      },
      "source": [
        "# Setup"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "JTzHVRgEcDUV"
      },
      "outputs": [],
      "source": [
        "!pip install graphistry[ai]\n",
        "!pip install dirty_cat\n",
        "!pip install umap-learn\n",
        "import umap\n",
        "import dirty_cat\n",
        "!pip install Biopython"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HK_8_7UB0mhx"
      },
      "source": [
        "# import /configure\n",
        "\n",
        "get a free api-key at https://www.graphistry.com/\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "id": "i42QrQ_ejC4h",
        "outputId": "e8ec41d6-1636-4d14-a2b7-9de2c295d67e"
      },
      "outputs": [
        {
          "data": {
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            },
            "text/plain": [
              "'0.29.2'"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "import pandas as pd\n",
        "import graphistry\n",
        "from time import time\n",
        "\n",
        "\n",
        "graphistry.register(api=3,protocol=\"https\", server=\"hub.graphistry.com\", username='dcolinmorgan', password='fXjJnkE3Gik6BWy') ## key id, secret key\n",
        "\n",
        "# graphistry.register(api=3,protocol=\"https\", server=\"hub.graphistry.com\", username='dcolinmorgan', password='***') ## key id, secret key\n",
        "graphistry.__version__"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_Y1VlCGy4_FR"
      },
      "source": [
        "# bio-ml dataset\n",
        "\n",
        "\n",
        "1.   [3 subjects x 10 time points](\n",
        "https://www.ebi.ac.uk/ena/browser/view/PRJNA544527)\n",
        "\n",
        "2.  [metadata](\n",
        "https://static-content.springer.com/esm/art%3A10.1038%2Fs41591-019-0559-3/MediaObjects/41591_2019_559_MOESM3_ESM.xlsx)\n",
        "\n",
        "3.   !wget https://raw.githubusercontent.com/dcolinmorgan/grph/main/ftp_PRJNA544527.txt\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CBN6Z_77Sduq"
      },
      "outputs": [],
      "source": [
        "!wget -nc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR922/006/SRR9224006/SRR9224006_1.fastq.gz\n",
        "!wget -nc ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR922/006/SRR9224006/SRR9224006_2.fastq.gz"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UjbF1ZfwULLt"
      },
      "outputs": [],
      "source": [
        "! gunzip SRR9224006_1.fastq.gz\n",
        "! gunzip SRR9224006_2.fastq.gz"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ljJnoNeLSnWq",
        "outputId": "64e52f66-1939-417f-ef5b-2ebc2306a133"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "@SRR9224006.1 7001174F:HVTFNBCXX161011:HVTFNBCXX:2:2206:18894:58151/1\n",
            "AAAAAAAACAAAATAATGGAAACAAAAAACATCTACTTCATCAGCGGCATTGATACAGATGCCGGAAAAAGCTATTGCACCGCCTGGTATGCCCGTGAGCT\n",
            "+\n",
            "DDDDDIIIIIIGIIIIIIIIIIIIHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIDHGIIIIGGHIIHHHIIIIIIIHHIIIIIIHIIIIIII\n",
            "@SRR9224006.2 7001174F:HVTFNBCXX161011:HVTFNBCXX:1:1108:17623:49640/1\n",
            "AAAAAAAACAACCCAATGCGATTCTGATCGCAATCTACATAAGTTACTACTGGTTATCTTCCCTGAAG\n",
            "+\n",
            "DDBDDIGIIHIIIGHGHFHIGHIIIEHIHIIIIIIIIIIIIEHHHGHIIIIIIIIIIIIIIEHHHHHH\n",
            "@SRR9224006.3 7001174F:HVTFNBCXX161011:HVTFNBCXX:1:1204:5345:82516/1\n",
            "AAAAAAAACAAGAGCTTTATTAAACACGTCTTGATCTTTTTTACACCTGCCGGAAATTCCATCGT\n"
          ]
        }
      ],
      "source": [
        "!head /content/SRR9224006_1.fastq"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UWDGFCWxSpv_"
      },
      "outputs": [],
      "source": [
        "from Bio import SeqIO\n",
        "import glob,os\n",
        "import pandas as pd\n",
        "B=pd.DataFrame()\n",
        "for i in glob.glob('/content/*.fastq'):\n",
        "    # j=os.path.basename(i)\n",
        "    fasta_sequences = SeqIO.parse(open(i),'fastq')\n",
        "    identifiers = []\n",
        "    sequences = []\n",
        "    for fasta in fasta_sequences:\n",
        "        name, sequence = fasta.id, str(fasta.seq)\n",
        "        identifiers.append(name)\n",
        "        sequences.append(sequence)\n",
        "\n",
        "    A=pd.DataFrame([identifiers,sequences]).T\n",
        "    A.columns=['ID','seq']\n",
        "    A.dropna(inplace=True)\n",
        "    B=B.append(A)\n",
        "    # A['ID']#=A.ID.str.split('-')[0:1]\n",
        "# B['ID']=B['ID'].str.split('-').str[0]+'_'+B['ID'].str.split('-').str[1]#.cat()\n",
        "B['ID']=B.ID.str.split('_length').str[0]\n",
        "B.index=B.ID"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fDvymRFjHPNB"
      },
      "source": [
        "# install [HUMAnN 3](https://huttenhower.sph.harvard.edu/humann), a method for efficiently and accurately profiling the abundance of microbial metabolic pathways and other molecular functions from metagenomic or metatranscriptomic sequencing data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "yIeeDXPBHN6D"
      },
      "outputs": [],
      "source": [
        "# !pip install humann --no-binary :all:\n",
        "!pip install metaphlan"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TnYmqPJ-bHPa"
      },
      "outputs": [],
      "source": [
        "### !humann_databases --download utility_mapping full /path/to/databases --update-config yes\n",
        "\n",
        "# !humann_test\n",
        "\n",
        "# !wget https://github.com/biobakery/humann/raw/master/examples/demo.fastq.gz\n",
        "# !humann -i demo.fastq.gz -o sample_results"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gtOG0QeoUoX6"
      },
      "source": [
        "### takes very long for running all samples\n",
        " (1day+ run on cluster)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XAtq7G_waxLd"
      },
      "outputs": [],
      "source": [
        "mkdir assemble epi_sam_out mpa4_out\n",
        "# !humann -i /content/All_MAGs/Sample_101_S75_bin_1.fa -o test_out\n",
        "%%bash\n",
        "seq=$(ls /content/*.fastq | cut -d / -f2| cut -d _ -f1)\n",
        "\n",
        "for i in $(eval \"echo \"$seq\" | cut -d _ -f1\")\n",
        "\n",
        "do\n",
        "metaphlan /content/${i}.fa --nproc 40 --input_type fasta -o /content/assemble/${i}/h4_out.txt -t rel_ab_w_read_stats\n",
        "done"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dHbwBIEn6Wxr"
      },
      "source": [
        "# umap and dbscan\n",
        "\n",
        "idea for metagenomic analysis based on [Quantifying Shared and Unique Gene Content across 17 Microbial Ecosystems\n",
        "](https://journals.asm.org/doi/full/10.1128/msystems.00118-23)\n",
        "\n",
        "(analyze all samples run on cluster)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RoIBLY3-670T"
      },
      "outputs": [],
      "source": [
        "# from sqlalchemy.util.compat import dataclass_fields\n",
        "!wget https://github.com/dcolinmorgan/grph/raw/main/PRJNA544527_mpa4out.txt\n",
        "\n",
        "data=pd.read_csv('/content/PRJNA544527_mpa4out.txt',sep='\\t',skiprows=1,index_col=0)\n",
        "data.index=data.reset_index().clade_name.str.split('|',expand=True)[6]\n",
        "data=data.reset_index().dropna(axis=0)\n",
        "data.index=data[6]\n",
        "data=data.drop(columns=6)\n",
        "\n",
        "!wget https://raw.githubusercontent.com/dcolinmorgan/grph/main/PRJNA544527-meta_inf.txt\n",
        "meta=pd.read_csv('/content/PRJNA544527-meta_inf.txt',sep='\\t',header=None)\n",
        "\n",
        "mm=pd.merge(data.T,meta[[3,5]],left_index=True,right_on=3)\n",
        "\n",
        "mm['id']=mm[5].str.split('-').str[0]\n",
        "mm['time']=mm[5].str.split('_').str[0].str.split('-').str[1]\n",
        "\n",
        "!wget https://static-content.springer.com/esm/art%3A10.1038%2Fs41591-019-0559-3/MediaObjects/41591_2019_559_MOESM3_ESM.xlsx\n",
        "metaa=pd.read_excel('/content/41591_2019_559_MOESM3_ESM.xlsx',sheet_name='SupTable2',skiprows=3)\n",
        "metaa=metaa[['Donor','Age','Sex','BMI']]\n",
        "\n",
        "Full_table=pd.merge(mm,metaa,left_on='id',right_on='Donor')\n",
        "Full_table=Full_table.drop(columns=[3,\t5,\t'id'])\n",
        "Full_table.time=pd.to_datetime(Full_table.time,unit='d')\n",
        "Full_table.time=Full_table.time.values.astype('datetime64[M]')\n",
        "\n",
        "data2=Full_table.melt(id_vars=['time','Donor','Age','Sex','BMI'])\n",
        "data2.to_csv('PRJNA544527_mpa4_annot_table.txt',sep='\\t')\n",
        "\n",
        "# final df stored here also\n",
        "# !wget https://raw.githubusercontent.com/dcolinmorgan/grph/main/PRJNA544527_mpa4_annot_table.txt\n",
        "# data2=pd.read_csv('PRJNA544527_mpa4_annot_table.txt',sep='\\t',index_col=0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "background_save": true
        },
        "id": "vBxk7SSfJAd8"
      },
      "outputs": [],
      "source": [
        "g = graphistry.nodes((data2[['variable','value']]))\n",
        "\n",
        "t=time()\n",
        "g2=g.featurize(memoize=False)\n",
        "g3=g2.umap(dbscan=True)\n",
        "print(\"\\n\"+str(time()-t))\n",
        "\n",
        "emb2=g3._node_embedding\n",
        "emb2['time']=data2.time\n",
        "emb2['_n']=data2.variable\n",
        "g22=graphistry.nodes(emb2.reset_index(),'index').edges(g3._edges,'_src_implicit','_dst_implicit').bind(point_x=\"x\",point_y=\"y\").settings(url_params={\"play\":0})\n",
        "g22.plot()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GCOYg71RDGPT"
      },
      "source": [
        "## protracted UMAP above, full below"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "background_save": true
        },
        "id": "6Zoo17Ui9zxh"
      },
      "outputs": [],
      "source": [
        "g = graphistry.nodes(data2.drop(columns='time'))\n",
        "\n",
        "t=time()\n",
        "g2=g.umap(dbscan=True)\n",
        "print(\"\\n\"+str(time()-t))\n",
        "\n",
        "g2.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "F_FNYu7hcKrR"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "_Y1VlCGy4_FR",
        "fDvymRFjHPNB"
      ],
      "machine_shape": "hm",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}