make_swangraph.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "93b9cc85-9ea8-4ad9-93ba-8f8c3359d4c2",
   "metadata": {},
   "source": [
    "# Create SwanGraph from C2C12 PacBio data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7701315a-6d40-44f7-8bce-63c9a38b58e5",
   "metadata": {},
   "source": [
    "## Load libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b2039baa-f6c2-4413-b859-722399b52446",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<frozen importlib._bootstrap>:219: RuntimeWarning: scipy._lib.messagestream.MessageStream size changed, may indicate binary incompatibility. Expected 56 from C header, got 64 from PyObject\n"
     ]
    }
   ],
   "source": [
    "import swan_vis as swan\n",
    "import pandas as pd\n",
    "import scanpy as sc\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52b3c069-13e5-47e1-87fc-4ba0c932c531",
   "metadata": {},
   "source": [
    "## Load data into Python\n",
    "We have data from C2C12 myoblasts and 3-day differentiated myotubes. Bulk long-read RNA-seq was performed using polyA primers and PacBio sequencing. To quantify transcript expression, reads were mapped to the mouse genome using Minimap2, then the reads were cleaned and pre-processed with TranscriptClean and [TALON](https://github.com/mortazavilab/TALON). \n",
    "\n",
    "In short, BAM files were downloaded from [this ENCODE data cart](https://www.encodeproject.org/carts/c2c12_bulk_pb/) and pre-processed with the walkthrough described [here](https://freese.gitbook.io/swan/tutorials/data_processing).\n",
    "\n",
    "We will use [Swan](https://freese.gitbook.io/swan/) to analyze and visualize the C2C12 long-read transcriptome."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "407be957-f819-4eae-9dbd-92cb6e63fefc",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Initialize a new SwanGraph\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "65d5091e-f8d5-44a2-80c8-7e2a617746cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "sg = swan.SwanGraph()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "73c7a449-1eb9-41c6-896f-cf9212dff0e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "annot_gtf = '../swan_data/gencode.vM21.primary_assembly.annotation_UCSC_names.gtf'\n",
    "data_gtf = '../swan_data/all_talon_observedOnly.gtf'\n",
    "ab_file = '../swan_data/all_talon_abundance_filtered.tsv'\n",
    "talon_db = '../swan_data/talon.db'\n",
    "pass_list = '../swan_data/all_pass_list.csv'\n",
    "meta = '../swan_data/swan_metadata.tsv'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eba87b3b-7334-48b6-bda3-2ada59ef36b0",
   "metadata": {},
   "source": [
    "## Adding a reference transcriptome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e1b3ff64-05ce-45bc-8541-0ec39565af2b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Adding annotation to the SwanGraph\n"
     ]
    }
   ],
   "source": [
    "# add an annotation transcriptome\n",
    "sg.add_annotation(annot_gtf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0a04798-490a-4ca2-a32e-813ab091f0fc",
   "metadata": {},
   "source": [
    "## Adding transcript models from a TALON database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "503cb785-db02-4810-a9d0-55fcee90c74e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Adding transcriptome to the SwanGraph\n"
     ]
    }
   ],
   "source": [
    "# add a dataset's transcriptome and abundance information to\n",
    "# the SwanGraph\n",
    "sg.add_transcriptome(talon_db, pass_list=pass_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21585b85-0f48-4984-8262-48ecd85e1ce9",
   "metadata": {},
   "source": [
    "## Adding datasets and their abundance\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6b25193b-8293-47a6-b4cb-f698cc1a5811",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Adding abundance for datasets ENCFF202MCY, ENCFF396UFT, ENCFF408MUF, ENCFF652XIT to SwanGraph.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/apps/python/3.8.0/lib/python3.8/site-packages/swan_vis/swangraph.py:313: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.\n",
      "  self.adata = anndata.AnnData(var=var, obs=obs, X=X)\n",
      "/opt/apps/python/3.8.0/lib/python3.8/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      "  warnings.warn(\"Transforming to str index.\", ImplicitModificationWarning)\n"
     ]
    }
   ],
   "source": [
    "# add each dataset's abundance information to the SwanGraph\n",
    "sg.add_abundance(ab_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc6fd8b6-faaa-4692-8243-8cc9d8adce19",
   "metadata": {},
   "source": [
    "## Adding metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a78e6205-7d42-4e23-aa66-efc807690970",
   "metadata": {},
   "outputs": [],
   "source": [
    "sg.add_metadata(meta)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77b6a8c7-0e7f-4240-9525-317e4a88a04a",
   "metadata": {},
   "source": [
    "## Saving and loading your SwanGraph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3b93dca0-da6b-4c58-8c52-c43aa9658e61",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving graph as ../swan_data/swan.p\n"
     ]
    }
   ],
   "source": [
    "# save the SwanGraph as a Python pickle file\n",
    "sg.save_graph('../swan_data/swan')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "55be14e9-d0c8-416e-b14e-b6b2ee122f93",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read in graph from ../swan_data/swan.p\n"
     ]
    }
   ],
   "source": [
    "# load up a saved SwanGraph from a pickle file\n",
    "sg = swan.read('../swan_data/swan.p')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "242fe4b0-7934-4a5a-a611-d88bc03b645f",
   "metadata": {},
   "source": [
    "## Differential gene expression "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ac9b0e5c-9d5b-4953-bb93-7d099a1b67f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_col = 'timepoint'\n",
    "obs_conditions = ['72hr', '0hr']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0e5c4130-e35f-4991-b700-16f2dfbf0b67",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "training location model: False\n",
      "training scale model: True\n",
      "iter   0: ll=6977754.095757\n",
      "iter   1: ll=6977754.095757, converged: 0.00% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter   2: ll=4184457.110859, converged: 75.25% (loc: 75.25%, scale update: True), in 114.07sec\n",
      "iter   3: ll=4184457.110859, converged: 75.25% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter   4: ll=3746292.546475, converged: 88.49% (loc: 88.49%, scale update: True), in 35.17sec\n",
      "iter   5: ll=3746292.546475, converged: 88.49% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter   6: ll=3739024.127695, converged: 96.88% (loc: 96.88%, scale update: True), in 19.12sec\n",
      "iter   7: ll=3739024.127695, converged: 96.88% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter   8: ll=3737746.290967, converged: 99.42% (loc: 99.42%, scale update: True), in 8.08sec\n",
      "iter   9: ll=3737746.290967, converged: 99.42% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter  10: ll=3737465.107693, converged: 99.87% (loc: 99.87%, scale update: True), in 4.22sec\n",
      "iter  11: ll=3737465.107693, converged: 99.87% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter  12: ll=3737440.844627, converged: 99.99% (loc: 99.99%, scale update: True), in 3.28sec\n",
      "iter  13: ll=3737440.844627, converged: 99.99% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter  14: ll=3737440.844627, converged: 100.00% (loc: 100.00%, scale update: True), in 3.66sec\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/apps/python/3.8.0/lib/python3.8/site-packages/dask/array/core.py:3141: RuntimeWarning: divide by zero encountered in divide\n",
      "  size = (limit / dtype.itemsize / largest_block) ** (1 / len(autos))\n"
     ]
    }
   ],
   "source": [
    "# perform a differential gene expression\n",
    "# Wald test on the provided metadata column and conditions\n",
    "test = sg.de_gene_test(obs_col, obs_conditions=obs_conditions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10ae5c84-c11d-4591-b138-1a00a447503f",
   "metadata": {},
   "source": [
    "## Differential transcript expression "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5d097b40-cbc2-450a-a578-797e6b2458be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "training location model: False\n",
      "training scale model: True\n",
      "iter   0: ll=19686469.746446\n",
      "iter   1: ll=19686469.746446, converged: 0.00% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter   2: ll=8719579.346120, converged: 77.21% (loc: 77.21%, scale update: True), in 309.03sec\n",
      "iter   3: ll=8719579.346120, converged: 77.21% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter   4: ll=6935936.328318, converged: 89.40% (loc: 89.40%, scale update: True), in 84.78sec\n",
      "iter   5: ll=6935936.328318, converged: 89.40% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter   6: ll=6901163.356511, converged: 96.14% (loc: 96.14%, scale update: True), in 48.71sec\n",
      "iter   7: ll=6901163.356511, converged: 96.14% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter   8: ll=6897957.642964, converged: 99.27% (loc: 99.27%, scale update: True), in 26.45sec\n",
      "iter   9: ll=6897957.642964, converged: 99.27% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter  10: ll=6897376.100704, converged: 99.85% (loc: 99.85%, scale update: True), in 13.16sec\n",
      "iter  11: ll=6897376.100704, converged: 99.85% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter  12: ll=6897335.704786, converged: 99.99% (loc: 99.99%, scale update: True), in 9.74sec\n",
      "iter  13: ll=6897335.704786, converged: 99.99% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter  14: ll=6897333.440464, converged: 100.00% (loc: 100.00%, scale update: True), in 7.55sec\n",
      "iter  15: ll=6897333.440464, converged: 100.00% (loc: 100.00%, scale update: False), in 0.00sec\n",
      "iter  16: ll=6897333.440464, converged: 100.00% (loc: 100.00%, scale update: True), in 4.49sec\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/apps/python/3.8.0/lib/python3.8/site-packages/dask/array/core.py:3141: RuntimeWarning: divide by zero encountered in divide\n",
      "  size = (limit / dtype.itemsize / largest_block) ** (1 / len(autos))\n"
     ]
    }
   ],
   "source": [
    "# perform a differential transcript expression\n",
    "# Wald test on the provided metadata column and conditions\n",
    "test = sg.de_transcript_test(obs_col, obs_conditions=obs_conditions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5193340-11d4-4730-ab75-045aa7b4cd7f",
   "metadata": {},
   "source": [
    "## Differential isoform expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9e115e54-f951-41a0-a1bc-69a6a50ef27d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Testing for DIE for each gene: 100%|██████████| 55868/55868 [16:55<00:00, 55.03it/s]\n"
     ]
    }
   ],
   "source": [
    "# find genes that exhibit DIE between HFFc6 and HepG2\n",
    "die_table = sg.die_gene_test(obs_col=obs_col,\n",
    "                             obs_conditions=obs_conditions,\n",
    "                             verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4d6ffc42-1a93-43a8-9bc1-65a10bc749ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>gid</th>\n",
       "      <th>p_val</th>\n",
       "      <th>dpi</th>\n",
       "      <th>adj_p_val</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ENSMUSG00000000058.6</td>\n",
       "      <td>6.010027e-01</td>\n",
       "      <td>2.666290</td>\n",
       "      <td>8.784116e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ENSMUSG00000000078.7</td>\n",
       "      <td>7.551110e-01</td>\n",
       "      <td>0.828447</td>\n",
       "      <td>9.830732e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ENSMUSG00000000088.7</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>0.459934</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ENSMUSG00000061689.15</td>\n",
       "      <td>3.539379e-19</td>\n",
       "      <td>66.746536</td>\n",
       "      <td>8.354474e-18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ENSMUSG00000000127.15</td>\n",
       "      <td>1.044654e-01</td>\n",
       "      <td>8.423532</td>\n",
       "      <td>2.900986e-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     gid         p_val        dpi     adj_p_val\n",
       "0   ENSMUSG00000000058.6  6.010027e-01   2.666290  8.784116e-01\n",
       "1   ENSMUSG00000000078.7  7.551110e-01   0.828447  9.830732e-01\n",
       "2   ENSMUSG00000000088.7  1.000000e+00   0.459934  1.000000e+00\n",
       "3  ENSMUSG00000061689.15  3.539379e-19  66.746536  8.354474e-18\n",
       "4  ENSMUSG00000000127.15  1.044654e-01   8.423532  2.900986e-01"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "die_table.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "da50ccd3-cb76-4e54-9f0b-b7e3860fe958",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving graph as ../swan_data/swan.p\n"
     ]
    }
   ],
   "source": [
    "# save the SwanGraph as a Python pickle file\n",
    "sg.save_graph('../swan_data/swan')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a36ee4f7-797f-4f8a-bdab-0b02311834e7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read in graph from ../swan_data/swan.p\n"
     ]
    }
   ],
   "source": [
    "# load up a saved SwanGraph from a pickle file\n",
    "sg = swan.read('../swan_data/swan.p')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59e90b45-c621-4084-b63c-c7c4f3a4f90a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}