Skip to content

Instantly share code, notes, and snippets.

@cfriedline

cfriedline/snp_sim

Created December 4, 2013 16:45
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save cfriedline/7790932 to your computer and use it in GitHub Desktop.
Save cfriedline/7790932 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
snp_sim
{
"metadata": {
"name": "snp_sim"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": false,
"input": "import numpy as np\nimport scipy as sp",
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "code",
"collapsed": false,
"input": "from Bio import Restriction as rst\nfrom Bio.Seq import Seq\nfrom Bio.Alphabet.IUPAC import IUPACAmbiguousDNA",
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "code",
"collapsed": false,
"input": "bases = [['A','T'],['C','G']]",
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 3
},
{
"cell_type": "code",
"collapsed": false,
"input": "def get_genome(size, gc):\n i = 0\n genome = \"\"\n while i < size:\n r = np.random.randint(0,100)\n b = bases[0]\n if r < gc:\n b = bases[1]\n genome += b[np.random.random_integers(0,1)]\n i += 1\n return genome",
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 4
},
{
"cell_type": "code",
"collapsed": false,
"input": "def get_gc(g):\n gc = 0\n for x in genome:\n if x == \"G\" or x == \"C\":\n gc += 1.0\n return gc/len(genome)",
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "code",
"collapsed": false,
"input": "genome = get_genome(1e6,60)\nget_gc(genome)",
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 6,
"text": "0.599752"
}
],
"prompt_number": 6
},
{
"cell_type": "code",
"collapsed": false,
"input": "amb = IUPACAmbiguousDNA()",
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 7
},
{
"cell_type": "code",
"collapsed": false,
"input": "g = Seq(genome, amb)",
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 8
},
{
"cell_type": "code",
"collapsed": false,
"input": "rb = rst.EcoRI + rst.MstI",
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 9
},
{
"cell_type": "code",
"collapsed": false,
"input": "rb",
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 10,
"text": "RestrictionBatch(['EcoRI', 'MstI'])"
}
],
"prompt_number": 10
},
{
"cell_type": "code",
"collapsed": false,
"input": "ana = rst.Analysis(rb, g, linear=True)",
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 11
},
{
"cell_type": "code",
"collapsed": false,
"input": "ana.with_sites()",
"language": "python",
"metadata": {},
"outputs": [
{
"ename": "AttributeError",
"evalue": "type object 'RestrictionType' has no attribute 'size'",
"output_type": "pyerr",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mAttributeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-12-21af6e87e288>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mana\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwith_sites\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;32m/data7/cfriedline/anaconda/envs/conda/lib/python2.7/site-packages/IPython/core/displayhook.pyc\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, result)\u001b[0m\n\u001b[0;32m 245\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstart_displayhook\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 246\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwrite_output_prompt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 247\u001b[1;33m \u001b[0mformat_dict\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmd_dict\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcompute_format_data\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 248\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwrite_format_data\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mformat_dict\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmd_dict\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 249\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mupdate_user_ns\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m/data7/cfriedline/anaconda/envs/conda/lib/python2.7/site-packages/IPython/core/displayhook.pyc\u001b[0m in \u001b[0;36mcompute_format_data\u001b[1;34m(self, result)\u001b[0m\n\u001b[0;32m 155\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 156\u001b[0m \"\"\"\n\u001b[1;32m--> 157\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshell\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdisplay_formatter\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 158\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 159\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mwrite_format_data\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mformat_dict\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmd_dict\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m/data7/cfriedline/anaconda/envs/conda/lib/python2.7/site-packages/IPython/core/formatters.pyc\u001b[0m in \u001b[0;36mformat\u001b[1;34m(self, obj, include, exclude)\u001b[0m\n\u001b[0;32m 150\u001b[0m \u001b[0mmd\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 151\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 152\u001b[1;33m \u001b[0mdata\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mformatter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 153\u001b[0m \u001b[1;32mexcept\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 154\u001b[0m \u001b[1;31m# FIXME: log the exception\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m/data7/cfriedline/anaconda/envs/conda/lib/python2.7/site-packages/IPython/core/formatters.pyc\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, obj)\u001b[0m\n\u001b[0;32m 479\u001b[0m \u001b[0mtype_pprinters\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtype_printers\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 480\u001b[0m deferred_pprinters=self.deferred_printers)\n\u001b[1;32m--> 481\u001b[1;33m \u001b[0mprinter\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpretty\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 482\u001b[0m \u001b[0mprinter\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflush\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 483\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mstream\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetvalue\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m/data7/cfriedline/anaconda/envs/conda/lib/python2.7/site-packages/IPython/lib/pretty.pyc\u001b[0m in \u001b[0;36mpretty\u001b[1;34m(self, obj)\u001b[0m\n\u001b[0;32m 345\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mcls\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtype_pprinters\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 346\u001b[0m \u001b[1;31m# printer registered in self.type_pprinters\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 347\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtype_pprinters\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mcls\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcycle\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 348\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 349\u001b[0m \u001b[1;31m# deferred printer\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m/data7/cfriedline/anaconda/envs/conda/lib/python2.7/site-packages/IPython/lib/pretty.pyc\u001b[0m in \u001b[0;36minner\u001b[1;34m(obj, p, cycle)\u001b[0m\n\u001b[0;32m 595\u001b[0m \u001b[0mp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtext\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m','\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 596\u001b[0m \u001b[0mp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbreakable\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 597\u001b[1;33m \u001b[0mp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpretty\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 598\u001b[0m \u001b[0mp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtext\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m': '\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 599\u001b[0m \u001b[0mp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpretty\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m/data7/cfriedline/anaconda/envs/conda/lib/python2.7/site-packages/IPython/lib/pretty.pyc\u001b[0m in \u001b[0;36mpretty\u001b[1;34m(self, obj)\u001b[0m\n\u001b[0;32m 331\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbegin_group\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 332\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 333\u001b[1;33m \u001b[0mobj_class\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'__class__'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mtype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 334\u001b[0m \u001b[1;31m# First try to find registered singleton printers for the type.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 335\u001b[0m \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m/data7/cfriedline/anaconda/envs/conda/lib/python2.7/site-packages/Bio/Restriction/Restriction.pyc\u001b[0m in \u001b[0;36m__len__\u001b[1;34m(cls)\u001b[0m\n\u001b[0;32m 322\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 323\u001b[0m length of the recognition site.\"\"\"\n\u001b[1;32m--> 324\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mcls\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 325\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 326\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__hash__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcls\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mAttributeError\u001b[0m: type object 'RestrictionType' has no attribute 'size'"
]
}
],
"prompt_number": 12
},
{
"cell_type": "code",
"collapsed": false,
"input": "import Bio\nprint Bio.__version__",
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": "1.62\n"
}
],
"prompt_number": 14
},
{
"cell_type": "code",
"collapsed": false,
"input": "import IPython\nprint IPython.__version__",
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": "1.1.0\n"
}
],
"prompt_number": 16
},
{
"cell_type": "code",
"collapsed": false,
"input": "import sys\nprint sys.version",
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": "2.7.6 |Anaconda 1.8.0 (64-bit)| (default, Nov 11 2013, 10:47:18) \n[GCC 4.1.2 20080704 (Red Hat 4.1.2-54)]\n"
}
],
"prompt_number": 17
},
{
"cell_type": "code",
"collapsed": false,
"input": "",
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment