celfie-simplified-demo.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CtkbWfJ9CAkI"
      },
      "source": [
        "## CelFiE Quick Start \n",
        "\n",
        "**By Nick Semenkovich \\<semenko@alum.mit.edu\\>**\n",
        "\n",
        "The original CelFiE project is by Christa Caggiano -- this simply reworks her code to be more practical and accept standard (.bed) inputs.\n",
        "\n",
        "In this example, we download and run CelFiE-simplified, to predict one sample's fractional tissue abundance -- a single mature neutrophil from Blueprint.\n",
        "\n",
        "Note: Note, CelFiE's TIM matrix was trained on Blueprint, so this is very overfit."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Aunzbkq-C3-Q"
      },
      "outputs": [],
      "source": [
        "# Install helpers\n",
        "! apt-get -qq install bedops bedtools\n",
        "! wget --no-config -q --no-clobber -O bigWigToWig https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64.v385/bigWigToWig\n",
        "! chmod +x bigWigToWig"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ToL63LyjDFLi",
        "outputId": "a52d6da5-30e2-4cf8-e100-16f83c8163d2"
      },
      "outputs": [],
      "source": [
        "# Download CelFiE\n",
        "! git clone https://github.com/semenko/celfie-simplified.git\n",
        "! git -C celfie-simplified pull\n",
        "print(\"Celfie git version: \", end = \"\")\n",
        "! git -C celfie-simplified rev-parse --short HEAD"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pVC6kkxYCK6C"
      },
      "outputs": [],
      "source": [
        "# Download S00K5EA1, a mature neutrophil\n",
        "# For details, see: http://dcc.blueprint-epigenome.eu/#/files\n",
        "\n",
        "# bs_cov = total read counts\n",
        "! wget --no-config -q --no-clobber -O S00K5EA1.bs_cov.bw http://ftp.ebi.ac.uk/pub/databases/blueprint/data/homo_sapiens/GRCh38/venous_blood/PB100713/mature_neutrophil/Bisulfite-Seq/CNAG/S00K5EA1.CPG_methylation_calls.bs_cov.GRCh38.20160531.bw\n",
        "\n",
        "# bs_call = methylation percentage (not #)\n",
        "! wget --no-config -q --no-clobber -O S00K5EA1.bs_call.bw http://ftp.ebi.ac.uk/pub/databases/blueprint/data/homo_sapiens/GRCh38/venous_blood/PB100713/mature_neutrophil/Bisulfite-Seq/CNAG/S00K5EA1.CPG_methylation_calls.bs_call.GRCh38.20160531.bw"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Y8RtK75jB-jC"
      },
      "outputs": [],
      "source": [
        "# Convert the two .bw files above into .beds\n",
        "! ./bigWigToWig S00K5EA1.bs_cov.bw stdout | wig2bed > S00K5EA1.bs_cov.bed\n",
        "! ./bigWigToWig S00K5EA1.bs_call.bw stdout | wig2bed > S00K5EA1.bs_call.bed\n",
        "\n",
        "# Combine into one .bed, containing:\n",
        "# chrom start stop %meth #depth\n",
        "! paste S00K5EA1.bs_call.bed S00K5EA1.bs_cov.bed | cut -f 1-3,5,10 > S00K5EA1.percent_meth.bed\n",
        "\n",
        "# Convert % methylation into counts, our required input format:\n",
        "# chrom start stop #meth #depth\n",
        "! awk 'BEGIN{OFS=\"\\t\"}{ print $1, $2, $3, int($4 * $5 + 0.5), int($5) }' S00K5EA1.percent_meth.bed > S00K5EA1.counts_meth.bed"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kmW00ONiMUQp",
        "outputId": "e65444a1-d7f6-4682-b1c0-8d5aefcbbb61"
      },
      "outputs": [],
      "source": [
        "# Let's make sure that looks valid\n",
        "! head -n 2 S00K5EA1.counts_meth.bed\n",
        "! wc -l S00K5EA1.counts_meth.bed\n",
        "\n",
        "# Optionally take a small sample for the git repo\n",
        "# ! bedtools sample -n 100000 -seed 42 -i S00K5EA1.counts_meth.bed | sort -k 1,1 -k2,2n > sample-neutrophil.bed"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "YfB1OJn-PIeb"
      },
      "outputs": [],
      "source": [
        "## Run CelFiE's sum_by_list -- basically intesects a .bed and sums two columns\n",
        "\n",
        "# The TIM matrix can't have a header for this\n",
        "! tail -n+2 celfie/tim_matrix.txt > celfie/tim_matrix.noheader.txt\n",
        "\n",
        "# For speed, pre-filter the .bed to regions intersecting the TIM matrix\n",
        "# -u prints only one entry of SAMPLE for each overlap with TIM_MATRIX\n",
        "# ! bedtools intersect -u -a S00K5EA1.counts_meth.bed -b celfie/tim_matrix.noheader.txt > S00K5EA1.counts_meth.filtered.bed\n",
        "\n",
        "! python celfie/TIMs/sum_by_list.py celfie/tim_matrix.noheader.txt S00K5EA1.counts_meth.filtered.bed S00K5EA1.summed.bed 1"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ItmHviVeP7ob",
        "outputId": "bb13e058-a08d-451f-fcca-664fbe90c981"
      },
      "outputs": [],
      "source": [
        "# Let's make sure that looks valid, too -- there should only be 1581 lines (one per TIM)\n",
        "! head -n 2 S00K5EA1.summed.bed\n",
        "! wc -l S00K5EA1.summed.bed"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fSti6ZmvMqRD",
        "outputId": "7e171941-d19f-47b8-a9d2-831d32768685"
      },
      "outputs": [],
      "source": [
        "# Now we need to stitch on a header\n",
        "! echo -e \"chrom\\tstart\\tend\\tsample_meth\\tsample_depth\" > S00K5EA1.summed.with_header.txt\n",
        "# We need to remove a ^M that sum_by_list.py leaves\n",
        "! cat S00K5EA1.summed.bed | tr -d $'\\r' >> S00K5EA1.summed.with_header.txt\n",
        "\n",
        "# Combine with TIM matrix\n",
        "! paste S00K5EA1.summed.with_header.txt celfie/tim_matrix.txt > S00K5EA1_celfie_input.txt\n",
        "\n",
        "# Finally, does that look OK?\n",
        "! head -n2 S00K5EA1_celfie_input.txt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "G7QpzaQtMYbP",
        "outputId": "6bedabaa-fa95-4d8a-9ddb-4528989f0177"
      },
      "outputs": [],
      "source": [
        "## Time to run CelFiE!\n",
        "! python celfie-simplified/celfie.py -u 0 S00K5EA1_celfie_input.txt S00K5EA1_output 1"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 363
        },
        "id": "MpYiaQC7OLOh",
        "outputId": "0c9f95ae-f594-4b0e-a634-b3e329346473"
      },
      "outputs": [],
      "source": [
        "# Let's plot it, using .ipynb code from CelFiE's demo\n",
        "import pandas as pd\n",
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "tissue_proportions = pd.read_csv(\"S00K5EA1_output/1_tissue_proportions.txt\", delimiter=\"\\t\")\n",
        "# Rename column 1\n",
        "tissue_proportions.rename(columns={\"Unnamed: 0\": \"samples\"}, inplace=True)\n",
        "# Melt to one entry per-line for sns plots\n",
        "tissue_proportions = tissue_proportions.melt(\"samples\",  var_name=\"tissue\", value_name=\"estimate\")\n",
        "sns.boxplot(x=\"tissue\", y=\"estimate\", data=tissue_proportions, palette=[\"#61c2d3\", \"#003f4b\"])\n",
        "plt.title('S00K5EA1 Neutrophil')\n",
        "plt.xticks(rotation=90)\n",
        "plt.ylabel(\"CelFiE estimate\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "OoHRZAwHTNjj",
        "outputId": "b763c80e-db67-481e-88a4-1aedf0043ffe"
      },
      "outputs": [],
      "source": [
        "! cat S00K5EA1_output/1_tissue_proportions.txt"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "quick-celfie-repro.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}