Sandbox_mm9.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import cmath\n",
      "import math\n",
      "import sys\n",
      "import numpy as np\n",
      "import glob \n",
      "import subprocess\n",
      "import re\n",
      "from matplotlib_venn import venn2\n",
      "import pandas as pd\n",
      "from collections import defaultdict\n",
      "from operator import itemgetter\n",
      "import matplotlib as mpl\n",
      "import matplotlib.pyplot as plt\n",
      "from optparse import OptionParser\n",
      "mpl.rcParams['savefig.dpi'] = 2 * mpl.rcParams['savefig.dpi']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "(1) Specify dataset - "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get raw data\n",
      "global outfilepath\n",
      "sampleName='MMhur'\n",
      "outfilepath=os.getcwd() + '/results/%s/'%sampleName"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 48
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "(2) Re-map the data (if desired, else pipeline will use existing Samfiles) - "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def runBowtie(fastqFiles):\n",
      "    # Useage: Short read mapping to reference (hg19).\n",
      "    # Input: Fastq files of replicates (not paired end).\n",
      "    # Output: Path to samfile for each read.\n",
      "    program = 'bowtie2'\n",
      "    mappedReads=[]\n",
      "    unMappedReads=[]\n",
      "    print \"Performing Bowtie...\"\n",
      "    # Parameters\n",
      "    # In -k mode, Bowtie 2 searches for up to N distinct, valid alignments for each read.\n",
      "    # N equals the integer specified with the -k parameter. \n",
      "    # That is, if -k 2 is specified, Bowtie 2 will search for at most 2 distinct alignments\n",
      "    # See http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml\n",
      "    for infastq in fastqFiles:\n",
      "\n",
      "        k='1'\n",
      "        # *** CHANGE INDEX ***\n",
      "        index = os.getcwd() + '/docs/mm9/mm9/mm9'\n",
      "        outfile = modifyName(infastq,\"mapped.sam\")\n",
      "        # *** CHANGE INDEX ***\n",
      "        unmapped = modifyName(infastq,\"notMappedTomm9.fastq\")\n",
      "        \n",
      "        print \"Input file:\"\n",
      "        print infastq \n",
      "        print 'Genome index:'\n",
      "        print index\n",
      "        print \"Output file (mapped):\"\n",
      "        print outfile\n",
      "        print \"Output file (unmapped)\"\n",
      "        print unmapped\n",
      "        proc = subprocess.Popen([program,'-x', index,'-k',k,'-U',infastq,'--un',unmapped,'-S',outfile,'2>>'],stdout=subprocess.PIPE,stderr=subprocess.PIPE)\n",
      "        out, err = proc.communicate()\n",
      "        result = out.decode()\n",
      "        error = err.decode()\n",
      "        print \"Result : \",result \n",
      "        print \"Error : \",error\n",
      "        mappedReads = mappedReads + [outfile]\n",
      "        unMappedReads = unMappedReads + [unmapped]\n",
      "    return (mappedReads,unMappedReads)\n",
      "\n",
      "# Run Bowtie\n",
      "print \"Run mapping\"  \n",
      "readsFormm9mapping=glob.glob(outfilepath+\"*_notMappedToRepeat.fastq\")\n",
      "mappedReads,unmappedReads=runBowtie(readsFormm9mapping)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Run mapping\n",
        "Performing Bowtie...\n"
       ]
      }
     ],
     "prompt_number": 49
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "(3) Set threshold param and run pipeline - "
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "*** STOP: CHECK TO MAKE SURE YOU ARE WRITING TO THE CORRECT DIRECTORY BEFORE STARTING PIPE! ***"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Hello\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Hello\n"
       ]
      }
     ],
     "prompt_number": 50
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print outfilepath"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "/arrayAhome/lmartin/CLIP/results/MMhur/\n"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run Bowtie\n",
      "print \"Run mapping\"  \n",
      "readsFormm9mapping=glob.glob(outfilepath+\"*_notMappedToRepeat.fastq\")\n",
      "mappedReads,unmappedReads=runBowtie(readsFormm9mapping)\n",
      "\n",
      "# Get hg19 mapped files and run Samtools without duplicate removal\n",
      "print \"Process mapped data\"  \n",
      "mappedReads=glob.glob(outfilepath+\"*_mapped.sam\")\n",
      "mappedBedFiles=runSamtools(mappedReads)\n",
      "    \n",
      "# Run filters\n",
      "blacklistedBedFiles=runBlacklistRegions(mappedBedFiles)\n",
      "maskedBedFiles=runRepeatMask(blacklistedBedFiles)\n",
      "\n",
      "# Seperate strands\n",
      "readsByStrand=seperateStrands(maskedBedFiles)\n",
      "\n",
      "# Modify strands\n",
      "negativeRTstop=isolate5prime(modifyNegativeStrand(readsByStrand[0])) # negative strand reads\n",
      "positiveRTstop=isolate5prime(readsByStrand[1]) # positive strand reads \n",
      "\n",
      "###################\n",
      "#### KEY PARAM ####\n",
      "threshold=8\n",
      "###################\n",
      "print \"Merge RT stops and run CLIPper\"    \n",
      "# Combine RT stops from replicates for each strand\n",
      "posMerged = outfilepath + sampleName + '_positivereads.mergedRT'\n",
      "negMerged = outfilepath + sampleName + '_negativereads.mergedRT'\n",
      "# Merge \n",
      "mergeRT(negativeRTstop,negMerged,threshold)\n",
      "mergeRT(positiveRTstop,posMerged,threshold)\n",
      "# Merge across strands\n",
      "negAndPosMerged = outfilepath + sampleName + '_threshold=%s'%threshold + '_allreads.mergedRT.bed'\n",
      "fileCat(negAndPosMerged,[posMerged,negMerged])\n",
      "\n",
      "# Run CLIPper\n",
      "CLIPPERio=runCLIPPER(negAndPosMerged)\n",
      "CLIPPERin=CLIPPERio[0]\n",
      "CLIPPERout=CLIPPERio[1]\n",
      "\n",
      "# Extract lowFDR reads based upon strand\n",
      "clipperStats=modCLIPPERout(CLIPPERin,CLIPPERout)\n",
      "CLIPPERlowFDR=clipperStats[0] # Low FDR reads returned filtred through CLIPper windows\n",
      "CLIPpeReadsPerCluster=clipperStats[1] # Number of reads per CLIPper cluster\n",
      "CLIPpeGeneList=clipperStats[2] # Gene names returned from the CLIPper file\n",
      "CLIPperOutBed=clipperStats[3] # CLIPper windows as a bed file\n",
      "\n",
      "# Make a bedGraph from all lowFDR CLIPper reads\n",
      "bedGraphCLIPout=makeBedGraph(CLIPPERlowFDR)\n",
      "CLIPPERlowFDRcenters=getBedCenterPoints(CLIPPERlowFDR)\n",
      "allLowFDRCentersBedGraph=makeBedGraph(CLIPPERlowFDRcenters)\n",
      "\n",
      "# --- Partition reads ---\n",
      "print \"Partition reads by type\"\n",
      "# Get names of lowFDR genes of each gene type\n",
      "pathToGeneLists=getLowFDRGeneTypes(CLIPpeGeneList)\n",
      "# Extract all lowFDR reads of each gene type\n",
      "pathToReadLists=getLowFDRReadTypes(CLIPPERlowFDR,pathToGeneLists)\n",
      "# BedGraph for protein coding\n",
      "proteinCodingReads = outfilepath+'clipGenes_proteinCoding_LowFDRreads.bed'\n",
      "proteinBedGraph=makeBedGraph(proteinCodingReads)\n",
      "\n",
      "\n",
      "'''\n",
      "filteredProteinCodingCenters=filterSnoRNAs(getBedCenterPoints(proteinCodingReads))\n",
      "filteredProteinCentersBedGraph=makeBedGraph(filteredProteinCodingCenters)\n",
      "# BedGraph for lincRNA\n",
      "lincRNAReads = outfilepath+'clipGenes_lincRNA_LowFDRreads.bed'\n",
      "filteredLincRNACenters=filterSnoRNAs(getBedCenterPoints(lincRNAReads))\n",
      "# Generate sorted gene lists\n",
      "pathToGeneLists=makeGeneLists(outfilepath,CLIPPERlowFDR)\n",
      "\n",
      "# --- Get intensity around cluster centers ---\n",
      "print \"Get binding intensity around cluster centers\"\n",
      "# Extract clipper input reads and cluster coordinates \n",
      "clipperIN=outfilepath+sampleName+\"_threshold=%s_allreads.mergedRT_CLIPPERin.bed\"%threshold\n",
      "CLIPclusters=outfilepath+sampleName+\"_threshold=%s_allreads.mergedRT_CLIP_clusters.bed\"%threshold\n",
      "# Make bedgraph file for clipper reads \n",
      "bedGraphCLIPin=makeBedGraph(clipperIN)\n",
      "# Extract the center coordinate of each clipper cluster \n",
      "centerCoordinates=makeClusterCenter(CLIPclusters) \n",
      "# Get histogram of read intensity around cluster center\n",
      "getClusterIntensity(bedGraphCLIPin,centerCoordinates)\n",
      "'''\n",
      "\n",
      "# --- Sort and print ---\n",
      "def makeGeneLists(outfilepath):\n",
      "    filesToSave=[]\n",
      "    readsPerGeneType=[]\n",
      "    geneTypeLabels=[]\n",
      "    for bedFile in glob.glob(outfilepath+\"*_LowFDRreads.bed\"):\n",
      "        head, tail = os.path.split(bedFile)\n",
      "        geneType=tail.strip().split('_')[1]\n",
      "        print geneType\n",
      "        geneTypeLabels=geneTypeLabels+[geneType]\n",
      "        # For protein coding genes ...\n",
      "        if geneType == \"proteinCoding\":\n",
      "            # *** No filtering *** \n",
      "            bedFile=outfilepath+'clipGenes_proteinCoding_LowFDRreads.bed'\n",
      "            geneNames=np.genfromtxt(bedFile,usecols=(9,),delimiter='\\t',dtype='string') # Note that column lookup is zero indexed\n",
      "            nameList=np.array([i.strip().split('_')[0] for i in geneNames],dtype='string')\n",
      "            # Count number of reads per gene name \n",
      "            try:\t\n",
      "                # Count the number of occourances of each gene name\n",
      "                count={}\n",
      "                for name in nameList:\n",
      "                    if count.has_key(name):\n",
      "                        count[name] += 1\n",
      "                    else:\n",
      "                        count[name] = 1\n",
      "                # Sort the dictionary by value\n",
      "                countSort=sorted(count.items(),key=lambda item: item[1])\n",
      "            except:\n",
      "                print \"No genes of this type.\"\n",
      "                \n",
      "            # Save the sorted gene names\n",
      "            tosave=outfilepath+'clipGenes_%s'%geneType\n",
      "            outfh = open(tosave, 'w')\n",
      "            for pair in countSort:\n",
      "                string=pair[0]+'\\t'+str(pair[1])+'\\n'\n",
      "                outfh.write(string)\n",
      "            outfh.close()\n",
      "\n",
      "makeGeneLists(outfilepath)\n",
      "\n",
      "'''\n",
      "# --- Sort reads by gene body ---\n",
      "print \"Intron and UTR analysis:\"\n",
      "fivePreads,notFivePreads,CDSreads,notCDSreads,threePreads,notThreePreads=extractUTRs(filteredProteinCodingCenters)\n",
      "UTRgenelists=[fivePreads,CDSreads,threePreads]\n",
      "saveSourceData(UTRgenelists)\n",
      "\n",
      "# --- Detailed mRNA body analysis --- \n",
      "print \"Gene body analysis:\"\n",
      "# Bed file with protein coding reads\n",
      "filteredProteinCoding = outfilepath+'clipGenes_proteinCoding_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved.bed'\n",
      "bedGraphProtein=makeBedGraph(filteredProteinCoding)\n",
      "# Generate gene-by-gene read coverage histogram \n",
      "makeAvgGraph(bedGraphProtein)\n",
      "\n",
      "# --- Generate coverage histograms for ncRNAs ---\n",
      "print \"Process non-coding RNAs:\"\n",
      "readPositionFiles=[]\n",
      "for geneTypeReads in pathToGeneLists:\n",
      "    # Do not perform this for protein coding or snoRNAs (performed using seperate gene reference file) \n",
      "    if 'snoRNA' not in geneTypeReads and 'proteinCoding' not in geneTypeReads:\n",
      "        readPos=getReadDist(geneTypeReads)\n",
      "        readPositionFiles=readPositionFiles+[readPos]\n",
      "partitionSnoRNAs(outfilepath)\n",
      "snoRNAfile = outfilepath + '/clipGenes_snoRNA_genePosition'\n",
      "fileList = [outfilepath+'/clipGenes_snoRNA_LowFDRreads_CDBOX_genePosition_snoRNA',outfilepath+'/clipGenes_snoRNA_LowFDRreads_HCAbox_genePosition_snoRNA',outfilepath+'/clipGenes_snoRNA_LowFDRreads_SCARNA_genePosition_snoRNA']\n",
      "fileCat(snoRNAfile,fileList)\n",
      "\n",
      "# --- Repeat RNA analysis ---\n",
      "print \"Reapat analysis\"\n",
      "negAndPosMerged = outfilepath + sampleName + 'repeatRNA_allreads.mergedRT'\n",
      "posMerged= outfilepath + sampleName + 'repeatRNA_positivereads.mergedRT'\n",
      "negMerged= outfilepath + sampleName + 'repeatRNA_negativereads.mergedRT'\n",
      "fileCat(negAndPosMerged,[posMerged,negMerged])\n",
      "repeatMapped=parseRepeatMapped(negAndPosMerged)\n",
      "\n",
      "# ---- Plots ---\n",
      "print \"Run plots\"\n",
      "runPlots(outfilepath,sampleName)\n",
      "'''"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Run mapping\n",
        "Performing Bowtie...\n",
        "Process mapped data\n",
        "Merge RT stops and run CLIPper"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Running CLIPper..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Partition reads by type"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Grabbing low FDR read by gene type..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Grep all reads for /arrayAhome/lmartin/CLIP/results/MMhur/clipGenes_proteinCoding\n",
        "proteinCoding"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 58,
       "text": [
        "'\\n# --- Sort reads by gene body ---\\nprint \"Intron and UTR analysis:\"\\nfivePreads,notFivePreads,CDSreads,notCDSreads,threePreads,notThreePreads=extractUTRs(filteredProteinCodingCenters)\\nUTRgenelists=[fivePreads,CDSreads,threePreads]\\nsaveSourceData(UTRgenelists)\\n\\n# --- Detailed mRNA body analysis --- \\nprint \"Gene body analysis:\"\\n# Bed file with protein coding reads\\nfilteredProteinCoding = outfilepath+\\'clipGenes_proteinCoding_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved.bed\\'\\nbedGraphProtein=makeBedGraph(filteredProteinCoding)\\n# Generate gene-by-gene read coverage histogram \\nmakeAvgGraph(bedGraphProtein)\\n\\n# --- Generate coverage histograms for ncRNAs ---\\nprint \"Process non-coding RNAs:\"\\nreadPositionFiles=[]\\nfor geneTypeReads in pathToGeneLists:\\n    # Do not perform this for protein coding or snoRNAs (performed using seperate gene reference file) \\n    if \\'snoRNA\\' not in geneTypeReads and \\'proteinCoding\\' not in geneTypeReads:\\n        readPos=getReadDist(geneTypeReads)\\n        readPositionFiles=readPositionFiles+[readPos]\\npartitionSnoRNAs(outfilepath)\\nsnoRNAfile = outfilepath + \\'/clipGenes_snoRNA_genePosition\\'\\nfileList = [outfilepath+\\'/clipGenes_snoRNA_LowFDRreads_CDBOX_genePosition_snoRNA\\',outfilepath+\\'/clipGenes_snoRNA_LowFDRreads_HCAbox_genePosition_snoRNA\\',outfilepath+\\'/clipGenes_snoRNA_LowFDRreads_SCARNA_genePosition_snoRNA\\']\\nfileCat(snoRNAfile,fileList)\\n\\n# --- Repeat RNA analysis ---\\nprint \"Reapat analysis\"\\nnegAndPosMerged = outfilepath + sampleName + \\'repeatRNA_allreads.mergedRT\\'\\nposMerged= outfilepath + sampleName + \\'repeatRNA_positivereads.mergedRT\\'\\nnegMerged= outfilepath + sampleName + \\'repeatRNA_negativereads.mergedRT\\'\\nfileCat(negAndPosMerged,[posMerged,negMerged])\\nrepeatMapped=parseRepeatMapped(negAndPosMerged)\\n\\n# ---- Plots ---\\nprint \"Run plots\"\\nrunPlots(outfilepath,sampleName)\\n'"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Funtions and modules - "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def makeGeneLists(outfilepath):\n",
      "    filesToSave=[]\n",
      "    readsPerGeneType=[]\n",
      "    geneTypeLabels=[]\n",
      "    for bedFile in glob.glob(outfilepath+\"*_LowFDRreads.bed\"):\n",
      "        head, tail = os.path.split(bedFile)\n",
      "        geneType=tail.strip().split('_')[1]\n",
      "        print geneType\n",
      "        geneTypeLabels=geneTypeLabels+[geneType]\n",
      "        # For protein coding genes ...\n",
      "        if geneType == \"proteinCoding\":\n",
      "            # *** No filtering *** \n",
      "            bedFile=outfilepath+'clipGenes_proteinCoding_LowFDRreads.bed'\n",
      "            geneNames=np.genfromtxt(bedFile,usecols=(9,),delimiter='\\t',dtype='string') # Note that column lookup is zero indexed\n",
      "            nameList=np.array([i.strip().split('_')[0] for i in geneNames],dtype='string')\n",
      "            # Count number of reads per gene name \n",
      "            try:\t\n",
      "                # Count the number of occourances of each gene name\n",
      "                count={}\n",
      "                for name in nameList:\n",
      "                    if count.has_key(name):\n",
      "                        count[name] += 1\n",
      "                    else:\n",
      "                        count[name] = 1\n",
      "                # Sort the dictionary by value\n",
      "                countSort=sorted(count.items(),key=lambda item: item[1])\n",
      "            except:\n",
      "                print \"No genes of this type.\"\n",
      "                \n",
      "            # Save the sorted gene names\n",
      "            tosave=outfilepath+'clipGenes_%s'%geneType\n",
      "            outfh = open(tosave, 'w')\n",
      "            for pair in countSort:\n",
      "                string=pair[0]+'\\t'+str(pair[1])+'\\n'\n",
      "                outfh.write(string)\n",
      "            outfh.close()\n",
      "\n",
      "makeGeneLists(outfilepath)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "proteinCoding\n"
       ]
      }
     ],
     "prompt_number": 55
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Modules - "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def compareLists(list1,list2,outname):\n",
      "    # Usage: Compare gene lists and output matches to the file. \n",
      "    # Input: Two gene lists.\n",
      "    # Output: Path file containing the matching genes.\n",
      "    # Set comparison, resulting in shared genes \n",
      "    f = open(list1, 'r')\n",
      "    g = open(list2, 'r')\n",
      "    content1 = set(f.readlines())\n",
      "    content2 = set(g.readlines())\n",
      "    commonGenes = content1 & content2\n",
      "    \n",
      "    # Write shared genes to an output file\n",
      "    geneCategory=outname.split('.')[1]\n",
      "    outputName=outfilepath+'clipGenes_'+geneCategory\n",
      "    outfh = open(outputName, 'w')\n",
      "    for gene in commonGenes:\n",
      "        outfh.write(gene)\n",
      "    outfh.close()\n",
      "    \n",
      "    return outputName\n",
      "    \n",
      "def getLowFDRGeneTypes(CLIPpeGeneList):\n",
      "    # Usage: Get all genes listed under each type, compare to CLIPer targets.\n",
      "    # Input: .bed file passed into CLIPper and the CLIPper windows file.\n",
      "    # Output: Path to file containing all CLIPper genes of each type.\n",
      "    readListByGeneType=[]\n",
      "    # Iterate through all gene types\n",
      "    # for geneType in os.listdir(os.getcwd() + '/docs/genes_types'):\n",
      "    for geneType in os.listdir(os.getcwd() + '/docs/mm9/genes_types'):\n",
      "        # Paths to gene lists and output directory\n",
      "        # genepath=os.getcwd() + '/docs/genes_types/' + geneType\n",
      "        genepath=os.getcwd() + '/docs/mm9/genes_types/' + geneType\n",
      "        # Compare list of lowFDR CLIP gene names and gene names by gene type\n",
      "        lowFDRreadlist=compareLists(CLIPpeGeneList,genepath,geneType)\n",
      "        # Update the list of paths to the resulting files\n",
      "        readListByGeneType=readListByGeneType+[lowFDRreadlist]\n",
      "    return readListByGeneType\n",
      "\n",
      "def getLowFDRReadTypes(CLIPPERlowFDR,pathToGeneLists):\n",
      "    # Usage: Given a list of genes, return all reads for the associated genes.\n",
      "    # Input: Gene list and the path to lowFDR read file.\n",
      "    # Output: List of reads assocaited with the given genes.\n",
      "    print \"Grabbing low FDR read by gene type...\"\n",
      "    lowFDRgenelist=[]\n",
      "    for path in pathToGeneLists:\n",
      "        # File path to which low FDR reads of each type will be written\n",
      "        outfile=path+'_LowFDRreads.bed'\n",
      "        \n",
      "        # Run file grep\n",
      "        print \"Grep all reads for %s\"%path\n",
      "        proc = subprocess.Popen('grep -F -f %s %s > %s'%(path,CLIPPERlowFDR,outfile),shell=True) # Run grep\n",
      "        return_code = proc.wait() # Wait for process to finish\n",
      "        lowFDRgenelist=lowFDRgenelist+[outfile]\n",
      "    \n",
      "    return lowFDRgenelist\n",
      "\n",
      "def filterSnoRNAs(proteinCodingReads):\n",
      "    # Usage: Filter snoRNA and miRNAs from protein coding reads.\n",
      "    # Input: .bed file with protein coding reads.\n",
      "    # Output: snoRNA and miR filtered .bed file.\n",
      "    program = 'intersectBed'\n",
      "    # *** MODIFIED 5' - 7 bp FOR DDX21; MAY BE A GENERAL EDIT ***\n",
      "    snoRNAindex = os.getcwd() + '/docs/snoRNA_reference/snoRNAmasker_formatted_5pExtend.bed'\n",
      "    miRNAindex = os.getcwd() + '/docs/miR_sort_clean.bed'\n",
      "    \n",
      "    # File name for low FDR reads\n",
      "    proteinWithoutsnoRNAs = proteinCodingReads.replace('.bed', '_snoRNAremoved.bed')\n",
      "    proteinWithoutmiRNAs = proteinWithoutsnoRNAs.replace('.bed', '_miRNAremoved.bed')\n",
      "    \n",
      "    # Intersect protein coding reads with snoRNA masker\n",
      "    outfh = open(proteinWithoutsnoRNAs, 'w')\n",
      "    proc = subprocess.Popen([program,'-a', proteinCodingReads,'-b', snoRNAindex,'-v','-s'],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    outfh.close()\n",
      "    \n",
      "    # Intersect protein coding reads with miR masker\n",
      "    outfh = open(proteinWithoutmiRNAs, 'w')\n",
      "    proc = subprocess.Popen([program,'-a', proteinWithoutsnoRNAs,'-b', miRNAindex,'-v','-s'],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    outfh.close()\n",
      "    return (proteinWithoutmiRNAs)\n",
      "\n",
      "def makeGeneLists(outfilepath,CLIPPERlowFDR):\n",
      "    # Usage: Make sorted gene lists for each RNA type, using specific reference file for snoRNAs and filtering miRNA, snoRNA from protein coding and lincRNAs.\n",
      "    # Input: Outfile path and low FDR reads.\n",
      "    # Output: Path to each file, number of reads assigned to each gene type, list of labels for the order gene types analyzed.\n",
      "    \n",
      "    filesToSave=[]\n",
      "    readsPerGeneType=[]\n",
      "    geneTypeLabels=[]\n",
      "    # Get all lowFDR files by each gene type\n",
      "    for bedFile in glob.glob(outfilepath+\"*_LowFDRreads.bed\"):\n",
      "        head, tail = os.path.split(bedFile)\n",
      "        geneType=tail.strip().split('_')[1]\n",
      "        print geneType\n",
      "        geneTypeLabels=geneTypeLabels+[geneType]\n",
      "        # For protein coding genes ...\n",
      "        if geneType == \"proteinCoding\":\n",
      "            # ... we will use the center coordinate file with snoRNAs filtered\n",
      "            bedFile=outfilepath+\"clipGenes_proteinCoding_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved.bed\"\n",
      "            geneNames=np.genfromtxt(bedFile,usecols=(3,),delimiter='\\t',dtype='string') # Note that column lookup is zero indexed\n",
      "            nameList=np.array([i.strip().split('_')[0] for i in geneNames],dtype='string')\n",
      "        # For lincRNAs genes (in which snoRNAs are contained) ...\n",
      "        elif geneType == \"lincRNA\":\n",
      "            # ... we will use the center coordinate file with snoRNAs filtered\n",
      "            bedFile=outfilepath+\"clipGenes_lincRNA_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved.bed\"\n",
      "            geneNames=np.genfromtxt(bedFile,usecols=(3,),delimiter='\\t',dtype='string') # Note that column lookup is zero indexed\n",
      "            nameList=np.array([i.strip().split('_')[0] for i in geneNames],dtype='string')\n",
      "        # For snoRNAs ...\n",
      "        elif geneType == \"snoRNA\":\n",
      "            # ... we will intersect centerpoints of all lowFDR reads with a reference file of well-annotated (accepted) snoRNAs\n",
      "            program = 'intersectBed'\n",
      "            snoRNAindex =  os.getcwd() + '/docs/snoRNA_reference/sno_coordinates_hg19_formatted.bed'\n",
      "            # This has the correct 6 column bedFile format\t\t\t\n",
      "            bedFile = outfilepath+'clipGenes_snoRNA_LowFDRreads.bed'\n",
      "            outfh = open(bedFile, 'w')\n",
      "            # This produces the primary bedGraph / bigWig for visualization\n",
      "            CLIPPERlowFDRcenters=getBedCenterPoints(CLIPPERlowFDR)\n",
      "            allLowFDRCentersBedGraph=makeBedGraph(CLIPPERlowFDRcenters)\n",
      "            # -u reports only one overlap for each read in -a, and -s reports strand\n",
      "            proc = subprocess.Popen([program,'-a', CLIPPERlowFDRcenters, '-b', snoRNAindex,'-wb'],stdout=outfh)\n",
      "            proc.communicate()\n",
      "            outfh.close()\t\n",
      "            snoRNAbedGraph=makeBedGraph(bedFile)\n",
      "            # Unique snoRNA ID is in column of the resulting file \n",
      "            nameList=np.genfromtxt(bedFile,usecols=(9,),delimiter='\\t',dtype='string')\n",
      "            # nameList=np.array([i.strip().split('.')[0] for i in geneNames],dtype='string')\n",
      "        else:\n",
      "            # All other lowFDR read files (extracted using Ensembl ID) have gene name in the ninth column \n",
      "            geneNames=np.genfromtxt(bedFile,usecols=(9,),delimiter='\\t',dtype='string')\n",
      "            nameList=np.array([i.strip().split('_')[0] for i in geneNames],dtype='string')\n",
      "        \n",
      "        # Count number of reads per gene name \n",
      "        # Count the number of occourances of each gene name\n",
      "        count={}\n",
      "        for name in nameList:\n",
      "            if count.has_key(name):\n",
      "                count[name] += 1\n",
      "            else:\n",
      "                count[name] = 1\n",
      "        # Sort the dictionary by value\n",
      "        countSort=sorted(count.items(),key=lambda item: item[1])\n",
      "        \n",
      "        # Save the sorted gene names\n",
      "        tosave=outfilepath+'clipGenes_%s'%geneType\n",
      "        outfh = open(tosave, 'w')\n",
      "        for pair in countSort:\n",
      "            string=pair[0]+'\\t'+str(pair[1])+'\\n'\n",
      "            outfh.write(string)\n",
      "        outfh.close()\n",
      "    \n",
      "        # Store files\n",
      "        filesToSave=filesToSave+[tosave]\n",
      "    \n",
      "        # Save list using \"SourceData\" handle for easy extraction \n",
      "        tosave=outfilepath+'SourceData_geneList_%s'%geneType\n",
      "        outfh = open(tosave, 'w')\n",
      "        for pair in countSort:\n",
      "            string=pair[0]+'\\t'+str(pair[1])+'\\n'\n",
      "            outfh.write(string)\n",
      "        outfh.close()\n",
      "    \n",
      "    return(filesToSave)\n",
      "\n",
      "print \"Filter and sort gene lists:\"\n",
      "\n",
      "'''\n",
      "# Get names of lowFDR genes of each gene type\n",
      "pathToGeneLists=getLowFDRGeneTypes(CLIPpeGeneList)\n",
      "# Extract all lowFDR reads of each gene type\n",
      "pathToReadLists=getLowFDRReadTypes(CLIPPERlowFDR,pathToGeneLists)\n",
      "# BedGraph for protein coding\n",
      "proteinCodingReads = outfilepath+'clipGenes_proteinCoding_LowFDRreads.bed'\n",
      "proteinBedGraph=makeBedGraph(proteinCodingReads)\n",
      "filteredProteinCodingCenters=filterSnoRNAs(getBedCenterPoints(proteinCodingReads))\n",
      "filteredProteinCentersBedGraph=makeBedGraph(filteredProteinCodingCenters)\n",
      "# BedGraph for lincRNA\n",
      "lincRNAReads = outfilepath+'clipGenes_lincRNA_LowFDRreads.bed'\n",
      "filteredLincRNACenters=filterSnoRNAs(getBedCenterPoints(lincRNAReads))\n",
      "# Generate sorted gene lists\n",
      "pathToGeneLists=makeGeneLists(outfilepath,CLIPPERlowFDR)\n",
      "'''"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Filter and sort gene lists:\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "\"\\n# Get names of lowFDR genes of each gene type\\npathToGeneLists=getLowFDRGeneTypes(CLIPpeGeneList)\\n# Extract all lowFDR reads of each gene type\\npathToReadLists=getLowFDRReadTypes(CLIPPERlowFDR,pathToGeneLists)\\n# BedGraph for protein coding\\nproteinCodingReads = outfilepath+'clipGenes_proteinCoding_LowFDRreads.bed'\\nproteinBedGraph=makeBedGraph(proteinCodingReads)\\nfilteredProteinCodingCenters=filterSnoRNAs(getBedCenterPoints(proteinCodingReads))\\nfilteredProteinCentersBedGraph=makeBedGraph(filteredProteinCodingCenters)\\n# BedGraph for lincRNA\\nlincRNAReads = outfilepath+'clipGenes_lincRNA_LowFDRreads.bed'\\nfilteredLincRNACenters=filterSnoRNAs(getBedCenterPoints(lincRNAReads))\\n# Generate sorted gene lists\\npathToGeneLists=makeGeneLists(outfilepath,CLIPPERlowFDR)\\n\""
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def makeClusterCenter(windowsFile):\n",
      "    # Usage: Generate a file of cluster centers.\n",
      "    # Input: Raw CLIPper output file.\n",
      "    # Output: File with coordinates for the center of each CLIPper cluster.\n",
      "    cleanBed = cleanBedFile(windowsFile)\n",
      "    centers=cleanBed.replace('.bed','.clusterCenter')\n",
      "    # Open new files\n",
      "    f = open(centers, 'w')\n",
      "    with open(cleanBed, 'r') as infile:\n",
      "        for line in infile:\n",
      "            # Each line from the CLIPper cluster file\n",
      "            elementList = line.strip().split('\\t')\n",
      "            # Lower coordinate in order to correct for differences in strand\n",
      "            minCoordinate=min(int(elementList[1]), int(elementList[2]))\n",
      "            diff=abs(int((int(elementList[1])-int(elementList[2]))/2))\n",
      "            # Write the center of each window\n",
      "            f.write(elementList[0] + '\\t' + str(minCoordinate+diff) + '\\t' + str(minCoordinate+diff+1) + '\\n')\n",
      "    f.close()\n",
      "    return centers\n",
      "\n",
      "def getClusterIntensity(bedGraph,centerCoordinates):\n",
      "    # Usage: Generate a matrix of read itensity values around CLIPper cluster center.\n",
      "    # Input: BedGraph and cluster center file.\n",
      "    # Output: Generates a matrix, which is passed into R.\n",
      "    program = os.getcwd() + '/bin/grep_chip-seq_intensity.pl'\n",
      "    program2 = 'wait'\n",
      "    proc = subprocess.Popen(['perl',program, centerCoordinates, bedGraph],)\n",
      "    proc.communicate()\n",
      "    print \"Waiting for Cluster Intensity file completion...\"\n",
      "    proc2 = subprocess.Popen(program2,shell=True)\n",
      "    proc2.communicate()\n",
      "\n",
      "'''\n",
      "# --- Get intensity around cluster centers ---\n",
      "print \"Get binding intensity around cluster centers\"\n",
      "# Extract clipper input reads and cluster coordinates \n",
      "clipperIN=outfilepath+sampleName+\"_threshold=%s_allreads.mergedRT_CLIPPERin.bed\"%threshold\n",
      "CLIPclusters=outfilepath+sampleName+\"_threshold=%s_allreads.mergedRT_CLIP_clusters.bed\"%threshold\n",
      "# Make bedgraph file for clipper reads \n",
      "bedGraphCLIPin=makeBedGraph(clipperIN)\n",
      "# Extract the center coordinate of each clipper cluster \n",
      "centerCoordinates=makeClusterCenter(CLIPclusters) \n",
      "# Get histogram of read intensity around cluster center\n",
      "getClusterIntensity(bedGraphCLIPin,centerCoordinates)\n",
      "'''"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Get binding intensity around cluster centers\n",
        "Waiting for Cluster Intensity file completion..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 168
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def extractUTRs(bedIn):\n",
      "    # Usage: Extract all UTR specific reads from the input file.\n",
      "    # Input: .bed file\n",
      "    # Output: Mutually exclusive partitions of the input file.\n",
      "    program = 'intersectBed'\n",
      "    fivePUTRBed = os.getcwd() + '/docs/5pUTRs_Ensbl_sort_clean_uniq.bed'\n",
      "    threePUTRBed = os.getcwd() + '/docs/3pUTRs_Ensbl_sort_clean_uniq.bed'\n",
      "    cdsBed = os.getcwd() + '/docs/Exons_Ensbl_sort_clean_uniq.bed'\n",
      "    \n",
      "    # Extract 5p reads and NOT 5p reads \n",
      "    fivePreads = bedIn.replace('.bed', '_5p.bed')\n",
      "    notFivePreads = bedIn.replace('.bed', '_NOT5p.bed')\n",
      "    outfh = open(fivePreads, 'w')\n",
      "    proc = subprocess.Popen([program, '-a', bedIn, '-b', fivePUTRBed,'-u','-s'],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    outfh.close()\n",
      "    outfh = open(notFivePreads, 'w')\n",
      "    proc = subprocess.Popen([program, '-a', bedIn, '-b', fivePUTRBed,'-v','-s'],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    outfh.close()\n",
      "    \n",
      "    # Extract 3p UTR reads and NOT 3pUTR\n",
      "    threePreads = bedIn.replace('.bed', '_3p.bed')\n",
      "    notThreePreads = bedIn.replace('.bed', '_NOT3p.bed')\n",
      "    outfh = open(threePreads, 'w')\n",
      "    proc = subprocess.Popen([program, '-a', notFivePreads, '-b', threePUTRBed,'-u','-s'],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    outfh.close()\n",
      "    outfh = open(notThreePreads, 'w')\n",
      "    proc = subprocess.Popen([program, '-a', notFivePreads, '-b', threePUTRBed,'-v','-s'],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    \n",
      "    # Extract CDS reads and NOT CDS reads \n",
      "    CDSreads = bedIn.replace('.bed', '_cds.bed')\n",
      "    notCDSreads = bedIn.replace('.bed', '_NOTcds.bed')\n",
      "    outfh = open(CDSreads, 'w')\n",
      "    proc = subprocess.Popen([program, '-a', notThreePreads, '-b', cdsBed,'-u','-s'],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    outfh.close()\n",
      "    outfh = open(notCDSreads, 'w')\n",
      "    proc = subprocess.Popen([program, '-a', notThreePreads, '-b', cdsBed,'-v','-s'],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    outfh.close()\n",
      "    \n",
      "    outfh.close()\n",
      "    return (fivePreads,notFivePreads,CDSreads,notCDSreads,threePreads,notThreePreads)\n",
      "\n",
      "def saveSourceData(UTRgenelists):\n",
      "    # Usage: Save source data for UTR analysis \n",
      "    # In: Path to lists of reads partitioned by type\n",
      "    # Out: None\n",
      "    \n",
      "    UTRlabelsTosave=['5p','cds','3p']\n",
      "    j=0\n",
      "    for UTRgenelist in UTRgenelists:\n",
      "        # Read all names\n",
      "        geneNames=np.genfromtxt(UTRgenelist,usecols=(3,),delimiter='\\t',dtype='string') # Note that column lookup is zero indexed\n",
      "        nameList=np.array([i.strip().split('_')[0] for i in geneNames],dtype='string')\n",
      "    \n",
      "        # Count reads per name\n",
      "        count={}\n",
      "        for name in nameList:\n",
      "            if count.has_key(name):\n",
      "                count[name] += 1\n",
      "            else:\n",
      "                count[name] = 1\n",
      "        # Sort the dictionary by value\n",
      "        countSort=sorted(count.items(),key=lambda item: item[1])\n",
      "    \n",
      "        # Save the sorted gene names\n",
      "        tosave=outfilepath+'SourceData_ReadsByUTR_%s' %UTRlabelsTosave[j]\n",
      "        outfh = open(tosave, 'w')\n",
      "        for pair in countSort:\n",
      "    \n",
      "            string=pair[0]+'\\t'+str(pair[1])+'\\n'\n",
      "            outfh.write(string)\n",
      "        outfh.close()\n",
      "        j += 1\n",
      "\n",
      "'''\n",
      "# --- Sort reads by gene body ---\n",
      "print \"Intron and UTR analysis:\"\n",
      "fivePreads,notFivePreads,CDSreads,notCDSreads,threePreads,notThreePreads=extractUTRs(filteredProteinCodingCenters)\n",
      "UTRgenelists=[fivePreads,CDSreads,threePreads]\n",
      "saveSourceData(UTRgenelists)\n",
      "'''"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Intron and UTR analysis:\n"
       ]
      }
     ],
     "prompt_number": 149
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def makeTab(bedGraph):\n",
      "    program = os.getcwd() + '/bin/bedGraph2tab.pl'\n",
      "    program2 = 'wait'\n",
      "    genesFile = os.getcwd() + '/docs/hg19_ensembl_genes.txt'\n",
      "    sizesFile = os.getcwd() + '/docs/hg19.sizes'\n",
      "    \n",
      "    outfile=bedGraph.replace('.bedgraph','.tab')\n",
      "    print \"Waiting for Tabfile completion...\"\n",
      "    proc = subprocess.Popen(['perl',program,genesFile,sizesFile,bedGraph,outfile],)\n",
      "    proc.communicate()\n",
      "    \n",
      "    proc2 = subprocess.Popen(program2,shell=True)\n",
      "    proc2.communicate()\n",
      "    \n",
      "    return outfile\n",
      "\n",
      "def makeAvgGraph(bedGraph):\n",
      "    # Usage: Generate a matrix of read itensity values across gene body.\n",
      "    # Input: BedGraph.\n",
      "    # Output: Generates two matricies, which are passed into R.\n",
      "    program= os.getcwd() + '/bin/averageGraph_scaled_tab.pl'\n",
      "    program2 = 'wait'\n",
      "    utrFile = os.getcwd() + '/docs/hg19_ensembl_UTR_annotation.txt'\n",
      "    print \"Make average graph...\"\n",
      "    tabFile=makeTab(bedGraph)\n",
      "    outhandle=tabFile.replace('.tab','_UTRs')\n",
      "    proc = subprocess.Popen(['perl',program,utrFile,tabFile,tabFile,outhandle],)\n",
      "    proc.communicate()\n",
      "    # Perl send this to a background process, so wait for completion before going ahead.\n",
      "    print \"Waiting for AverageGraph completion...\"\n",
      "    proc2 = subprocess.Popen(program2,shell=True)\n",
      "    proc2.communicate()\n",
      "    \n",
      "'''\n",
      "# --- Detailed mRNA body analysis --- \n",
      "print \"Gene body analysis:\"\n",
      "# Bed file with protein coding reads\n",
      "filteredProteinCoding = outfilepath+'clipGenes_proteinCoding_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved.bed'\n",
      "bedGraphProtein=makeBedGraph(filteredProteinCoding)\n",
      "# Generate gene-by-gene read coverage histogram \n",
      "makeAvgGraph(bedGraphProtein)\n",
      "'''"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Gene body analysis:\n",
        "Make average graph..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Waiting for Tabfile completion...\n",
        "Waiting for AverageGraph completion..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 150
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def getReadDist(geneList):\n",
      "\t# Usage: Extract raw reads and associated gene start/stop associated with each Ensembl gene ID\n",
      "\t# Input: Gene list\n",
      "\t# Output: Bed file with read IDs and assocaited start/stop for the gene\n",
      "\tgeneStartStopRepo = os.getcwd() + '/docs/all_genes.txt'\n",
      "\t# Reads associated with each gene type\n",
      "\tgeneReadList=geneList+'_LowFDRreads.bed'\n",
      "\t# Ensure that we read genes from the filtered lncRNAs\n",
      "\tif 'lincRNA' in geneList:\n",
      "\t\tgeneReadList=outfilepath+'/clipGenes_lincRNA_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved.bed'\n",
      "\treadsByGenePosition=geneList+'_genePosition'\n",
      "\t# Open the output file\n",
      "\toutfh = open(readsByGenePosition, 'w')\n",
      "\twith open(geneList, 'r') as infile:\n",
      "\t\tfor ensemblID in infile:\n",
      "\t\t\t# Gene name\n",
      "\t\t\tensemblID=ensemblID.strip().split('\\t')[0]\n",
      "\t\t\t# Extract basic gene data (from the first of the list, if there are repeats)\n",
      "\t\t\tstorePosition=grep(ensemblID,geneStartStopRepo)\n",
      "\t\t\tif storePosition:\n",
      "\t\t\t\t# Parse information about the gene\n",
      "\t\t\t\tID,chrom,geneStart,geneEnd=storePosition[0].strip().split('\\t')\n",
      "\t\t\t\t# Then extract all raw reads with the gene ID and output data for each read\n",
      "\t\t\t\tresults=grep(ensemblID,geneReadList)\n",
      "\t\t\t\t# Check to make sure reads were recovered\n",
      "\t\t\t\tif results:\n",
      "\t\t\t\t\tfor read in results:\n",
      "\t\t\t\t\t\treadData=read.strip().split('\\t')\n",
      "\t\t\t\t\t\toutfh.write('\\t'.join((ID,geneStart,geneEnd,readData[1],readData[2],'\\n')))\n",
      "\treturn readsByGenePosition\n",
      "\n",
      "def grep(pattern,filein):\n",
      "\t# Usage: Open a file and search all lines for a pattern.\n",
      "\t# Input: Pattern to search (gene name) and file name.\n",
      "\t# Output: List of lines in a file that have the pattern.\n",
      "\tr = []\n",
      "\tfilein_open = open(filein, 'r')\n",
      "\tfor line in filein_open:\n",
      "\t \tif re.search(pattern,line):\n",
      "\t \t\tr.append(line)\n",
      "\tfilein_open.close()\n",
      "\treturn r\n",
      "\n",
      "def partitionSnoRNAs(outfilepath):\n",
      "\t# Usage: Plot abundance of reads that map to each snoRNA type and coverage histograms for each.\n",
      "\t# Input: Output file path.\n",
      "\t# Output: Figure.\n",
      "\n",
      "\tprogram='intersectBed'\n",
      "\tsnoRNAindex = os.getcwd() + '/docs/snoRNA_reference/sno_coordinates_hg19_formatted.bed'\n",
      "\t# File name for low FDR snoRNA reads\n",
      "\tsnoRNAreads = outfilepath+'clipGenes_snoRNA_LowFDRreads.bed'\n",
      "\t# snoRNA read types \n",
      "\thcBOX_READS = outfilepath+'clipGenes_snoRNA_LowFDRreads_HCAbox.bed'\n",
      "\tSCA_READS = outfilepath+'clipGenes_snoRNA_LowFDRreads_SCARNA.bed'\n",
      "\tCDBOX_READS = outfilepath+'clipGenes_snoRNA_LowFDRreads_CDBOX.bed'\n",
      "\toutfh1 = open(hcBOX_READS, 'w')\n",
      "\toutfh2 = open(SCA_READS, 'w')\n",
      "\toutfh3 = open(CDBOX_READS,'w')\n",
      "\n",
      "\t# Read file into independent files for each snoRNA type\n",
      "\tscaRNAs=0\n",
      "\thcaBox=0\n",
      "\tcdBox=0\n",
      "\twith open(snoRNAreads, 'r') as infile:\n",
      "\t\tfor line in infile:\t\n",
      "\t\t\telementList = line.strip().split('\\t')\n",
      "\t\t\tsnoType = elementList[10]\n",
      "\t\t\telementList.append('\\n')\n",
      "\t\t\tif snoType == 'C':\n",
      "\t\t\t\toutfh3.write('\\t'.join(elementList))\n",
      "\t\t\t\tcdBox += 1\n",
      "\t\t\telif snoType == 'H':\n",
      "\t\t\t\toutfh1.write('\\t'.join(elementList))\n",
      "\t\t\t\thcaBox += 1\n",
      "\t\t\telif snoType == 's':\n",
      "\t\t\t\toutfh2.write('\\t'.join(elementList))\n",
      "\t\t\t\tscaRNAs += 1\n",
      "\t# Close the input files before moving forward and counting reads\n",
      "\toutfh1.close()\n",
      "\toutfh2.close()\n",
      "\toutfh3.close()\n",
      "\n",
      "\ttotalCount=float(scaRNAs+hcaBox+cdBox)\n",
      "\t\n",
      "\t# Check to make sure snoRNA reads recovered\n",
      "\tif totalCount > 0:\n",
      "\n",
      "\t\t# Make a dictionary of gene start,stop for each unique snoRNA ID\n",
      "\t\td={}\n",
      "\t\tsnoStartStopRepo =  os.getcwd() + '/docs/snoRNA_reference/sno_coordinates_hg19_formatted.bed'\n",
      "\t\twith open(snoStartStopRepo, 'r') as infile:\n",
      "\t\t\tfor gene in infile:\t\n",
      "\t\t\t\tchrom,geneStart,geneEnd,name,snoType,strand=gene.strip().split('\\t')\n",
      "\t\t\t\td[name]=[geneStart,geneEnd]\n",
      "\t\t\tinfile.close()\n",
      "\t\t\n",
      "\t\t# Read each snoRNA file\n",
      "\t\tsnoRNAfiles=[CDBOX_READS,hcBOX_READS,SCA_READS]\n",
      "\t\t# Store for position files\n",
      "\t\tpositionFiles=[]\n",
      "\t\tfor snoFile in snoRNAfiles:\n",
      "\t\t\t# Create a list of the output files\n",
      "\t\t\treadsByGenePosition=snoFile.replace('.bed', '_genePosition_snoRNA')\n",
      "\t\t\tpositionFiles=positionFiles+[readsByGenePosition]\n",
      "\t\t\t# Open the output file\n",
      "\t\t\toutfh = open(readsByGenePosition, 'w')\n",
      "\t\t\t# Read through the snoRNA file\n",
      "\t\t\tcount=0\n",
      "\t\t\twith open(snoFile, 'r') as infile:\n",
      "\t\t\t\tfor read in infile:\n",
      "\t\t\t\t\t# Extract the elements of each read\n",
      "\t\t\t\t\telementList=read.strip().split('\\t')\n",
      "\t\t\t\t\t# Extract the snoRNAID\n",
      "\t\t\t\t\tsnoID=elementList[9].strip()\n",
      "\t\t\t\t\t# Extract information for each gene\n",
      "\t\t\t\t\tgeneStartStop=d[snoID]\n",
      "\t\t\t\t\t# Write to the output file: start and end of snoRNA based upon reference file, and start and end of snoRNA read\n",
      "\t\t\t\t\toutfh.write('\\t'.join((snoID,geneStartStop[0],geneStartStop[1],elementList[1],elementList[2],'\\n')))\n",
      "\t\t\t\t\tcount += 1\n",
      "\t\t\tinfile.close()\n",
      "\t\t\toutfh.close()\n",
      "\n",
      "'''\n",
      "# --- Generate coverage histograms for ncRNAs ---\n",
      "print \"Process non-coding RNAs:\"\n",
      "readPositionFiles=[]\n",
      "for geneTypeReads in pathToGeneLists:\n",
      "    # Do not perform this for protein coding or snoRNAs (performed using seperate gene reference file) \n",
      "    if 'snoRNA' not in geneTypeReads and 'proteinCoding' not in geneTypeReads:\n",
      "        readPos=getReadDist(geneTypeReads)\n",
      "        readPositionFiles=readPositionFiles+[readPos]\n",
      "partitionSnoRNAs(outfilepath)\n",
      "snoRNAfile = outfilepath + '/clipGenes_snoRNA_genePosition'\n",
      "fileList = [outfilepath+'/clipGenes_snoRNA_LowFDRreads_CDBOX_genePosition_snoRNA',outfilepath+'/clipGenes_snoRNA_LowFDRreads_HCAbox_genePosition_snoRNA',outfilepath+'/clipGenes_snoRNA_LowFDRreads_SCARNA_genePosition_snoRNA']\n",
      "fileCat(snoRNAfile,fileList)\n",
      "'''"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Process non-coding RNAs:\n"
       ]
      }
     ],
     "prompt_number": 111
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def parseRepeatMapped(RTmerged):\n",
      "\t# Usage: Parse repeat bed file for each \n",
      "\t# In: bedFile with merged RT stops\n",
      "\t# Out: List of lists (each has RT stops assocaited with a single gene)\n",
      "\t\n",
      "\t# Load start coordinates\n",
      "\tstart=np.loadtxt(RTmerged,usecols=(1,),delimiter='\\t',dtype=int)\n",
      "\tstart[start < 0]=0 # Ensure no negative start coordinates and normalize \n",
      "\trtStops=start+15 # Regenerate the RT stop position: RTmerged is +/- 15 bases from the RT stop\n",
      "\trnu_u1=rtStops[rtStops < 165]\n",
      "\trnu_u2=rtStops[(rtStops > 164) & (rtStops < 353)]\n",
      "\trnu_u4_1=rtStops[(rtStops > 352) & (rtStops < 497)]\n",
      "\trnu_u4_atac=rtStops[(rtStops > 496) & (rtStops < 627)]\n",
      "\trnu_u5=rtStops[(rtStops > 626) & (rtStops < 744)]\n",
      "\trnu_u6=rtStops[(rtStops > 743) & (rtStops < 850)]\n",
      "\trnu_u6_atac=rtStops[(rtStops > 849) & (rtStops < 975)]\n",
      "\trnu_u7=rtStops[(rtStops > 974) & (rtStops < 1038)]\n",
      "\trnu_u11=rtStops[(rtStops > 1037) & (rtStops < 1173)]\n",
      "\trnu_u12=rtStops[(rtStops > 1172) & (rtStops < 1322)]\n",
      "\trn7_sl2=rtStops[(rtStops > 1321) & (rtStops < 1621)]\n",
      "\trn7_sk=rtStops[(rtStops > 1620) & (rtStops < 1953)]\n",
      "\trny_1=rtStops[(rtStops > 1952) & (rtStops < 2066)]\n",
      "\trny_3=rtStops[(rtStops > 2065) & (rtStops < 2168)]\n",
      "\trny_4=rtStops[(rtStops > 2167) & (rtStops < 2264)]\n",
      "\trny_5=rtStops[(rtStops > 2263) & (rtStops < 2348)]\n",
      "\tu3a=rtStops[(rtStops > 2347) & (rtStops < 2565)]\n",
      "\trna_5s=rtStops[(rtStops > 2564) & (rtStops < 2686)]\n",
      "\t\n",
      "\t# Define the end of the rRNA regions\n",
      "\tglobal rRNAstart\n",
      "\trRNAstart=2686\n",
      "\trDNA=rtStops[(rtStops > rRNAstart)]\n",
      "\treadLists=[rnu_u1,rnu_u2,rnu_u4_1,rnu_u4_atac,rnu_u5,rnu_u6,rnu_u6_atac,rnu_u7,rnu_u11,rnu_u12,rn7_sl2,rn7_sk,rny_1,rny_3,rny_4,rny_5,u3a,rna_5s,rDNA]\n",
      "\t\n",
      "\t# Save the files \n",
      "\tlabels=['U1','U2','U4_1','U4_ATAC','U5','U6','U6_ATAC','U7','U11','U12','7SL-SRP','7SK','Y_1','Y_3','Y_4','Y_5','U3A','5s','rDNA']\n",
      "\tfor i in range(len(readLists)):\n",
      "\t\trepeatData=outfilepath+'SourceData_RepeatRNA_%s'%labels[i]\n",
      "\t\tnp.savetxt(repeatData,readLists[i],delimiter='\\t',fmt=\"%s\")\n",
      "\n",
      "\treturn readLists\n",
      "\n",
      "'''\n",
      "# --- Repeat RNA analysis ---\n",
      "print \"Reapat analysis\"\n",
      "negAndPosMerged = outfilepath + sampleName + 'repeatRNA_allreads.mergedRT'\n",
      "repeatMapped=parseRepeatMapped(negAndPosMerged)\n",
      "'''"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "IOError",
       "evalue": "[Errno 2] No such file or directory: '/arrayAhome/lmartin/CLIP/results/ddx3wtT2/ddx3wtT2repeatRNA_allreads.mergedRT'",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mIOError\u001b[0m                                   Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-151-034e8a3511d7>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     44\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[1;34m\"Reapat analysis\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     45\u001b[0m \u001b[0mnegAndPosMerged\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0moutfilepath\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0msampleName\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m'repeatRNA_allreads.mergedRT'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 46\u001b[1;33m \u001b[0mrepeatMapped\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mparseRepeatMapped\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnegAndPosMerged\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[1;32m<ipython-input-151-034e8a3511d7>\u001b[0m in \u001b[0;36mparseRepeatMapped\u001b[1;34m(RTmerged)\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m         \u001b[1;31m# Load start coordinates\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m         \u001b[0mstart\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloadtxt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mRTmerged\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0musecols\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdelimiter\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'\\t'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mint\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      8\u001b[0m         \u001b[0mstart\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mstart\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m \u001b[1;31m# Ensure no negative start coordinates and normalize\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      9\u001b[0m         \u001b[0mrtStops\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mstart\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m15\u001b[0m \u001b[1;31m# Regenerate the RT stop position: RTmerged is +/- 15 bases from the RT stop\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;32m/usr/local/lib/python2.7/site-packages/numpy/lib/npyio.pyc\u001b[0m in \u001b[0;36mloadtxt\u001b[1;34m(fname, dtype, comments, delimiter, converters, skiprows, usecols, unpack, ndmin)\u001b[0m\n\u001b[0;32m    713\u001b[0m                 \u001b[0mfh\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0miter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbz2\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mBZ2File\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfname\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    714\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 715\u001b[1;33m                 \u001b[0mfh\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0miter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'U'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    716\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    717\u001b[0m             \u001b[0mfh\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0miter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfname\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mIOError\u001b[0m: [Errno 2] No such file or directory: '/arrayAhome/lmartin/CLIP/results/ddx3wtT2/ddx3wtT2repeatRNA_allreads.mergedRT'"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Reapat analysis\n"
       ]
      }
     ],
     "prompt_number": 151
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def lineCount(filename):\n",
      "\t# Usage: Get and return line count for a file.\n",
      "\t# Input: File \n",
      "\t# Outline: Lines\n",
      "\tprint \"Counting lines for ...\"\n",
      "\tprint filename\n",
      "\ti=0\n",
      "\twith open(filename) as f:\n",
      "\t\tfor i,l in enumerate(f):\n",
      "\t\t\tpass\n",
      "\tprint \"Lines counted:\"\n",
      "\tprint i\n",
      "\treturn i+1\n",
      "\n",
      "def plot_ReadAccounting(outfilepath,sampleName):\n",
      "\t# Usage: Make hbar chart of read count at specific steps in pipeline\n",
      "\t# Input: File path and sample name\n",
      "\t# Output: None\n",
      "\n",
      "\t# Paths to files for which we want to count reads\n",
      "\trawRead1=outfilepath+sampleName+'_R1.fastq'\n",
      "\trawRead2=outfilepath+sampleName+'_R2.fastq'\n",
      "\treads3pTrim=[outfilepath+sampleName+'_R1_3ptrimmed.fastq',outfilepath+sampleName+'_R2_3ptrimmed.fastq']\n",
      "\treadsFilter=[outfilepath+sampleName+'_R1_3ptrimmed_filter.fastq',outfilepath+sampleName+'_R2_3ptrimmed_filter.fastq']\n",
      "\treadsNoDupes=[outfilepath+sampleName+'_R1_3ptrimmed_filter_nodupe.fastq',outfilepath+sampleName+'_R2_3ptrimmed_filter_nodupe.fastq']\n",
      "\treadsMappedReapeat=[outfilepath+sampleName+'_R1_3ptrimmed_filter_nodupe_5ptrimmed_mappedToRepeatRNA.bed',outfilepath+sampleName+'_R2_3ptrimmed_filter_nodupe_5ptrimmed_mappedToRepeatRNA.bed']\n",
      "\treadsMappedHg19=[outfilepath+sampleName+'_R1_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes.bed',outfilepath+sampleName+'_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes.bed']\n",
      "\treadsMappedBlacklist=[outfilepath+sampleName+'_R1_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes.bed',outfilepath+sampleName+'_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes.bed']\n",
      "\treadsMappedRepeatMask=[outfilepath+sampleName+'_R1_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat.bed',outfilepath+sampleName+'_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat.bed']\n",
      "\tclipperIN=outfilepath+sampleName+'_threshold=%s_allreads.mergedRT_CLIPPERin.bed'%threshold\n",
      "\tclipperOUT=outfilepath+sampleName+'_threshold=%s_allreads.mergedRT_CLIP_clusters_lowFDRreads.bed'%threshold\n",
      "\t# File name labels for plotting\n",
      "\tfileNames=['Raw (R1)','Raw (R2)','3p Trim (R1)','3p Trim (R2)','Filter (R1)','Filter (R2)','No dupes (R1)','No dupes (R2)','RepeatMapped(R1)','RepeatMaped(R2)','Hg19Mapped (R1)','Hg19Mapped(R2)','Blacklist (R1)','Blacklist (R2)','RepeatMask(R1)','RepeatMask(R2)','ClipperIn','ClipperOut']\n",
      "\t# Determine read count for specified files\n",
      "\tfilesToCount=[rawRead1,rawRead2,reads3pTrim[0],reads3pTrim[1],readsFilter[0],readsFilter[1],readsNoDupes[0],readsNoDupes[1],readsMappedReapeat[0],readsMappedReapeat[1],readsMappedHg19[0],readsMappedHg19[1],readsMappedBlacklist[0],readsMappedBlacklist[1],readsMappedRepeatMask[0],readsMappedRepeatMask[1],clipperIN,clipperOUT]\n",
      "\tcounts=[]\n",
      "\tcounter=0\n",
      "\tfor fileString in filesToCount:\n",
      "\t\ttemp=lineCount(fileString)\n",
      "\t\tif counter < 8:\n",
      "\t\t\ttemp=temp/4  # Ensure that the read count for all .fastq files is corrected (1 read = 4 lines)\n",
      "\t\tcounts=counts+[temp]\n",
      "\t\tcounter += 1\n",
      "\t# Plot\n",
      "\t#fig1=plt.figure(1)\n",
      "\tplt.subplot(2,3,1)\n",
      "\tind = np.arange(len(counts)) + 0.5\n",
      "\tplt.barh(ind,list(reversed(counts)),align='center',color='blue')\n",
      "\tplt.xlabel('Read count per file',fontsize=5)\n",
      "\tplt.tick_params(axis='xticks',labelsize=5) \n",
      "\tplt.yticks(ind,list(reversed(fileNames)),fontsize=5)\n",
      "\tplt.tick_params(axis='yticks',labelsize=5) \n",
      "\t# Remove tick lines \n",
      "\tax=plt.gca()\n",
      "\tfor line in ax.get_yticklines():\n",
      "\t    line.set_markersize(0)\n",
      "\tplt.title('Read counts per step',fontsize=10)\n",
      "\t# Save plot data \n",
      "\treadDF=pd.DataFrame()\n",
      "\treadDF['File_name']=fileNames\n",
      "\treadDF['Reads_per_file']=counts\n",
      "\toutfilepath=outfilepath + 'SourceData_plot_ReadAccounting'\n",
      "\treadDF.to_csv(outfilepath)\n",
      "\n",
      "def plot_BoundGeneTypes(outfilepath,sampleName):\n",
      "\t# Usage: Quantify bound genes of select type\n",
      "\t# Input: Gene lists to plot, file path\n",
      "\t# Output: None\n",
      "\n",
      "\t# Gene lists\n",
      "\tgeneListToPlot=[outfilepath+'clipGenes_proteinCoding',outfilepath+'clipGenes_pseudogenes',outfilepath+'clipGenes_snoRNA',outfilepath+'clipGenes_lincRNA',outfilepath+'clipGenes_miR']\n",
      "\t# Count the number of genes in each list\n",
      "\tcounts=[]\n",
      "\tlabels=[]\n",
      "\tfor fileString in geneListToPlot:\n",
      "\t\thead, tail = os.path.split(fileString)\n",
      "\t\tgeneTypeName = tail.split(\"_\")[1]\n",
      "\t\tlabels=labels+[geneTypeName]\n",
      "\t\ttemp=lineCount(fileString)\n",
      "\t\tcounts=counts+[temp]\n",
      "\t# Sort the counts \n",
      "\tlabels=np.array(labels)\n",
      "\tcounts=np.array(counts)\n",
      "\tsortedIndex=[i[0] for i in sorted(enumerate(counts),key=lambda x:x[1])]\n",
      "\t# Plot\n",
      "\t#fig1=plt.figure(1)\n",
      "\tplt.subplot(2,3,6)\n",
      "\tind = np.arange(len(counts)) + 0.5\n",
      "\tplt.bar(ind,counts[sortedIndex],align='center',color='blue')\n",
      "\tlocs,pltlabels = plt.xticks(ind,labels[sortedIndex],fontsize=5)\n",
      "\tplt.setp(pltlabels, rotation=90, fontsize=5)\n",
      "\tplt.tick_params(axis='xticks',labelsize=5) \n",
      "\tax=plt.gca()\n",
      "\tfor line in ax.get_xticklines():\n",
      "\t\tline.set_markersize(0)\n",
      "\tplt.ylabel('Gene count',fontsize=5)\n",
      "\tplt.tick_params(axis='yticks',labelsize=5)\n",
      "\tplt.title('Read counts per step',fontsize=10)\n",
      "\tplt.title('Bound gene type and count',fontsize=10)\n",
      "\t# Save plot data \n",
      "\treadDF=pd.DataFrame()\n",
      "\treadDF['Gene_type']=labels[sortedIndex]\n",
      "\treadDF['Number_bound_genes']=counts[sortedIndex]\n",
      "\toutfilepath=outfilepath + 'SourceData_plot_BoundGeneTypes'\n",
      "\treadDF.to_csv(outfilepath)\n",
      "\n",
      "def plot_ReadsPerCluster(outfilepath,sampleName):\n",
      "\t# Usage: Make histogram of RT stops per cluster\n",
      "\t# Input: File path and sample name\n",
      "\t# Output: None\n",
      "\n",
      "\t# Extract file containing the number of reads per cluster\n",
      "\tCLIPpeReadsPerCluster=outfilepath+sampleName+'_threshold=%s_allreads.mergedRT_CLIP_clusters.readsPerCluster'%threshold\n",
      "\tCLIPpeReadsPerClusterList=[]\n",
      "\twith open(CLIPpeReadsPerCluster) as fin:\n",
      "\t\tCLIPpeReadsPerClusterList = [int(line.strip()) for line in fin]\n",
      "\t# Make hisogram, with each bin (except last) including lower value, excluding upper\n",
      "\tbins=range(min(CLIPpeReadsPerClusterList)-10,max(CLIPpeReadsPerClusterList)+10,10)\n",
      "\thist,bins=np.histogram(CLIPpeReadsPerClusterList,bins=bins)\n",
      "\t# Plot\n",
      "\t#fig1=plt.figure(1)\n",
      "\tplt.subplot(2,3,2) \n",
      "\txlimit=100 # The maximum number of reads per cluster visualized \n",
      "\twidth=0.7*(bins[1]-bins[0]) # Bar width is 70% of bin size\n",
      "\tcenter=(bins[:-1] + bins[1:])/2 # The center position of each bin\n",
      "\tplt.bar(center, hist,align='center',width=width)\n",
      "\tplt.xlabel('Reads per cluster (bin=10)',fontsize=5)\n",
      "\tplt.ylabel('Frequency (RT stop count)',fontsize=5)\n",
      "\tplt.tick_params(axis='yticks',labelsize=5) \n",
      "\tplt.xlim(0,xlimit)\n",
      "\tplt.title('Reads per cluster',fontsize=10)\n",
      "\n",
      "def plot_ClusterSizes(outfilepath,sampleName):\n",
      "\t# Usage: Make boxplot of cluster lengths\n",
      "\t# Input: File path\n",
      "\t# Output: None\n",
      "\n",
      "\t# Cluster file returned from clipper\n",
      "\tclipClusters=outfilepath+sampleName+\"_threshold=%s_allreads.mergedRT_CLIP_clusters\"%threshold\n",
      "\t# Obtain lengths of each cluster \n",
      "\tclusterLengths=[]\n",
      "\twith open(clipClusters) as fin:\n",
      "\t\tfor line in fin:\n",
      "\t\t\ttry:\n",
      "\t\t\t\telements=line.strip().split('\\t')\n",
      "\t\t\t\tstart=int(elements[1])\n",
      "\t\t\t\tend=int(elements[2])\n",
      "\t\t\t\tlengthOfCluster=math.fabs(start-end)\n",
      "\t\t\t\tclusterLengths=clusterLengths+[lengthOfCluster]\n",
      "\t\t\texcept:\n",
      "\t\t\t\tprint \"Error in line syntax.\"\n",
      "\t\t\t\tprint line\n",
      "\n",
      "\t# Plot\n",
      "\t#fig1=plt.figure(1)\n",
      "\tplt.subplot(2,3,3)\n",
      "\tylimit=500 # Maximum cluster length to show \n",
      "\tplt.boxplot(clusterLengths)\n",
      "\tplt.tick_params(axis='x',labelbottom='off') \n",
      "\tax=plt.gca()\n",
      "\tfor line in ax.get_xticklines():\n",
      "\t\tline.set_markersize(0)\n",
      "\tplt.ylabel('Cluster length (bases)',fontsize=5)\n",
      "\tplt.tick_params(axis='yticks',labelsize=5)\n",
      "\tplt.ylim(0,ylimit)\n",
      "\tplt.title('Cluster size',fontsize=10)\n",
      "\n",
      "def plot_clusterBindingIntensity(outfilepath,sampleName):\n",
      "\t# Usage: Read distribution around cluster center point\n",
      "\t# Input: File path\n",
      "\t# Output: None \n",
      "\n",
      "\t# Get heatmap of read intensity around cluster centers\t\n",
      "\tclusterCenterHeatmap=outfilepath+sampleName+'_threshold=%s_allreads.mergedRT_CLIP_clusters_cleaned_sorted.clusterCenter_heatmap.txt'%threshold\n",
      "\theatmapCols=40\n",
      "\t# Skip the header row (with column numbers) and the first column (cluster coordinates)\n",
      "\theatmapData=np.loadtxt(clusterCenterHeatmap,skiprows=1,dtype='float',usecols=range(1,heatmapCols+1))\n",
      "\t# Compute the sum of each row, where each row is centered on the coordiantes of each cluster center\n",
      "\tclusterSums=heatmapData.sum(axis=1)\n",
      "\t# Sort by sum of read abundance per cluster\n",
      "\tsortedIndex=[i[0] for i in sorted(enumerate(clusterSums),key=lambda x:x[1])]\n",
      "\tsortedHeatMap=heatmapData[sortedIndex,:]\n",
      "\t# Plot\n",
      "\t#fig1=plt.figure(1)\n",
      "\tplt.subplot(2,3,4)\n",
      "\tylimit=sortedHeatMap.shape[0]\n",
      "\t# The final value in the matrix (largest, in this case) is plotted at the top of the heatmap\n",
      "\tp=plt.pcolormesh(sortedHeatMap,cmap='Blues')\n",
      "\tplt.colorbar(p)\n",
      "\tplt.xticks(range(sortedHeatMap.shape[1]))\n",
      "\tplt.tick_params(axis='x',labelbottom='off') \n",
      "\tplt.xlabel('Cluster position',fontsize=5)\n",
      "\tplt.ylim(0,ylimit)\n",
      "\tplt.ylabel('Cluster number',fontsize=5)\n",
      "\tplt.title('Read distribution',fontsize=10)\n",
      "\n",
      "def plot_readsBymRNAregion(outfilepath,sampleName):\n",
      "\tutrData=['5p','cds','3p'] \n",
      "\treadCounts=[]\n",
      "\tfor utr in utrData:\n",
      "\t\t# Get read data from source files (sorted gene list) and check against raw .bed file\n",
      "\t\tsource=outfilepath + 'SourceData_ReadsByUTR_%s'%utr\n",
      "\t\tgeneCount=np.genfromtxt(source,usecols=(1,),delimiter='\\t',dtype='int')\n",
      "\t\tbedSource=outfilepath + 'clipGenes_proteinCoding_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved_%s.bed'%utr\n",
      "\t\tbedFile=np.genfromtxt(bedSource,usecols=(3,),delimiter='\\t',dtype='string')\n",
      "\t\treadCounts=readCounts+[bedFile.size]\n",
      "\t# Plot\n",
      "\t# Pie chart based upon read based accounting\n",
      "\tallProteinCoding=outfilepath + 'clipGenes_proteinCoding_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved.bed'\n",
      "\tallProteinCoding=np.genfromtxt(allProteinCoding,usecols=(3,),delimiter='\\t',dtype='string')\n",
      "\ttotalReads=allProteinCoding.size\n",
      "\tnonIntronic=sum(readCounts) # Reads that fall into 5p UTR, CDS, 3p UTR\n",
      "\tvector=[float(totalReads-nonIntronic)/totalReads,float(readCounts[0])/totalReads,float(readCounts[1])/totalReads,float(readCounts[2])/totalReads]\n",
      "\tax=plt.subplot(2,3,5)\n",
      "\tpie_wedges=ax.pie(vector,labels=[\"Intronic\",\"5p UTR\",\"CDS\",\"3pUTR\"],labeldistance=1.1,autopct='%1.1f%%')\n",
      "\tplt.rcParams['font.size']=10\n",
      "\tfor wedge in pie_wedges[0]:\n",
      "\t\twedge.set_edgecolor('black')\n",
      "\t\twedge.set_lw(1)\n",
      "\n",
      "def plot_mRNAgeneBodyDist(outfilepath,sampleName):\n",
      "\t# Bed file with protein coding reads\n",
      "\tfilteredProteinCoding = outfilepath+'clipGenes_proteinCoding_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved.bed'\n",
      "\t# Generate gene-by-gene read coverage histogram \n",
      "\taverageGraph=outfilepath+'clipGenes_proteinCoding_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved_cleaned_sorted_UTRs_scaled_cds200_abt0_averageGraph.txt'\n",
      "\t# Number of columns \n",
      "\tavgGraphCols=600\n",
      "\tavgGraphData=np.loadtxt(averageGraph,skiprows=1,dtype='float',usecols=range(1,avgGraphCols+1))\n",
      "\tfig2=plt.figure(2)\n",
      "\tplt.subplot2grid((2,3),(0,0),colspan=3)\n",
      "\tylimit=max(avgGraphData[1,:])*1.1\n",
      "\tplt.plot(avgGraphData[1,:],color='blue',linewidth='2')\n",
      "\tplt.ylim(0,ylimit)\n",
      "\tplt.vlines(200,0,ylimit,linestyles='dashed')\n",
      "\tplt.vlines(400,0,ylimit,linestyles='dashed')\n",
      "\tplt.tick_params(axis='x',labelbottom='off') \n",
      "\tplt.xlabel('mRNA gene body (5pUTR, CDS, 3pUTR)')\n",
      "\tplt.ylabel('Abundance')\n",
      "\tplt.tick_params(axis='y',labelsize=5) \n",
      "\tplt.title('CLIP signal across average mRNA transcript.',fontsize=10)\n",
      "\n",
      "def convertENBLids(inNames):\n",
      "\t# Usage: Converl ENST to ENSG (unique ID) using ENSEMBL annotation file\n",
      "\t# Input: List of ENST IDs\n",
      "\t# Output: List of ENSG IDs\n",
      "\tgenesFile = os.getcwd() + '/docs/hg19_ensembl_genes.txt'\n",
      "\tensemblIDfile=np.genfromtxt(genesFile,usecols=(1,12,),delimiter='\\t',dtype='string') # Note that column lookup is zero indexed\n",
      "\ttemp=[]\n",
      "\tfor name in inNames:\n",
      "\t\toutName=ensemblIDfile[ensemblIDfile[:,0]==name,1]\n",
      "\t\ttemp=temp+[outName]\n",
      "\ttemp=np.array(temp)\n",
      "\treturn temp\n",
      "\n",
      "def extractGenesByUTR(utrType,treatMatrixData,geneNames):\n",
      "\t# Usage: Extract genes with exclusive binding \n",
      "\t# Input: UTR type, data, names\n",
      "\t# Output: Extracted and sorted data, names\n",
      "\t# Extract specific regions\n",
      "\tif utrType=='5p':\n",
      "\t\tindexer=(treatMatrixData[:,range(200,600)].sum(axis=1) == 0) & (treatMatrixData[:,range(0,200)].sum(axis=1) > 0)\n",
      "\telif utrType=='3p':\n",
      "\t\tindexer=(treatMatrixData[:,range(0,400)].sum(axis=1) == 0) & (treatMatrixData[:,range(400,600)].sum(axis=1) > 0)\n",
      "\telse:\n",
      "\t\tindexer=(treatMatrixData[:,range(0,200)].sum(axis=1) == 0) & (treatMatrixData[:,range(400,600)].sum(axis=1) == 0) & (treatMatrixData[:,range(200,400)].sum(axis=1) > 0)\n",
      "\t# Extract data    \n",
      "\textractedData=treatMatrixData[indexer,:]\n",
      "\textractedNames=geneNames[indexer]\n",
      "\t# Get total abundance / gene\n",
      "\tgeneSums=extractedData.sum(axis=1)  \n",
      "\textractedDataSort=extractedData[np.argsort(geneSums),:]\n",
      "\textractedNamesSort=extractedNames[np.argsort(geneSums)]\n",
      "\treturn (extractedDataSort,extractedNamesSort)\n",
      "\n",
      "def plot_geneBodyPartition(outfilepath,sampleName):\n",
      "\t# Usage: \n",
      "\t# In: Outfilepath\n",
      "\t# Out: Source data, heatmaps\n",
      "\n",
      "\t# Binding across gene body (5p UTR, CDS, 3p UTR) of protein coding genes\n",
      "\ttreatMatrixCols=600\n",
      "\ttreatMatrix=outfilepath+'clipGenes_proteinCoding_LowFDRreads_centerCoord_snoRNAremoved_miRNAremoved_cleaned_sorted_UTRs_scaled_cds200_abt0_treatmatrix.txt'\n",
      "\ttreatMatrixData=np.genfromtxt(treatMatrix,skip_header=1,usecols=range(1,treatMatrixCols+1),delimiter='\\t',dtype='float')\n",
      "\tgeneNames=np.loadtxt(treatMatrix,dtype='string',skiprows=1,usecols=(0,),delimiter='\\t')\n",
      "\t# Convert to ENSG IDs\n",
      "\tgeneNames=convertENBLids(geneNames)\n",
      "\t# Ensure that genes are in the CLIPper protein coding list (treat analysis can return some spurious genes)\n",
      "\tmasterList = outfilepath+'clipGenes_proteinCoding'\n",
      "\tmasterNames = np.genfromtxt(masterList,usecols=(0,),delimiter='\\t',dtype='string')\n",
      "\tindexer=[]\n",
      "\tfor geneName in geneNames:\n",
      "\t\tif geneName in masterNames:\n",
      "\t\t\tindexer=indexer+[1]\n",
      "\t\telse:\n",
      "\t\t\tindexer=indexer+[0]\n",
      "\tindexer=np.array(indexer,dtype=bool)\n",
      "\tgeneNames=geneNames[indexer]\n",
      "\ttreatMatrixData=treatMatrixData[indexer,:]\n",
      "\t# Isolate genes and treatmatrix data by region\n",
      "\ttry:\n",
      "\t\tdata_5p,names_5p=extractGenesByUTR('5p',treatMatrixData,geneNames)\n",
      "\texcept:\n",
      "\t\tprint \"No 5' exclusive reads.\"\n",
      "\t\tdata_5p=[]\n",
      "\t\tnames_5p=[]\n",
      "\ttry:\n",
      "\t\tdata_cds,names_cds=extractGenesByUTR('cds',treatMatrixData,geneNames)\n",
      "\texcept:\n",
      "\t\tprint \"No cds exclusive reads.\"\n",
      "\t\tdata_cds=[]\n",
      "\t\tnames_cds=[]\n",
      "\ttry:   \n",
      "\t\tdata_3p,names_3p=extractGenesByUTR('3p',treatMatrixData,geneNames)\n",
      "\texcept:\n",
      "\t\tprint \"No CDS' exclusive reads.\"\n",
      "\t\tdata_3p=[]\n",
      "\t\tnames_3p=[]\n",
      "\n",
      "\t# All protein coding genes\n",
      "\tmasterList = outfilepath+'clipGenes_proteinCoding'\n",
      "\tmasterNames = np.genfromtxt(masterList,usecols=(0,),delimiter='\\t',dtype='string')\n",
      "\t# All 5p UTR, CDS, 3p UTR genes\n",
      "\tallTreatMatrixGenes=list(geneNames[:,0])\n",
      "\t# Intronic gene exctraction\n",
      "\ttosave=outfilepath + 'SourceData_geneBodyPartition_ExclusiveIntronic' \n",
      "\tintronicBoundGenes=list(set(masterNames)-set(allTreatMatrixGenes))\n",
      "\tnp.savetxt(tosave,np.array(intronicBoundGenes),fmt=\"%s\")\n",
      "\n",
      "\t# Check relative to previously identified genes\n",
      "\tutrBinding_Names_Fig1=['5p','cds','3p'] \n",
      "\tutrBinding_Names_Fig2=[names_5p[:,0],names_cds[:,0],names_3p[:,0]]\n",
      "\tutrBinding_Data_Fig2=[data_5p,data_cds,data_3p]\n",
      "\tdiffs=[]\n",
      "\n",
      "\tfor i in range(0,3):\n",
      "\t\t# Identified genes with 5p UTR, CDS, 3p UTR binding from read-based paritioning\n",
      "\t\tutr=utrBinding_Names_Fig1[i]\n",
      "\t\tsource=outfilepath + 'SourceData_ReadsByUTR_%s'%utr\n",
      "\t\ttemp1e=np.genfromtxt(source,usecols=(0,),delimiter='\\t',dtype=str)\n",
      "\t\t# Identified genes with 5p UTR, CDS, 3p UTR from gene body analysis\n",
      "\t\ttemp2b=np.unique(utrBinding_Names_Fig2[i]) \n",
      "\t\t# Non-overlaps are due to reads that are ambiguously mapped into different 5p UTR, CDS, 3p UTR region on isoforms\n",
      "\t\tdiff=list(set(temp2b)-set(temp1e))\n",
      "\t\tdiffs.append(diff)\n",
      "\n",
      "\tnamesFilter=[]\n",
      "\tdataFilter=[]\n",
      "\tgeneTypes=['5p','cds','3p'] \n",
      "\ttitles=['5pUTR','cds','3pUTR'] \n",
      "\tdepth=50\n",
      "\tfor i in range(0,3):\n",
      "\t\t# Names of bound genes\n",
      "\t\tnames=utrBinding_Names_Fig2[i]\n",
      "\t\tdata=utrBinding_Data_Fig2[i]\n",
      "\t\t# Indexer\n",
      "\t\tindexer=np.array([geneName not in diffs[i] for geneName in names])\n",
      "\t\ttempNames=names[indexer]\n",
      "\t\tnamesFilter.append(tempNames)\n",
      "\t\ttempData=data[indexer,:]\n",
      "\t\tdataFilter.append(tempData)\n",
      "\t\t# Save data\n",
      "\t\ttosave=outfilepath + 'SourceData_geneBodyPartition_Exclusive%s'%titles[i] \n",
      "\t\tnp.savetxt(tosave,np.array(tempData),fmt=\"%s\")\n",
      "\t\t# Plot\n",
      "\t\tplt.subplot2grid((2,3),(1,i),colspan=1)\n",
      "\t\tp=plt.pcolormesh(dataFilter[i][-depth:-1,:],cmap='Blues')\n",
      "\t\tplt.title(titles[i],fontsize=10)\n",
      "\t\tplt.vlines(200,0,depth,linestyles='dashed')\n",
      "\t\tplt.vlines(400,0,depth,linestyles='dashed')\n",
      "\t\tplt.tick_params(axis='x',labelbottom='off') \n",
      "\t\tplt.tick_params(axis='y',labelleft='off') \n",
      "\t\tplt.ylim(0,depth)\n",
      "\t\tplt.ylabel('Ranked genes (highest on bottom)',fontsize=5)\n",
      "\t\tplt.xticks(visible=False)\n",
      "\t\tplt.yticks(visible=False)\n",
      "\t\tplt.title('%s specific genes: %s'%(titles[i],np.unique(namesFilter[i]).size),fontsize=7.5)\n",
      "\n",
      "def plot_repeatRNA(outfilepath,sampleName):\n",
      "\trepeatGenomeBuild=os.getcwd() + '/docs/repeat/repeatRNA.fa'\n",
      "\trepeat_genome=np.genfromtxt(repeatGenomeBuild,dtype='string')\n",
      "\trepeat_genome_bases=repeat_genome[1]\n",
      "\trepeat_genome_size=len(repeat_genome[1])\n",
      "\t# Repeat index positions\n",
      "\trepeatAnnotation=os.getcwd() + '/docs/repeat/Hs_repeatIndex_positions.txt'\n",
      "\trepeatAnnotDF=pd.DataFrame(pd.read_table(repeatAnnotation,header=None))\n",
      "\trepeatAnnotDF.columns=['Name','Length','IndexStart','IndexEnd']\n",
      "\t# Python list extraction is not end index inclusive; to extract sequence, use end + 1.\n",
      "\trepeatAnnotDF['End_for_extraction']=repeatAnnotDF['IndexEnd']+1 \n",
      "\n",
      "\t# Get all repeat RNA files\n",
      "\trepeatRNA_RTstops=glob.glob(outfilepath+'SourceData_RepeatRNA_*')\n",
      "\tplotDim=math.ceil(math.sqrt(len(repeatRNA_RTstops)))\n",
      "\n",
      "\tfor i in range(len(repeatRNA_RTstops)):\n",
      "\t\t# RT stop file\n",
      "\t\tfilePath=repeatRNA_RTstops[i]\n",
      "\t\t# Extract RNA ID\n",
      "\t\tRNAid=filePath.split('RepeatRNA_')[1]\n",
      "\t\t# Get start,end for repeat RNA and extract the sequence\n",
      "\t\tcoords=np.array(repeatAnnotDF[repeatAnnotDF['Name']==RNAid][['IndexStart','End_for_extraction']])\n",
      "\t\tsequence=np.array(list(repeat_genome_bases[coords[0][0]:coords[0][1]]))\n",
      "\t\t# Create data frame for storage \n",
      "\t\tstorageDF=pd.DataFrame()\n",
      "\t\t# Get RT stops\n",
      "\t\tRTpositions=np.loadtxt(filePath,delimiter='\\t',dtype=int)\n",
      "\t\t# Histogram\n",
      "\t\tbins=range(coords[0][0],coords[0][1]+1,1) # Make sure bins are end coordinate inclusive\n",
      "\t\thist,bins=np.histogram(RTpositions,bins=bins)\n",
      "\t\twidth=0.7*(bins[1]-bins[0])\n",
      "\t\tcenter=(bins[:-1] + bins[1:])/2\n",
      "\t\t# Noramalize histogram frequencies by the total numer of RT stops \n",
      "\t\thistPlot=np.array(hist,dtype=float)\n",
      "\t\thistPlot=np.array(histPlot/float(len(RTpositions)),dtype=float)\n",
      "\t\t# Plot\n",
      "\t\tplt.subplot(plotDim,plotDim,i+1)\n",
      "\t\tplt.bar(center,histPlot,align='center',width=width,color='blue',alpha=0.45)\n",
      "\t\tplt.tick_params(axis='x',labelsize=2.5) \n",
      "\t\tplt.tick_params(axis='y',labelsize=2.5)  \n",
      "\t\tplt.title('RT stops for %s: %s'%(RNAid,len(RTpositions)),fontsize=5)\n",
      "\t\tplt.xlim(coords[0][0],coords[0][1])\n",
      "\t\t# Extract sequence and record counts\n",
      "\t\tstorageDF['Sequence']=sequence\n",
      "\t\treadsPerBase=np.array(list(hist))\n",
      "\t\treadsPerBaseNorm=np.array(list(histPlot))\n",
      "\t\tcolName=sampleName+'RT_stops'\n",
      "\t\tstorageDF[colName]=readsPerBase\n",
      "\t\tcolNameNorm=sampleName+'RT_stops_norm'\n",
      "\t\tstorageDF[colNameNorm]=readsPerBaseNorm\n",
      "\t\t# Layout\n",
      "\t\tplt.tight_layout()\n",
      "\n",
      "def plot_rDNA(outfilepath,sampleName):\n",
      "\trepeatGenomeBuild=os.getcwd() + '/docs/repeat/repeatRNA.fa'\n",
      "\trepeat_genome=np.genfromtxt(repeatGenomeBuild,dtype='string')\n",
      "\trepeat_genome_bases=repeat_genome[1]\n",
      "\trepeat_genome_size=len(repeat_genome[1])\n",
      "\t# Repeat index positions\n",
      "\trepeatAnnotation=os.getcwd() + '/docs/repeat/Hs_repeatIndex_positions.txt'\n",
      "\trepeatAnnotDF=pd.DataFrame(pd.read_table(repeatAnnotation,header=None))\n",
      "\trepeatAnnotDF.columns=['Name','Length','IndexStart','IndexEnd']\n",
      "\t# Python list extraction is not end index inclusive; to extract sequence, use end + 1.\n",
      "\trepeatAnnotDF['End_for_extraction']=repeatAnnotDF['IndexEnd']+1 \n",
      "\n",
      "\tax=plt.subplot2grid((3,3),(0,0),colspan=3)\n",
      "\tRNAid='rDNA'\n",
      "\trepeatRNA_RTstops=glob.glob(outfilepath+'SourceData_RepeatRNA_*')\n",
      "\tfilePath = [s for s in repeatRNA_RTstops if RNAid in s][0]\n",
      "\tcoords=np.array(repeatAnnotDF[repeatAnnotDF['Name']==RNAid][['IndexStart','End_for_extraction']])\n",
      "\trRNAstart=coords[0][0]\n",
      "\tsequence=np.array(list(repeat_genome_bases[coords[0][0]:coords[0][1]]))\n",
      "\t# Create data frame for storage \n",
      "\tstorageDF=pd.DataFrame()\n",
      "\t# Get RT stops\n",
      "\tRTpositions=np.loadtxt(filePath,delimiter='\\t',dtype=int)\n",
      "\t# Histogram\n",
      "\tbins=range(coords[0][0],coords[0][1]+1,1) # Make sure bins are end coordinate inclusive\n",
      "\thist,bins=np.histogram(RTpositions,bins=bins)\n",
      "\twidth=0.7*(bins[1]-bins[0])\n",
      "\tcenter=(bins[:-1] + bins[1:])/2\n",
      "\t# Noramalize histogram frequencies by the total numer of RT stops \n",
      "\thistPlot=np.array(hist,dtype=float)\n",
      "\thistPlot=np.array(histPlot/float(len(RTpositions)),dtype=float)\n",
      "\t# Plot\n",
      "\tplt.bar(center,histPlot,align='center',width=width,color='blue',alpha=0.45)\n",
      "\tplt.tick_params(axis='x',labelsize=2.5) \n",
      "\tplt.tick_params(axis='y',labelsize=2.5)  \n",
      "\tplt.title('RT stops for %s: %s'%(RNAid,len(RTpositions)),fontsize=5)\n",
      "\tplt.xlim(coords[0][0],coords[0][1])\n",
      "\t# Features of rDNA with respect to start of the bowtie index (index=0)\n",
      "\tstart18s=3657\n",
      "\tend18s=5527\n",
      "\tstart5s=6623\n",
      "\tend5s=6779\n",
      "\tstart28s=7935\n",
      "\tend28s=12969\n",
      "\t# Overlay regions on plot \n",
      "\tplt.axvspan(start18s+rRNAstart,end18s+rRNAstart,facecolor='g',alpha=0.5)\n",
      "\tplt.axvspan(start5s+rRNAstart,end5s+rRNAstart,facecolor='r',alpha=0.5)\n",
      "\tplt.axvspan(start28s+rRNAstart,end28s+rRNAstart,facecolor='b',alpha=0.5)\n",
      "\n",
      "\t# Generate histogram for transcribed region\n",
      "\tplt.subplot2grid((3,3),(1,0),colspan=3)\n",
      "\t# Generate histogram\n",
      "\tdatarDNAOnly=RTpositions-rRNAstart\n",
      "\tbins=range((coords[0][0]-rRNAstart),(coords[0][1]-rRNAstart+1),1) # Make sure bins are end coordinate inclusive\n",
      "\thist,bins=np.histogram(datarDNAOnly,bins=bins)\n",
      "\twidth=0.7*(bins[1]-bins[0])\n",
      "\tcenter=(bins[:-1] + bins[1:])/2\n",
      "\t# Noramalize histogram frequencies by the total numer of RT stops \n",
      "\thistPlot=np.array(hist,dtype=float)\n",
      "\thistPlot=np.array(histPlot/float(len(RTpositions)),dtype=float)\n",
      "\t# Plot\n",
      "\tplt.bar(center,histPlot,align='center',width=width,color='blue',alpha=0.45)\n",
      "\tplt.tick_params(axis='x',labelsize=2.5) \n",
      "\tplt.tick_params(axis='y',labelsize=2.5)  \n",
      "\tplt.xlabel('rRNA locus position (bin=1 base)',fontsize=5)\n",
      "\tplt.ylabel('Normalized RT stop / bin',fontsize=2.5)\n",
      "\t# Overlay regions\n",
      "\tplt.axvspan(start18s,end18s,facecolor='g',alpha=0.5)\n",
      "\tplt.axvspan(start5s,end5s,facecolor='r',alpha=0.5)\n",
      "\tplt.axvspan(start28s,end28s,facecolor='b',alpha=0.5)\n",
      "\t# Set x-axis to end position of rRNA locus\n",
      "\trRNAend=13314\n",
      "\tplt.xlim(0,rRNAend)\n",
      "\n",
      "\t# Individual regions \n",
      "\tplt.subplot2grid((3,3),(2,0),colspan=1)\n",
      "\tplt.bar(center,histPlot,align='center',width=width,color='green',alpha=0.75)\n",
      "\tplt.xlim(start18s,end18s)\n",
      "\tplt.xlabel('rRNA locus position (bin=1 base)',fontsize=5)\n",
      "\tplt.ylabel('Normalized RT stop / bin',fontsize=2.5)\n",
      "\tplt.tick_params(axis='x',labelsize=5) \n",
      "\tplt.tick_params(axis='y',labelsize=5) \n",
      "\tplt.title('18s Region',fontsize=10)\n",
      "\tplt.subplot2grid((3,3),(2,1),colspan=1)\n",
      "\tplt.bar(center,histPlot,align='center',width=width,color='red',alpha=0.75)\n",
      "\tplt.xlim(start5s,end5s)\n",
      "\tplt.xlabel('rRNA locus position (bin=1 base)',fontsize=5)\n",
      "\tplt.tick_params(axis='x',labelsize=5) \n",
      "\tplt.tick_params(axis='y',labelsize=5) \n",
      "\tplt.title('5.8s Region',fontsize=10)\n",
      "\tplt.subplot2grid((3,3),(2,2),colspan=1)\n",
      "\tplt.bar(center,histPlot,align='center',width=width,color='blue',alpha=0.75)\n",
      "\tplt.xlim(start28s,end28s)\n",
      "\tplt.xlabel('rRNA locus position (bin=1 base)',fontsize=5)\n",
      "\tplt.tick_params(axis='x',labelsize=5) \n",
      "\tplt.tick_params(axis='y',labelsize=5)  \n",
      "\tplt.title('28s Region',fontsize=10)\n",
      "\tplt.tight_layout()\n",
      "\n",
      "def plot_snoRNAdetail(outfilepath,sampleName):\n",
      "\t# Read through source files\n",
      "\thcBOX_READS = outfilepath+'clipGenes_snoRNA_LowFDRreads_HCAbox.bed'\n",
      "\tSCA_READS = outfilepath+'clipGenes_snoRNA_LowFDRreads_SCARNA.bed'\n",
      "\tCDBOX_READS = outfilepath+'clipGenes_snoRNA_LowFDRreads_CDBOX.bed'\n",
      "\treadList=[SCA_READS,hcBOX_READS,CDBOX_READS]\n",
      "\tcounts=[]\n",
      "\tfor readFile in readList:\n",
      "\t\tcounts.append(lineCount(readFile))\n",
      "\t\n",
      "\ttotalCount=float(counts[0]+counts[1]+counts[2])\n",
      "\n",
      "\tax=plt.subplot(2,2,1)\n",
      "\t# Check to make sure snoRNA reads recovered\n",
      "\tif totalCount > 0:\n",
      "\t\tinputVector=[float(counts[0])/totalCount,float(counts[1])/totalCount,float(counts[2])/totalCount]\n",
      "\t\tpie_wedges=ax.pie(inputVector,labels=[\"scaRNAs\",\"hBox\",\"CDbox\"],labeldistance=1.1,autopct='%1.1f%%')\n",
      "\t\tplt.rcParams['font.size']=7.5\n",
      "\t\tfor wedge in pie_wedges[0]:\n",
      "\t\t\twedge.set_edgecolor('black')\n",
      "\t\t\twedge.set_lw(1)\n",
      "\t\t\tplt.tight_layout()\n",
      "\n",
      "\t# Position files\n",
      "\tpositionFiles=glob.glob(outfilepath+'*_genePosition_snoRNA')\n",
      "\t# Generate coverage histogram for each snoRNA type \n",
      "\tsnoNames=[\"scaRNA\",\"C/D box\",\"H/ACA box\"]\n",
      "\tfor j in range(0,len(positionFiles)):\n",
      "\t\t# Skip the first column, which has name of gene, and order is gene start / end and cluster start / end\n",
      "\t\tdata=np.genfromtxt(positionFiles[j],usecols=(1,2,3,4,5),delimiter='\\t',dtype='float')\n",
      "\t\t# Add exception handling for case in which no snoRNA reads come through\n",
      "\t\ttry:\n",
      "\t\t\t# Ratio of read distance from start of gene to full gene length\n",
      "\t\t\ttemp1=np.subtract(np.transpose(data[:,0]),np.transpose(data[:,2]))\n",
      "\t\t\ttemp2=np.subtract(np.transpose(data[:,0]),np.transpose(data[:,1]))\n",
      "\t\t\tpercentage=np.absolute(np.divide(temp1,temp2))\n",
      "\t\t\t# Correct any erronous values that are greater than one\n",
      "\t\t\tpercentage[percentage > 1]=1\n",
      "\t\texcept:\n",
      "\t\t\tprint \"No snoRNA reads of this type\"\n",
      "\t\t\tpercentage=np.zeros(1)\n",
      "\t\t# Plot histograms\n",
      "\t\tplt.subplot(2,2,j+2)\n",
      "\t\tbins=np.arange(0,1,0.01)\n",
      "\t\thist,bins=np.histogram(percentage,bins=bins)\n",
      "\t\thist=np.array(hist/float(percentage.size),dtype=float)\n",
      "\t\twidth=0.7*(bins[1]-bins[0])\n",
      "\t\tcenter=(bins[:-1] + bins[1:])/2\n",
      "\t\tplt.bar(center, hist,align='center',width=width,color='blue',alpha=0.75)\n",
      "\t\tplt.tick_params(axis='x',labelsize=5) \n",
      "\t\tplt.tick_params(axis='y',labelsize=5)  \n",
      "\t\tplt.xlabel('Fraction of gene body',fontsize=5)\n",
      "\t\t# For spacing the figures, only add ylabel to two bottom subplots\n",
      "\t\tif j+2 != 2:\n",
      "\t\t\tplt.ylabel('Normalized RT stop / bin',fontsize=5)\n",
      "\t\tplt.title(snoNames[j] + ' RT stops: %s'%str(percentage.size),fontsize=7.5)\n",
      "\t\tj += 1\n",
      "\t\tplt.tight_layout()\n",
      "\n",
      "def plot_ncRNAs(outfilepath,sampleName):\n",
      "\t# Extract all files that list read start, stop as well as start,stop coordiantes for the gene.\n",
      "\tgenePositionFiles=glob.glob(outfilepath+\"*_genePosition\")\n",
      "\tplotDim=math.ceil(math.sqrt(len(genePositionFiles)))\n",
      "\tj=1\n",
      "\tfor i in range(0,len(genePositionFiles)):\n",
      "\t\t# Get the file name\n",
      "\t\thead, tail = os.path.split(genePositionFiles[i])\n",
      "\t\tgeneTypeName = tail.split(\"_\")[1]\n",
      "\t\t# Exclude certain gene names from the analysis; be sure to exclude snoRNAs since those are indepdently dealt with \n",
      "\t\tif geneTypeName != 'IG' and geneTypeName != 'LRG' and geneTypeName != 'proteinCoding' and geneTypeName != 'TR' and geneTypeName != 'snRNA' and geneTypeName != 'rRNA':\n",
      "\t\t\ttry:\t\n",
      "\t\t\t\t# Skip the first column, which has name of gene, and order is gene start / end and cluster start / end\n",
      "\t\t\t\tdata=np.loadtxt(genePositionFiles[i],dtype='float',usecols=range(1,5),delimiter='\\t')\n",
      "\t\t\t\tclusterCenterPoint=np.zeros([data.shape[0],1])\n",
      "\t\t\t\tclusterCenterPoint=np.mean(data[:,[2,3]],axis=1)\n",
      "\t\t\t\tclusterCenterPoint=np.around(clusterCenterPoint,0)\n",
      "\t\t\t\t# Subtract cluster center from the gene start coordinate, and divide by gene length\n",
      "\t\t\t\ttemp1=np.subtract(np.transpose(data[:,[0]]),clusterCenterPoint) # Must correct dimentions for row-wise subtraction \n",
      "\t\t\t\ttemp2=np.transpose(np.subtract(data[:,[0]],data[:,[1]]))\n",
      "\t\t\t\tpercentage=np.absolute(np.divide(temp1,temp2))\n",
      "\t\t\t\tpercentage[percentage > 1]=1 \n",
      "\t\t\texcept:\n",
      "\t\t\t\tprint \"No data to read.\"\n",
      "\t\t\t\tpercentage=np.zeros(1)\t\n",
      "\t\t\t# Save source data for each ncRNA\n",
      "\t\t\tallReadCounts=outfilepath+'SourceData_ncRNAs_%s' %geneTypeName\n",
      "\t\t\t# Save array as a column vector\n",
      "\t\t\tnp.savetxt(allReadCounts,percentage,delimiter='\\t',fmt=\"%s\")\n",
      "\t\t\t# Make subplot\n",
      "\t\t\tplt.subplot(3,4,j)\n",
      "\t\t\tbins=np.arange(0,1,0.01)\n",
      "\t\t\thist,bins=np.histogram(percentage,bins=bins)\n",
      "\t\t\t# Noramlize histogram by total number of RT stops for the given ncRNA\n",
      "\t\t\thist=np.array(hist/float(percentage.size),dtype=float)\n",
      "\t\t\twidth=0.7*(bins[1]-bins[0])\n",
      "\t\t\tcenter=(bins[:-1] + bins[1:])/2\n",
      "\t\t\tplt.xlabel('Fraction of ncRNA gene body',fontsize=5)\n",
      "\t\t\tplt.ylabel('Normalized RT stop / bin',fontsize=5)\n",
      "\t\t\tplt.bar(center,hist,align='center',width=width,color='blue',alpha=0.75)\n",
      "\t\t\tplt.tick_params(axis='x',labelbottom='off') \n",
      "\t\t\tplt.tick_params(axis='y',labelsize=5)  \n",
      "\t\t\tplt.title(geneTypeName+' RT stops: %s'%str(percentage.size),fontsize=5)\n",
      "\t\t\tplt.xlim(0,1)\n",
      "\t\t\tj += 1\n",
      "\t\t\tplt.tight_layout()\n",
      "\n",
      "def runPlots(outfilepath,sampleName):\n",
      "\n",
      "\tfig1=plt.figure(1)\n",
      "\tplot_ReadAccounting(outfilepath,sampleName)\n",
      "\tplot_ReadsPerCluster(outfilepath,sampleName)\n",
      "\tplot_ClusterSizes(outfilepath,sampleName)\n",
      "\tplot_clusterBindingIntensity(outfilepath,sampleName)\n",
      "\tplot_readsBymRNAregion(outfilepath,sampleName)\n",
      "\tplot_BoundGeneTypes(outfilepath,sampleName)\n",
      "\tfig1.tight_layout()\n",
      "\tfig1.savefig(outfilepath+'Figure1.png',format='png',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\tfig1.savefig(outfilepath+'Figure1.pdf',format='pdf',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\n",
      "\tfig2=plt.figure(2)\n",
      "\tplot_mRNAgeneBodyDist(outfilepath,sampleName)\n",
      "\tplot_geneBodyPartition(outfilepath,sampleName)\n",
      "\tfig2.tight_layout()\n",
      "\tfig2.savefig(outfilepath+'Figure2.png',format='png',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\tfig2.savefig(outfilepath+'Figure2.pdf',format='pdf',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\n",
      "\tfig3=plt.figure(3)\n",
      "\tplot_repeatRNA(outfilepath,sampleName)\n",
      "\tfig3.tight_layout()\n",
      "\tfig3.savefig(outfilepath+'Figure3.png',format='png',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\tfig3.savefig(outfilepath+'Figure3.pdf',format='pdf',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\n",
      "\tfig4=plt.figure(4)\n",
      "\tplot_rDNA(outfilepath,sampleName)\n",
      "\tfig4.tight_layout()\n",
      "\tfig4.savefig(outfilepath+'Figure4.png',format='png',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\tfig4.savefig(outfilepath+'Figure4.pdf',format='pdf',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\n",
      "\tfig5=plt.figure(5)\n",
      "\tplot_snoRNAdetail(outfilepath,sampleName)\n",
      "\tfig5.tight_layout()\n",
      "\tfig5.savefig(outfilepath+'Figure5.png',format='png',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\tfig5.savefig(outfilepath+'Figure5.pdf',format='pdf',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\n",
      "\tfig6=plt.figure(6)\n",
      "\tplot_ncRNAs(outfilepath,sampleName)\n",
      "\tfig6.tight_layout()\n",
      "\tfig6.savefig(outfilepath+'Figure6.png',format='png',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\tfig6.savefig(outfilepath+'Figure6.pdf',format='pdf',bbox_inches='tight',dpi=150,pad_inches=0.5)\n",
      "\n",
      "# ---- Plots ---\n",
      "print \"Run plots\"\n",
      "runPlots(outfilepath,sampleName)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Run plots\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R1.fastq\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "52697023\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R2.fastq\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "36234779\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R1_3ptrimmed.fastq\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "47992015\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R2_3ptrimmed.fastq\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "32654179\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R1_3ptrimmed_filter.fastq\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "45402363\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R2_3ptrimmed_filter.fastq\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "29651319\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R1_3ptrimmed_filter_nodupe.fastq\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "22604951\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R2_3ptrimmed_filter_nodupe.fastq\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "16514195\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R1_3ptrimmed_filter_nodupe_5ptrimmed_mappedToRepeatRNA.bed\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "2214451\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R2_3ptrimmed_filter_nodupe_5ptrimmed_mappedToRepeatRNA.bed\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1870185\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R1_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes.bed\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1141491\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes.bed\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "651851\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R1_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes.bed\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1141491\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes.bed\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "651851\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R1_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat.bed\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "465092\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat.bed\n",
        "Lines counted:\n",
        "235052\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_threshold=8_allreads.mergedRT_CLIPPERin.bed\n",
        "Lines counted:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "40627\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/Rec_threshold=8_allreads.mergedRT_CLIP_clusters_lowFDRreads.bed\n",
        "Lines counted:\n",
        "49443\n",
        "Error in line syntax."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "track name=\"/arrayAhome/lmartin/CLIP/results/Rec/Rec_threshold=8_allreads.mergedRT_CLIP_clusters\" visibility=2 colorByStrand=\"0,0,0 0,0,0\"\n",
        "\n",
        "Counting lines for ..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/clipGenes_proteinCoding\n",
        "Lines counted:\n",
        "229\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/clipGenes_pseudogenes\n",
        "Lines counted:\n",
        "56\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/clipGenes_snoRNA\n",
        "Lines counted:\n",
        "19\n",
        "Counting lines for ...\n",
        "/arrayAhome/lmartin/CLIP/results/Rec/clipGenes_lincRNA\n",
        "Lines counted:"
       ]
      }
     ],
     "prompt_number": "*"
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Funtions - "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def runSamtools(samfiles):\n",
      "    # Useage: Samfile processing.\n",
      "    # Input: Sam files from Bowtie mapping.\n",
      "    # Output: Duplicate removed, sorted bedFiles.\n",
      "    program = 'samtools'\n",
      "    program2 = 'bamToBed'\n",
      "    outBedFiles=[]\n",
      "    \n",
      "    for samfile in samfiles:\n",
      "            \n",
      "        # Convert to bamfile\n",
      "        bamfile = samfile.replace('.sam', '.bam')  \n",
      "        proc = subprocess.Popen( [program, 'view', '-bS', '-o', bamfile, samfile])\n",
      "        proc.communicate()\n",
      "        \n",
      "        # Sort the bamfile and note that samtools sort adds the .bam handle\n",
      "        bamfile_sort = bamfile.replace('.bam', '_sorted') \n",
      "        proc2 = subprocess.Popen( [program, 'sort', bamfile, bamfile_sort])\n",
      "        proc2.communicate()\n",
      "        \n",
      "        # Convert to bedFile\n",
      "        bedFile = bamfile_sort.replace('_sorted', '_withDupes.bed')\n",
      "        outfh = open(bedFile, 'w')\n",
      "        proc3 = subprocess.Popen( [program2,'-i', bamfile_sort+'.bam'],stdout=outfh)\n",
      "        proc3.communicate()\n",
      "        \n",
      "        outBedFiles=outBedFiles+[bedFile]\n",
      "        \n",
      "    return outBedFiles\n",
      "\n",
      "def runRepeatMask(mappedReads):\n",
      "    # Usage: Remove repeat regions from bedfile following mapping.\n",
      "    # Input: .bed file after mapping (duplicates removed by samtools) and blastlist regions removed.\n",
      "    # Output: Bedfile with repeat regions removed.\n",
      "    program = 'intersectBed'\n",
      "    # *** Chance repeat masker ***\n",
      "    # repeatregions = os.getcwd() + '/docs/repeat_masker.bed'\n",
      "    repeatregions = os.getcwd() + '/docs/mm9/mm9_repeatMasker.bed'\n",
      "    masked=[]\n",
      "    try:\n",
      "        for bedIn in mappedReads:\n",
      "            # File name for low FDR reads\n",
      "            noRepeat = bedIn.replace('.bed', '_noRepeat.bed')\n",
      "            # Intersect input reads with the repeat region, and return those that do not intersect \n",
      "            outfh = open(noRepeat, 'w')\n",
      "            proc = subprocess.Popen([program, '-a', bedIn, '-b', repeatregions, '-v'],stdout=outfh)\n",
      "            proc.communicate()\n",
      "            outfh.close()\n",
      "            masked=masked+[noRepeat]\n",
      "        return (masked)\n",
      "    except:\n",
      "        print \"Problem with repeat masking.\"\n",
      "        \n",
      "def runBlacklistRegions(mappedReads):\n",
      "\t# Usage: Remove blacklisted regions from bedfile following mapping.\n",
      "\t# Input: .bed file after mapping (duplicates removed by samtools).\n",
      "\t# Output: Bedfile with blacklisted regions removed.\n",
      "\tprogram = 'intersectBed'\n",
      "\tblacklistregions = os.getcwd() + '/docs/wgEncodeDukeMapabilityRegionsExcludable.bed'\n",
      "\tblackListed=[]\n",
      "\ttry:\n",
      "\t\tfor bedIn in mappedReads:\n",
      "\t\t\t# File name for low FDR reads\n",
      "\t\t\tnoBlacklist = bedIn.replace('.bed', '_noBlacklist.bed')\n",
      "\t\t\t# Intersect input reads with the blacklist region, and return those that do not intersect \n",
      "\t\t\toutfh = open(noBlacklist, 'w')\n",
      "\t\t\tproc = subprocess.Popen([program, '-a', bedIn, '-b', blacklistregions, '-v'],stdout=outfh)\n",
      "\t\t\tproc.communicate()\n",
      "\t\t\toutfh.close()\n",
      "\t\t\tblackListed=blackListed+[noBlacklist]\n",
      "\t\treturn (blackListed)\n",
      "\texcept:\n",
      "\t\tprint \"Problem with blacklist.\"\n",
      "\n",
      "def modifyName(filepath,newTag):\n",
      "\t# Useage: Modifies the filepath name. \n",
      "\t# Input: File path of format <path>/<name>.fastq and a string to add to the name.\n",
      "\t# Output: Returns the modified path of type <old path>_<new modifier>.fastq\n",
      "    head, tail = os.path.split(filepath)\n",
      "    oldname = tail.split('.')[0]\n",
      "    newName = head+\"/\"+oldname+\"_\"+newTag\n",
      "    return newName\n",
      "        \n",
      "def seperateStrands(mappedReads):\n",
      "\t# Useage: Seperate positive and negative strands.\n",
      "\t# Input: Paths to two bed files from Samtools.\n",
      "\t# Output: Paths to bed files isolated by strand.\n",
      "    # Create list for storing file names\n",
      "    negativeStrand=[]\n",
      "    positiveStrand=[]\n",
      "    # For each file in the input list\n",
      "    for mapFile in mappedReads:\n",
      "        # Open the file\n",
      "        with open(mapFile, 'r') as infile:\n",
      "            # Create new file handles\n",
      "            neg_strand=modifyName(mapFile,'neg.bed')\n",
      "            pos_strand=modifyName(mapFile,'pos.bed')\t\n",
      "            # Open new files\n",
      "            neg = open(neg_strand, 'w')\n",
      "            pos = open(pos_strand, 'w')\n",
      "            negativeStrand=negativeStrand+[neg_strand]\n",
      "            positiveStrand=positiveStrand+[pos_strand]\n",
      "            # Read one line at the time to memory, and write to outputs.\n",
      "            for line in infile:\t\n",
      "                # Sort read based upon strand, which is field six (index=5) \n",
      "                if str(line.strip().split('\\t')[5]) == '-':\n",
      "                    neg.write(line)\n",
      "                elif str(line.strip().split('\\t')[5]) == '+':\n",
      "                    pos.write(line)\n",
      "    return (negativeStrand,positiveStrand)\n",
      "        \n",
      "def modifyNegativeStrand(negativeStrandReads):\n",
      "    # Useage: For negative stranded reads, ensure 5' position (RT stop) is listed first.\n",
      "    # Input: Bed file paths to all negative stranded.\n",
      "    # Output: Paths to modified bed files.\n",
      "    negativeStrandEdit=[]\n",
      "    # For each file in the input list\n",
      "    for negativeRead in negativeStrandReads:\n",
      "        # Outpit file name\n",
      "        neg_strand_edited = modifyName(negativeRead,'edit.bed')\n",
      "        negativeStrandEdit=negativeStrandEdit+[neg_strand_edited]\n",
      "        # Open new files\n",
      "        neg_edit = open(neg_strand_edited, 'w')\n",
      "        with open(negativeRead, 'r') as infile:\n",
      "            for line in infile:\t\n",
      "                chrom,start,end,name,quality,strand=line.strip().split('\\t')\n",
      "                # For negative stranded reads, invert so that 5' position is listed first (3' position is an arbitrary 30 bases beyond)\n",
      "                neg_edit.write('\\t' .join((chrom,end,str(int(end)+30),name,quality,strand,'\\n')))\n",
      "    return negativeStrandEdit\n",
      "        \n",
      "def isolate5prime(strandedReads):\n",
      "\t# Useage: Isolate only the Chr, 5' position (RT stop), and strand.\n",
      "\t# Input: Bed file paths to strand seperated reads.\n",
      "\t# Output: Paths to 5' isolated reads.\n",
      "    # For each file in the input list\n",
      "    RTstops=[]\n",
      "    for reads in strandedReads:\n",
      "        # Outpit file name\n",
      "        RTstop = modifyName(reads,'RTstop.bed')\n",
      "        # Open new files\n",
      "        f = open(RTstop, 'w')\n",
      "        with open(reads, 'r') as infile:\n",
      "            RTstops=RTstops+[RTstop]\n",
      "            for line in infile:\t\n",
      "                chrom,start,end,name,quality,strand=line.strip().split('\\t')\n",
      "                \n",
      "                f.write('\\t' .join((chrom,start,strand,'\\n')))\n",
      "\n",
      "    return RTstops\n",
      "\n",
      "def fileCat(destinationFile,fileList):\n",
      "\t# Useage: Concatenate two files.\n",
      "\t# Input: Output file path, as well as a list input files.\n",
      "\t# Output: Nothing (outfile name is specified in input).\n",
      "    f = open(destinationFile, \"w\")\n",
      "    for tempfile in fileList:\n",
      "        # Read each file into the destrination file\n",
      "        readfile = open(tempfile, \"r\")\n",
      "        f.write(readfile.read())\n",
      "        readfile.close()\n",
      "    f.close()\n",
      "\n",
      "def mergeRT(RTstops,outfilename,threshold):\n",
      "\t# Useage: Merge RT stops between replicates.\n",
      "\t# Input: Paths to RT stop files (stranded reads) and output filename.\n",
      "\t# Output: Nothing (outfile name is specified in input)\n",
      "\n",
      "    # Create object for storing dictionaries\n",
      "    store=[0,0]\n",
      "    # Flag for storing each dictionary\n",
      "    i=0\n",
      "    # Iterare through each file in input\n",
      "    for RT in RTstops:\n",
      "        # Create a dictionary with each value initialized to zero\n",
      "        d = defaultdict(int)\n",
      "        # Open the file \n",
      "        with open(RT, 'r') as infile:\t\t\n",
      "            for line in infile:\t\n",
      "                if line in d:\n",
      "                    d[line] += 1\n",
      "                else:\n",
      "                    # Make each RT stop a unique key in the dictionary\n",
      "                    d[line] += 1\n",
      "        # Store the dictionary\n",
      "        store[i]=d\n",
      "        i += 1\n",
      "        \n",
      "    rt_rep1=[k for k in store[0].iteritems()]\n",
      "    rt_rep2=[k for k in store[1].iteritems()]\n",
      "\n",
      "    # Open output files\n",
      "    f = open(outfilename, 'w')\n",
      "\n",
      "    # Iderate through each RT stop in dictionary 1\t\n",
      "    for key in store[0]:\n",
      "        # Check if same RT stop is in dictionary 2\n",
      "        name=store[1].get(key,None)\n",
      "\n",
      "        # If so, then the read is preserved in both files\n",
      "        if name:\n",
      "            \n",
      "            # Get total number of RT stops\n",
      "            totalRT = store[0][key] + store[1][key]\n",
      "\n",
      "            # Record if the count if it exceeds the threshold\n",
      "            if totalRT > threshold:\n",
      "                \n",
      "                chrom,start,strand=key.strip().split('\\t')\n",
      "                # Make sure the start of the read is greater than 15 bases from end of the chrom\n",
      "                if int(start)>15:\n",
      "                    # Create a new read centered +/- 15 bases around the RT stop coordinate\n",
      "                    read='\\t' .join((chrom,str(int(start)-15),str(int(start)+15),'CLIPread','255',strand,'\\n'))\n",
      "                else:\n",
      "                    read='\\t' .join((chrom,str(int(start)),str(int(start)+15),'CLIPread','255',strand,'\\n'))\n",
      "                # Write the read to the output, and repeat for the total number of instances the RT stop appears\n",
      "                f.write(read*(store[0][key]+store[1][key]))\n",
      "                \n",
      "def runCLIPPER(RTclusterfile):\n",
      "    # Useage: Process the mergedRT file and pass through CLIPper FDR script.\n",
      "    # Input: Merged RT file.\n",
      "    # Output: CLIPper input (.bed) file and output file.\n",
      "    program = 'bedToBam'\n",
      "    genomeFile = os.getcwd()+'/docs/human.hg19.genome'\n",
      "    program2 = 'samtools'\n",
      "    program3 = 'bamToBed'\n",
      "    program4 = 'clipper'\n",
      "    print \"Running CLIPper...\"\n",
      "    \n",
      "    # Create and open bamfile\n",
      "    bamfile = RTclusterfile.replace('.bed', '.bam')  \n",
      "    outfh = open(bamfile, 'w')\n",
      "    proc = subprocess.Popen([program, '-i', RTclusterfile,'-g',genomeFile],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    \n",
      "    bamfile_sort = bamfile.replace('.bam', '.srt')\n",
      "    proc2 = subprocess.Popen([program2, 'sort', bamfile, bamfile_sort])\n",
      "    proc2.communicate()\n",
      "    \n",
      "    bamfile_sorted=bamfile_sort+'.bam'\n",
      "    mapStats = bamfile_sorted.replace('.srt.bam', '.mapStats.txt') \n",
      "    outfh = open(mapStats, 'w')\n",
      "    proc3 = subprocess.Popen([program2, 'flagstat', bamfile_sorted],stdout=outfh)\n",
      "    proc3.communicate()\n",
      "    \n",
      "    proc4 = subprocess.Popen([program2, 'index', bamfile_sorted])\n",
      "    proc4.communicate()\n",
      "    \n",
      "    CLIPPERin = bamfile_sorted.replace('.srt.bam', '_CLIPPERin.bed') \n",
      "    outfh = open(CLIPPERin, 'w')\n",
      "    proc5 = subprocess.Popen([program3, '-i', bamfile_sorted],stdout=outfh)\n",
      "    proc5.communicate()\n",
      "    \n",
      "    CLIPPERout = CLIPPERin.replace('_CLIPPERin.bed', '_CLIP_clusters') \n",
      "    # *** Change CLIPper call ***\n",
      "    # proc6 = subprocess.Popen([program4, '--bam', bamfile_sorted,'-shg19','--outfile=%s'%CLIPPERout],)\n",
      "    proc6 = subprocess.Popen([program4, '--bam', bamfile_sorted,'-smm9','--outfile=%s'%CLIPPERout],)\n",
      "    proc6.communicate()\n",
      "    outfh.close()\n",
      "    \n",
      "    return (CLIPPERin,CLIPPERout)\n",
      "\n",
      "def modCLIPPERout(CLIPPERin,CLIPPERout):\n",
      "    # Usage: Process the CLIPper output and isolate lowFDR reads based upon CLIPper windows.\n",
      "    # Input: .bed file passed into CLIPper and the CLIPper windows file.\n",
      "    # Output: Low FDR reads recovered using the CLIPer windows file, genes per cluster, gene list of CLIPper clusters, and CLIPper windows as .bed.\n",
      "    program = 'intersectBed'\n",
      "    # Output format from CLIPper will be: <same_name>_CLIP_clusters\n",
      "    CLIPperOutBed=CLIPPERout+'.bed'\n",
      "    CLIPpeReadsPerCluster=CLIPPERout+'.readsPerCluster'\n",
      "    CLIPpeGeneList=CLIPPERout+'.geneNames'\n",
      "    \n",
      "    # Open new files\n",
      "    f = open(CLIPperOutBed, 'w')\n",
      "    g = open(CLIPpeReadsPerCluster, 'w')\n",
      "    h = open(CLIPpeGeneList, 'w')\n",
      "    with open(CLIPPERout, 'r') as infile:\n",
      "        # For each CLIPper window.\n",
      "        for line in infile:\t\n",
      "            try:\n",
      "                # Version of CLIPper used here includes a header that cannot be parsed. Handle this.\n",
      "                chrom,start,end,name,stats,strand,start_2,end_2 = line.strip().split('\\t')\n",
      "                # Old CLIPPER: Ensembl genes are parsed with <name>_<cluster>_<count>\n",
      "                readPerCluster=name.strip().split('_')[2]\n",
      "                geneName=name.strip().split('_')[0]\n",
      "                # Re-write the CLIPper windows file\n",
      "                f.write('\\t' .join((chrom,start,end,name,stats,strand,'\\n')))\n",
      "                g.write((readPerCluster+'\\n'))\n",
      "                h.write((geneName+'\\n'))\n",
      "            except:\n",
      "                print \"\"\n",
      "\n",
      "    f.close()\n",
      "    g.close()\n",
      "    h.close()\n",
      "    \n",
      "    # File name for low FDR reads\n",
      "    CLIPPERlowFDR = CLIPperOutBed.replace('.bed', '_lowFDRreads.bed')\n",
      "    # Intersect input reads with the CLIPper windows\n",
      "    outfh = open(CLIPPERlowFDR, 'w')\n",
      "    # Better docs: http://bedtools.readthedocs.org/en/latest/content/tools/intersect.html\n",
      "    # -wa should report full alignment from a overlap.\n",
      "    proc = subprocess.Popen([program,'-a', CLIPPERin, '-b', CLIPperOutBed,'-wa','-wb','-s'],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    outfh.close()\n",
      "    \n",
      "    return (CLIPPERlowFDR,CLIPpeReadsPerCluster,CLIPpeGeneList,CLIPperOutBed)\n",
      "\n",
      "def cleanBedFile(inBed):\n",
      "    # Usage: Sort and recover only first 6 fields from a bed file.\n",
      "    # Input: BedFile.\n",
      "    # Output: Sorted bedFile with correct number of fields.\n",
      "    program='sortBed'\n",
      "    \n",
      "    # Make sure bedfile only has 5 fields\n",
      "    CLIPperOutBed=inBed.replace('.bed','_cleaned.bed')\t\n",
      "    sortedBed=CLIPperOutBed.replace('_cleaned.bed','_cleaned_sorted.bed')\n",
      "    \n",
      "    # Open new files\n",
      "    f = open(CLIPperOutBed, 'w')\n",
      "    with open(inBed, 'r') as infile:\n",
      "        for line in infile:\t\n",
      "            elementList = line.strip().split('\\t')\n",
      "            # Re-write the CLIPper windows file\n",
      "            f.write('\\t' .join((elementList[0],elementList[1],elementList[2],elementList[3],elementList[4],elementList[5],'\\n')))\n",
      "    f.close()\n",
      "    \n",
      "    # Sort the resulting bedFile\n",
      "    outfh = open(sortedBed, 'w')\n",
      "    proc = subprocess.Popen([program, '-i', CLIPperOutBed],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    outfh.close()\n",
      "    \n",
      "    return sortedBed\n",
      "\n",
      "def makeBedGraph(lowFDRreads):\n",
      "    # Usage: From a bedFile, generate a bedGraph and bigWig.\n",
      "    # Input: BedFile.\n",
      "    # Output: BedGraph file.\n",
      "    program = 'genomeCoverageBed'\n",
      "    program2 = os.getcwd() + '/bin/bedGraphToBigWig'\n",
      "    # Note include this reference to account for the fact that non-traditional chroms (e.g., chrUn_gl000220) can appear in the data.\n",
      "    # Without this, bedGraph creation will throw an error.\n",
      "    \n",
      "    # *** Change sizes reference ***\n",
      "    #sizesFile = os.getcwd()+'/docs/human.hg19.genome'\n",
      "    #sizesFile2 = os.getcwd()+'/docs/hg19.sizes'\n",
      "    sizesFile = os.getcwd()+'/docs/mm9/mm9.sizes' # Note that a larger sizes file is needed for human at this tep : '/docs/human.hg19.genome'\n",
      "    sizesFile2 = os.getcwd()+'/docs/mm9/mm9.sizes'\n",
      "    \n",
      "    cleanBed = cleanBedFile(lowFDRreads)\n",
      "    outname = cleanBed.replace('.bed','.bedgraph')\n",
      "    outname2 = cleanBed.replace('.bed','.bw')\n",
      "    \n",
      "    outfh = open(outname, 'w')\n",
      "    proc = subprocess.Popen([program, '-bg','-split','-i',cleanBed,'-g',sizesFile],stdout=outfh)\n",
      "    proc.communicate()\n",
      "    \n",
      "    outfh2 = open(outname2, 'w')\n",
      "    proc2 = subprocess.Popen([program2,outname,sizesFile,outname2],stdout=subprocess.PIPE)\n",
      "    proc2.communicate()\n",
      "    return outname\n",
      "\n",
      "def getBedCenterPoints(inBed):\n",
      "    # Usage: Obtain ceter coordiantes of bedFile\n",
      "    # Input: BedFile.\n",
      "    # Output: Center coodinates returned.\n",
      "    \n",
      "    # Make sure bedfile only has 5 fields\n",
      "    outBed=inBed.replace('.bed','_centerCoord.bed')\t\n",
      "    # Open new files\n",
      "    f = open(outBed, 'w')\n",
      "    with open(inBed, 'r') as infile:\n",
      "        for line in infile:\t\n",
      "            elementList = line.strip().split('\\t')\n",
      "            # Re-write the CLIPper windows file\n",
      "            f.write('\\t' .join((elementList[0],str(int(elementList[1])+15),str(int(elementList[1])+16),elementList[9],elementList[4],elementList[5],'\\n')))\n",
      "    f.close()\n",
      "    return outBed"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Scratch space - "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "### Check samtools output ###\n",
      "# - Sort files - \n",
      "# cat ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_nodupes.bed | sort > temp_mapped_R1_nodupes.bed\n",
      "# cat ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes.bed | sort > temp_mapped_R1_withdupes.bed\n",
      "# - Result -\n",
      "# 10182670 temp_mapped_R1_withdupes.bed\n",
      "# 6912195 temp_mapped_R1_nodupes.bed\n",
      "# => Therefore, ~ 2x effect with duplicate removal by samtools. We may want to remove this.\n",
      "\n",
      "\n",
      "# - Samtools dupes versus no dupes - \n",
      "# 5736213 ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat.bed\n",
      "# 3967167 ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_nodupes_noBlacklist_noRepeat.bed\n",
      "\n",
      "# - Uniq lines without Samtools rmdupes (therefore ~ 1M RT stop are duplicated) - \n",
      "# cat ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat.bed | awk -F '\\t' '{print $1, \"\\t\", $2, \"\\t\", $6}' | sort | uniq | wc -l \n",
      "# = 3892401\n",
      "# - Uniq lines with Samtools rmdupes (therefore nearly all RT stop are uniq already) - \n",
      "# cat ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_nodupes_noBlacklist_noRepeat.bed | awk -F '\\t' '{print $1, \"\\t\", $2, \"\\t\", $6}' | sort | uniq | wc -l\n",
      "# = 3811599\n",
      "\n",
      "\n",
      "# Examine sorted RT, pos RT stop - \n",
      "\n",
      "# cat ddx21v3_R1_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat_pos_RTstop.bed | sort > temp_R1_posRT\n",
      "# cat ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat_pos_RTstop.bed | sort > temp_R2_posRT\n",
      "# cat ddx21v3_positivereads.mergedRT | sort > temp_posMerge\n",
      "\n",
      "\n",
      "# - Check neg edit - \n",
      "# head ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat_neg.bed\n",
      "# chr1\t12132\t12187\t11449178-2\t0\t-\n",
      "\n",
      "# head ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat_neg_edit.bed\n",
      "# chr1\t12187\t12217\t11449178-2\t0\t-\n",
      "\n",
      "# Therefore, it works: adds 30 the flipped 5' coordinate.\n",
      "\n",
      "\n",
      "# - Check total RT stops - \n",
      "# 5736213 ddx21v3_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedToRepeat_mapped_withDupes_noBlacklist_noRepeat.bed\n",
      "2945362+2790851\n",
      "\n",
      "\n",
      "'''\n",
      "# - First RT stop +/- 15 from the RT stop and total number is repeats - \n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "chr10\t101996648\t101996678\tCLIPread\t255\t+\t\n",
      "'''\n",
      "\n",
      "'''\n",
      "# - temp_R1_posRT (remember to add 15 : 101996648+15 = 101996663 - \n",
      "chr10\t101996663\t+\t\n",
      "chr10\t101996663\t+\t\n",
      "chr10\t101996663\t+\t\n",
      "chr10\t101996663\t+\t\n",
      "chr10\t101996663\t+\t\n",
      "chr10\t101996663\t+\n",
      "'''\n",
      "\n",
      "'''\n",
      "# - temp_R2_posRT (remember to add 15) - \n",
      "chr10\t101996663\t+\t\n",
      "chr10\t101996663\t+\t\n",
      "chr10\t101996663\t+\t\n",
      "chr10\t101996663\t+\t\n",
      "chr10\t101996663\t+\t\n",
      "'''\n",
      "\n",
      "### Examine data from old useage of intersectBed ###\n",
      "\n",
      "# 1) For ddx3wtAT2 data, first cluster is:\n",
      "# chr10\t100004483\t100004505\tENSG00000166024_9_4\t0.00365984682734\t+\t100004491\t100004495\n",
      "\n",
      "# 2) CLIPper in data that overlaps\n",
      "# chr10\t100004473\t100004503\tCLIPread\t255\t+\n",
      "# chr10\t100004473\t100004503\tCLIPread\t255\t+\n",
      "# chr10\t100004485\t100004515\tCLIPread\t255\t+\n",
      "# chr10\t100004485\t100004515\tCLIPread\t255\t+\n",
      "\n",
      "# 3) Reads recorded\n",
      "# chr10\t100004483\t100004503\tCLIPread\t255\t+\tchr10\t100004483\t100004505\tENSG00000166024_9_4\t0.00365984682734\t+\n",
      "# chr10\t100004483\t100004503\tCLIPread\t255\t+\tchr10\t100004483\t100004505\tENSG00000166024_9_4\t0.00365984682734\t+\n",
      "# chr10\t100004485\t100004505\tCLIPread\t255\t+\tchr10\t100004483\t100004505\tENSG00000166024_9_4\t0.00365984682734\t+\n",
      "# chr10\t100004485\t100004505\tCLIPread\t255\t+\tchr10\t100004483\t100004505\tENSG00000166024_9_4\t0.00365984682734\t+\n",
      "\n",
      "# Therefore, it clearly stops at the CLIP cluster windows. This is a problem on the 5' end, as we loose RT stop information.\n",
      "# Based upon intersectBed docs, -wa should address this."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "### Examine data from new useage of intersectBed ###\n",
      "\n",
      "# First cluster in sorted clusters file: ddx21_clipper_clusters\n",
      "# chr10\t102114184\t102114208\tENSG00000099194_1_15\t1.86517468137e-14\t+\t102114192\t102114196\n",
      "\n",
      "# INFILE: ddx21_clipper_in\n",
      "'''\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "'''\n",
      "\n",
      "# OUTFILE: ddx21_clipper_outreads\n",
      "'''\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "chr10\t102114177\t102114207\tCLIPread\t255\t+\n",
      "'''\n",
      "\n",
      "# Therefore, modifying -wa fixes the thresholding issue."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}