OP-WORKFLOW-CAGEscan-short-reads-v2.0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OP-WORKFLOW-CAGEscan-short-reads-v2.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This document is an example of how to process a C1 CAGE library with a Jupyter notebook from raw reads to single molecule count. All the steps are described in the [tutorial](https://github.com/Population-Transcriptomics/C1-CAGE-preview/blob/master/tutorial.md) section of this repository. In the following section we assume that:\n",
    "- you are familiar with python programming and Jupyter (IPython notebook)\n",
    "- the softwares mentioned in the [prerequisite](https://github.com/Population-Transcriptomics/C1-CAGE-preview/blob/master/prerequisite.md) section are all installed and referenced in your PATH\n",
    "- you've already downloaded the available [C1 CAGE library](http://dx.doi.org/10.5281/zenodo.48478) and extracted it in your current directory\n",
    "- The reference genome is already indexed with BWA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our hands this notebook worked without trouble on a machine running Debian GNU/Linux 8. We noticed that the behavior of tagdust2 in single-end mode was different on Mac OSX. In short, the order of the reads1 is changed after extraction on Mac OSX which is a problem because syncpairs expect the order of reads1 and reads2 to be the same. One way to overcome this issue is sort reads1 and reads2 separately after the extraction then syncpairs will work properly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import subprocess, os, csv, signal, pysam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "remove_extension = lambda x: x.split('.')[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Declare the function that deals with inputs and outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_args(read1, read2, ref_genome, output_folders):\n",
    "    '''Set the input and output path for a given pair of reads'''\n",
    "    r1_shortname = remove_extension(os.path.basename(read1))\n",
    "\n",
    "    args = {  \n",
    "        'r1_input': read1,\n",
    "        'r2_input': read2,\n",
    "        'ref_genome': ref_genome,\n",
    "    }\n",
    "    \n",
    "    output_paths = {folder: os.path.join('output', folder, r1_shortname) for folder in output_folders}\n",
    "    \n",
    "    return dict(args, **output_paths)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the required softwares are not in the PATH you can manually set their location here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tagdust2_path = 'tagdust'\n",
    "bwa_path = 'bwa'\n",
    "samtools_path = 'samtools'\n",
    "paired_bam_to_bed12_path = 'pairedBamToBed12'\n",
    "umicountFP_path = 'umicountFP'\n",
    "syncpairs_path = 'syncpairs'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Path to the reference genome you want to align your reads against"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ref_genome = './GRCh38.fa'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "softwares = {    \n",
    "    'bwa': bwa_path,\n",
    "    'tagdust': tagdust2_path,\n",
    "    'syncpairs': syncpairs_path,\n",
    "    'samtools': samtools_path,\n",
    "    'pairedBamToBed12': paired_bam_to_bed12_path,\n",
    "    'umicountFP': umicountFP_path}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The name of the output folders for each command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "output_folders = [ 'tagdust_r1', 'unzip_r2'                    # Demultiplexed R1, unziped R2\n",
    "                 , 'extracted_r1', 'extracted_r2'              # Synced R1 and R2\n",
    "                 , 'cleaned_reads', 'cleaned_r1', 'cleaned_r2' # rRNA reads removed\n",
    "                 , 'r1_sai', 'r2_sai', 'sampe'                 # Intermediate files from BWA\n",
    "                 , 'genome_mapped', 'properly_paired'          # Final output in BAM format\n",
    "                 , 'cagescan_pairs', 'cagescan_fragments'      # Final output in BED12 format\n",
    "                 ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for folder in output_folders:\n",
    "    os.makedirs(os.path.join('output', folder))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The actual command to run. See the [tutorial](https://github.com/Population-Transcriptomics/C1-CAGE-preview/blob/master/tutorial.md) section for more details about each command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cmds = [\n",
    "    \n",
    "    '{tagdust} -t8 -o {tagdust_r1} -1 F:NNNNNNNN -2 S:TATAGGG -3 R:N {r1_input}',\n",
    "    \n",
    "    'gunzip -c {r2_input} > {unzip_r2}.fq',\n",
    "        \n",
    "    '{syncpairs} {tagdust_r1}.fq {unzip_r2}.fq {extracted_r1}.fq {extracted_r2}.fq',\n",
    "    \n",
    "    '{tagdust} -arch SimpleArchitecture.txt -ref ercc_and_human_rRNA_and_tagdust.fa -o {cleaned_reads} {extracted_r1}.fq {extracted_r2}.fq',\n",
    "    \n",
    "    'cp {cleaned_reads}_READ1.fq {cleaned_r1}.fq',\n",
    "    \n",
    "    'cp {cleaned_reads}_READ2.fq {cleaned_r2}.fq',\n",
    "    \n",
    "    '{bwa} aln -t8 {ref_genome} {cleaned_r1}.fq > {r1_sai}.sai',\n",
    "    \n",
    "    '{bwa} aln -t8 {ref_genome} {cleaned_r2}.fq > {r2_sai}.sai',\n",
    "    \n",
    "    '{bwa} sampe -a 2000000 {ref_genome} {r1_sai}.sai {r2_sai}.sai {cleaned_r1}.fq {cleaned_r2}.fq > {sampe}.sam',\n",
    "    \n",
    "    '{samtools} view -uSo - {sampe}.sam | {samtools} sort - -o {genome_mapped}.bam',\n",
    "    \n",
    "    '{samtools} view -f 0x0002 -F 0x0100 -uo - {genome_mapped}.bam | {samtools} sort -n - -o {properly_paired}.bam',\n",
    "    \n",
    "    '{pairedBamToBed12} -i {properly_paired}.bam > {cagescan_pairs}.bed',\n",
    "    \n",
    "    '{umicountFP} -f {cagescan_pairs}.bed > {cagescan_fragments}.bed'\n",
    "    \n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the reads. Here we assume that the reads are in the current directory, in a folder named following the MiSeq run id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fastq_path = './150519_M00528_0125_000000000-ACUAB-links/'\n",
    "root, folders, files = os.walk(fastq_path).__next__()\n",
    "\n",
    "files = [f for f in files if not f.startswith('.')] #remove hidden files if there exist\n",
    "reads1 = sorted([os.path.join(root, f) for f in files if 'R1' in f])\n",
    "reads2 = sorted([os.path.join(root, f) for f in files if 'R2' in f])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the commands for all the pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for read1, read2 in zip(reads1, reads2):\n",
    "    args = get_args(read1, read2, ref_genome, output_folders)\n",
    "    args = dict(args, **softwares)\n",
    "    \n",
    "    for cmd in cmds:\n",
    "#        print(cmd.format(**args))\n",
    "        subprocess.call(cmd.format(**args), preexec_fn=lambda: signal.signal(signal.SIGPIPE, signal.SIG_DFL), shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the level1 file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "root, folders, files = os.walk('./output/genome_mapped/').__next__()\n",
    "files = [os.path.join(root, f) for f in files if f.endswith('bam')]\n",
    "level1 = 'level1.py -o output/mylevel1file.l1.osc.gz -f 0x0042 -F 0x0104 --fingerprint {files}'.format(files=' '.join(files))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "subprocess.call(level1, shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate logs (triplet)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we generate four summary files that will be used for [QC](https://github.com/Population-Transcriptomics/C1-CAGE-preview/blob/master/QC.md) and place them in the 'output' directory. \n",
    "\n",
    "1.  mapped.log: The number of mapped reads per cell\n",
    "2.  extracted.log: The number of remaining reads after filtering for ribosomal DNA and unreadable UMIs\n",
    "3.  filtered.log: The detailed number of ribosomal DNA extracted per cell\n",
    "4.  transcript_count.log: The exact number of unique transcript per cell\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_cmd       = \"grep 'total input reads'      {tagdust_r1}_logfile.txt    | cut -f 2\"\n",
    "extracted_cmd   = \"grep 'successfully extracted' {tagdust_r1}_logfile.txt    | cut -f 2\"\n",
    "cleaned_cmd     = \"grep 'successfully extracted' {cleaned_reads}_logfile.txt | cut -f 2\"\n",
    "rdna_cmd        = \"grep 'ribosomal'              {cleaned_reads}_logfile.txt | cut -f 2\"\n",
    "spikes_cmd      = \"grep 'ERCC'                   {cleaned_reads}_logfile.txt | cut -f 2 | paste -s -d+ - | bc\"\n",
    "mapped_cmd      = \"{samtools} view -u -f 0x40 {genome_mapped}.bam | {samtools} flagstat - | grep 'mapped (' | cut -f 1 -d ' '\"\n",
    "properpairs_cmd = \"{samtools} view -u -f 0x42 {genome_mapped}.bam | {samtools} flagstat - | grep 'mapped (' | cut -f 1 -d ' '\"\n",
    "counts_cmd      = \"wc -l {cagescan_fragments}.bed | cut -f 1 -d ' '\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#remove _R1 from the file's name\n",
    "custom_rename = lambda x: x.replace('_R1', '')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total, cleaned, extracted, rdna, spikes, mapped, properpairs, counts = ([], [], [], [], [], [], [], [])\n",
    "\n",
    "def run_qc(command, dest, keyword):\n",
    "    output = subprocess.check_output(command.format(**args), shell=True).strip().decode()\n",
    "    dest.append([keyword, custom_rename(r1_shortname), output])\n",
    "\n",
    "for read1 in reads1:\n",
    "    r1_shortname = remove_extension(os.path.basename(read1))\n",
    "    \n",
    "    args = {'tagdust_r1':         os.path.join('output', 'tagdust_r1',         r1_shortname),\n",
    "            'genome_mapped':      os.path.join('output', 'genome_mapped',      r1_shortname),\n",
    "            'cagescan_fragments': os.path.join('output', 'cagescan_fragments', r1_shortname),\n",
    "            'cleaned_reads':      os.path.join('output', 'cleaned_reads',      r1_shortname),\n",
    "            'samtools': samtools_path}\n",
    "    \n",
    "    run_qc(total_cmd,       total,       'total')\n",
    "    run_qc(extracted_cmd,   extracted,   'extracted')\n",
    "    run_qc(cleaned_cmd,     cleaned,     'cleaned')\n",
    "    run_qc(rdna_cmd,        rdna,        'rdna')\n",
    "    run_qc(spikes_cmd,      spikes,      'spikes')\n",
    "    run_qc(mapped_cmd,      mapped,      'mapped')\n",
    "    run_qc(properpairs_cmd, properpairs, 'properpairs')\n",
    "    run_qc(counts_cmd,      counts,      'counts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def write_logs(file_dest, results):\n",
    "    with open(file_dest, 'w') as handler:\n",
    "        writer = csv.writer(handler, delimiter='\\t')\n",
    "        writer.writerows(results)\n",
    "\n",
    "write_logs('output/total.log',       total)\n",
    "write_logs('output/cleaned.log',     cleaned)\n",
    "write_logs('output/extracted.log',   extracted)\n",
    "write_logs('output/rdna.log',        rdna)\n",
    "write_logs('output/spikes.log',      spikes)\n",
    "write_logs('output/mapped.log',      mapped)\n",
    "write_logs('output/properpairs.log', properpairs)\n",
    "write_logs('output/counts.log',      counts)"
   ]
  }
 ],
 "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.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}