rnaseqc_gc_bias_generator.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pybedtools\n",
      "import gffutils\n",
      "from Bio import SeqIO"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 38
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gffutils.__version__"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "'0.7'"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hg19_base = \"/nas3/yeolab/Genome/ensembl/gtf/gencode.v17.annotation\"\n",
      "hg19_gff = hg19_base + \".gtf\"\n",
      "hg19_db = hg19_gff + \".db.old\"\n",
      "db = gffutils.FeatureDB(hg19_db)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "transcripts = db.features_of_type('transcript')\n",
      "transcript_lst = []\n",
      "for transcript in transcripts:\n",
      "    transcript_lst.append(pybedtools.create_interval_from_list([transcript.chrom, \n",
      "                                                                str(transcript.start), \n",
      "                                                                str(transcript.stop), \n",
      "                                                                transcript.attributes['transcript_id'], \n",
      "                                                                \"0\", \n",
      "                                                                transcript.strand]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = pybedtools.BedTool(transcript_lst)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sequences = x.sequence(\"/nas3/yeolab/Genome/ucsc/hg19/chromosomes/all.fa\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 37
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lncRNA = SeqIO.parse(open(\"/nas3/yeolab/Genome/ensembl/gtf/gencode.v17.lncRNA_transcripts.fa\"), \"fasta\")\n",
      "pc_RNA = SeqIO.parse(open(\"/nas3/yeolab/Genome/ensembl/gtf/gencode.v17.pc_transcripts.fa\"), \"fasta\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with open(\"gencode.v17.gc.txt\", 'w') as outfile:\n",
      "    for rna in lncRNA:\n",
      "        outfile.write(\"\\t\".join([rna.id.split(\"|\")[0], str((rna.seq.count(\"A\") + rna.seq.count(\"G\")) / float(len(rna.seq)))]) + \"\\n\")\n",
      "    for rna in pc_RNA:\n",
      "        outfile.write(\"\\t\".join([rna.id.split(\"|\")[0], str((rna.seq.count(\"A\") + rna.seq.count(\"G\")) / float(len(rna.seq)))]) + \"\\n\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 57
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}