test-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dda2ce31-2259-4dd0-9388-22c27573e3a7",
   "metadata": {},
   "source": [
    "# answers 00: _Moriarty's game of chaos_\n",
    "\n",
    "Welcome to MCB112, Python, and Jupyter Notebook!\n",
    "\n",
    "The week 0 (\"w00\") pset is an example that gives you a walkthrough of what an MCB112 pset looks like, and a feeling for what we expect an answer to look like. Starting with w01, you'll be turning in your pset answers in the form of a single Jupyter Notebook page (a `.ipynb` file) like this one. Jupyter Notebook allows you to interleave Python code with text, written in Markdown format, a simple markup language for text.\n",
    "\n",
    "This notebook will be way more wordy than you need to be in an answer though, because I'm going to explain a bunch of stuff in a lot of detail, for people who have never seen Python before. I will set these long explanations off in their own _explanations_ sections so they don't interrupt the flow _too_ much.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6cb24c5-bfc7-4740-8128-f0c5ba915a55",
   "metadata": {},
   "source": [
    "## our imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "9146a832-d097-4897-a5e1-30eabda3280b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                    \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline                     "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fb2f8af-f688-4c87-b18e-2acb30d41b33",
   "metadata": {},
   "source": [
    "## _explanation_\n",
    "\n",
    "This example uses the [NumPy](https://numpy.org/doc/stable/index.html) module for scientific numeric computing, and the [matplotlib](https://matplotlib.org/) module for data visualization. The w01 pset will use only basic python, and you won't need to know or use NumPy or matplotlib until w02.\n",
    "\n",
    "We will also get to the [Pandas](https://pandas.pydata.org/) data analysis module in w03, and some of the [SciPy](https://scipy.org/) module later in the course.\n",
    "\n",
    "For each pset, we'll state which modules you can use. We don't want you to learn to code by googling and finding a bunch of random Python packages. At least not yet. There's a time and place for that, but get your fundamentals down first.\n",
    "\n",
    "The `%matplotlib inline` sets your notebook page so it will render matplotlib figures as (static) images embedded in the notebook page. In fancier notebooks, you can use `%matplotlib widget` (among other options) to use _interactive_ plots. We're not going to do that, at least not for now."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7349be22-041e-446c-ab3c-815cb73a3b77",
   "metadata": {},
   "source": [
    "## 1. \"Write Python code with functions for...\"\n",
    "\n",
    "The pset conveniently suggested one good way of writing modular and reusable code. We can guess that we're going to have sequences that we read from the Genbank files and random sequences that we generate ourself, so we want to separate running the chaos game algorithm on an input sequence from the different ways we might obtain that sequence. It's also useful to separate _generating_ the data from _plotting_ the data, to make our plotting code easier to read, and easier to steal from in the future when we want to make some other kind of 2D scatterplot."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c723741-025b-439c-af1e-9f13d6a1d6e3",
   "metadata": {},
   "source": [
    "###  a) reading one DNA sequence from a Genbank format file\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1b98e28d-1ca6-4e9a-ad2a-d6f37708866b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_genbank_file(seqfile):\n",
    "    \"\"\"Read one sequence from a Genbank-format file.\n",
    "\n",
    "    Input:  <seqfile>: path to a Genbank format file containing one DNA sequence.\n",
    "    Output: <seq>:     sequence, as a list of upper-case characters (typically A|C|G|T)\n",
    "\n",
    "    If the sequence contains characters other than A|C|G|T, such as degeneracy codes like N|R|Y, these will still be included.\n",
    "    \"\"\"\n",
    "    in_sequence = False\n",
    "    seq         = []\n",
    "    with open(seqfile) as f:\n",
    "        for line in f:\n",
    "            if   line.startswith('//'):     break\n",
    "            elif in_sequence:               seq.extend( [ x.upper() for x in line if x.isalpha() ] )\n",
    "            elif line.startswith('ORIGIN'): in_sequence = True\n",
    "    return seq\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4a43250-12df-4359-a56a-635519736699",
   "metadata": {},
   "source": [
    "### _explanation_\n",
    "\n",
    "If you're new to Python, the tricksiest line in this code was the line with `seq.extend()` on it. Let's dwell on that for a second because we can learn a fair amount of Python all at once.\n",
    "\n",
    "`seq` is a **list**. You can see that because of how we initialize it (`seq = []`). We can look at the documentation for Python lists and see that we can add to them either one item at a time with `append()` or we can add a list of several items at once with `extend()`. Variables in Python (lists, and integers, and floats, and dicts, and whatever) are _objects_ with _methods_: so when we do something like `seq.extend()` we are saying \"invoke the `extend()` method in the `seq` object\", and the stuff in the parentheses is the argument to `extend()`: it will be interpreted as a list. \n",
    "\n",
    "Here, it's a list of the A|C|G|T characters on one line of the GenBank file. You can see us open the file (`open(seqfile) as f`) so that we have an open **file object** `f`, and then we're reading one line at a time from the file in a loop (`for line in f`) in a way where the variable `line` is a **string** containing something like \"1240 ACGTACGATCAG ACGACTGATCTCA\", because lines in Genbank files start with coordinates and also have spaces. We're waiting in the file until the sequence part is about to start (after a line that says 'ORIGIN'), and we stop when there's a line that says '//', by checking `line` with a string method `.startswith()` that tests for exact matching to a query string at the start of `line`. So that's what the `in_sequence` flag is doing: it gets set to `True` when we're in the sequence part of the file.\n",
    "\n",
    "OK, back to the argument to `seq.extend()`, which needs to be a list, but `line` is a string, not a list. We also want to skip the coordinates and the spaces in the line, and only extend `seq` with a list of A|C|G|T characters on `line`. We can iterate over each character in `line` with a for loop: `for x in line`. We could write a few lines of code to initialize a new list, iterate over each character x in `line`, check whether the character is A|C|G|T before we include it, and convert the character to upper case (if it's a lower-case a|c|g|t), but we can do all of this very compactly using something called a **list comprehension**. If we did `[x for x in line]` we would get a list of the characters in `line`. The first `x` there is \"put this thing in the new list...\" and the `for x in line` is \"... for each successive x in line\". We can put anything we want into the new list as we do the `for x in line`, including some sort of operation on `x`: here, we're using the `x.upper()` method to convert `x` to upper case. Finally, Python's list comprehension syntax allows us to append an `if` test for each iteration of the for loop: here, `if x.isalpha()`, where `isalpha()` is a method with no argument that returns True if `x` is an alphabetic character: this is the piece of our code that is skipping spaces and numbers.\n",
    "\n",
    "Whew. Got that?\n",
    "\n",
    "No, it would be easy to still be confused. How should you figure this out as you're learning? The technique I advise to use throughout the course (and indeed throughout life) is, whatever you're trying to understand, create the smallest possible example where you know what the answer is supposed to be, and make sure you can do that first. Then move outwards to more and more complex things, once you've nailed the small things. And in the case of Python: use a lot of print statements when you do it.\n",
    "\n",
    "So here's an example; let's peel apart that list comprehension in a code block of its own in our Markdown, where we can play.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "7c53cf44-8eb3-4d23-9be6-ca4eb96e42e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'str'>\n",
      "<class 'list'>\n",
      "['1', '2', '0', ' ', 'A', 'C', 'T', ' ', 'g', 'a', 't']\n",
      "['hey!', 'hey!', 'hey!', 'hey!', 'hey!', 'hey!', 'hey!', 'hey!', 'hey!', 'hey!', 'hey!']\n",
      "['1', '2', '0', ' ', 'A', 'C', 'T', ' ', 'G', 'A', 'T']\n",
      "['A', 'C', 'T', 'G', 'A', 'T']\n",
      "['A', 'C', 'T', 'G', 'A', 'T']\n"
     ]
    }
   ],
   "source": [
    "line = '120 ACT gat'                # Forget all the Genbank file stuff. Make sure we can turn this one little line into list('ACTGAT').\n",
    "print( type(line) )                 # The built-in method type() is your friend for figuring out what things are: here, line is a `str` class\n",
    "L    = [x for x in line]            # Simple list comprehension. \n",
    "print( type(L) )                    # The result L is a list of the characters in `line`\n",
    "print (L)      \n",
    "L    = [ 'hey!' for x in line ]     # Nothing sacred about that first x. We can put anything into the new list as we iterate through `line`\n",
    "print (L)\n",
    "L    = [ x.upper() for x in line ]  # but usually it would be something to do with the `x` we're iterating over in line, of course\n",
    "print (L)\n",
    "L    = [ x.upper() for x in line if x.isalpha() ]   # and finally here's what we did to make the argument we passed to `seq.extend()`\n",
    "print (L)\n",
    "L    = list('ACTGAT')               # if it weren't for the numbers and spaces - converting a string to a list is easy with `list()`\n",
    "print (L)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a00f19a2-c35d-47a0-a15f-74f0e5c92796",
   "metadata": {},
   "source": [
    "A notebook cell like that is a (mostly) safe space where you can try different things, run the cell, and see the result quickly. \n",
    "\n",
    "(I say \"mostly\" safe because if you change variables that you're using somewhere else in the notebook, you can break things.)\n",
    "\n",
    "Now back to our regular programming..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1cdabb0-c56b-4d9a-b566-4c7844229645",
   "metadata": {},
   "source": [
    "### b) doing the chaos game algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "2c65234f-a7b0-4afe-a4ba-4caca1fca68c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def chaos_game(seq):\n",
    "    \"\"\"Given a DNA sequence, produce a chaos game representation, as a 2D array of coordinates\n",
    "\n",
    "    Input:  <seq> : a list or array of A|C|G|T characters\n",
    "    Output: 2D numpy array of shape (L x 2): array of 2D coords for each of L points\n",
    "\n",
    "    Assumes that the <seq> list only contains A|C|G|T characters.\n",
    "    \"\"\"\n",
    "    vertices = { 'A' : [0, 0], 'C': [0,1], 'G': [1,1], 'T': [1,0] }\n",
    "    L = len(seq)\n",
    "\n",
    "    data = np.zeros( (L,2) )\n",
    "    data[0] = vertices[seq[0]]\n",
    "    for i in range(1,L):\n",
    "        data[i] = (data[i-1] + vertices[seq[i]]) / 2\n",
    "    return data\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2240fb4d-caaa-4163-a006-25873a49e815",
   "metadata": {},
   "source": [
    "### _explanation_\n",
    "\n",
    "One thing to notice here is that `data[i]` is row `i` of a 2D matrix, so `data[i]` is itself a vector of length 2. Meanwhile we've created `vertices` to be a dict, indexed by the four DNA characters: `vertices['A']` is a vector of length 2, the coords of the 'A' vertex. So when we do `data[i] + vertices[seq[i]]`, we are doing **vector addition**, an example of something that NumPy excels at both in syntax and in internal speed. Add some `type()` and `print()` statements in there to watch it work - or again, create an isolated small example in a new cell:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "7fbe2b54-b504-48d6-bc9a-aebcefeaa6cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'numpy.ndarray'>\n",
      "(3, 2)\n",
      "<class 'numpy.ndarray'>\n",
      "(2,)\n",
      "<class 'numpy.ndarray'>\n",
      "(2,)\n",
      "[0.7 0.7]\n"
     ]
    }
   ],
   "source": [
    "D = np.array( [ [0.1, 0.1 ],       # Make a small example of the Lx2 array\n",
    "                [0.2, 0.2 ], \n",
    "                [0.3, 0.3 ] ] )\n",
    "vertex = np.array([0.5, 0.5])     # and an example vertex\n",
    "print(type(D))                    # `D` is a numpy array\n",
    "print(D.shape)                    #   ... of \"shape\" 3x2\n",
    "print(type(D[1]))                 # `D[1]` is one row of the array.. but still just another numpy array\n",
    "print(D[1].shape)                 #   ... now of shape 1x2, but since it's a vector, numpy interprets it as `2,`\n",
    "print(type(vertex))               # `vertex` is also a numpy array\n",
    "print(vertex.shape)               #   ... also of shape `2,`\n",
    "r = D[1] + vertex                 # the two vectors have the same shape, so we can add them, elementwise\n",
    "print (r)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36f35a32-95ee-43d1-a6f9-4137e476cf55",
   "metadata": {},
   "source": [
    "And let's also have a closer look at `vertices`. Are you confused about Python `dict` versus `list` objects? I sure was at first. The key idea is the following:\n",
    "\n",
    "A `list` is an ordered set of `n` elements (of any type) that you index by **position** `0..n-1` in the list: `L[i]` gives you the i'th element in `L`.\n",
    "\n",
    "A `dict` is an unordered set of elements (of any type) that you index by a **string**: `D['A']` gives you the element in `D` that you indexed with key `A`.\n",
    "\n",
    "So I used a dict for `vertices` because I wanted to say `vertices['C']`, looking up the coords for the `C` vertex using `C` itself as the index. \n",
    "\n",
    "Again, play with creating and manipulating lists and dicts in an isolated cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "16cff3bb-5634-4eda-9986-49ab3f9aea30",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\n",
      "postdoc\n",
      "[10, 20, 30, 40]\n",
      "{'Holmes': 'postdoc', 'Watson': 'student', 'Wiggins': 'Baker Street Irregular'}\n"
     ]
    }
   ],
   "source": [
    "L = [ 10, 20, 30 ]\t\t\t\t\t\t         \t# lists: initialization\n",
    "D = { 'Holmes': 'postdoc', 'Watson': 'student' }  \t# dicts: initialization\n",
    "print(L[1])\n",
    "print(D['Holmes'])\n",
    "L.append(40)                              # append to the end of a list\n",
    "print(L)\n",
    "D['Wiggins'] = 'Baker Street Irregular'   # add element to a dict\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d4a3de5-f690-43ea-95b2-99040e7ccd41",
   "metadata": {},
   "source": [
    "and again back to our regular programming..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ae7b51b-8697-4be8-86a1-85a921587e3f",
   "metadata": {},
   "source": [
    "### c) plotting the chaos game data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "5a7cb2f8-0598-4030-aa55-1f25683a0844",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_chaos(data):\n",
    "    \"\"\"Create a scatterplot figure from an Lx2 array of chaos game points\n",
    "\n",
    "    Input <data>: Lx2 array of coords for L points\n",
    "    Result:       an inline figure in our notebook\n",
    "    \"\"\"\n",
    "    fig, ax = plt.subplots()                                                # Creates Figure, Axes objects to plot in\n",
    "    ax.scatter(data[:,0], data[:,1], s=0.1, marker='.', color='black')      # scatter plot with black points of size 0.1\n",
    "\n",
    "    ax.text(0, 1, 'C', horizontalalignment='right')                         # here we're labeling A,C,G,T vertices\n",
    "    ax.text(1, 1, 'G')\n",
    "    ax.text(0, 0, 'A', verticalalignment='top', horizontalalignment='right')\n",
    "    ax.text(1, 0, 'T', verticalalignment='top')\n",
    "\n",
    "    ax.axis('equal')      # make sure the plot is shown as a square, not a rectangle\n",
    "    ax.set_axis_off()     # turn off everything associated with the axes: tick marks, labels, spines\n",
    "\n",
    "    fig.tight_layout()    # tighten margins, reduce unnecessary white space\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e184c726-83e2-48c9-b4e2-b7486f7f52f7",
   "metadata": {},
   "source": [
    "### _explanation_\n",
    "\n",
    "One thing to notice here: `data[:,0]` is what's called a **slice** of an array. matplotlib (and Pandas) are very column-oriented in how they think about 2D data tables. The scatter plot function `ax.scatter()` expects to get a list/array of `n` x coordinates and a corresponding list of `n` y coordinates for the points it will plot; but our data are in an Lx2 array. `data[:,0]` is the syntax for \"give me all of the rows\" (that's the ':') and give me just column 0 (that's the '0'): i.e., extract the X coords as a single array. `data[:,1]` are the y coords. You know the drill: look up array slicing, make a small example in an isolated cell, and play with lots of print statements.\n",
    "\n",
    "matplotlib is big and complicated and I won't say much about it for now; best to start with their tutorial, to understand the syntax. But here's one thing to know up front. There are two main ways of using matplotlib. In one way, you basically assume that you're just making a simple single figure, and you directly call matplotlib functions with `plt.` prefixes, like `plt.plot()`. The other way is \"object oriented\", where you create and manipulate Figure and Axes objects. That's what I'm doing here: creating a Figure and Axes object for my figure with `fig, ax = plt.subplots`, and then using object-oriented methods calls on those objects with things like `ax.scatter()` to create and design my plot. I recommend following this pattern and using the object-oriented style, becuse it's easier to control in something like a jupyter notebook where you have lots of figures at once. I _don't_ recommend mixing the two styles, which is what you will get if you learn matplotlib by googling StackOverflow.\n",
    "\n",
    "OK, enough with the code explanations. We can proceed with the biology now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76661518-d086-4e36-b72d-9e956a8da7e3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}