get_data_example_notebook.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "collapsed_sections": [
        "-jhshEcOx3Ze"
      ],
      "authorship_tag": "ABX9TyMSi21woKYMVogfY3E2XbXH"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "**In this notebook you will see how to go from fastq sequencing files to Unspliced and Spliced count matrices**\n",
        "\n",
        "We use the loom file output to store the U and S matrices for meK-Means inference."
      ],
      "metadata": {
        "id": "txlOjYhZxoty"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## **Install Packages**"
      ],
      "metadata": {
        "id": "-jhshEcOx3Ze"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip3 install loompy --quiet"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "cwHLCYylyn1U",
        "outputId": "06a824de-ec71-4914-829c-3e2fa63e167c"
      },
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m4.8/4.8 MB\u001b[0m \u001b[31m14.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25h  Preparing metadata (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Building wheel for loompy (setup.py) ... \u001b[?25l\u001b[?25hdone\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install kb-python==0.27.2 --quiet"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "95UBhb-x0wmX",
        "outputId": "9a91f9b9-ddfc-4b5f-d3c7-3ccb6fc22a36"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m7.5/7.5 MB\u001b[0m \u001b[31m2.0 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25h  Preparing metadata (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m122.1/122.1 kB\u001b[0m \u001b[31m12.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m45.2/45.2 MB\u001b[0m \u001b[31m19.4 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m2.1/2.1 MB\u001b[0m \u001b[31m71.5 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m21.9/21.9 MB\u001b[0m \u001b[31m54.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m56.8/56.8 kB\u001b[0m \u001b[31m6.8 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25h  Preparing metadata (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m85.7/85.7 kB\u001b[0m \u001b[31m10.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m79.8/79.8 kB\u001b[0m \u001b[31m9.0 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25h  Building wheel for kb-python (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Building wheel for session-info (setup.py) ... \u001b[?25l\u001b[?25hdone\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "mV-VU3ioxUNB",
        "outputId": "ff67fa89-aacf-4051-d3ab-32abae579128"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.10/dist-packages/loompy/bus_file.py:68: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
            "  def twobit_to_dna(twobit: int, size: int) -> str:\n",
            "/usr/local/lib/python3.10/dist-packages/loompy/bus_file.py:85: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
            "  def dna_to_twobit(dna: str) -> int:\n",
            "/usr/local/lib/python3.10/dist-packages/loompy/bus_file.py:102: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
            "  def twobit_1hamming(twobit: int, size: int) -> List[int]:\n"
          ]
        }
      ],
      "source": [
        "import loompy as lp\n",
        "\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!wget --output-document sratoolkit.tar.gz https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-ubuntu64.tar.gz"
      ],
      "metadata": {
        "id": "brx4yhJbykUB"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!tar -vxzf sratoolkit.tar.gz"
      ],
      "metadata": {
        "id": "ke-NDg910Hnq"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "sFmk08sL0SRo"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "OIKuIveR0SUM"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## **Download fastqs**\n",
        "\n",
        "> We will be using the scMixology dataset with three human lung adenocarcinoma cell lines HCC827, H1975 and H2228 mixed together and sequenced with 10xv2 scRNAseq reagents (see the sc_10x data on the [scMixology Github](https://github.com/LuyiTian/sc_mixology))."
      ],
      "metadata": {
        "id": "YlKbnAysyktK"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#3 cell line mix GSM3022245\n",
        "#Only download one for now because of Colab disk limits\n",
        "\n",
        "!./sratoolkit.3.1.0-ubuntu64/bin/prefetch SRR6782109 --max-size 6000000000 -O ./ && ./sratoolkit.3.1.0-ubuntu64/bin/fasterq-dump --include-technical --split-files SRR6782109\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nh81DqFnynOO",
        "outputId": "0ad439fa-02a0-4f88-d30e-29149b1e10a0"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "2024-04-04T17:41:03 prefetch.3.1.0: Current preference is set to retrieve SRA Normalized Format files with full base quality scores.\n",
            "2024-04-04T17:41:04 prefetch.3.1.0: 1) Downloading 'SRR6782109'...\n",
            "2024-04-04T17:41:04 prefetch.3.1.0: SRA Normalized Format file is being retrieved, if this is different from your preference, it may be due to current file availability.\n",
            "2024-04-04T17:41:04 prefetch.3.1.0:  Downloading via HTTPS...\n",
            "2024-04-04T17:47:17 prefetch.3.1.0:  HTTPS download succeed\n",
            "2024-04-04T17:47:43 prefetch.3.1.0:  'SRR6782109' is valid\n",
            "2024-04-04T17:47:43 prefetch.3.1.0: 1) 'SRR6782109' was downloaded successfully\n",
            "2024-04-04T17:47:43 prefetch.3.1.0: 'SRR6782109' has 0 unresolved dependencies\n",
            "spots read      : 109,178,700\n",
            "reads read      : 218,357,400\n",
            "reads written   : 218,357,400\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "p4rHAfFBapQY"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "btQYGAshzkBU"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "OuV3E_kBzkDp"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## **Run kb-python to get U/S Counts**\n",
        "\n",
        "> We will get loom file output from the kb-python run which stores these count matrices, and we will use the cells in the filtered barcodes output."
      ],
      "metadata": {
        "id": "6qh-IAqk1Ww1"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "kIZa3DCm1i6v"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Since this is human data we need the human reference kallisto index and intron/exon annotations (GRCh38-2020).\n",
        "\n"
      ],
      "metadata": {
        "id": "AOTwGdN-1loV"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "SA2DexRsCWvC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Read in saved index and intron/exon files from kb ref\n",
        "\n",
        "#If wget does not work, use the download link here https://doi.org/10.22002/6wyra-tar37 for the refdata-gex-GRCh38-2020-A.tar.gz\n",
        "!wget --content-disposition https://data.caltech.edu/records/6wyra-tar37/files/refdata-gex-GRCh38-2020-A.tar.gz?download=1\n",
        "!tar -vxzf refdata-gex-GRCh38-2020-A.tar.gz"
      ],
      "metadata": {
        "id": "GLrI5HL1Bhe_",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "61d9f827-9ef8-46fb-cfd3-771d740edea5"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "--2024-04-04 18:08:33--  https://data.caltech.edu/records/6wyra-tar37/files/refdata-gex-GRCh38-2020-A.tar.gz?download=1\n",
            "Resolving data.caltech.edu (data.caltech.edu)... 35.155.11.48\n",
            "Connecting to data.caltech.edu (data.caltech.edu)|35.155.11.48|:443... connected.\n",
            "HTTP request sent, awaiting response... 302 FOUND\n",
            "Location: https://s3.us-west-2.amazonaws.com/caltechdata/ad/fe/e87c-c65f-4034-ae56-2b57af502588/data?response-content-type=application%2Foctet-stream&response-content-disposition=attachment%3B%20filename%3Drefdata-gex-GRCh38-2020-A.tar.gz&X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=AKIARCVIVNNAP7NNDVEA%2F20240404%2Fus-west-2%2Fs3%2Faws4_request&X-Amz-Date=20240404T180834Z&X-Amz-Expires=60&X-Amz-SignedHeaders=host&X-Amz-Signature=26bf1cf1422719e595adaca4cf05fc13d2d6d3851c1a6dc16df8121e0d8350af [following]\n",
            "--2024-04-04 18:08:34--  https://s3.us-west-2.amazonaws.com/caltechdata/ad/fe/e87c-c65f-4034-ae56-2b57af502588/data?response-content-type=application%2Foctet-stream&response-content-disposition=attachment%3B%20filename%3Drefdata-gex-GRCh38-2020-A.tar.gz&X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=AKIARCVIVNNAP7NNDVEA%2F20240404%2Fus-west-2%2Fs3%2Faws4_request&X-Amz-Date=20240404T180834Z&X-Amz-Expires=60&X-Amz-SignedHeaders=host&X-Amz-Signature=26bf1cf1422719e595adaca4cf05fc13d2d6d3851c1a6dc16df8121e0d8350af\n",
            "Resolving s3.us-west-2.amazonaws.com (s3.us-west-2.amazonaws.com)... 52.92.237.144, 52.218.153.24, 52.92.196.248, ...\n",
            "Connecting to s3.us-west-2.amazonaws.com (s3.us-west-2.amazonaws.com)|52.92.237.144|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 31441161559 (29G) [application/octet-stream]\n",
            "Saving to: ‘refdata-gex-GRCh38-2020-A.tar.gz’\n",
            "\n",
            "refdata-gex-GRCh38- 100%[===================>]  29.28G  19.0MB/s    in 29m 57s \n",
            "\n",
            "2024-04-04 18:38:32 (16.7 MB/s) - ‘refdata-gex-GRCh38-2020-A.tar.gz’ saved [31441161559/31441161559]\n",
            "\n",
            "refdata-gex-GRCh38-2020-A/\n",
            "refdata-gex-GRCh38-2020-A/cdna_t2c.txt\n",
            "refdata-gex-GRCh38-2020-A/index.idx\n",
            "tar: refdata-gex-GRCh38-2020-A/index.idx: Wrote only 6656 of 10240 bytes\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# #Equivalent mouse references here:\n",
        "# !wget --content-disposition https://data.caltech.edu/records/1dd7a-cc411/files/refdata-gex-mm10-2020-A.tar.gz?download=1\n",
        "# !tar -vxzf refdata-gex-mm10-2020-A.tar.gz"
      ],
      "metadata": {
        "id": "JqUzBeaY81yH"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!mkdir scMix"
      ],
      "metadata": {
        "id": "ZMV1pDzQ1pje"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "! kb count --verbose \\\n",
        "-i ./refdata-gex-GRCh38-2020-A/index.idx \\\n",
        "-g ./refdata-gex-GRCh38-2020-A/t2g_grch38.txt \\\n",
        "-x 10xv2 \\\n",
        "-o ./scMix \\\n",
        "-t 2 \\\n",
        "-c1 ./refdata-gex-GRCh38-2020-A/cdna_t2c.txt \\\n",
        "-c2 ./refdata-gex-GRCh38-2020-A/intron_t2c.txt \\\n",
        "--workflow lamanno --filter bustools --overwrite --loom \\\n",
        "./SRR6782109_1.fastq \\\n",
        "./SRR6782109_2.fastq\n",
        "\n",
        "#There are only two threads, t, on Colab"
      ],
      "metadata": {
        "id": "TjTd74HD1fzT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Commented out code below shows how to create the index and txt files above\n",
        "(Will need a new index for kb-python >=0.28)"
      ],
      "metadata": {
        "id": "MLuMoNGYMXJN"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# #Make kallisto index, -->runs out of memory on Colab\n",
        "\n",
        "# #Download genome information\n",
        "\n",
        "# # This does\n",
        "# # 1.   Make a kallisto index of introns and exons.\n",
        "# # 2.   Generate U/S counts for each cell across the reference transcriptome.\n",
        "\n",
        "\n",
        "# !wget --content-disposition https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz\n",
        "# #refdata-gex-mm10-2020-A for mouse\n",
        "\n",
        "# !tar -zxvf refdata-gex-GRCh38-2020-A.tar.gz\n",
        "\n",
        "# !kb ref -i ./refdata-gex-GRCh38-2020-A/index.idx \\\n",
        "# -g ./refdata-gex-GRCh38-2020-A/t2g_grch38.txt \\\n",
        "# -f1 ./refdata-gex-GRCh38-2020-A/cdna.fa \\\n",
        "# -f2 ./refdata-gex-GRCh38-2020-A/intron.fa \\\n",
        "# -c1 ./refdata-gex-GRCh38-2020-A/cdna_t2c.txt \\\n",
        "# -c2 ./refdata-gex-GRCh38-2020-A/intron_t2c.txt \\\n",
        "# --workflow lamanno \\\n",
        "# ./refdata-gex-GRCh38-2020-A/fasta/genome.fa \\\n",
        "# ./refdata-gex-GRCh38-2020-A/genes/genes.gtf"
      ],
      "metadata": {
        "id": "GmgMdBgfWF9C"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "pVRwRhWYWGC-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Here is the final loom file that we will use for analysis**"
      ],
      "metadata": {
        "id": "_YyzpaLXWH5Y"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "ds = lp.connect('./scMix/counts_filtered/adata.loom')"
      ],
      "metadata": {
        "id": "4R-8St67WKWr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "unspliced and spliced counts are stored in the layers"
      ],
      "metadata": {
        "id": "XUfk26VdWK8a"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#matrices are genexcell\n",
        "print(ds.layers['unspliced'][:,:].shape) #or 'spliced'"
      ],
      "metadata": {
        "id": "cyLOJQz8WOJV"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "gene names and cell barcodes are stored under 'gene_name' and 'barcode' respectively"
      ],
      "metadata": {
        "id": "WXRMNVmMWQAX"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "ds.ra['gene_name'][0:10]"
      ],
      "metadata": {
        "id": "HLCKJqaWWeFz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "ds.ca['barcode'][0:10]"
      ],
      "metadata": {
        "id": "fJzPCUyX0WOW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "ds.close()"
      ],
      "metadata": {
        "id": "nCtH6-sbfze7"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}