ViralGenome_analysis.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sys\n",
      "import os\n",
      "import numpy as np\n",
      "import matplotlib\n",
      "import matplotlib.pyplot as plt \n",
      "import pandas as pd\n",
      "import glob \n",
      "import math\n",
      "import re\n",
      "from collections import defaultdict\n",
      "from matplotlib.font_manager import FontProperties\n",
      "import subprocess\n",
      "import scipy\n",
      "from scipy.stats.stats import pearsonr \n",
      "%matplotlib inline\n",
      "matplotlib.rcParams['savefig.dpi'] = 3 * matplotlib.rcParams['savefig.dpi']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get raw data\n",
      "global outfilepath\n",
      "sampleName='ddx3wt'\n",
      "outfilepath=os.getcwd() + '/results/%s/'%sampleName"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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 runBowtie(fastqFiles,viralGenome):\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",
      "        index = os.getcwd() + '/docs/%s/%s'%(viralGenome,viralGenome)\n",
      "        outfile = modifyName(infastq,\"mappedTo%s.sam\"%viralGenome)\n",
      "        unmapped = modifyName(infastq,\"notMappedTo%s.fastq\"%viralGenome)\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,'-a','-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",
      "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",
      "# Run bowtie\n",
      "print \"Mapping to viral index.\"  \n",
      "trimmedReads5p=glob.glob(outfilepath+'*5ptrimmed.fastq')\n",
      "mappedReads,unMappedReads=runBowtie(trimmedReads5p,'jfh1')\n",
      "\n",
      "# Run Samtools\n",
      "print \"Process mapped data\"  \n",
      "mappedBedFiles=runSamtools(mappedReads)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Mapping to viral index.\n",
        "Performing Bowtie...\n",
        "Input file:\n",
        "/arrayAhome/lmartin/CLIP/results/hcvwt/hcvwt_R2_3ptrimmed_filter_nodupe_5ptrimmed.fastq\n",
        "Genome index:\n",
        "/arrayAhome/lmartin/CLIP/docs/jfh1/jfh1\n",
        "Output file (mapped):\n",
        "/arrayAhome/lmartin/CLIP/results/hcvwt/hcvwt_R2_3ptrimmed_filter_nodupe_5ptrimmed_mappedTojfh1.sam\n",
        "Output file (unmapped)\n",
        "/arrayAhome/lmartin/CLIP/results/hcvwt/hcvwt_R2_3ptrimmed_filter_nodupe_5ptrimmed_notMappedTojfh1.fastq\n",
        "Result : "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \n",
        "Error :  Use of uninitialized value $l in scalar chomp at /seq/bowtie2-2.1.0/bowtie2 line 437, <BT> line 4858856.\n",
        "Use of uninitialized value $l in substitution (s///) at /seq/bowtie2-2.1.0/bowtie2 line 438, <BT> line 4858856.\n",
        "Use of uninitialized value $l in print at /seq/bowtie2-2.1.0/bowtie2 line 442, <BT> line 4858856.\n",
        "\n",
        "Input file:\n",
        "/arrayAhome/lmartin/CLIP/results/hcvwt/hcvwt_R1_3ptrimmed_filter_nodupe_5ptrimmed.fastq\n",
        "Genome index:\n",
        "/arrayAhome/lmartin/CLIP/docs/jfh1/jfh1\n",
        "Output file (mapped):\n",
        "/arrayAhome/lmartin/CLIP/results/hcvwt/hcvwt_R1_3ptrimmed_filter_nodupe_5ptrimmed_mappedTojfh1.sam\n",
        "Output file (unmapped)\n",
        "/arrayAhome/lmartin/CLIP/results/hcvwt/hcvwt_R1_3ptrimmed_filter_nodupe_5ptrimmed_notMappedTojfh1.fastq\n",
        "Result : "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \n",
        "Error :  \n",
        "Process mapped data\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Index\n",
      "outfilepathjfh1=os.getcwd() + '/docs/jfh1/jfh1.fa'\n",
      "jfh1_genome=np.genfromtxt(outfilepathjfh1,dtype='string')\n",
      "jfh1_genome_size=len(jfh1_genome[1])\n",
      "\n",
      "# Bed files for mapped reads\n",
      "sampleName='jfh1'\n",
      "outfilepathjfh1=os.getcwd() + '/results/'+sampleName+'/'\n",
      "mappedBedFiles1=glob.glob(outfilepathjfh1+'*mappedToJFH1.bed')\n",
      "sampleName='jfh1v2'\n",
      "outfilepathjfh1=os.getcwd() + '/results/'+sampleName+'/'\n",
      "mappedBedFiles2=glob.glob(outfilepathjfh1+'*mappedToJFH1.bed')\n",
      "sampleName='hcvwt'\n",
      "outfilepathjfh1=os.getcwd() + '/results/'+sampleName+'/'\n",
      "mappedBedFiles3=glob.glob(outfilepathjfh1+'*mappedToJFH1.bed')\n",
      "mappedBedFiles=mappedBedFiles1+mappedBedFiles2+mappedBedFiles3\n",
      "\n",
      "# Global plotting stats\n",
      "bins=np.array(range(0,jfh1_genome_size,5)) # Bin the genome by 5 bp\n",
      "width=0.7*(bins[1]-bins[0]) # Bar width for plotting\n",
      "center=(bins[:-1] + bins[1:])/2 # Center position of each bin\n",
      "\n",
      "# Compute histogram for each\n",
      "coverageHistograms=[]\n",
      "for bedFile in mappedBedFiles:\n",
      "    data=np.genfromtxt(bedFile,usecols=(1,),delimiter='\\t',dtype='int') # 5' coordinate (RT stop)\n",
      "    hist,bins=np.histogram(data,bins=bins) # Compute histogram\n",
      "    hist=hist/float(len(data)) # Normalize to total number of reads mapped to viral genome\n",
      "    coverageHistograms.append(hist)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Histograms to filter\n",
      "HistToFilter=[0,1,2,3]\n",
      "Background=[4,5]\n",
      "filteredHist=[]\n",
      "# For each histogram\n",
      "for i in HistToFilter:\n",
      "    # Copy the values into temp object\n",
      "    tempHist=coverageHistograms[i]\n",
      "    # For each index in the histogram\n",
      "    for index in range(0,len(tempHist)):\n",
      "        # Check if background histograms (no infection) have counts\n",
      "        if coverageHistograms[Background[0]][index] > 0 or coverageHistograms[Background[1]][index] > 0:\n",
      "            # If so, then do not count\n",
      "            tempHist[index]=0\n",
      "    filteredHist.append(tempHist)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "IndexError",
       "evalue": "list index out of range",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mIndexError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-12-207c43447f61>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mHistToFilter\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m     \u001b[1;31m# Copy the values into temp object\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 8\u001b[1;33m     \u001b[0mtempHist\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcoverageHistograms\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      9\u001b[0m     \u001b[1;31m# For each index in the histogram\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     10\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0mindex\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtempHist\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mIndexError\u001b[0m: list index out of range"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Scatter plot\n",
      "geneType='HCV genome'\n",
      "sampleName1='JFH-1 Rep1'\n",
      "sampleName2='JFH-1 Rep2'\n",
      " \n",
      "# Compute based upon histograms bins\n",
      "bins=np.array(range(0,jfh1_genome_size,5)) # Bin the genome by 5 bp\n",
      "width=0.7*(bins[1]-bins[0]) # Bar width for plotting\n",
      "center=(bins[:-1] + bins[1:])/2 # Center position of each bin\n",
      "\n",
      "# Annotation\n",
      "UTR5end=340\n",
      "UTR3start=9280\n",
      "\n",
      "# Isolate UTR data from new JFH-1 histograms\n",
      "boolarray=[element < UTR5end or element > UTR3start for element in center]\n",
      "boolarray=np.array(boolarray,dtype=bool)\n",
      "center=center[boolarray]\n",
      "sample1=filteredHist[2][boolarray]\n",
      "sample2=filteredHist[3][boolarray]\n",
      "\n",
      "# Fit_fn is now a function which takes in x and returns an estimate for y\n",
      "fit = np.polyfit(sample1,sample2,1)\n",
      "fit_fn = np.poly1d(fit)\n",
      "rS=scipy.stats.linregress(sample1,sample2)\n",
      "print rS\n",
      "\n",
      "fig3=plt.figure(3)\n",
      "\n",
      "# Iterate through cdenter index of each bin\n",
      "plt.subplot(1,2,1)\n",
      "for i in range(0,len(center)):\n",
      "    # Get position at the point\n",
      "    position=center[i]\n",
      "    # Annotate\n",
      "    if position < UTR5end:\n",
      "        color='blue'\n",
      "        plt.scatter(sample1[i],sample2[i],color=color,alpha=0.35) \n",
      "    elif position > UTR3start:\n",
      "        color='red'\n",
      "        plt.scatter(sample1[i],sample2[i],color=color,alpha=0.35)\n",
      "plt.plot(sample1,fit_fn(sample1),alpha=0.5,color='black')\n",
      "plt.xlabel(\"Normalized gene read count: %s\"%sampleName1,fontsize=5)\n",
      "plt.ylabel(\"Normalized gene read count: %s\"%sampleName2,fontsize=5)\n",
      "plt.tick_params(axis='xticks',labelsize=5) \n",
      "plt.tick_params(axis='yticks',labelsize=5)\n",
      "plt.xticks(rotation=90)\n",
      "maxVal=0.04\n",
      "plt.xlim(-0.1*maxVal,maxVal)\n",
      "plt.ylim(-0.1*maxVal,maxVal)\n",
      "\n",
      "plt.subplot(1,2,2)\n",
      "for i in range(0,len(center)):\n",
      "    # Get position at the point\n",
      "    position=center[i]\n",
      "    # Annotate\n",
      "    if position < UTR5end:\n",
      "        color='blue'\n",
      "        plt.scatter(sample1[i],sample2[i],color=color,alpha=0.35) \n",
      "    elif position > UTR3start:\n",
      "        color='red'\n",
      "        plt.scatter(sample1[i],sample2[i],color=color,alpha=0.35)\n",
      "plt.plot(sample1,fit_fn(sample1),alpha=0.5,color='black')\n",
      "plt.xlabel(\"Normalized gene read count: %s\"%sampleName1,fontsize=5)\n",
      "plt.ylabel(\"Normalized gene read count: %s\"%sampleName2,fontsize=5)\n",
      "plt.tick_params(axis='xticks',labelsize=5) \n",
      "plt.tick_params(axis='yticks',labelsize=5)\n",
      "plt.xticks(rotation=90)\n",
      "maxVal=0.01\n",
      "plt.xlim(-0.1*maxVal,maxVal)\n",
      "plt.ylim(-0.1*maxVal,maxVal)\n",
      "plt.tight_layout()\t\n",
      "\n",
      "project='pcbp'\n",
      "outfilepath=os.getcwd() + '/metaanalysis/%s/'%project\n",
      "fig3.savefig(outfilepath+'Fig5a.png',format='png',bbox_inches='tight',dpi=200,pad_inches=0.5)\n",
      "fig3.savefig(outfilepath+'FIg5a.pdf',format='pdf',bbox_inches='tight',dpi=200,pad_inches=0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig4=plt.figure(4)\n",
      "\n",
      "bins=np.array(range(0,jfh1_genome_size,5)) # Bin the genome by 5 bp\n",
      "width=0.7*(bins[1]-bins[0]) # Bar width for plotting\n",
      "center=(bins[:-1] + bins[1:])/2 # Center position of each bin\n",
      "\n",
      "# 5' UTR specific\n",
      "plt.subplot(1,2,1)\n",
      "color='blue'\n",
      "for i in range(0,4):    \n",
      "    plt.bar(center,filteredHist[i],align='center',width=width,color=color,alpha=0.3)\n",
      "    plt.tick_params(axis='x',labelsize=5)\n",
      "    plt.tick_params(axis='y',labelsize=5)\n",
      "plt.title(\"Coverage histogram for JFH-1 5' UTR.\",fontsize=7.5)\n",
      "plt.xlim([-10,400])\n",
      "plt.ylim([0,0.01])\n",
      "\n",
      "# 3' UTR specific\n",
      "plt.subplot(1,2,2)\n",
      "color='red'\n",
      "for i in range(0,4):    \n",
      "    plt.bar(center,filteredHist[i],align='center',width=width,color=color,alpha=0.3)\n",
      "    plt.tick_params(axis='x',labelsize=5)\n",
      "    plt.tick_params(axis='y',labelsize=5)\n",
      "plt.title(\"Coverage histogram for JFH-1 3' UTR.\",fontsize=7.5)\n",
      "plt.xlim([9280,jfh1_genome_size+10])\n",
      "plt.ylim([0,0.05])\n",
      "\n",
      "project='pcbp'\n",
      "outfilepath=os.getcwd() + '/metaanalysis/%s/'%project\n",
      "fig4.savefig(outfilepath+'Fig5b.png',format='png',bbox_inches='tight',dpi=200,pad_inches=0.5)\n",
      "fig4.savefig(outfilepath+'Fig5b.pdf',format='pdf',bbox_inches='tight',dpi=200,pad_inches=0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig5=plt.figure(5)\n",
      "\n",
      "# Compute based upon histograms bins\n",
      "bins=np.array(range(0,jfh1_genome_size,5)) # Bin the genome by 5 bp\n",
      "width=0.7*(bins[1]-bins[0]) # Bar width for plotting\n",
      "center=(bins[:-1] + bins[1:])/2 # Center position of each bin\n",
      "\n",
      "# Scatter plot\n",
      "geneType='HCV genome'\n",
      "sampleName1='JFH-1 Rep1'\n",
      "sampleName2='JFH-1 Rep2'\n",
      "\n",
      "# Annotation\n",
      "UTR5end=340\n",
      "UTR3start=9280\n",
      " \n",
      "# Fit_fn is now a function which takes in x and returns an estimate for y\n",
      "fit = np.polyfit(filteredHist[2],filteredHist[3],1)\n",
      "fit_fn = np.poly1d(fit)\n",
      "rS=scipy.stats.linregress(filteredHist[2],filteredHist[3])\n",
      "print rS\n",
      "\n",
      "# Iterate through cdenter index of each bin\n",
      "plt.subplot(1,2,1)\n",
      "for i in range(0,len(center)):\n",
      "    # Get position at the point\n",
      "    position=center[i]\n",
      "    # Annotate\n",
      "    if position < UTR5end:\n",
      "        color='blue'\n",
      "    elif position > UTR3start:\n",
      "        color='red'\n",
      "    else:\n",
      "        color='cyan'\n",
      "    # Extract points from newer JFH-1 data histograms\n",
      "    plt.scatter(filteredHist[2][i],filteredHist[3][i],color=color,alpha=0.35)\n",
      "plt.plot(filteredHist[2],fit_fn(filteredHist[2]),alpha=0.5,color='black')\n",
      "plt.xlabel(\"Normalized gene read count: %s\"%sampleName1,fontsize=5)\n",
      "plt.ylabel(\"Normalized gene read count: %s\"%sampleName2,fontsize=5)\n",
      "plt.title(\"Comparison for %s\"%geneType)\n",
      "plt.tick_params(axis='xticks',labelsize=5) \n",
      "plt.tick_params(axis='yticks',labelsize=5)\n",
      "plt.xticks(rotation=90)\n",
      "maxVal=0.04\n",
      "plt.xlim(-0.1*maxVal,maxVal)\n",
      "plt.ylim(-0.1*maxVal,maxVal)\n",
      "\n",
      "plt.subplot(1,2,2)\n",
      "for i in range(0,len(center)):\n",
      "    # Get position at the point\n",
      "    position=center[i]\n",
      "    # Annotate\n",
      "    if position < UTR5end:\n",
      "        color='blue'\n",
      "    elif position > UTR3start:\n",
      "        color='red'\n",
      "    else:\n",
      "        color='cyan'\n",
      "    # Extract points from newer JFH-1 data histograms\n",
      "    plt.scatter(filteredHist[2][i],filteredHist[3][i],color=color,alpha=0.35)\n",
      "plt.plot(filteredHist[2],fit_fn(filteredHist[2]),alpha=0.5,color='black')\n",
      "plt.title(\"Zoom\")\n",
      "plt.tick_params(axis='xticks',labelsize=5) \n",
      "plt.tick_params(axis='yticks',labelsize=5)\n",
      "plt.xticks(rotation=90)\n",
      "maxVal=0.01\n",
      "plt.xlim(-0.1*maxVal,maxVal)\n",
      "plt.ylim(-0.1*maxVal,maxVal)\n",
      "plt.tight_layout()\t\n",
      "\n",
      "outfilepath=os.getcwd() + '/metaanalysis/pcbp/'\n",
      "fig5.savefig(outfilepath+'HCV_fullGenomeCoor_Supp.png',format='png',bbox_inches='tight',dpi=200,pad_inches=0.5)\n",
      "fig5.savefig(outfilepath+'HCV_fullGenomeCoor_Supp.pdf',format='pdf',bbox_inches='tight',dpi=200,pad_inches=0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# SUPPLEMENT\n",
      "\n",
      "fig6=plt.figure(6)\n",
      "\n",
      "color='cyan'\n",
      "\n",
      "plt.subplot2grid((2,3),(0,0),colspan=3)\n",
      "for i in range(0,4): # Full genome \n",
      "    plt.bar(center,filteredHist[i],align='center',width=width,color=colors,alpha=0.3)\n",
      "    plt.tick_params(axis='x',labelsize=5)\n",
      "    plt.tick_params(axis='y',labelsize=5)\n",
      "plt.title(\"Coverage histogram for JFH-1.\",fontsize=7.5)\n",
      "plt.ylim([0,0.04])\n",
      "\n",
      "plt.subplot2grid((2,3),(1,0),colspan=1)\n",
      "for i in range(0,4):    \n",
      "    plt.bar(center,filteredHist[i],align='center',width=width,color=color,alpha=0.3)\n",
      "    plt.tick_params(axis='x',labelsize=5)\n",
      "    plt.tick_params(axis='y',labelsize=5)\n",
      "plt.xlim([400,2500])\n",
      "plt.ylim([0,0.01])\n",
      "\n",
      "plt.subplot2grid((2,3),(1,1),colspan=1)\n",
      "for i in range(0,4):    \n",
      "    plt.bar(center,filteredHist[i],align='center',width=width,color=color,alpha=0.3)\n",
      "    plt.tick_params(axis='x',labelsize=5)\n",
      "    plt.tick_params(axis='y',labelsize=5)\n",
      "plt.xlim([6000,7000])\n",
      "plt.ylim([0,0.01])\n",
      "\n",
      "plt.subplot2grid((2,3),(1,2),colspan=1)\n",
      "for i in range(0,4):    \n",
      "    plt.bar(center,filteredHist[i],align='center',width=width,color=color,alpha=0.3)\n",
      "    plt.tick_params(axis='x',labelsize=5)\n",
      "    plt.tick_params(axis='y',labelsize=5)\n",
      "plt.xlim([8000,9280])\n",
      "plt.ylim([0,0.01])\n",
      "plt.tight_layout()\n",
      "\n",
      "outfilepath=os.getcwd() + '/metaanalysis/pcbp/'\n",
      "fig6.savefig(outfilepath+'HCV_fullGenomeCoor_hist_Supp.png',format='png',bbox_inches='tight',dpi=200,pad_inches=0.5)\n",
      "fig6.savefig(outfilepath+'HCV_fullGenomeCoor_hist_Supp.pdf',format='pdf',bbox_inches='tight',dpi=200,pad_inches=0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Data frame with all HCV data\n",
      "storageDFHCV=pd.DataFrame()\n",
      "storageDFHCV['Genome_Position']=center\n",
      "storageDFHCV['JFH-1_CLIP_Rep_1']=filteredHist[0]\n",
      "storageDFHCV['JFH-1_CLIP_Rep_2']=filteredHist[1]\n",
      "storageDFHCV['JFH-1_CLIP_Rep_3']=filteredHist[2]\n",
      "storageDFHCV['JFH-1_CLIP_Rep_4']=filteredHist[3]\n",
      "# Save source data\n",
      "saveName='JFH-1_CLIP_CoverageData'\n",
      "project='pcbp'\n",
      "outfilepath=os.getcwd() + '/metaanalysis/%s/%s'%(project,saveName)\n",
      "storageDFHCV.to_csv(outfilepath)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}