tutorials/preprocess.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Preparing a dataset for Basenji training involves a series of design choices.\n",
    "\n",
    "The input you bring to the pipeline is:\n",
    "* BigWig coverage tracks\n",
    "* Genome FASTA file\n",
    "\n",
    "First, make sure you have an hg19 FASTA file visible. If you have it already, put a symbolic link into the data directory. Otherwise, I have a machine learning friendly simplified version you can download in the next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100 2943M  100 2943M    0     0  28.9M      0  0:01:41  0:01:41 --:--:-- 30.3M\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100   729  100   729    0     0   3973      0 --:--:-- --:--:-- --:--:--  3983\n"
     ]
    }
   ],
   "source": [
    "import os, subprocess\n",
    "\n",
    "if not os.path.isfile('data/hg19.ml.fa'):\n",
    "    subprocess.call('curl -o data/hg19.ml.fa https://storage.googleapis.com/basenji_tutorial_data/hg19.ml.fa', shell=True)\n",
    "    subprocess.call('curl -o data/hg19.ml.fa.fai https://storage.googleapis.com/basenji_tutorial_data/hg19.ml.fa.fai', shell=True)                "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's grab a few CAGE datasets from FANTOM5 related to heart biology.\n",
    "\n",
    "These data were processed by\n",
    "1. Aligning with Bowtie2 with very sensitive alignment parameters.\n",
    "2. Distributing multi-mapping reads and estimating genomic coverage with [bam_cov.py](https://github.com/calico/basenji/blob/master/bin/bam_cov.py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100 1563M  100 1563M    0     0  31.9M      0  0:00:48  0:00:48 --:--:-- 27.3M\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100 1236M  100 1236M    0     0  36.2M      0  0:00:34  0:00:34 --:--:-- 27.5M\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100 3484M  100 3484M    0     0  30.1M      0  0:01:55  0:01:55 --:--:-- 26.6M\n"
     ]
    }
   ],
   "source": [
    "if not os.path.isfile('data/CNhs11760.bw'):\n",
    "    subprocess.call('curl -o data/CNhs11760.bw https://storage.googleapis.com/basenji_tutorial_data/CNhs11760.bw', shell=True)\n",
    "    subprocess.call('curl -o data/CNhs12843.bw https://storage.googleapis.com/basenji_tutorial_data/CNhs12843.bw', shell=True)\n",
    "    subprocess.call('curl -o data/CNhs12856.bw https://storage.googleapis.com/basenji_tutorial_data/CNhs12856.bw', shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we'll write out these BigWig files and labels to a samples table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lines = [['index','identifier','file','clip','sum_stat','description']]\n",
    "lines.append(['0', 'CNhs11760', 'data/CNhs11760.bw', '384', 'sum', 'aorta'])\n",
    "lines.append(['1', 'CNhs12843', 'data/CNhs12843.bw', '384', 'sum', 'artery'])\n",
    "lines.append(['2', 'CNhs12856', 'data/CNhs12856.bw', '384', 'sum', 'pulmonic_valve'])\n",
    "\n",
    "samples_out = open('data/heart_wigs.txt', 'w')\n",
    "for line in lines:\n",
    "    print('\\t'.join(line), file=samples_out)\n",
    "samples_out.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we want to choose genomic sequences to form batches for stochastic gradient descent, divide them into training/validation/test sets, and construct TFRecords to provide to downstream programs.\n",
    "\n",
    "The script [basenji_data.py](https://github.com/calico/basenji/blob/master/bin/basenji_data.py) implements this procedure.\n",
    "\n",
    "The most relevant options here are:\n",
    "\n",
    "| Option/Argument | Value | Note |\n",
    "|:---|:---|:---|\n",
    "| -s | 0.1 | Down-sample the genome to 10% to speed things up here. |\n",
    "| -g | data/unmap_macro.bed | Dodge large-scale unmappable regions like assembly gaps. |\n",
    "| -l | 131072 | Sequence length. |\n",
    "| --local | True | Run locally, as opposed to on my SLURM scheduler. |\n",
    "| -o | data/heart_l131k | Output directory |\n",
    "| -p | 8 | Uses multiple concourrent processes to read/write. |\n",
    "| -t | .1 | Hold out 10% sequences for testing. |\n",
    "| -v | .1 | Hold out 10% sequences for validation. |\n",
    "| -w | 128 | Pool the nucleotide-resolution values to 128 bp bins. |\n",
    "| fasta_file| data/hg19.ml.fa | FASTA file to extract sequences from. |\n",
    "| targets_file | data/heart_wigs.txt | Target samples table with BigWig paths. |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "stride_train 1 converted to 131072.000000\n",
      "stride_test 1 converted to 131072.000000\n",
      "Remove output directory data/heart_l131k or use --restart option.\n"
     ]
    }
   ],
   "source": [
    "! basenji_data.py -s .1 -g data/unmap_macro.bed -l 131072 --local -o data/heart_l131k -p 8 -t .1 -v .1 -w 128 data/hg19.ml.fa data/heart_wigs.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, data/heart_l131k contains relevant data for training.\n",
    "\n",
    "*contigs.bed* contains the original large contiguous regions from which training sequences were taken (possibly strided).\n",
    "*sequences.bed*\n",
    "\n",
    "contains the train/valid/test sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    179 test\n",
      "   1499 train\n",
      "    180 valid\n"
     ]
    }
   ],
   "source": [
    "! cut -f4 data/heart_l131k/sequences.bed | sort | uniq -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr2\t140425791\t140556863\ttrain\n",
      "chr16\t27143973\t27275045\ttrain\n",
      "chr14\t72972403\t73103475\ttrain\n"
     ]
    }
   ],
   "source": [
    "! head -n3 data/heart_l131k/sequences.bed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr11\t116271608\t116402680\tvalid\n",
      "chr2\t26239628\t26370700\tvalid\n",
      "chr6\t102241368\t102372440\tvalid\n"
     ]
    }
   ],
   "source": [
    "! grep valid data/heart_l131k/sequences.bed | head -n3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr14\t76007464\t76138536\ttest\n",
      "chr13\t99391330\t99522402\ttest\n",
      "chr15\t52280212\t52411284\ttest\n"
     ]
    }
   ],
   "source": [
    "! grep test data/heart_l131k/sequences.bed | head -n3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r----- 1 drk Domain Users 10610356 Feb 12  2021 data/heart_l131k/tfrecords/test-0.tfr\n",
      "-rw-r----- 1 drk Domain Users 15168641 Feb 12  2021 data/heart_l131k/tfrecords/train-0.tfr\n",
      "-rw-r----- 1 drk Domain Users 15153054 Feb 12  2021 data/heart_l131k/tfrecords/train-1.tfr\n",
      "-rw-r----- 1 drk Domain Users 15169161 Feb 12  2021 data/heart_l131k/tfrecords/train-2.tfr\n",
      "-rw-r----- 1 drk Domain Users 15182385 Feb 12  2021 data/heart_l131k/tfrecords/train-3.tfr\n",
      "-rw-r----- 1 drk Domain Users 15151675 Feb 12  2021 data/heart_l131k/tfrecords/train-4.tfr\n",
      "-rw-r----- 1 drk Domain Users 12977773 Feb 12  2021 data/heart_l131k/tfrecords/train-5.tfr\n",
      "-rw-r----- 1 drk Domain Users 10662202 Feb 12  2021 data/heart_l131k/tfrecords/valid-0.tfr\n"
     ]
    }
   ],
   "source": [
    "! ls -l data/heart_l131k/tfrecords/*.tfr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}