Scasat_Pre-process_rep1_BAMPE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ATACseq pipeline Human Batch-1"
   ]
  },
  {
   "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": [
    "__Prerequisites__\n",
    "The following has to be installed before hand in order to run the pipeline. This pipeline runs in the distributed server with multiple threads.\n",
    "\n",
    "- [bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)\n",
    "- [samtools](http://www.htslib.org)\n",
    "- [macs2](https://github.com/taoliu/MACS)\n",
    "- [picard](http://broadinstitute.github.io/picard/)\n",
    "- [bedtools](http://bedtools.readthedocs.io/en/latest/)\n",
    "- [trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic)\n",
    "- [bdg2bw](https://github.com/ManchesterBioinference/Scasat/blob/master/bdg2bw_local)\n",
    "\n",
    "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\n",
    "\n",
    "__Acknowledgement:__\n",
    "\n",
    "- Portion of the code was modified from C1-CAGE processing pipeline by Charles Plessy <plessy@riken.jp>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Description\n",
    "\n",
    "In this notebook we process the single cell ATAC-seq in Batch-1. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data processing summary\n",
    "\n",
    "We used `trimmomatic` to trim the adapter sequences. Paired-end reads were aligned to `hg19` 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."
   ]
  },
  {
   "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",
    "One can 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",
    "### 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",
    "- bdgbw is used to convert the bedGraph files to bigWig file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing python modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "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": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Setting the output folder where all the results would be stored\n",
    "outputFolder = '../../my-single-cell-share/scratch-backup-dpsf/Syed_scATAC_and_scRNA_comparison/Comparison_Rep1_Improved_MACS2/output/'\n",
    "\n",
    "# inputFolder : Folder name with all the fastq files\n",
    "intputFolder = '../Conor_Rogerson_ATAC_seq_rep1/fastqs/'\n",
    "\n",
    "# Blacklisted sequences\n",
    "blackListFile = '../../my-dpsf-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 = '../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": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bowtie_path = '../Downloads/bowtie2-2.2.6/bowtie2'\n",
    "samtools_path = 'samtools'\n",
    "macs2_path = 'macs2'\n",
    "bdg2bw_path = '../my-hydra-share/Packages/bdg2bw'\n",
    "picard_path = '../picard/dist/picard.jar'\n",
    "intersectBed_path = 'intersectBed'\n",
    "MarkDuplicate_path = '../my-dpsf-share/Packages/picard-tools-1.115/MarkDuplicates.jar'\n",
    "Trimmomatic_path = '../my-dpsf-share/Packages/Trimmomatic-0.35/trimmomatic-0.35.jar'\n",
    "fastqc_path ='../my-dpsf-share/Packages/FastQC/fastqc'\n",
    "Trimmomatic_Adapter = '../my-dpsf-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": 4,
   "metadata": {},
   "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": 5,
   "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'              \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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "for folder in output_folders:\n",
    "    if not os.path.exists(os.path.join(outputFolder, folder)):\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": 6,
   "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": 7,
   "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]\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",
    "  \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": 8,
   "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:/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",
    "    \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",
    "    '{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} -p 0.0001 -g hs -f BAMPE --nomodel --nolambda -B --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": 9,
   "metadata": {},
   "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": "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": 11,
   "metadata": {},
   "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": 12,
   "metadata": {},
   "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": 13,
   "metadata": {
    "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": 12,
   "metadata": {},
   "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_cleaned.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_cleaned.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": 13,
   "metadata": {},
   "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('summits_shifted.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": {},
   "source": [
    "## Generating read summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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_rep1_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": {
  "anaconda-cloud": {},
  "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": 1
}