examples/H1ESC_Nanog_gkmsvm/Nanog_GkmExplain_Generate_Data.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Nanog GkmExplain Generate Data.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/kundajelab/tfmodisco/blob/master/examples/H1ESC_Nanog_gkmsvm/Nanog_GkmExplain_Generate_Data.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "clFJEqUmWACq"
      },
      "source": [
        "### Download all the requisite data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NkDxJNiFWJUT",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "3fe88097-3fd4-470d-c527-2a112a1a1198"
      },
      "source": [
        "!apt-get install bedtools\n",
        "!git clone https://github.com/kundajelab/lsgkm.git lsgkm\n",
        "%cd lsgkm/src\n",
        "!make\n",
        "%cd ../..\n",
        "\n",
        "#Download ENCODE-processed peak files to get the foreground and background\n",
        "! [[ -e conservative_peaks.bed.gz ]] || wget https://www.encodeproject.org/files/ENCFF148PBJ/@@download/ENCFF148PBJ.bed.gz -O conservative_peaks.bed.gz\n",
        "! [[ -e optimal_peaks.bed.gz ]] || wget https://www.encodeproject.org/files/ENCFF379EPK/@@download/ENCFF379EPK.bed.gz -O optimal_peaks.bed.gz\n",
        "# This DNAse dataset was obtained from ENCODE (accessions: ENCSR000EMU, ENCSR000EMU_ENCSR794OFW)\n",
        "# and reprocessed using the Kundaje lab's ATAC/DNAse processing pipeline (https://github.com/kundajelab/atac_dnase_pipelines)\n",
        "# by Daniel Kim.\n",
        "! [[ -e bg_dnase.bed.gz ]] || wget https://raw.githubusercontent.com/AvantiShri/model_storage/8947701/gkmexplain/ENCSR000EMU_ENCSR794OFW.H1_Cells.UW_Stam.DNase-seq_rep1-pr.IDR0.1.narrowPeak.gz -O bg_dnase.bed.gz\n",
        "  \n",
        "#Get hg19 fasta by download 2bit and then converting to fa\n",
        "! [[ -f hg19.2bit ]] || wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.2bit -O hg19.2bit  \n",
        "! [[ -f twoBitToFa ]] || wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa -O twoBitToFa\n",
        "!chmod a+x twoBitToFa\n",
        "! [[ -f hg19.genome.fa ]] || ./twoBitToFa hg19.2bit hg19.genome.fa"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Reading package lists... Done\n",
            "Building dependency tree       \n",
            "Reading state information... Done\n",
            "The following NEW packages will be installed:\n",
            "  bedtools\n",
            "0 upgraded, 1 newly installed, 0 to remove and 13 not upgraded.\n",
            "Need to get 577 kB of archives.\n",
            "After this operation, 2,040 kB of additional disk space will be used.\n",
            "Get:1 http://archive.ubuntu.com/ubuntu bionic/universe amd64 bedtools amd64 2.26.0+dfsg-5 [577 kB]\n",
            "Fetched 577 kB in 1s (488 kB/s)\n",
            "Selecting previously unselected package bedtools.\n",
            "(Reading database ... 149406 files and directories currently installed.)\n",
            "Preparing to unpack .../bedtools_2.26.0+dfsg-5_amd64.deb ...\n",
            "Unpacking bedtools (2.26.0+dfsg-5) ...\n",
            "Setting up bedtools (2.26.0+dfsg-5) ...\n",
            "Cloning into 'lsgkm'...\n",
            "remote: Enumerating objects: 59, done.\u001b[K\n",
            "remote: Counting objects: 100% (59/59), done.\u001b[K\n",
            "remote: Compressing objects: 100% (44/44), done.\u001b[K\n",
            "remote: Total 382 (delta 34), reused 32 (delta 15), pack-reused 323\u001b[K\n",
            "Receiving objects: 100% (382/382), 585.86 KiB | 2.33 MiB/s, done.\n",
            "Resolving deltas: 100% (247/247), done.\n",
            "/content/lsgkm/src\n",
            "g++ -Wall -Wconversion -O3 -fPIC -c libsvm.cpp\n",
            "g++ -Wall -Wconversion -O3 -fPIC -c libsvm_gkm.c\n",
            "g++ -Wall -Wconversion -O3 -fPIC gkmtrain.c libsvm.o libsvm_gkm.o -o gkmtrain -lm -lpthread\n",
            "g++ -Wall -Wconversion -O3 -fPIC gkmpredict.c libsvm.o libsvm_gkm.o -o gkmpredict -lm -lpthread\n",
            "g++ -Wall -Wconversion -O3 -fPIC gkmexplain.c libsvm.o libsvm_gkm.o -o gkmexplain -lm -lpthread\n",
            "/content\n",
            "--2021-02-27 10:54:20--  https://www.encodeproject.org/files/ENCFF148PBJ/@@download/ENCFF148PBJ.bed.gz\n",
            "Resolving www.encodeproject.org (www.encodeproject.org)... 34.211.244.144\n",
            "Connecting to www.encodeproject.org (www.encodeproject.org)|34.211.244.144|:443... connected.\n",
            "HTTP request sent, awaiting response... 307 Temporary Redirect\n",
            "Location: https://encode-public.s3.amazonaws.com/2017/12/08/aa3b4cf0-2a7a-403e-9a02-a12130198638/ENCFF148PBJ.bed.gz?response-content-disposition=attachment%3B%20filename%3DENCFF148PBJ.bed.gz&AWSAccessKeyId=ASIATGZNGCNXVPBOCGPS&Signature=Gj7GooQPReRWoK1MJyueGe6sZ3k%3D&x-amz-security-token=IQoJb3JpZ2luX2VjEEIaCXVzLXdlc3QtMiJHMEUCIGBA7aWLyN1o7hwXi9LsVcAOQ9bxJPx9uucjcNYCfFdrAiEArAnse8FDrdc8y3%2F7zEOe57FOGar9tEMEbJURozgFLB8qtAMIWxAAGgwyMjA3NDg3MTQ4NjMiDJgf%2BFpO5OshvOO5NyqRA232H7TmeYodqzstAHKxwkDrGJqSgqGbHNyozOYuAxjiuS%2BpOn8fZaj3X%2BTLOhBzWMWvIDmlJcQUTfFCoqZovhnu%2FVpwAvpGOLoa5%2Be4m2oT0E3lwFqOqxyAfTU8CvzkBTmkjsBPIIEdqNpnclDRVsYLko1%2FOuNbmRuZ3sToc1ieTfS5WAm7QALthc2j4KNK0sGvy%2BnK2Fv0gVxH6aUPtYRC1jiSZ0Sb6A0MhigaJos8ceg6Guq6dTEHF9Z9zcwmYUW227hQweqPCnarHo3l%2F5j6jW3TWRV3YMoufq3gHdHjjdrNn4KiHnjWu8LLv6q8fhdzTraKJ4XRjE0xhuAjHqC%2BgYW91bb%2FwkHCZaWGQgrRijT%2BPy0804z%2BfeSW4B%2Ft8fxxh3SldIfq7ImvCScPK64xnbUMuNW1aOq50vr7As92vVNjUVcGfI2g2ipRywFJkxiWGshfEfTZ18wh21L9MCX58HdMQkBsdOlwvl7IxBTqHPt2ZCztvPPHW%2FSe1AOfF7iU4rin9xYCzwnkSASMRhQCMPKo6IEGOusBM4odUjqylUYoEt80ImQ9%2BbOlQGQvFV4kb2bNrO2C6eFP2NAMh2%2FrYWP19zYNk%2BwXkfa5%2BqKb6iseMZQv3n6EaSlm5fenaJgpCOiegml%2F9LiHBCvUsVzZkdv3nA6WJgymXVy29RM%2B5YUwmEtAJhe%2FoME1%2B6KkJQgVKI4M8sz2VMr5QjEL79EI8Dy%2BLRHO5an6jRpC5jijrnom3l8arXHxq3OafyjBRz7FRAPIRWKVm7FJvuQ8ACvglg5RYKI%2Bt%2FT2x8LrAyF9jI%2BgFcYQbTQzjTyJQb2v0f3OTHvHX4QcL0NIdz0dJEbw%2Bu81Rw%3D%3D&Expires=1614552860 [following]\n",
            "--2021-02-27 10:54:20--  https://encode-public.s3.amazonaws.com/2017/12/08/aa3b4cf0-2a7a-403e-9a02-a12130198638/ENCFF148PBJ.bed.gz?response-content-disposition=attachment%3B%20filename%3DENCFF148PBJ.bed.gz&AWSAccessKeyId=ASIATGZNGCNXVPBOCGPS&Signature=Gj7GooQPReRWoK1MJyueGe6sZ3k%3D&x-amz-security-token=IQoJb3JpZ2luX2VjEEIaCXVzLXdlc3QtMiJHMEUCIGBA7aWLyN1o7hwXi9LsVcAOQ9bxJPx9uucjcNYCfFdrAiEArAnse8FDrdc8y3%2F7zEOe57FOGar9tEMEbJURozgFLB8qtAMIWxAAGgwyMjA3NDg3MTQ4NjMiDJgf%2BFpO5OshvOO5NyqRA232H7TmeYodqzstAHKxwkDrGJqSgqGbHNyozOYuAxjiuS%2BpOn8fZaj3X%2BTLOhBzWMWvIDmlJcQUTfFCoqZovhnu%2FVpwAvpGOLoa5%2Be4m2oT0E3lwFqOqxyAfTU8CvzkBTmkjsBPIIEdqNpnclDRVsYLko1%2FOuNbmRuZ3sToc1ieTfS5WAm7QALthc2j4KNK0sGvy%2BnK2Fv0gVxH6aUPtYRC1jiSZ0Sb6A0MhigaJos8ceg6Guq6dTEHF9Z9zcwmYUW227hQweqPCnarHo3l%2F5j6jW3TWRV3YMoufq3gHdHjjdrNn4KiHnjWu8LLv6q8fhdzTraKJ4XRjE0xhuAjHqC%2BgYW91bb%2FwkHCZaWGQgrRijT%2BPy0804z%2BfeSW4B%2Ft8fxxh3SldIfq7ImvCScPK64xnbUMuNW1aOq50vr7As92vVNjUVcGfI2g2ipRywFJkxiWGshfEfTZ18wh21L9MCX58HdMQkBsdOlwvl7IxBTqHPt2ZCztvPPHW%2FSe1AOfF7iU4rin9xYCzwnkSASMRhQCMPKo6IEGOusBM4odUjqylUYoEt80ImQ9%2BbOlQGQvFV4kb2bNrO2C6eFP2NAMh2%2FrYWP19zYNk%2BwXkfa5%2BqKb6iseMZQv3n6EaSlm5fenaJgpCOiegml%2F9LiHBCvUsVzZkdv3nA6WJgymXVy29RM%2B5YUwmEtAJhe%2FoME1%2B6KkJQgVKI4M8sz2VMr5QjEL79EI8Dy%2BLRHO5an6jRpC5jijrnom3l8arXHxq3OafyjBRz7FRAPIRWKVm7FJvuQ8ACvglg5RYKI%2Bt%2FT2x8LrAyF9jI%2BgFcYQbTQzjTyJQb2v0f3OTHvHX4QcL0NIdz0dJEbw%2Bu81Rw%3D%3D&Expires=1614552860\n",
            "Resolving encode-public.s3.amazonaws.com (encode-public.s3.amazonaws.com)... 52.218.137.187\n",
            "Connecting to encode-public.s3.amazonaws.com (encode-public.s3.amazonaws.com)|52.218.137.187|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 99222 (97K) [binary/octet-stream]\n",
            "Saving to: ‘conservative_peaks.bed.gz’\n",
            "\n",
            "conservative_peaks. 100%[===================>]  96.90K  --.-KB/s    in 0.1s    \n",
            "\n",
            "2021-02-27 10:54:21 (884 KB/s) - ‘conservative_peaks.bed.gz’ saved [99222/99222]\n",
            "\n",
            "--2021-02-27 10:54:21--  https://www.encodeproject.org/files/ENCFF379EPK/@@download/ENCFF379EPK.bed.gz\n",
            "Resolving www.encodeproject.org (www.encodeproject.org)... 34.211.244.144\n",
            "Connecting to www.encodeproject.org (www.encodeproject.org)|34.211.244.144|:443... connected.\n",
            "HTTP request sent, awaiting response... 307 Temporary Redirect\n",
            "Location: https://encode-public.s3.amazonaws.com/2017/12/08/b2d61850-1fa9-483a-a27e-cf9da64aed58/ENCFF379EPK.bed.gz?response-content-disposition=attachment%3B%20filename%3DENCFF379EPK.bed.gz&AWSAccessKeyId=ASIATGZNGCNX6NQC4OXW&Signature=ypb9TjX9gqrUs3LVUcY2civBjxk%3D&x-amz-security-token=IQoJb3JpZ2luX2VjEEMaCXVzLXdlc3QtMiJHMEUCIGyyclDJbElI89Sq9MSi9EZ7RsEY8Z0mJ4gbw2EGtOx9AiEAku64f6gHwX2FDzCOKDhJbULCg%2BkwnyZh6NwdHYCCE38qtAMIXBAAGgwyMjA3NDg3MTQ4NjMiDNfv572k8pZyCqDmKCqRA2JWvzQTJOn3F36hlPe8CTYfKCMNbCGcDMC5BNQDy3ImjhHetQvN87SPYawIsz0wySLSyV4gXHIdIUZ62uJCupza0PhS2oqlHctM7cuh50UIFgqNaqzEoo%2BiokN3MKN7I5h7vGdutcgsykxZePpHqvEj98GLhz1onzhV0OVv4ANJyqauVaU5SEoVYkIwM%2BlwtH1HhS8SomUoygOVUIYGUT3RTSOsMwKklzPD%2B0E1lx%2BN2cR9APXOyjc0qU7Dxz0wIiu8WX0LR9Z%2FGY8Nj%2FuZ8e27FTf4ZTXpLDmJE9h7juH1PANz5H0KZ1ORbycASBGd8NVn7NhOUHbMT6990Gor9fm7JjDCjMrJ0g9TqGDlUDlVCWGB%2FmAjkpY3lqRNeBqzM8ilA4IgdnbT0M4XYKu6M4Z%2B6AYyACNsYGVhXhNtWNlJ%2FiNaR3v4EfFz5XuIzKIFm3L2iPTalLiKXpjvQEl71RIh3C3Y4Sn3plIdmfZwb2mQ31hjYMQ5oRf80805uem1jymk9g%2BtFqdD2kDtTiEjQ2hPMJnE6IEGOusBETjCkybcKy09ZCfjs%2FOEZ%2FBR66TLXCMQ5ZlZ8v9zjDU9eRxKKSA8isGFX%2FtogjWtwhVbDDH%2Bl%2Bx45nFZh311KGzt8Rps4OEeWJvhAvWu6KmlBdq%2FxOPoOoQr9m8vJ6AefXnCBzTI%2F7Z02p%2Fe0eogjY%2BU53zSNLFPS%2B9pnq7zvkGOfQmFYXFgB%2FOY50jKtIQJcp9poPsIgnJ88wqi0AvnKHXfVeLQWFKg%2BnebVa1lx8uvTixZ8pvtMvxUUFvOiko3hrg%2Ft%2FMs%2BXfjtC1vE2tu0In4vMRVC7kfRZ2QcTtW3u%2FKYmuzMEl4twzkTg%3D%3D&Expires=1614552861 [following]\n",
            "--2021-02-27 10:54:21--  https://encode-public.s3.amazonaws.com/2017/12/08/b2d61850-1fa9-483a-a27e-cf9da64aed58/ENCFF379EPK.bed.gz?response-content-disposition=attachment%3B%20filename%3DENCFF379EPK.bed.gz&AWSAccessKeyId=ASIATGZNGCNX6NQC4OXW&Signature=ypb9TjX9gqrUs3LVUcY2civBjxk%3D&x-amz-security-token=IQoJb3JpZ2luX2VjEEMaCXVzLXdlc3QtMiJHMEUCIGyyclDJbElI89Sq9MSi9EZ7RsEY8Z0mJ4gbw2EGtOx9AiEAku64f6gHwX2FDzCOKDhJbULCg%2BkwnyZh6NwdHYCCE38qtAMIXBAAGgwyMjA3NDg3MTQ4NjMiDNfv572k8pZyCqDmKCqRA2JWvzQTJOn3F36hlPe8CTYfKCMNbCGcDMC5BNQDy3ImjhHetQvN87SPYawIsz0wySLSyV4gXHIdIUZ62uJCupza0PhS2oqlHctM7cuh50UIFgqNaqzEoo%2BiokN3MKN7I5h7vGdutcgsykxZePpHqvEj98GLhz1onzhV0OVv4ANJyqauVaU5SEoVYkIwM%2BlwtH1HhS8SomUoygOVUIYGUT3RTSOsMwKklzPD%2B0E1lx%2BN2cR9APXOyjc0qU7Dxz0wIiu8WX0LR9Z%2FGY8Nj%2FuZ8e27FTf4ZTXpLDmJE9h7juH1PANz5H0KZ1ORbycASBGd8NVn7NhOUHbMT6990Gor9fm7JjDCjMrJ0g9TqGDlUDlVCWGB%2FmAjkpY3lqRNeBqzM8ilA4IgdnbT0M4XYKu6M4Z%2B6AYyACNsYGVhXhNtWNlJ%2FiNaR3v4EfFz5XuIzKIFm3L2iPTalLiKXpjvQEl71RIh3C3Y4Sn3plIdmfZwb2mQ31hjYMQ5oRf80805uem1jymk9g%2BtFqdD2kDtTiEjQ2hPMJnE6IEGOusBETjCkybcKy09ZCfjs%2FOEZ%2FBR66TLXCMQ5ZlZ8v9zjDU9eRxKKSA8isGFX%2FtogjWtwhVbDDH%2Bl%2Bx45nFZh311KGzt8Rps4OEeWJvhAvWu6KmlBdq%2FxOPoOoQr9m8vJ6AefXnCBzTI%2F7Z02p%2Fe0eogjY%2BU53zSNLFPS%2B9pnq7zvkGOfQmFYXFgB%2FOY50jKtIQJcp9poPsIgnJ88wqi0AvnKHXfVeLQWFKg%2BnebVa1lx8uvTixZ8pvtMvxUUFvOiko3hrg%2Ft%2FMs%2BXfjtC1vE2tu0In4vMRVC7kfRZ2QcTtW3u%2FKYmuzMEl4twzkTg%3D%3D&Expires=1614552861\n",
            "Resolving encode-public.s3.amazonaws.com (encode-public.s3.amazonaws.com)... 52.218.193.154\n",
            "Connecting to encode-public.s3.amazonaws.com (encode-public.s3.amazonaws.com)|52.218.193.154|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 146369 (143K) [binary/octet-stream]\n",
            "Saving to: ‘optimal_peaks.bed.gz’\n",
            "\n",
            "optimal_peaks.bed.g 100%[===================>] 142.94K  --.-KB/s    in 0.1s    \n",
            "\n",
            "2021-02-27 10:54:21 (1.23 MB/s) - ‘optimal_peaks.bed.gz’ saved [146369/146369]\n",
            "\n",
            "--2021-02-27 10:54:21--  https://raw.githubusercontent.com/AvantiShri/model_storage/8947701/gkmexplain/ENCSR000EMU_ENCSR794OFW.H1_Cells.UW_Stam.DNase-seq_rep1-pr.IDR0.1.narrowPeak.gz\n",
            "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...\n",
            "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 2519952 (2.4M) [application/octet-stream]\n",
            "Saving to: ‘bg_dnase.bed.gz’\n",
            "\n",
            "bg_dnase.bed.gz     100%[===================>]   2.40M  --.-KB/s    in 0.07s   \n",
            "\n",
            "2021-02-27 10:54:22 (33.2 MB/s) - ‘bg_dnase.bed.gz’ saved [2519952/2519952]\n",
            "\n",
            "--2021-02-27 10:54:22--  http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.2bit\n",
            "Resolving hgdownload.cse.ucsc.edu (hgdownload.cse.ucsc.edu)... 128.114.119.163\n",
            "Connecting to hgdownload.cse.ucsc.edu (hgdownload.cse.ucsc.edu)|128.114.119.163|:80... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 816241703 (778M)\n",
            "Saving to: ‘hg19.2bit’\n",
            "\n",
            "hg19.2bit           100%[===================>] 778.43M   169MB/s    in 4.7s    \n",
            "\n",
            "2021-02-27 10:54:27 (165 MB/s) - ‘hg19.2bit’ saved [816241703/816241703]\n",
            "\n",
            "--2021-02-27 10:54:27--  http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa\n",
            "Resolving hgdownload.soe.ucsc.edu (hgdownload.soe.ucsc.edu)... 128.114.119.163\n",
            "Connecting to hgdownload.soe.ucsc.edu (hgdownload.soe.ucsc.edu)|128.114.119.163|:80... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 9326144 (8.9M)\n",
            "Saving to: ‘twoBitToFa’\n",
            "\n",
            "twoBitToFa          100%[===================>]   8.89M  56.1MB/s    in 0.2s    \n",
            "\n",
            "2021-02-27 10:54:27 (56.1 MB/s) - ‘twoBitToFa’ saved [9326144/9326144]\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pmXn9gkPaMi6"
      },
      "source": [
        "### Process the downloaded data to get the fasta sequences for pos and neg set"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_bYqIeK4WdxH"
      },
      "source": [
        "#positive set is 200bp around the consiervative IDR otpimal peaks\n",
        "! zcat conservative_peaks.bed.gz | perl -lane 'print($F[0].\"\\t\".($F[1]+$F[9]-100).\"\\t\".($F[1]+$F[9]+100))' | gzip -c > positive_set_full.bed.gz\n",
        "#negative set is 200bp around summit of H1 accessible peaks that don't overlap 1kb within any optimal or conservative peak\n",
        "! zcat conservative_peaks.bed.gz optimal_peaks.bed.gz | perl -lane 'print($F[0].\"\\t\".($F[1]+$F[9]-500).\"\\t\".($F[1]+$F[9]+500))' | gzip -c > 1kb_around_optimal_or_conservative_peaks.bed.gz\n",
        "! zcat bg_dnase.bed.gz | perl -lane 'print($F[0].\"\\t\".($F[1]+$F[9]-100).\"\\t\".($F[1]+$F[9]+100))' | gzip -c > prefiltering_neg_set.bed.gz\n",
        "! bedtools intersect -a prefiltering_neg_set.bed.gz -b 1kb_around_optimal_or_conservative_peaks.bed.gz -v -wa | gzip -c > neg_set_full.bed.gz"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oNNw4gOLcD9F"
      },
      "source": [
        "#subsample the negative set to have approx. the same number of regions as the pos set\n",
        "! zcat neg_set_full.bed.gz | perl -lane 'if ($.%20==1) {print $_}' | gzip -c > subsampled_neg_set.bed.gz"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kEGc09DucJZ0"
      },
      "source": [
        "#use chr1 and 2 for the test set                                                \n",
        "! zcat positive_set_full.bed.gz | egrep -w 'chr1|chr2' | gzip -c > positives_test_set.bed.gz\n",
        "! zcat positive_set_full.bed.gz | egrep -w -v 'chr1|chr2' | gzip -c > positives_train_set.bed.gz\n",
        "! zcat subsampled_neg_set.bed.gz | egrep -w 'chr1|chr2' | gzip -c > negatives_test_set.bed.gz\n",
        "! zcat subsampled_neg_set.bed.gz | egrep -w -v 'chr1|chr2' | gzip -c > negatives_train_set.bed.gz"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mfMgskXDcMUD",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "88165321-de05-48b3-f762-4a102abdb2c5"
      },
      "source": [
        "#Extract the underlying fasta regions\n",
        "! bedtools getfasta -fi hg19.genome.fa -bed positives_train_set.bed.gz > positives_train.fa\n",
        "! bedtools getfasta -fi hg19.genome.fa -bed positives_test_set.bed.gz > positives_test.fa\n",
        "! bedtools getfasta -fi hg19.genome.fa -bed negatives_train_set.bed.gz > negatives_train.fa\n",
        "! bedtools getfasta -fi hg19.genome.fa -bed negatives_test_set.bed.gz > negatives_test.fa"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "index file hg19.genome.fa.fai not found, generating...\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aFziRxc9dgZq"
      },
      "source": [
        "### Train the model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Jb2xThYMcZgc",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "dfa63452-5a30-4fb4-a10c-b4a193f448e0"
      },
      "source": [
        "#To save time, we can download the pre-trained model\n",
        "! [[ -e lsgkm_defaultsettings_t2.model.txt.gz ]] || wget https://raw.githubusercontent.com/AvantiShri/model_storage/5dcfc2b/gkmexplain/lsgkm_defaultsettings_t2.model.txt.gz -O lsgkm_defaultsettings_t2.model.txt.gz\n",
        "! zcat lsgkm_defaultsettings_t2.model.txt.gz > lsgkm_defaultsettings_t2.model.txt\n",
        "\n",
        "#To train the model from scratch, run the lines below:\n",
        "##Model is trained using kernel=2, which is the standard gkm kernel (no position weighting)\n",
        "## I used the standard gkm kernel so that the method from Gandhi et al. (2014)\n",
        "## would be applicable\n",
        "! [[ -e lsgkm_defaultsettings_t2.model.txt ]] || lsgkm/src/gkmtrain -T 16 -t 2 positives_train.fa negatives_train.fa lsgkm_defaultsettings_t2\n",
        "\n",
        "#Make predictions to assess performance\n",
        "! [[ -e preds_test_positives.txt ]] || lsgkm/src/gkmpredict -T 16 positives_test.fa lsgkm_defaultsettings_t2.model.txt preds_test_positives.txt\n",
        "! [[ -e preds_test_negatives.txt ]] || lsgkm/src/gkmpredict -T 16 negatives_test.fa lsgkm_defaultsettings_t2.model.txt preds_test_negatives.txt\n",
        "\n",
        "from sklearn.metrics import roc_auc_score                                                                                                                    \n",
        "pos_preds = [float(x.rstrip().split(\"\\t\")[1])                                   \n",
        "             for x in open(\"preds_test_positives.txt\")]                         \n",
        "neg_preds = [float(x.rstrip().split(\"\\t\")[1])                                   \n",
        "             for x in open(\"preds_test_negatives.txt\")]                         \n",
        "print(roc_auc_score(y_true=[1 for x in pos_preds]+[0 for x in neg_preds],       \n",
        "                    y_score = pos_preds+neg_preds))"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "--2021-02-27 10:54:51--  https://raw.githubusercontent.com/AvantiShri/model_storage/5dcfc2b/gkmexplain/lsgkm_defaultsettings_t2.model.txt.gz\n",
            "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...\n",
            "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 655201 (640K) [application/octet-stream]\n",
            "Saving to: ‘lsgkm_defaultsettings_t2.model.txt.gz’\n",
            "\n",
            "lsgkm_defaultsettin 100%[===================>] 639.84K  --.-KB/s    in 0.03s   \n",
            "\n",
            "2021-02-27 10:54:51 (18.9 MB/s) - ‘lsgkm_defaultsettings_t2.model.txt.gz’ saved [655201/655201]\n",
            "\n",
            "INFO 2021-02-27 10:54:51: Number of threads is set to 16\n",
            "INFO 2021-02-27 10:54:51: load model lsgkm_defaultsettings_t2.model.txt\n",
            "INFO 2021-02-27 10:54:52: reading... 1000/8873\n",
            "INFO 2021-02-27 10:54:53: reading... 2000/8873\n",
            "INFO 2021-02-27 10:54:53: reading... 3000/8873\n",
            "INFO 2021-02-27 10:54:54: reading... 4000/8873\n",
            "INFO 2021-02-27 10:54:54: reading... 5000/8873\n",
            "INFO 2021-02-27 10:54:55: reading... 6000/8873\n",
            "INFO 2021-02-27 10:54:55: reading... 7000/8873\n",
            "INFO 2021-02-27 10:54:56: reading... 8000/8873\n",
            "INFO 2021-02-27 10:54:57: write prediction result to preds_test_positives.txt\n",
            "INFO 2021-02-27 10:55:01: 100 scored\n",
            "INFO 2021-02-27 10:55:05: 200 scored\n",
            "INFO 2021-02-27 10:55:10: 300 scored\n",
            "INFO 2021-02-27 10:55:14: 400 scored\n",
            "INFO 2021-02-27 10:55:19: 500 scored\n",
            "INFO 2021-02-27 10:55:23: 600 scored\n",
            "INFO 2021-02-27 10:55:28: 700 scored\n",
            "INFO 2021-02-27 10:55:32: 800 scored\n",
            "INFO 2021-02-27 10:55:36: 900 scored\n",
            "INFO 2021-02-27 10:55:39: 960 scored\n",
            "INFO 2021-02-27 10:55:39: Number of threads is set to 16\n",
            "INFO 2021-02-27 10:55:39: load model lsgkm_defaultsettings_t2.model.txt\n",
            "INFO 2021-02-27 10:55:40: reading... 1000/8873\n",
            "INFO 2021-02-27 10:55:40: reading... 2000/8873\n",
            "INFO 2021-02-27 10:55:41: reading... 3000/8873\n",
            "INFO 2021-02-27 10:55:42: reading... 4000/8873\n",
            "INFO 2021-02-27 10:55:42: reading... 5000/8873\n",
            "INFO 2021-02-27 10:55:43: reading... 6000/8873\n",
            "INFO 2021-02-27 10:55:43: reading... 7000/8873\n",
            "INFO 2021-02-27 10:55:44: reading... 8000/8873\n",
            "INFO 2021-02-27 10:55:45: write prediction result to preds_test_negatives.txt\n",
            "INFO 2021-02-27 10:55:49: 100 scored\n",
            "INFO 2021-02-27 10:55:53: 200 scored\n",
            "INFO 2021-02-27 10:55:58: 300 scored\n",
            "INFO 2021-02-27 10:56:02: 400 scored\n",
            "INFO 2021-02-27 10:56:06: 500 scored\n",
            "INFO 2021-02-27 10:56:11: 600 scored\n",
            "INFO 2021-02-27 10:56:15: 700 scored\n",
            "INFO 2021-02-27 10:56:20: 800 scored\n",
            "INFO 2021-02-27 10:56:24: 900 scored\n",
            "INFO 2021-02-27 10:56:27: 964 scored\n",
            "0.83056978302213\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XwEKM0cquHD3"
      },
      "source": [
        "### Generate dinuc shuffled sequences for computing null distribution of importance scores"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RcXPVrB1uNUO",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "77f053a0-09d1-42ac-b917-b411593be11c"
      },
      "source": [
        "!pip install deeplift\n",
        "from deeplift.dinuc_shuffle import dinuc_shuffle\n",
        "\n",
        "import numpy as np\n",
        "import random\n",
        "np.random.seed(1234)\n",
        "random.seed(1234)\n",
        "\n",
        "num_dinuc_shuffled_seqs = 500\n",
        "#Generate the dinucleotide shuffled sequences and write to a file\n",
        "fasta_seqs_no_N = [x.rstrip() for (i,x) in enumerate(open(\"positives_test.fa\"))\n",
        "                   if (i%2==1 and ('N' not in x))]\n",
        "open(\"dnshuff_seqs.fa\", 'w').write(\n",
        " \"\\n\".join([\">seq\"+str(i)+\"\\n\"+dinuc_shuffle(\n",
        "            str(np.random.choice(fasta_seqs_no_N)))\n",
        "            for i in range(num_dinuc_shuffled_seqs)]))\n",
        "\n",
        "#We can also download the pre-generated file\n",
        "! [[ -e dnshuff_seqs.fa.gz ]] || wget https://raw.githubusercontent.com/AvantiShri/model_storage/aae0902/gkmexplain/dnshuff_seqs.fa.gz -O dnshuff_seqs.fa.gz\n",
        "! zcat dnshuff_seqs.fa.gz > dnshuff_seqs.fa"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Collecting deeplift\n",
            "  Downloading https://files.pythonhosted.org/packages/d2/48/e8c4a331664c32682d6f7f55f1148f59224e32cbf4f22c90f3f961eb5a40/deeplift-0.6.13.0.tar.gz\n",
            "Requirement already satisfied: numpy>=1.9 in /usr/local/lib/python3.7/dist-packages (from deeplift) (1.19.5)\n",
            "Building wheels for collected packages: deeplift\n",
            "  Building wheel for deeplift (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Created wheel for deeplift: filename=deeplift-0.6.13.0-cp37-none-any.whl size=36447 sha256=c56e49565a3d5f89ba24a9e5adb3f8fff57b4c834df85a1da36aeb0ec907eb0e\n",
            "  Stored in directory: /root/.cache/pip/wheels/39/a2/1b/a2eac3afbfedc4fb40e213ec4f8d97d26598325f187ae0dc5a\n",
            "Successfully built deeplift\n",
            "Installing collected packages: deeplift\n",
            "Successfully installed deeplift-0.6.13.0\n",
            "--2021-02-27 10:56:32--  https://raw.githubusercontent.com/AvantiShri/model_storage/aae0902/gkmexplain/dnshuff_seqs.fa.gz\n",
            "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...\n",
            "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 31875 (31K) [application/octet-stream]\n",
            "Saving to: ‘dnshuff_seqs.fa.gz’\n",
            "\n",
            "dnshuff_seqs.fa.gz  100%[===================>]  31.13K  --.-KB/s    in 0.002s  \n",
            "\n",
            "2021-02-27 10:56:32 (16.1 MB/s) - ‘dnshuff_seqs.fa.gz’ saved [31875/31875]\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kh0TPYGMr3_u"
      },
      "source": [
        "### Compute hypothetical importance scores"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Cb4Jedzvr3B7",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "df6d1ec0-d9d1-4d1d-da45-e64ff29bb649"
      },
      "source": [
        "#The actual importance scores can be derived from the hypothetical importance\n",
        "# scores by doing an elementwise multiplication of the hypothetical importance\n",
        "# with the one-hot encoded sequence.\n",
        "! [[ -e gkmexplain_positives_hypimpscores.txt ]] || lsgkm/src/gkmexplain -m 1 positives_test.fa lsgkm_defaultsettings_t2.model.txt gkmexplain_positives_hypimpscores.txt\n",
        "! [[ -e gkmexplain_dnshuff_hypimpscores.txt ]] || lsgkm/src/gkmexplain -m 1 dnshuff_seqs.fa lsgkm_defaultsettings_t2.model.txt gkmexplain_dnshuff_hypimpscores.txt"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "INFO 2021-02-27 10:56:33: Number of threads is set to 1\n",
            "INFO 2021-02-27 10:56:33: load model lsgkm_defaultsettings_t2.model.txt\n",
            "INFO 2021-02-27 10:56:33: reading... 1000/8873\n",
            "INFO 2021-02-27 10:56:33: reading... 2000/8873\n",
            "INFO 2021-02-27 10:56:33: reading... 3000/8873\n",
            "INFO 2021-02-27 10:56:33: reading... 4000/8873\n",
            "INFO 2021-02-27 10:56:34: reading... 5000/8873\n",
            "INFO 2021-02-27 10:56:34: reading... 6000/8873\n",
            "INFO 2021-02-27 10:56:34: reading... 7000/8873\n",
            "INFO 2021-02-27 10:56:34: reading... 8000/8873\n",
            "INFO 2021-02-27 10:56:34: write prediction result to gkmexplain_positives_hypimpscores.txt\n",
            "INFO 2021-02-27 10:57:47: 100 scored\n",
            "INFO 2021-02-27 10:58:58: 200 scored\n",
            "INFO 2021-02-27 11:00:10: 300 scored\n",
            "INFO 2021-02-27 11:01:21: 400 scored\n",
            "INFO 2021-02-27 11:02:32: 500 scored\n",
            "INFO 2021-02-27 11:03:41: 600 scored\n",
            "INFO 2021-02-27 11:04:50: 700 scored\n",
            "INFO 2021-02-27 11:05:59: 800 scored\n",
            "INFO 2021-02-27 11:07:07: 900 scored\n",
            "INFO 2021-02-27 11:07:48: 960 scored\n",
            "INFO 2021-02-27 11:07:48: Number of threads is set to 1\n",
            "INFO 2021-02-27 11:07:48: load model lsgkm_defaultsettings_t2.model.txt\n",
            "INFO 2021-02-27 11:07:48: reading... 1000/8873\n",
            "INFO 2021-02-27 11:07:48: reading... 2000/8873\n",
            "INFO 2021-02-27 11:07:48: reading... 3000/8873\n",
            "INFO 2021-02-27 11:07:49: reading... 4000/8873\n",
            "INFO 2021-02-27 11:07:49: reading... 5000/8873\n",
            "INFO 2021-02-27 11:07:49: reading... 6000/8873\n",
            "INFO 2021-02-27 11:07:49: reading... 7000/8873\n",
            "INFO 2021-02-27 11:07:49: reading... 8000/8873\n",
            "INFO 2021-02-27 11:07:49: write prediction result to gkmexplain_dnshuff_hypimpscores.txt\n",
            "INFO 2021-02-27 11:08:57: 100 scored\n",
            "INFO 2021-02-27 11:10:04: 200 scored\n",
            "INFO 2021-02-27 11:11:12: 300 scored\n",
            "INFO 2021-02-27 11:12:19: 400 scored\n",
            "INFO 2021-02-27 11:13:26: 500 scored\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yWoIBvlBXP7p"
      },
      "source": [
        ""
      ],
      "execution_count": 8,
      "outputs": []
    }
  ]
}