ATACseq_rep2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ATACseq pipeline Human rep2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This jupyter-notebook automates the ATACseq pre-processing pipeline as well as the basic analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following section we assume that:\n",
    "- The reference genome has to be already indexed with bowtie\n",
    "- The tutorials are introduced with the example file mentioned in the tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Description\n",
    "\n",
    "In this pipeline we used data from GEO series GSE65360. This dataset is baed on the paper Buenrostro JD, Wu B, Litzenburger UM, Ruff D et al. Single-cell chromatin accessibility reveals principles of regulatory variation. Nature 2015 Jul 23;523(7561):486-90. PMID: [26083756](https://www.ncbi.nlm.nih.gov/pubmed/26083756)\n",
    "\n",
    "They profiled single cell epigenomes, scATAC-seq assay, across 8 cell types and 4 targeted cell manipulations. The complete data set contains a total of 1,632 assayed wells.\n",
    "\n",
    "For their experiment they used four different platforms\n",
    "- Illumina MiSeq (Homo sapiens)\n",
    "- Illumina HiSeq 2500 (Homo sapiens)\n",
    "- Illumina NextSeq 500 (Homo sapiens)\n",
    "- Illumina NextSeq 500 (Mus musculus)\n",
    "\n",
    "The fastq files have varying length of sequencing reads. An overview of different cell types they used are as follows\n",
    "\n",
    "__Human__\n",
    "- singles-H1ESC: The H1 human embryonic stem cells were obtained from WiCell and grown in chemically defined mTeSR1 medium (STEMCELL Technologies) on Matrigel-coated plates.\n",
    "- singles-GM: GM12878 lymphoblastoid cells (a kind gift from the Snyder lab) were grown in RPMI 1640 with 2 mmol/L L-glutamine, 15% FBS and 1% Pen/Strep\n",
    "- singles-GM-TNFa6h: GM12878 lymphoblastoid cells (a kind gift from the Snyder lab) were grown in RPMI 1640 with 2 mmol/L L-glutamine, 15% FBS and 1% Pen/Strep. GM12878 cells were treated with 25 ng/mL rhTNFα (eBiosciences) for 6 h\n",
    "- singles-K562: K562 (ATCC) chronic myeloid leukemia cells were maintained in Iscove’s modified Dulbecco’s medium (IMDM) containing 10% FBS (HyClone, Thermo Scientific) and 1% Penicillin Streptomycin (Pen/Strep).\n",
    "- singles-K562-CDKi: K562 (ATCC) chronic myeloid leukemia cells were maintained in Iscove’s modified Dulbecco’s medium (IMDM) containing 10% FBS (HyClone, Thermo Scientific) and 1% Penicillin Streptomycin (Pen/Strep). K562 cells were treated with 1 µM CDK4/6 inhibitor (PD 0332991, Pfizer) for 24 h. [Acute myeloid leukemia (AML), also known as acute myelogenous leukemia or acute nonlymphocytic leukemia (ANLL), is a cancer of the myeloid line of blood cells, characterized by the rapid growth of abnormal white blood cells that accumulate in the bone marrow and interfere with the production of normal blood cells. AML is the most common acute leukemia affecting adults, and its incidence increases with age. Although AML is a relatively rare disease, accounting for roughly 1.2% of cancer deaths in the United States,[1] its incidence is expected to increase as the population ages]\n",
    "- singles-K562-Imat1hr: K562 (ATCC) chronic myeloid leukemia cells were maintained in Iscove’s modified Dulbecco’s medium (IMDM) containing 10% FBS (HyClone, Thermo Scientific) and 1% Penicillin Streptomycin (Pen/Strep). K562 cells were treated with 1 µM Imatinib (Gleevec, Novartis) for 1 h.\n",
    "- singles-K562-JNKi: K562 (ATCC) chronic myeloid leukemia cells were maintained in Iscove’s modified Dulbecco’s medium (IMDM) containing 10% FBS (HyClone, Thermo Scientific) and 1% Penicillin Streptomycin (Pen/Strep). K562 cells were treated with 10 µM JNK inhibitor VIII (CAS 894804-07-0, Calbiochem) for 24 h\n",
    "- singles-TF1: The erythroleukemia cell line TF-1 (kind gift from the Majeti lab, Stanford) was maintained in RPMI 1640 with 10% FBS, 1% Pen/Strep and 2 ng/mL rhGM-CSF (Peprotech).[Acute erythroid leukemia or Di Guglielmo syndrome is a rare form of acute myeloid leukemia (less than 5% of AML cases[1]) where the myeloproliferation is of erythroblastic precursors. It is defined as type \"M6\"]\n",
    "- singles-HL60: The promyelocytic leukemia cells HL-60 were grown in IMDM containing 20% FBS and 1% Pen/Strep. [Acute promyelocytic leukemia (APML, APL) is the M3 subtype of acute myelogenous leukemia (AML), a cancer of the white blood cells.[1] In APL, there is an abnormal accumulation of immature granulocytes called promyelocytes. The disease is characterized by a chromosomal translocation involving the retinoic acid receptor alpha (RARα or RARA) gene and is distinguished from other forms of AML by its responsiveness to all-trans retinoic acid (ATRA; also known as tretinoin) therapy]\n",
    "- singles-BJ: Human BJ fibroblasts (ATCC) were maintained in Eagle’s minimum essential medium (EMEM, ATCC) supplemented with 10% FBS and 1% Pen/Strep.\n",
    "\n",
    "__Mouse__\n",
    "- singles-EML: The mouse EML hematopoietic cells were grown in IMDM containing 20% horse serum, 2 mM L-Glutamine, 1% Pen/Strep and 100 ng/mL recombinant SCF (Peprotech). [Hematopoietic stem cells (HSCs) or hemocytoblasts are the stem cells that give rise to all the other blood cells through the process of haematopoiesis]\n",
    "- singles-mESC: V6.5 mouse embryonic stem cells were grown in Glasgow Minimum Essential Medium (GMEM) supplemented with 10% FBS, 2 mM glutamine, 1% Pen/Strep, 1 mM sodium pyruvate, 1000 units/ml Leukemia Inhibitory Factor (LIF, Millipore), 1X Minimum Essential Medium NonEssential Amino Acids (MEM NEAA, Invitrogen) and 50 uM β-Mercaptoethanol. Cells were maintained on gelatin-coated dishes without feeders\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data processing summary\n",
    "\n",
    "`trimmomatic` was supposed to be used to trim the adapter sequences. However using it causes macs2 not to find peaks. So we commented the line. Paired-end reads were aligned to `hg19` or `mm10` using `BOWTIE2`. We used the `-X2000` parameter in `BOWTIE2` allowing fragments of up to 2 Kb to align. This considers that the read pairs within 2000 bp are concordant mapping. Duplicates were removed using `PICARD` tools `Markduplicate`. Reads are then filtered for alignment quality of >Q30 and were required to be properly paired. Reads mapping to the mitochondria, unmapped contigs and chromosome Y are not considered.\n",
    "\n",
    "We used MACS2 to call ATAC-seq peaks. The parameters that we used for MACS2 are (--nomodel --nolambda --keep-dup all --call-summits). Peaks were filtered using the consensus excludable [ENCODE blacklist](http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/). We have saved it here in our local disk. Using the filtered peak set, peak summits were extended +/-250 bps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explanation of different parameters for the tools\n",
    "\n",
    "For convenience of the user an explanation of the parameters used with different tools are given below\n",
    "\n",
    "### Bowtie options\n",
    "- Pass the trimmed fastq files to bowtie2 for mapping.\n",
    "- Bowtie2 mapping with default parameters using 16 cores, except pass a -X 2000, \n",
    "- which treat read pairs that are within 2000 bp as concordant mapping.\n",
    "- This only affects concordant rate, it won't affect mapping rate.\n",
    "- -X 2000 is suggested by the Greenleaf paper. That's the reason I used this.\n",
    "\n",
    "### Samtools option\n",
    "- Convert the sam file, including sam header, to bam file, and remove non-mapped reads.\n",
    "- The following command will do these in one go.\n",
    "- -S option indicates the input is a sam file\n",
    "- -b option asked samtools to write output as a bam file\n",
    "- -h option asked samtools to print header, it is very important to include this \"-h\"\n",
    "- when converting between sam <-> bam files.\n",
    "- -F 4 option removes non-mapped reads\n",
    "- -f2 only include reads with all bits set in INT set in FLAG [0]\n",
    "- -q only include reads with mapping quality >= INT [0]\n",
    "\n",
    "### Picard tools\n",
    "- Picard tools is used to remove the duplicates and to estimate the library size with looking into the column ESTIMATE_LIBRARY_SIZE\n",
    "\n",
    "### Samtools to remove mitochondrial, unmapped contigs and chrY\n",
    "Reads mapping to the mitochondria, unmapped contigs and chromosome Y were removed and not considered. We need to cleanse BAM files of chrM, and unassembled \"random\" contigs before running ATAC-seq analysis\n",
    "\n",
    "### Removing the Blacklists\n",
    "Functional genomics experiments based on next-gen sequencing (e.g. ChIP-seq, MNase-seq, DNase-seq, FAIRE-seq) that measure biochemical activity of various elements in the genome often produce artifact signal in certain regions of the genome. It is important to keep track of and filter artifact regions that tend to show artificially high signal (excessive unstructured anomalous reads mapping). Below is a list of comprehensive empirical blacklists identified by the ENCODE and modENCODE consortia. Note that these blacklists were empirically derived from large compendia of data using a combination of automated heuristics and manual curation. These blacklists are applicable to functional genomic data based on short-read sequencing (20-100bp reads). These are not directly applicable to RNA-seq or any other transcriptome data types. The blacklisted regions typically appear u niquely mappable so simple mappability filters do not remove them. These regions are often found at specific types of repeats such as centromeres, telomeres and satellite repeats. It is especially important to remove these regions that computing measures of similarity such as Pearson correlation between genome-wide tracks that are especially affected by outliers.\n",
    "\n",
    "download the `hg19 blacklist` into the data folder by `wget` [hg19 Blacklist](http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeDacMapabilityConsensusExcludable.bed.gz)\n",
    "\n",
    "### Perform peak calling using macs2 (2 2.1.0.20151222)\n",
    "Perform peak calling and generate bedGraph using 50 bp reads centred on 5 prime cutting end, within a 50-bp smoothing window for bedGraph signal track. The parameters for macs2 are as follows\n",
    "\n",
    "- `-t/--treatment` FILENAME: This is the only REQUIRED parameter for MACS. File can be in any supported format specified by –format option. Check –format for detail. If you have more than one alignment files, you can specify them as `-t A B C`. MACS will pool up all these files together.\n",
    "- `-n/--name`: The name string of the experiment. MACS will use this string NAME to create output files like ‘NAME_peaks.xls’, ‘NAME_negative_peaks.xls’, ‘NAME_peaks.bed’ , ‘NAME_summits.bed’, ‘NAME_model.r’ and so on. So please avoid any confliction between these filenames and your existing files.\n",
    "- `-g/--gsize` <br><br> PLEASE assign this parameter to fit your needs! <br> <br> It’s the mappable genome size or effective genome size which is defined as the genome size which can be sequenced. Because of the repetitive features on the chromsomes, the actual mappable genome size will be smaller than the original size, about 90% or 70% of the genome size. The default hs – 2.7e9 is recommended for UCSC human hg18 assembly. Here are all precompiled parameters for effective genome size: <br> __hs__:\t2.7e9 <br> __mm__:\t1.87e9 <br> __ce__:\t9e7 <br> __dm__: 1.2e8\n",
    "- `-f/--format FORMAT`: Format of tag file, can be “ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE” or “BAMPE”. Default is “AUTO” which will allow MACS to decide the format automatically. “AUTO” is also usefule when you combine different formats of files.\n",
    "- `--nomodel`: While on, MACS will bypass building the shifting model.\n",
    "- `--shift` : <br> Note, this is NOT the legacy –shiftsize option which is replaced by –extsize! You can set an arbitrary shift in bp here. Please Use discretion while setting it other than default value (0). When –nomodel is set, MACS will use this value to move cutting ends (5’) then apply –extsize from 5’ to 3’ direction to extend them to fragments. When this value is negative, ends will be moved toward 3’->5’ direction, otherwise 5’->3’ direction. Recommended to keep it as default 0 for ChIP-Seq datasets, or -1 * half of EXTSIZE together with –extsize option for detecting enriched cutting loci such as certain DNAseI-Seq datasets. Note, you can’t set values other than 0 if format is BAMPE for paired-end data. Default is 0. <br><br> Here are some examples for combining –shift and –extsize: <br><br>1. To find enriched cutting sites such as some DNAse-Seq datasets. In this case, all 5’ ends of sequenced reads should be extended in both direction to smooth the pileup signals. If the wanted smoothing window is 200bps, then use ‘–nomodel –shift -100 –extsize 200’.<br><br>2. For certain nucleosome-seq data, we need to pileup the centers of nucleosomes using a half-nucleosome size for wavelet analysis (e.g. NPS algorithm). Since the DNA wrapped on nucleosome is about 147bps, this option can be used: ‘–nomodel –shift 37 –extsize 73’.\n",
    "- `--extsize`: While ‘–nomodel’ is set, MACS uses this parameter to extend reads in 5’->3’ direction to fix-sized fragments. For example, if the size of binding region for your transcription factor is 200 bp, and you want to bypass the model building by MACS, this parameter can be set as 200. This option is only valid when –nomodel is set or when MACS fails to build model and –fix-bimodal is on.\n",
    "- `-B/--bdg`: If this flag is on, MACS will store the fragment pileup, control lambda, -log10(pvalue) and -log10(qvalue) scores in bedGraph files. The bedGraph files will be stored in current directory named NAME+’_treat_pileup.bdg’ for treatment data, NAME+’_control_lambda.bdg’ for local lambda values from control, NAME+’_treat_pvalue.bdg’ for Poisson pvalue scores (in -log10(pvalue) form), and NAME+’_treat_qvalue.bdg’ for q-value scores from Benjamini–Hochberg–Yekutieli procedure [False discovery rate](http://en.wikipedia.org/wiki/False_discovery_rate#Dependent_tests).\n",
    "- `--SPMR`: Normalises the coverage plot values by millions of reads. \n",
    "- `--call-summits`: MACS will now reanalyze the shape of signal profile (p or q-score depending on cutoff setting) to deconvolve subpeaks within each peak called from general procedure. It’s highly recommended to detect adjacent binding events. While used, the output subpeaks of a big peak region will have the same peak boundaries, and different scores and peak summit positions.\n",
    "\n",
    "- <span style=\"color:green\"> __ To call peaks and generate bedGraph using a different smoothing window, change --shift 0.5*window --extsize window  for example, --shift 75 --extsize 150 will use 150 bp reads centred on 5 prime cutting end to call peaks, and use a 150-bp smoothing window to generate bedGraphgreen __ </span>\n",
    "\n",
    "### macs2 _Output Files_ Format\n",
    "1. NAME_peaks.xls is a tabular file which contains information about called peaks. You can open it in excel and sort/filter using excel functions. Information include:\n",
    "    * chromosome name\n",
    "    * start position of peak\n",
    "    * end position of peak\n",
    "    * length of peak region\n",
    "    * absolute peak summit position\n",
    "    * pileup height at peak summit, -log10(pvalue) for the peak summit (e.g. pvalue =1e-10, then this value should be 10)\n",
    "    * fold enrichment for this peak summit against random Poisson distribution with local lambda, -log10(qvalue) at peak summit\n",
    "2. NAME_peaks.narrowPeak is BED6+4 format file which contains the peak locations together with peak summit, pvalue and qvalue. You can load it to UCSC genome browser. Definition of some specific columns are:\n",
    "    * 5th: integer score for display\n",
    "    * 7th: fold-change\n",
    "    * 8th: -log10(pvalue)\n",
    "    * 9th: -log10(qvalue)\n",
    "    * 10th: relative summit position to peak start\n",
    "    The file can be loaded directly to UCSC genome browser. Remove the beginning track line if you want to analyze it by other tools.\n",
    "\n",
    "\n",
    "\n",
    "### Convert bedGraph to bigWig\n",
    "- To call peaks and generate bedGraph using a different smoothing window, change --shift 0.5*window --extsize window\n",
    "- For example, --shift 75 --extsize 150 will use 150 bp reads centred on 5 prime cutting end to call peaks, and use a 150-bp smoothing window to generate bedGraph\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing python modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import subprocess, os, sys, signal, pip\n",
    "from os import listdir\n",
    "from os.path import isfile, join\n",
    "import multiprocessing\n",
    "import threading\n",
    "import time\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting the folders and tools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we configure the __Folders__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Setting the output folder where all the results would be stored\n",
    "outputFolder = '/home/baker/my-dpsf-scratch-share/Syed_scATAC_and_scRNA_comparison/Comparison_Rep2/output/'\n",
    "\n",
    "# inputFolder : Folder name with all the fastq files\n",
    "intputFolder = '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs'\n",
    "\n",
    "# Blacklisted sequences\n",
    "blackListFile = '/home/baker/my-hydra-share/Packages/Blacklist_regions/consensusBlacklist.bed'\n",
    "\n",
    "# Setting up the genome folder\n",
    "#ref_genome = '/home/baker/my-mm10-index-share/Mus_musculus.GRCm38.71.fa'\n",
    "ref_genome = '/home/baker/my-edBritton-Bowtie2-genome/genome'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the softwares\n",
    "\n",
    "We configure the __software__ parameters here. If the required softwares are not in the PATH you can manually set their location here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bowtie_path = '/home/baker/Downloads/bowtie2-2.2.6/bowtie2'\n",
    "samtools_path = 'samtools'\n",
    "macs2_path = 'macs2'\n",
    "bdg2bw_path = '/home/baker/my-hydra-share/Packages/bdg2bw'\n",
    "picard_path = '/home/baker/picard/dist/picard.jar'\n",
    "intersectBed_path = 'intersectBed'\n",
    "MarkDuplicate_path = '/home/baker/my-hydra-share/Packages/picard-tools-1.115/MarkDuplicates.jar'\n",
    "Trimmomatic_path = '/home/baker/my-hydra-share/Packages/Trimmomatic-0.35/trimmomatic-0.35.jar'\n",
    "fastqc_path ='/home/baker/my-hydra-share/Packages/FastQC/fastqc'\n",
    "Trimmomatic_Adapter = '/home/baker/my-hydra-share/Packages/Trimmomatic-0.35/adapters/NexteraPE-PE.fa:2:30:10:1:true TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:25'\n",
    "awk_cmd_path = '''awk 'BEGIN {OFS = \"\"\"\\t\"\"\"} ; {print $1, $2 - 250, $3 + 250, $4, $5}' '''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "softwares = {    \n",
    "    'bowtie2': bowtie_path,\n",
    "    'macs2': macs2_path,\n",
    "    'bdg2bw': bdg2bw_path,\n",
    "    'samtools': samtools_path,\n",
    "    'picard': picard_path,\n",
    "    'intersectBed': intersectBed_path,\n",
    "    'MarkDuplicate': MarkDuplicate_path,\n",
    "    'trimmomatic': Trimmomatic_path,\n",
    "    'fastqc': fastqc_path,\n",
    "    'Trimmomatic_Adapter': Trimmomatic_Adapter,\n",
    "    'awk_cmd' : awk_cmd_path\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating the output folders\n",
    "We now create the output folders where all the processed files will be stored"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "output_folders = [ 'Bowtie_files', 'Blacklist_removed' # Trimmo, Bowtie, Blacklist Removed Files\n",
    "                 , 'Trimmomatic_Files', 'Fastqc_SE_Files' \n",
    "                 , 'Duplicates_removed', 'Macs2_files'              # Synced R1 and R2\n",
    "                 , 'Merged_BAM', 'Merged_Macs2_files'\n",
    "                 , 'Bam_Files_Filtered_on_Library_Size', 'Merged_Filtered_BAM'\n",
    "                 ,  'Merged_Filtered_Macs2_files', 'Filtered_Macs2_files'\n",
    "                 ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for folder in output_folders:\n",
    "    os.makedirs(os.path.join(outputFolder, folder))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Custom functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we define __custom__ functions that we will use for file processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "remove_extension = lambda x: x.split('.')[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following functions deals with the inputs and the outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_args(read1, read2, ref_genome, blackListFile, 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",
    "    #sample_name = r1_shortname.split('_')[0] + '_'+ r1_shortname.split('_')[1]\n",
    "    sample_name = r1_shortname.split('_')[0]\n",
    "\n",
    "    args = {  \n",
    "        'r1_input': read1,\n",
    "        'r2_input': read2,\n",
    "        'ref_genome': ref_genome,\n",
    "        'blackListFile': blackListFile,\n",
    "        'sample_name': sample_name,\n",
    "    }\n",
    "    \n",
    "    # output_paths = {folder: os.path.join(outputFolder, folder, r1_shortname) for folder in output_folders}\n",
    "    output_paths = {folder: os.path.join(outputFolder, folder, sample_name) for folder in output_folders}\n",
    "    \n",
    "    return dict(args, **output_paths)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the pipeline\n",
    "\n",
    "In order to run the pipeline we first define the commands. A brief explanation of the tools are given below\n",
    "\n",
    "- `trimmomatic`: To trim the adapter and low quality reads\n",
    "- `bowti2`: To map the reads with genome. Here we are mapping with Human genome. The genome is defined accordingly in the `{ref_genome}` folder\n",
    "- `samtools`: In this samtools view command we retain only the properly paired reads and reads with quality above 30\n",
    "- `intersectBed`: We remove the blacklists\n",
    "- `MarkDuplicate`: With this picard tool we are removing the duplicates [where duplicate reads are defined as originating from a single fragment of DNA which can arise during sample preparation e.g. library construction using PCR]. MarkDuplicate just marks the duplicate sequences but to make sure that it removes the duplcate from the bam file the parameter `REMOVE_DUPLICATES` is set to `TRUE`\n",
    "- `samtools idxstats`: samtools idxstas are used to generate the mapping quality of the reads\n",
    "- `macs2 callpeak`: Here we call the peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cmds = [\n",
    "    \n",
    "     #'java -jar {trimmomatic} PE -phred33 {r1_input} {r2_input} {Trimmomatic_Files}_r1_paired.fq {Trimmomatic_Files}_r1_unpaired.fq {Trimmomatic_Files}_r2_paired.fq {Trimmomatic_Files}_r2_upaired.fq ILLUMINACLIP:{Trimmomatic_Adapter}',\n",
    "    \n",
    "     #'{bowtie2} -p 12 -X 2000 --dovetail -x {ref_genome} -1 {Trimmomatic_Files}_r1_paired.fq -2 {Trimmomatic_Files}_r2_paired.fq -S {Bowtie_files}.sam 2> {Bowtie_files}_allignment.log',\n",
    "\n",
    "    '{bowtie2} -p 12 -X 2000 --dovetail -x {ref_genome} -1 {r1_input} -2 {r2_input} -S {Bowtie_files}.sam 2> {Bowtie_files}_allignment.log',\n",
    "    \n",
    "    '{samtools} view -SbhF 4 -f2 -q30 {Bowtie_files}.sam > {Bowtie_files}.bam' ,\n",
    "    \n",
    "    'intersectBed -v -abam {Bowtie_files}.bam -b {blackListFile} > {Blacklist_removed}_noexclu.bam ' ,\n",
    "    \n",
    "    '{samtools} sort {Blacklist_removed}_noexclu.bam {Blacklist_removed}_noexclu_sorted' ,\n",
    "    \n",
    "    'java -Xmx2g -jar {MarkDuplicate}  INPUT={Blacklist_removed}_noexclu_sorted.bam OUTPUT={Duplicates_removed}_nodup.bam METRICS_FILE={Duplicates_removed}_nodup_stats REMOVE_DUPLICATES=True' ,\n",
    "    \n",
    "    '{samtools} sort {Duplicates_removed}_nodup.bam {Duplicates_removed}_nodup_sorted' ,\n",
    "    \n",
    "    '{samtools} index {Duplicates_removed}_nodup_sorted.bam',\n",
    "            \n",
    "    '{samtools} view -b {Duplicates_removed}_nodup_sorted.bam chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrX  > {Duplicates_removed}_nodup_sorted_cleaned.bam' ,\n",
    "        \n",
    "    '{samtools} index {Duplicates_removed}_nodup_sorted_cleaned.bam',\n",
    "    \n",
    "    '{samtools} idxstats {Duplicates_removed}_nodup_sorted_cleaned.bam > {Duplicates_removed}_nodup_sorted_cleaned_chrom_stat.txt' ,\n",
    "    \n",
    "    'macs2 callpeak -t {Duplicates_removed}_nodup_sorted_cleaned.bam -n {Macs2_files} -q 0.01 -g hs -f BAM --nomodel --nolambda --shift -75 --extsize 150 -B --SPMR --keep-dup all --call-summits',\n",
    "    \n",
    "    '{awk_cmd}  {Macs2_files}_summits.bed > {Macs2_files}_summits_shifted.bed',\n",
    "        \n",
    "    '{bdg2bw} {Macs2_files}_treat_pileup.bdg /home/baker/my-hydra-share/Packages/Homer_4.6/bin/hg19.chrom.sizes',\n",
    "        \n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the reads from the `inputFolder`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for root, folders, files in os.walk(intputFolder):\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": [
    "Print the read names. This is to make sure that the correct files are picked up. This can be run just once as it is expected to produce large output. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E01_S17_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E02_S9_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E03_S1_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E04_S65_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E05_S57_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E06_S49_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E07_S18_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E08_S10_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E09_S2_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E10_S66_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E11_S58_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E12_S50_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E13_S19_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E14_S11_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E15_S3_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E16_S67_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E17_S59_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E18_S51_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E19_S20_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E20_S12_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E21_S4_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E22_S68_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E23_S60_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E24_S52_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E25_S5_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E26_S13_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E27_S21_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E28_S53_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E29_S61_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E30_S69_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E31_S6_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E32_S14_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E33_S22_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E34_S54_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E35_S62_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E36_S70_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E37_S7_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E38_S15_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E39_S23_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E40_S55_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E41_S63_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E42_S71_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E43_S8_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E44_S16_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E45_S24_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E46_S56_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E47_S64_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E48_S72_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E49_S25_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E50_S33_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E51_S41_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E52_S73_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E53_S81_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E54_S89_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E55_S26_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E56_S34_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E57_S42_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E58_S74_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E59_S82_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E60_S90_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E61_S27_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E62_S35_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E63_S43_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E64_S75_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E65_S83_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E66_S91_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E67_S28_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E68_S36_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E69_S44_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E70_S76_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E71_S84_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E72_S92_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E73_S45_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E74_S37_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E75_S29_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E76_S93_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E77_S85_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E78_S77_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E79_S46_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E80_S38_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E81_S30_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E82_S94_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E83_S86_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E84_S78_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E85_S47_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E86_S39_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E87_S31_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E88_S95_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E89_S87_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E90_S79_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E91_S48_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E92_S40_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E93_S32_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E94_S96_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E95_S88_R1_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E96_S80_R1_001.fastq.gz'] ['/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E01_S17_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E02_S9_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E03_S1_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E04_S65_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E05_S57_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E06_S49_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E07_S18_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E08_S10_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E09_S2_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E10_S66_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E11_S58_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E12_S50_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E13_S19_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E14_S11_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E15_S3_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E16_S67_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E17_S59_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E18_S51_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E19_S20_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E20_S12_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E21_S4_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E22_S68_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E23_S60_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E24_S52_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E25_S5_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E26_S13_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E27_S21_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E28_S53_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E29_S61_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E30_S69_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E31_S6_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E32_S14_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E33_S22_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E34_S54_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E35_S62_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E36_S70_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E37_S7_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E38_S15_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E39_S23_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E40_S55_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E41_S63_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E42_S71_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E43_S8_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E44_S16_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E45_S24_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E46_S56_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E47_S64_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E48_S72_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E49_S25_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E50_S33_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E51_S41_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E52_S73_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E53_S81_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E54_S89_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E55_S26_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E56_S34_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E57_S42_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E58_S74_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E59_S82_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E60_S90_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E61_S27_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E62_S35_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E63_S43_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E64_S75_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E65_S83_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E66_S91_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E67_S28_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E68_S36_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E69_S44_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E70_S76_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E71_S84_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E72_S92_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E73_S45_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E74_S37_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E75_S29_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E76_S93_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E77_S85_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E78_S77_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E79_S46_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E80_S38_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E81_S30_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E82_S94_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E83_S86_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E84_S78_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E85_S47_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E86_S39_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E87_S31_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E88_S95_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E89_S87_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E90_S79_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E91_S48_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E92_S40_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E93_S32_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E94_S96_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E95_S88_R2_001.fastq.gz', '/home/baker/Conor_Rogerson_ATAC_seq_rep2/fastqs/E96_S80_R2_001.fastq.gz']\n"
     ]
    }
   ],
   "source": [
    "print (reads1, reads2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The platform is now set to run the commands. Below we run the command sequentially. However, it is always good to first check that the right commands are executed. To test that one can uncomment the print command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def run_cmd(cmds, args):\n",
    "    for cmd in cmds: \n",
    "        #print(cmd.format(**args))\n",
    "        subprocess.call(cmd.format(**args), shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the number of cores in the node"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total Cores: 8\n"
     ]
    }
   ],
   "source": [
    "print(\"Total Cores:\",multiprocessing.cpu_count())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sequential run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for read1, read2 in zip(reads1, reads2):\n",
    "    args = get_args(read1, read2, ref_genome, blackListFile, output_folders)\n",
    "    args = dict(args, **softwares)\n",
    "    run_cmd(cmds,args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering low Library sized BAM files\n",
    "We filter out the cells for which the library size as estimated by picard is $<LIBRARY\\_SIZE\\_THRESHOLD$. We look at the stats generated during the duplicate check stat and decide based on the parameter value `ESTIMATED_LIBRARY_SIZE`. By default $LIBRARY\\_SIZE\\_THRESHOLD = 10000$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filtering the BAM files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "LIBRARY_SIZE_THRESHOLD = 10000\n",
    "\n",
    "NoDupStatPath = outputFolder + 'Duplicates_removed'\n",
    "FilteredFilePath = outputFolder + 'Bam_Files_Filtered_on_Library_Size/'\n",
    "for root, folders, files in os.walk(NoDupStatPath):\n",
    "    files = [os.path.join(root, f) for f in files if (f.endswith('nodup_stats'))]\n",
    "LibrarySize = 0\n",
    "for f in files:\n",
    "    with open(f,'rt') as DupStatFile:\n",
    "        FLAG = 0\n",
    "        for line in DupStatFile:\n",
    "            if 'ESTIMATED_LIBRARY_SIZE' in line:\n",
    "                FLAG = 1\n",
    "            elif FLAG == 1:\n",
    "                LibrarySize = float(line.rsplit(None, 1)[-1])\n",
    "                FLAG = 0\n",
    "        if LibrarySize > LIBRARY_SIZE_THRESHOLD:        \n",
    "            TempCellName = os.path.basename(f).split('_')\n",
    "            CellName = TempCellName[0]\n",
    "            BamFiles = [join(NoDupStatPath, f) for f in listdir(NoDupStatPath) if (f.startswith(CellName) & f.endswith('nodup_sorted.bam'))]\n",
    "            CopyFilteredCmd = 'cp ' + BamFiles[0] + ' ' + FilteredFilePath\n",
    "            subprocess.call(CopyFilteredCmd, shell=True)\n",
    "            BamFiles = [join(NoDupStatPath, f) for f in listdir(NoDupStatPath) if (f.startswith(CellName) & f.endswith('nodup_sorted.bam.bai'))]\n",
    "            CopyFilteredCmd = 'cp ' + BamFiles[0] + ' ' + FilteredFilePath\n",
    "            subprocess.call(CopyFilteredCmd, shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filtering the macs2 files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "LIBRARY_SIZE_THRESHOLD = 10000\n",
    "\n",
    "NoDupStatPath = outputFolder + 'Duplicates_removed'\n",
    "FilteredFilePath = outputFolder + 'Filtered_Macs2_files/'\n",
    "Macs2Files = outputFolder + 'Macs2_files'\n",
    "for root, folders, files in os.walk(NoDupStatPath):\n",
    "    files = [os.path.join(root, f) for f in files if (f.endswith('nodup_stats'))]\n",
    "LibrarySize = 0\n",
    "for f in files:\n",
    "    with open(f,'rt') as DupStatFile:\n",
    "        FLAG = 0\n",
    "        for line in DupStatFile:\n",
    "            if 'ESTIMATED_LIBRARY_SIZE' in line:\n",
    "                FLAG = 1\n",
    "            elif FLAG == 1:\n",
    "                LibrarySize = float(line.rsplit(None, 1)[-1])\n",
    "                FLAG = 0\n",
    "        if LibrarySize > LIBRARY_SIZE_THRESHOLD:        \n",
    "            TempCellName = os.path.basename(f).split('_')\n",
    "            CellName = TempCellName[0]\n",
    "            MacsFiles = [join(Macs2Files, f) for f in listdir(Macs2Files) if (f.startswith(CellName) & f.endswith('control_lambda.bdg'))]\n",
    "            CopyFilteredCmd = 'cp ' + MacsFiles[0] + ' ' + FilteredFilePath\n",
    "            subprocess.call(CopyFilteredCmd, shell=True)\n",
    "            MacsFiles = [join(Macs2Files, f) for f in listdir(Macs2Files) if (f.startswith(CellName) & f.endswith('peaks.narrowPeak'))]\n",
    "            CopyFilteredCmd = 'cp ' + MacsFiles[0] + ' ' + FilteredFilePath\n",
    "            subprocess.call(CopyFilteredCmd, shell=True)\n",
    "            MacsFiles = [join(Macs2Files, f) for f in listdir(Macs2Files) if (f.startswith(CellName) & f.endswith('peaks.xls'))]\n",
    "            CopyFilteredCmd = 'cp ' + MacsFiles[0] + ' ' + FilteredFilePath\n",
    "            subprocess.call(CopyFilteredCmd, shell=True)\n",
    "            MacsFiles = [join(Macs2Files, f) for f in listdir(Macs2Files) if (f.startswith(CellName) & f.endswith('summits.bed'))]\n",
    "            CopyFilteredCmd = 'cp ' + MacsFiles[0] + ' ' + FilteredFilePath            \n",
    "            subprocess.call(CopyFilteredCmd, shell=True)\n",
    "            MacsFiles = [join(Macs2Files, f) for f in listdir(Macs2Files) if (f.startswith(CellName) & f.endswith('treat_pileup.bdg'))]\n",
    "            CopyFilteredCmd = 'cp ' + MacsFiles[0] + ' ' + FilteredFilePath            \n",
    "            subprocess.call(CopyFilteredCmd, shell=True)\n",
    "            MacsFiles = [join(Macs2Files, f) for f in listdir(Macs2Files) if (f.startswith(CellName) & f.endswith('treat_pileup.bw'))]\n",
    "            CopyFilteredCmd = 'cp ' + MacsFiles[0] + ' ' + FilteredFilePath            \n",
    "            subprocess.call(CopyFilteredCmd, shell=True)          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Generating mapping summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "LogFilePath = outputFolder + 'Bowtie_files'\n",
    "CellMapSummary = 'Cell_ID Number_of_input_reads Unmapped_reads_number Percent_Unmapped Uniquely_mapped_reads_number Percent_of_uniquely_mapped \\n'\n",
    "for root, folder, files in os.walk(LogFilePath):\n",
    "    files = [os.path.join(root, f) for f in files if (f.endswith('_allignment.log'))]\n",
    "    for f in files:\n",
    "        with open(f,'rt') as LogFile:\n",
    "            for line in LogFile:\n",
    "                if 'reads; of these:' in line:\n",
    "                    CellMapSummary = CellMapSummary + f.split('/')[-1].split('_')[0]+ '_' +f.split('/')[-1].split('_')[1]+ ' ' + line.split(' reads')[0].strip()\n",
    "                if 'aligned concordantly exactly 1 time' in line:\n",
    "                    CellMapSummary = CellMapSummary + ' ' + line.split('(')[0].strip(' ').strip()\n",
    "                if 'aligned concordantly exactly 1 time' in line:                    \n",
    "                    CellMapSummary = CellMapSummary + ' ' + line.split('(')[1].split('%')[0].strip()\n",
    "                if ') aligned concordantly 0 times' in line:\n",
    "                    CellMapSummary = CellMapSummary +  ' ' + line.split('(')[0].strip()\n",
    "                if ') aligned concordantly 0 times' in line:\n",
    "                    CellMapSummary = CellMapSummary + ' ' + line.split('(')[1].split('%')[0].strip()\n",
    "            CellMapSummary = CellMapSummary + '\\n'\n",
    "\n",
    "fileToWrite = outputFolder + \"MappingSummary_rep2_scATAC_withoutTrimmo.csv\"\n",
    "target = open(fileToWrite, 'w')                    \n",
    "target.write(CellMapSummary)\n",
    "target.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}