dataset-notebooks/10x_neuron_1k_v3chem_python/10x_neuron_1k_v3chem.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys, collections, os, argparse\n",
    "%matplotlib inline  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Download the 10x Dataset `1k Brain Cells from an E18 Mouse (v3 chemistry)`\n",
    "\n",
    "10x datasets are available at\n",
    "https://support.10xgenomics.com/single-cell-gene-expression/datasets\n",
    "\n",
    "The page for the `1k Brain Cells from an E18 Mouse (v3 chemistry)` dataset is\n",
    "https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0/neuron_1k_v3\n",
    "\n",
    "But the FASTQ files (7.4GB) can be downloaded with `wget` directly (without giving them email info) from http://cf.10xgenomics.com/samples/cell-exp/3.0.0/neuron_1k_v3/neuron_1k_v3_fastqs.tar\n",
    "\n",
    "In the cell below we check if the dataset file `neuron_1k_v3_fastqs.tar` already exists. If not we download the dataset to the same directory as this notebook\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2018-12-03 17:29:29--  http://cf.10xgenomics.com/samples/cell-exp/3.0.0/neuron_1k_v3/neuron_1k_v3_fastqs.tar\n",
      "Resolving cf.10xgenomics.com (cf.10xgenomics.com)... 13.35.99.77, 13.35.99.80, 13.35.99.113, ...\n",
      "Connecting to cf.10xgenomics.com (cf.10xgenomics.com)|13.35.99.77|:80... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 7086786560 (6.6G) [application/x-tar]\n",
      "Saving to: ‘neuron_1k_v3_fastqs.tar’\n",
      "\n",
      "100%[====================================>] 7,086,786,560 85.3MB/s   in 78s    \n",
      "\n",
      "2018-12-03 17:30:52 (86.8 MB/s) - ‘neuron_1k_v3_fastqs.tar’ saved [7086786560/7086786560]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Check if the file was downloaded already before doing wget:\n",
    "if not (os.path.isfile('./neuron_1k_v3_fastqs.tar')): \n",
    "    # the `!` means we're running a command line statement (rather than python)  \n",
    "    !wget http://cf.10xgenomics.com/samples/cell-exp/3.0.0/neuron_1k_v3/neuron_1k_v3_fastqs.tar\n",
    "else: print('Dataset already downloaded!')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because this dataset was run on two lanes, we need to uncompress the FASTQ files and concatenate them before using kallisto bus. If we had a single file kallisto could take gz files as is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "neuron_1k_v3_fastqs/\n",
      "neuron_1k_v3_fastqs/neuron_1k_v3_S1_L002_I1_001.fastq.gz\n",
      "neuron_1k_v3_fastqs/neuron_1k_v3_S1_L001_R2_001.fastq.gz\n",
      "neuron_1k_v3_fastqs/neuron_1k_v3_S1_L001_R1_001.fastq.gz\n",
      "neuron_1k_v3_fastqs/neuron_1k_v3_S1_L002_R1_001.fastq.gz\n",
      "neuron_1k_v3_fastqs/neuron_1k_v3_S1_L001_I1_001.fastq.gz\n",
      "neuron_1k_v3_fastqs/neuron_1k_v3_S1_L002_R2_001.fastq.gz\n"
     ]
    }
   ],
   "source": [
    "# now we untar the fastq files on neuron_1k_v3_fastqs folder\n",
    "!tar -xvf ./neuron_1k_v3_fastqs.tar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Buiding the kallisto index\n",
    "\n",
    "First make sure that kallisto is installed and the version is greater than 0.45\n",
    "\n",
    "If it's not installed, see instructions at https://pachterlab.github.io/kallisto/download"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kallisto, version 0.45.0\n"
     ]
    }
   ],
   "source": [
    "!kallisto version "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we build the kallisto index for the dataset. \n",
    "The index is built from the published reference transcriptome for each organism.\n",
    "Building the index takes a few minutes and needs to be done only once for each organism.\n",
    "\n",
    "### Download reference transcriptome from ensembl\n",
    "In order to do that we first download the mouse transcriptome from ensembl, you can see the reference genomes they have at https://uswest.ensembl.org/info/data/ftp/index.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2018-12-03 17:31:00--  ftp://ftp.ensembl.org/pub/release-94/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz\n",
      "           => ‘Mus_musculus.GRCm38.cdna.all.fa.gz’\n",
      "Resolving ftp.ensembl.org (ftp.ensembl.org)... 193.62.193.8\n",
      "Connecting to ftp.ensembl.org (ftp.ensembl.org)|193.62.193.8|:21... connected.\n",
      "Logging in as anonymous ... Logged in!\n",
      "==> SYST ... done.    ==> PWD ... done.\n",
      "==> TYPE I ... done.  ==> CWD (1) /pub/release-94/fasta/mus_musculus/cdna ... done.\n",
      "==> SIZE Mus_musculus.GRCm38.cdna.all.fa.gz ... 50809568\n",
      "==> PASV ... done.    ==> RETR Mus_musculus.GRCm38.cdna.all.fa.gz ... done.\n",
      "Length: 50809568 (48M) (unauthoritative)\n",
      "\n",
      "100%[======================================>] 50,809,568  3.24MB/s   in 15s    \n",
      "\n",
      "2018-12-03 17:31:22 (3.26 MB/s) - ‘Mus_musculus.GRCm38.cdna.all.fa.gz’ saved [50809568]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Check if the file was downloaded already before doing wget:\n",
    "if not (os.path.isfile('Mus_musculus.GRCm38.cdna.all.fa.gz')): \n",
    "    # the `!` means we're running a command line statement (rather than python)  \n",
    "    !wget ftp://ftp.ensembl.org/pub/release-94/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz\n",
    "else: print('Mouse transcriptome already downloaded!')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "[build] loading fasta file Mus_musculus.GRCm38.cdna.all.fa.gz\n",
      "[build] k-mer length: 31\n",
      "[build] warning: clipped off poly-A tail (longer than 10)\n",
      "        from 600 target sequences\n",
      "[build] warning: replaced 3 non-ACGUT characters in the input sequence\n",
      "        with pseudorandom nucleotides\n",
      "[build] counting k-mers ... done.\n",
      "[build] building target de Bruijn graph ...  done \n",
      "[build] creating equivalence classes ...  done\n",
      "[build] target de Bruijn graph has 711215 contigs and contains 98989067 k-mers \n",
      "\n"
     ]
    }
   ],
   "source": [
    "### Now we can build the index\n",
    "if not (os.path.isfile('mouse_transcripts.idx')): \n",
    "    !kallisto index -i mouse_transcripts.idx Mus_musculus.GRCm38.cdna.all.fa.gz\n",
    "else: print ('Mouse transcript index already exist!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparing transcript_to_gene.tsv file process the single cell data with kallisto bus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depending on which transcriptome you used, you will need to create a file translating transcripts to genes. This notebook assumes the file is in `transcript_to_gene.tsv`, for ensembl transcriptomes these can be generated using biomart.\n",
    "\n",
    "The general format of `transcript_to_gene.tsv` is\n",
    "\n",
    "```\n",
    "ENST00000632684.1\tENSG00000282431.1\n",
    "ENST00000434970.2\tENSG00000237235.2\n",
    "ENST00000448914.1\tENSG00000228985.1\n",
    "ENST00000415118.1\tENSG00000223997.1\n",
    "ENST00000631435.1\tENSG00000282253.1\n",
    "...\n",
    "```\n",
    "\n",
    "To create the `transcript_to_gene.tsv` we fetch and parse the mouse GTF file from ensembl.\n",
    "\n",
    "The reference GTF files are available at https://uswest.ensembl.org/info/data/ftp/index.html\n",
    "\n",
    "The mouse ones which we use are at ftp://ftp.ensembl.org/pub/release-94/gtf/mus_musculus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2018-12-03 17:36:33--  ftp://ftp.ensembl.org/pub/release-94/gtf/mus_musculus/Mus_musculus.GRCm38.94.gtf.gz\n",
      "           => ‘Mus_musculus.GRCm38.94.gtf.gz’\n",
      "Resolving ftp.ensembl.org (ftp.ensembl.org)... 193.62.193.8\n",
      "Connecting to ftp.ensembl.org (ftp.ensembl.org)|193.62.193.8|:21... connected.\n",
      "Logging in as anonymous ... Logged in!\n",
      "==> SYST ... done.    ==> PWD ... done.\n",
      "==> TYPE I ... done.  ==> CWD (1) /pub/release-94/gtf/mus_musculus ... done.\n",
      "==> SIZE Mus_musculus.GRCm38.94.gtf.gz ... 29397390\n",
      "==> PASV ... done.    ==> RETR Mus_musculus.GRCm38.94.gtf.gz ... done.\n",
      "Length: 29397390 (28M) (unauthoritative)\n",
      "\n",
      "100%[======================================>] 29,397,390  2.01MB/s   in 15s    \n",
      "\n",
      "2018-12-03 17:36:55 (1.90 MB/s) - ‘Mus_musculus.GRCm38.94.gtf.gz’ saved [29397390]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Check if the file was downloaded already before doing wget:\n",
    "if not (os.path.isfile('Mus_musculus.GRCm38.94.gtf.gz') or os.path.isfile('Mus_musculus.GRCm38.94.gtf')): \n",
    "    # the `!` means we're running a command line statement (rather than python)  \n",
    "    !wget ftp://ftp.ensembl.org/pub/release-94/gtf/mus_musculus/Mus_musculus.GRCm38.94.gtf.gz\n",
    "else: print('Mouse transcriptome already downloaded!')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Unzip the file\n",
    "!gunzip ./Mus_musculus.GRCm38.94.gtf.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create transcript_to_gene.tsv\n",
    "\n",
    "Now we can use the cells below to parse the GTF file and keep only the transcript mapping as a tsv file in the format below.\n",
    "```\n",
    "ENST00000632684.1\tENSG00000282431.1\n",
    "ENST00000434970.2\tENSG00000237235.2\n",
    "ENST00000448914.1\tENSG00000228985.1\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_transcript_list(input, use_name = False, use_version = True):\n",
    "    r = {}\n",
    "    for line in input:\n",
    "        if len(line) == 0 or line[0] == '#':\n",
    "            continue\n",
    "        l = line.strip().split('\\t')\n",
    "        if l[2] == 'transcript':\n",
    "            info = l[8]\n",
    "            d = {}\n",
    "            for x in info.split('; '):\n",
    "                x = x.strip()\n",
    "                p = x.find(' ')\n",
    "                if p == -1:\n",
    "                    continue\n",
    "                k = x[:p]\n",
    "                p = x.find('\"',p)\n",
    "                p2 = x.find('\"',p+1)\n",
    "                v = x[p+1:p2]\n",
    "                d[k] = v\n",
    "\n",
    "\n",
    "            if 'transcript_id' not in d or 'gene_id' not in d:\n",
    "                continue\n",
    "\n",
    "            tid = d['transcript_id']\n",
    "            gid = d['gene_id']\n",
    "            if use_version:\n",
    "                if 'transcript_version' not in d or 'gene_version' not in d:\n",
    "                    continue\n",
    "\n",
    "                tid += '.' + d['transcript_version']\n",
    "                gid += '.' + d['gene_version']\n",
    "            gname = None\n",
    "            if use_name:\n",
    "                if 'gene_name' not in d:\n",
    "                    continue\n",
    "                gname = d['gene_name']\n",
    "\n",
    "            if tid in r:\n",
    "                continue\n",
    "\n",
    "            r[tid] = (gid, gname)\n",
    "    return r\n",
    "\n",
    "\n",
    "\n",
    "def print_output(output, r, use_name = True):\n",
    "    for tid in r:\n",
    "        if use_name:\n",
    "            output.write(\"%s\\t%s\\t%s\\n\"%(tid, r[tid][0], r[tid][1]))\n",
    "        else:\n",
    "            output.write(\"%s\\t%s\\n\"%(tid, r[tid][0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created transcript_to_gene.tsv file\n"
     ]
    }
   ],
   "source": [
    "with open('./Mus_musculus.GRCm38.94.gtf') as file:\n",
    "    r = create_transcript_list(file, use_name = False, use_version = True)\n",
    "with open('transcript_to_gene.tsv', \"w+\") as output:\n",
    "    print_output(output, r, use_name = False)\n",
    "print('Created transcript_to_gene.tsv file')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run kallisto bus\n",
    "kallisto bus supports several single cell sequencing technologies, as you can see below. We'll be using 10xv3            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "List of supported single cell technologies\n",
      "\n",
      "short name       description\n",
      "----------       -----------\n",
      "10xv1            10x version 1 chemistry\n",
      "10xv2            10x version 2 chemistry\n",
      "10xv3            10x version 3 chemistry\n",
      "CELSeq           CEL-Seq\n",
      "CELSeq2          CEL-Seq version 2\n",
      "DropSeq          DropSeq\n",
      "inDrops          inDrops\n",
      "SCRBSeq          SCRB-Seq\n",
      "SureCell         SureCell for ddSEQ\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!kallisto bus --list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "[index] k-mer length: 31\n",
      "[index] number of targets: 115,270\n",
      "[index] number of k-mers: 98,989,067\n",
      "[index] number of equivalence classes: 419,171\n",
      "[quant] will process sample 1: ./neuron_1k_v3_fastqs/neuron_1k_v3_S1_L001_R1_001.fastq.gz\n",
      "                               ./neuron_1k_v3_fastqs/neuron_1k_v3_S1_L001_R2_001.fastq.gz\n",
      "[quant] will process sample 2: ./neuron_1k_v3_fastqs/neuron_1k_v3_S1_L002_R1_001.fastq.gz\n",
      "                               ./neuron_1k_v3_fastqs/neuron_1k_v3_S1_L002_R2_001.fastq.gz\n",
      "[quant] finding pseudoalignments for the reads ... done\n",
      "[quant] processed 92,902,231 reads, 58,058,974 reads pseudoaligned\n"
     ]
    }
   ],
   "source": [
    "!kallisto453 bus -i mouse_transcripts.idx  -o out_1k_mouse_brain_v3 -x 10xv3 -t 4 \\\n",
    "./neuron_1k_v3_fastqs/neuron_1k_v3_S1_L001_R1_001.fastq.gz \\\n",
    "./neuron_1k_v3_fastqs/neuron_1k_v3_S1_L001_R2_001.fastq.gz \\\n",
    "./neuron_1k_v3_fastqs/neuron_1k_v3_S1_L002_R1_001.fastq.gz \\\n",
    "./neuron_1k_v3_fastqs/neuron_1k_v3_S1_L002_R2_001.fastq.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `matrix.ec` file\n",
    "\n",
    "The `matrix.ec` is generated by kallisto and connects the equivalence class ids to sets of transcripts. The format looks like\n",
    "~~~\n",
    "0\t0\n",
    "1\t1\n",
    "2\t2\n",
    "3\t3\n",
    "4\t4\n",
    "...\n",
    "\n",
    "884398\t26558,53383,53384,69915,69931,85319,109252,125730\n",
    "884399\t7750,35941,114698,119265\n",
    "884400\t9585,70083,92571,138545,138546\n",
    "884401\t90512,90513,134202,159456\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load transcript to gene file\n",
    "tr2g = {}\n",
    "trlist = []\n",
    "with open('./transcript_to_gene.tsv') as f:\n",
    "    for line in f:\n",
    "        l = line.split()\n",
    "        tr2g[l[0]] = l[1]\n",
    "        trlist.append(l[0])\n",
    "\n",
    "genes = list(set(tr2g[t] for t in tr2g))\n",
    "\n",
    "# load equivalence classes\n",
    "ecs = {}\n",
    "with open('./out_1k_mouse_brain_v3/matrix.ec') as f:\n",
    "    for line in f:\n",
    "        l = line.split()\n",
    "        ec = int(l[0])\n",
    "        trs = [int(x) for x in l[1].split(',')]\n",
    "        ecs[ec] = trs\n",
    "        \n",
    "def ec2g(ec):\n",
    "    if ec in ecs:\n",
    "        return list(set(tr2g[trlist[t]] for t in ecs[ec]))        \n",
    "    else:\n",
    "        return []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Processing the BUS file\n",
    "\n",
    "For these notebooks we will work with the text file that `BUStools` produces, rather than the raw `BUS` file. \n",
    "To install `BUStools` see https://github.com/BUStools/bustools\n",
    "\n",
    "We discard any barcodes that don't have more 10 UMIs \n",
    "\n",
    "To produce the text file, starting with the `output.bus` file produced by kallisto, we first sort it on bustools:\n",
    "```\n",
    "bustools sort -o output.sorted output.bus\n",
    "```\n",
    "Then we convert it to txt:\n",
    "```\n",
    "bustools text -o output.sorted.txt output.sorted\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read in 58058974 number of busrecords\n",
      "All sorted\n"
     ]
    }
   ],
   "source": [
    "#sort bus file\n",
    "!bustools sort -o ./out_1k_mouse_brain_v3/output.sorted ./out_1k_mouse_brain_v3/output.bus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read in 27302856 number of busrecords\n"
     ]
    }
   ],
   "source": [
    "# convert the sorted busfile to txt\n",
    "!bustools text -o ./out_1k_mouse_brain_v3/output.sorted.txt ./out_1k_mouse_brain_v3/output.sorted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot the bus file results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "from collections import defaultdict\n",
    "\n",
    "# precompute because this is constant per ec\n",
    "ec2g = {ec:frozenset(tr2g[trlist[t]] for t in ecs[ec]) for ec in ecs}\n",
    "\n",
    "# first pass: collect gene sets\n",
    "bcu_gs = dict()\n",
    "\n",
    "with open('./out_1k_mouse_brain_v3/output.sorted.txt') as f:\n",
    "    rdr = csv.reader(f, delimiter='\\t')\n",
    "    for bar,umi,ec,_ in rdr:\n",
    "        gs = ec2g[int(ec)]\n",
    "\n",
    "        if (bar,umi) in bcu_gs:\n",
    "            bcu_gs[bar,umi].intersection_update(gs)\n",
    "        else:\n",
    "            bcu_gs[bar,umi] = set(gs)\n",
    "\n",
    "# second pass: compute gene counts\n",
    "cell_gene = defaultdict(lambda: defaultdict(float))\n",
    "\n",
    "for (bar,umi),gs in bcu_gs.items():\n",
    "    for g in gs:\n",
    "        cell_gene[bar][g] += 1.0 / len(gs)\n",
    "\n",
    "# finally: filter out barcodes below threshold\n",
    "cell_gene = {bar:cell_gene[bar] for bar in cell_gene\n",
    "              if sum(cell_gene[bar].values()) >= 10.0}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "barcode_hist = collections.defaultdict(int)\n",
    "for barcode in cell_gene:\n",
    "    cg = cell_gene[barcode]\n",
    "    s = len([cg[g] for g in cg])\n",
    "    barcode_hist[barcode] += s\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download the 10x whitelist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2018-12-03 17:59:45--  https://github.com/10XGenomics/supernova/blob/master/tenkit/lib/python/tenkit/barcodes/737K-august-2016.txt\n",
      "Resolving github.com (github.com)... 192.30.255.112, 192.30.255.113\n",
      "Connecting to github.com (github.com)|192.30.255.112|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: unspecified [text/html]\n",
      "Saving to: ‘737K-august-2016.txt’\n",
      "\n",
      "    [ <=>                                   ] 56,530      --.-K/s   in 0.06s   \n",
      "\n",
      "2018-12-03 17:59:50 (991 KB/s) - ‘737K-august-2016.txt’ saved [56530]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget https://github.com/10XGenomics/supernova/blob/master/tenkit/lib/python/tenkit/barcodes/737K-august-2016.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "whitelist = set(x.strip() for x in open('737K-august-2016.txt'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1364\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAD9dJREFUeJzt3WuM3NV5x/HvU3NLSRTjsLVcjLtGQalQpQBdIRBR1UIgJCCgEkIg1DgNlaVepKSplJryKlJfQFslIWpVYgVatyLElEBBkJZSh6iqVDmxy/3iYIhpsAAbCiFppTZOnr6Ys2Zwd5nZ3ZmdmWe/H2m0/9vsnLNn9rdnnv9/ZiMzkSRNvp8ZdQMkSYNhoEtSEQa6JBVhoEtSEQa6JBVhoEtSEQa6JBVhoEtSEQa6JBVx1HI+2IknnpjT09PL+ZCSNPF27979amZO9TpuWQN9enqaXbt2LedDStLEi4gX+jnOkoskFWGgS1IRBrokFWGgS1IRBrokFWGgS1IRBrokFWGgS1IRBrokFbGs7xRdiukt9x9e3nfDxSNsiSSNJ2foklSEgS5JRRjoklSEgS5JRRjoklSEgS5JRRjoklREX9ehR8Q+4IfAT4BDmTkTEWuA7cA0sA+4MjNfH04zJUm9LGSG/muZeXpmzrT1LcCOzDwV2NHWJUkjspSSy2XAtra8Dbh86c2RJC1Wv4GewD9FxO6I2Ny2rc3Ml9ryy8Daue4YEZsjYldE7Dp48OASmytJmk+/n+XyoczcHxE/BzwYEc9078zMjIic646ZuRXYCjAzMzPnMZKkpetrhp6Z+9vXA8DdwFnAKxGxDqB9PTCsRkqSeusZ6BFxfES8Z3YZuBB4ArgX2NQO2wTcM6xGSpJ666fksha4OyJmj/9qZv5jRHwHuCMirgVeAK4cXjMlSb30DPTMfB744BzbXwPOH0ajJEkL5ztFJakIA12SijDQJakIA12SijDQJakIA12SijDQJakIA12SijDQJakIA12SijDQJakIA12SijDQJakIA12SijDQJakIA12SijDQJakIA12Siujnf4qOnekt9x9e3nfDxSNsiSSND2foklSEgS5JRRjoklSEgS5JRRjoklSEgS5JRRjoklSEgS5JRRjoklSEgS5JRRjoklRE34EeEasi4uGIuK+tb4yInRGxNyK2R8Qxw2umJKmXhczQPwU83bV+I/CFzHw/8Dpw7SAbJklamL4CPSLWAxcDX2nrAZwH3NkO2QZcPowGSpL60+8M/YvAZ4GftvX3AW9k5qG2/iJw0oDbJklagJ6BHhGXAAcyc/diHiAiNkfErojYdfDgwcV8C0lSH/qZoZ8LXBoR+4Cv0Sm13ASsjojZf5CxHtg/150zc2tmzmTmzNTU1ACaLEmaS89Az8zrMnN9Zk4DVwHfzMxrgIeAK9phm4B7htZKSVJPS7kO/Q+Bz0TEXjo19VsG0yRJ0mIs6H+KZua3gG+15eeBswbfJEnSYvhOUUkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCJ6BnpEHBcR346IRyPiyYj4XNu+MSJ2RsTeiNgeEccMv7mSpPn0M0P/H+C8zPwgcDpwUUScDdwIfCEz3w+8Dlw7vGZKknrpGejZ8aO2enS7JXAecGfbvg24fCgtlCT1pa8aekSsiohHgAPAg8BzwBuZeagd8iJw0nCaKEnqx1H9HJSZPwFOj4jVwN3AL/b7ABGxGdgMsGHDhsW0UdIYm95y/+HlfTdcPMKWaEFXuWTmG8BDwDnA6oiY/YOwHtg/z322ZuZMZs5MTU0tqbGSpPn1c5XLVJuZExHvAi4AnqYT7Fe0wzYB9wyrkZKk3vopuawDtkXEKjp/AO7IzPsi4ingaxHxx8DDwC1DbKckqYeegZ6ZjwFnzLH9eeCsYTRKkrRwvlNUkoow0CWpCANdkoow0CWpCANdkoow0CWpCANdkoow0CWpCANdkoow0CWpiL4+PlfS8vCjaLUUztAlqQgDXZKKMNAlqYiJr6Fbc1RVPre1UM7QJakIA12SijDQJamIia+hS5NoGPVxa+5yhi5JRRjoklSEgS5JRVhDl7D+rBqcoUtSEQa6JBVhoEtSEdbQpQnWXfuXnKFLUhEGuiQVYaBLUhE9a+gRcTLwN8BaIIGtmXlTRKwBtgPTwD7gysx8fXhN7c1riTXO5qt3L2cdfBx/R8axTZOqnxn6IeAPMvM04GzgdyPiNGALsCMzTwV2tHVJ0oj0DPTMfCkz/70t/xB4GjgJuAzY1g7bBlw+rEZKknpbUA09IqaBM4CdwNrMfKnteplOSUaSNCJ9X4ceEe8Gvg58OjPfjIjD+zIzIyLnud9mYDPAhg0bltZaSQOruS+ldu317+Oprxl6RBxNJ8xvy8y72uZXImJd278OODDXfTNza2bOZObM1NTUINosSZpDz0CPzlT8FuDpzPx81657gU1teRNwz+CbJ0nqVz8ll3OB3wAej4hH2rY/Am4A7oiIa4EXgCuH00RJUj96Bnpm/isQ8+w+f7DNkSZDv/XnYdS7h2G+7z+s68KXsz8r6dp23ykqSUUY6JJUxIr4+NyV+vJLg7XSL9Vb6f2fBM7QJakIA12SijDQJamIFVFD72Y9XcM0jnXmcWjTOLRhJXCGLklFGOiSVISBLklFrLgaulaGfs6VVK7rrpS+Deo8WJVza87QJakIA12SijDQJakIa+gqYxz+NVs1lWvxFTlDl6QiDHRJKsJAl6QiytbQrf1pLsN4Xvhce8s4/ixW0jkRZ+iSVISBLklFGOiSVETZGvooraSa3aj1U7Mdx7quFsb3GPTHGbokFWGgS1IRBrokFbGia+jz1eUq1ta0ONbfx9tyjs9C6++jqNc7Q5ekIgx0SSrCQJekInrW0CPiVuAS4EBm/lLbtgbYDkwD+4ArM/P14TVzci30f1uu9Pr9Sv9foFo+S3kejevvbD8z9L8GLjpi2xZgR2aeCuxo65KkEeoZ6Jn5L8B/HrH5MmBbW94GXD7gdkmSFmixNfS1mflSW34ZWDug9kiSFmnJ16FnZkZEzrc/IjYDmwE2bNiw1IeTVJjnR5ZmsTP0VyJiHUD7emC+AzNza2bOZObM1NTUIh9OktTLYgP9XmBTW94E3DOY5kiSFqtnoEfE7cC/AR+IiBcj4lrgBuCCiHgW+HBblySNUM8aemZePc+u8wfclrExrteYzhr39kkr1ajPAfhOUUkqwkCXpCIMdEkqYkV/Hrp6m69e/061wvmO6/f+khbHGbokFWGgS1IRllx6OLI0sJTLBMf5csNBtm2+coplFk2CSX6eOkOXpCIMdEkqwkCXpCKsoY+ZYV/mN0n1wXE+5yDNGqffKWfoklSEgS5JRRjoklSENfRCKr/NvkIfpGFzhi5JRRjoklSEgS5JRVhDX4JhXyc9qrqxn8UiTSZn6JJUhIEuSUUY6JJUhDX0BRq3+rJ1bUmznKFLUhEGuiQVYaBLUhHW0Edk2LVva+vSyuMMXZKKMNAlqQgDXZKKWFKgR8RFEbEnIvZGxJZBNUqStHCLDvSIWAX8BfBR4DTg6og4bVANkyQtzFJm6GcBezPz+cz8X+BrwGWDaZYkaaGWEugnAd/vWn+xbZMkjcDQr0OPiM3A5rb6o4jYM+zHXEYnAq+OuhFDYL8mR8U+QbF+xY3A0vr0C/0ctJRA3w+c3LW+vm17m8zcCmxdwuOMrYjYlZkzo27HoNmvyVGxT1CzX8vRp6WUXL4DnBoRGyPiGOAq4N7BNEuStFCLnqFn5qGI+D3gAWAVcGtmPjmwlkmSFmRJNfTM/AbwjQG1ZRKVLCVhvyZJxT5BzX4NvU+RmcN+DEnSMvCt/5JUhIF+hIg4OSIeioinIuLJiPhU274mIh6MiGfb1xPa9oiIL7WPP3gsIs7s+l6b2vHPRsSmUfWpqz2rIuLhiLivrW+MiJ2t7dvbyW0i4ti2vrftn+76Hte17Xsi4iOj6clbImJ1RNwZEc9ExNMRcc6kj1VE/H577j0REbdHxHGTOFYRcWtEHIiIJ7q2DWxsIuKXI+Lxdp8vRUSMsF9/2p6Dj0XE3RGxumvfnOMQ83x0ynxj3ZfM9NZ1A9YBZ7bl9wDfpfPRBn8CbGnbtwA3tuWPAf8ABHA2sLNtXwM8376e0JZPGHHfPgN8Fbivrd8BXNWWbwZ+uy3/DnBzW74K2N6WTwMeBY4FNgLPAatG3KdtwG+15WOA1ZM8VnTenPc94F1dY/SJSRwr4FeAM4EnurYNbGyAb7djo933oyPs14XAUW35xq5+zTkO7fYccEp73j4KnNY15v9vrPtq2yietJN0A+4BLgD2AOvatnXAnrb8ZeDqruP3tP1XA1/u2v6240bQj/XADuA84L72S/Bq15PwHOCBtvwAcE5bPqodF8B1wHVd3/PwcSPq03tb+MUR2yd2rHjrHdhr2s/+PuAjkzpWwPQRwTeQsWn7nuna/rbjlrtfR+z7deC2tjznOHSPYfdx7/R72c/Nkss7aC9fzwB2Amsz86W262VgbVue7yMQxu2jEb4IfBb4aVt/H/BGZh5q693tO9z2tv8H7fhx69NG4CDwV62U9JWIOJ4JHqvM3A/8GfAfwEt0fva7mfyxmjWosTmpLR+5fRx8ks4rBlh4v97p97InA30eEfFu4OvApzPzze592fnTOTGXB0XEJcCBzNw96rYM2FF0Xvr+ZWaeAfwXnZfxh03gWJ1A50PuNgI/DxwPXDTSRg3JpI1NPyLieuAQcNsoHt9An0NEHE0nzG/LzLva5lciYl3bvw440LbP9xEIfX00wjI5F7g0IvbR+VTM84CbgNURMftehO72HW572/9e4DXGq0/Qmb28mJk72/qddAJ+ksfqw8D3MvNgZv4YuIvO+E36WM0a1Njsb8tHbh+ZiPgEcAlwTftjBQvv12vMP9Y9GehHaGfKbwGezszPd+26F5g9w76JTm19dvvH21n6s4EftJeUDwAXRsQJbdZ1Ydu27DLzusxcn5nTdE6cfTMzrwEeAq5ohx3Zp9m+XtGOz7b9qnZlxUbgVDonpkYiM18Gvh8RH2ibzgeeYoLHik6p5eyI+Nn2XJzt00SPVZeBjE3b92ZEnN1+Th/v+l7LLiIuolPSvDQz/7tr13zjMOdHp7Sxm2+se1vukyTjfgM+ROdl4GPAI+32MTq1rR3As8A/A2va8UHnH308BzwOzHR9r08Ce9vtN0fdt9amX+Wtq1xOaU+uvcDfAce27ce19b1t/yld97++9XUPy3RVQY/+nA7sauP193SuhJjosQI+BzwDPAH8LZ0rJCZurIDb6ZwH+DGdV1PXDnJsgJn2M3oO+HOOODm+zP3aS6cmPpsZN/cah5Yr3237ru/aPudY93PznaKSVIQlF0kqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCIMdEkqwkCXpCL+D/pkqcXHY71DAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bcv = [x for b,x in barcode_hist.items() if x > 600 and x < 12000]\n",
    "_ = plt.hist(bcv,bins=100)\n",
    "print(len(bcv))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}