Lab 1. Python, Jupyter, and characterizing nucleotide composition.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Biol526H Computational Genetics, Fall 2016\n",
    "## Lab 1: Python, Jupyter and characterizing nucleotide composition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The relationship of recitations, computer labs, and assignments\n",
    "This computer lab, like the others that are scheduled for recitations, is structured as a tutorial.  Complete as much as you can of the tutorial during the recitation session.  I suggest going through the tutorial individually, but discussing it in groups of 2-3 as you go.\n",
    "\n",
    "Each lab will be have an accompanying assignment. For this first lab, the assignment consists of only one problem.  Look ahead, and take advantage of the recitation time to ask questions so you can complete the assignment later. \n",
    "\n",
    "When you are done, you should upload your answers (typically as an .ipynb notebook, explained below) to your dropbox on Sakai. The answers are due before the start of the following lab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Our toolkit\n",
    "For these labs, we will be using [Python](https://docs.python.org/3.1/tutorial/index.html), which is free, available on all major platforms, powerful, high-level, well documented, relatively easy to get started, and widely used in the scientific community. Python is an interpreted language. You can enter commands interactively, or you can store commands in a file or as a function and execute them within another block of code.\n",
    "\n",
    "We will be taking advantage of functionality from [Biopython](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc2), a large library of Python objects and methods that simplify working with various bioinformatic datatypes, formats, libraries and tools. \n",
    "\n",
    "Finally, we will be doing our work within [Jupyter notebooks](http://jupyter-notebook-beginner-guide.readthedocs.io/en/latest/what_is_jupyter.html).  Jupyter provides an interactive environment for documentation, coding and analysis and gives us files that are both human readable and machine executable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Installation and getting started\n",
    "Follow the instructions here to obtain the [Anaconda distribution](http://jupyter.readthedocs.io/en/latest/install.html) with both Python 3 and Jupyter. The page also provides instructions on using pip to install Juypter if you already have Python 3 installed.\n",
    "\n",
    "Once you have it installed, you can follow along the online \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting started\n",
    "Once you have Anaconda installed, \n",
    "\n",
    "How to get help. How to write a comment.  How to save and reopen.\n",
    "\n",
    "Simple commands.\n",
    "\n",
    "Data types.  Lists, control, logical operators, etc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Characterizing a genome\n",
    "\n",
    "*Haemophilus influenzae* is a bacterium first described in 1892 by Dr. Robert Pfeiffer during an influenza pandemic and responsible for a variety of clinical diseases. Its genome was the first to be sequenced and assembled in a free-living organism. It contains about 1.8 million\n",
    "base pairs and is estimated to have 1,740 genes.\n",
    "\n",
    "Today we will just be exploring the genomic DNA sequence. To what else a genome record typically includes, go to [Genbank][http://www.ncbi.nlm.nih.gov] and search for accession number NC000907.  Take a moment to explore the item found in the nucleotide database and look over its contents, in the standard Genbank format, in the graphics display, and in the (barebones) FASTA text format, which we will be working with today. \n",
    "\n",
    "[Obtain the data.\n",
    "\n",
    "Packaged by textbook at: \n",
    "http://www.computational-genomics.net/case_studies/haemophilus_demo.html. \n",
    "\n",
    "Save it to a directory that you can easily find and unpack it. Now set your path to the local directory on your computer containing the files.]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-7-95a90014cf66>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-7-95a90014cf66>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    Hflu = getgenbank(`NC_000907',`SequenceOnly',true);\u001b[0m\n\u001b[0m                      ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "# Get sequence from Genbank\n",
    "Hflu = getgenbank(`NC_000907',`SequenceOnly',true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-8-c81f14d254a0>, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-8-c81f14d254a0>\"\u001b[0;36m, line \u001b[0;32m2\u001b[0m\n\u001b[0;31m    whos Hflu\u001b[0m\n\u001b[0m            ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "# What kind of variable is it, and how long is the sequence?\n",
    "whos Hflu\n",
    "# What is the composition of nucleotides?\n",
    "basecount(Hflu)\n",
    "# "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an aside, take a look at the function you executed to see how you can build interactive commands into more sophisticated functions that can later be reused."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "type basecount\n",
    "# Note that others symbols occur in the sequence besides A, C, G and T. They\n",
    "correspond to various kinds of ambiguity codes:\n",
    "# N: aNy base\n",
    "# K: G or T\n",
    "# R: A or G (puRine)\n",
    "# Y: C or T (pYrimidine)\n",
    "# M: A or C (aMino)\n",
    "# S: C or G\n",
    "# W: A or T\n",
    "# basecount(Hflu,`Others',`Full')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The DNA sequence you have just looked at is that of the 5'-3' strand. Thanks to Watson-Crick base pairing, you should be able to guess the composition of the other strand already. Can you? Verify whether you are correct by examining the reverse complement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-9-6f45c75d2376>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-9-6f45c75d2376>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    > HfluRC = seqrcomplement(Hflu);\u001b[0m\n\u001b[0m    ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "HfluRC = seqrcomplement(Hflu);\n",
    "basecount(HfluRC,`Others',`Full')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Local fluctuations in the frequencies of nucleotides can be informative about features such as gene composition or horizontally transferred segments of DNA. You can often see patterns in the local base composition with a sliding window of variable size. Take a look using windows of size 10,000 bp. This may be very slow, so give it some time to think, and if it is taking too long, type `Ctrl-C'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-10-dfbdffd28585>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-10-dfbdffd28585>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    ntdensity1(Hflu,`window',10000)\u001b[0m\n\u001b[0m                    ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "ntdensity1(Hflu,`window',10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A different sequence and a different approach\n",
    "Let us now do some basic compositional analysis of a somewhat shorter DNA\n",
    "sequence. The data come from a large-insert clone from the nuclear genome of the plant *Mimulus guttatus*.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "> Mim = getgenbank(`AC182567',`SequenceOnly',true);\n",
    "> whos Mim\n",
    "> basecount(Mim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-11-be48e2a373f2>, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-11-be48e2a373f2>\"\u001b[0;36m, line \u001b[0;32m2\u001b[0m\n\u001b[0;31m    begin = (0:140)*1000+1;\u001b[0m\n\u001b[0m              ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "# Make two vectors that will be used to break the string into 1 Kb chunks\n",
    "begin = (0:140)*1000+1;\n",
    "finish = (10:150)*1000+1;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Break the DNA sequence into 1 Kb pieces\n",
    "> w = zeros(0);\n",
    "for i = 1:size(begin,2);\n",
    "w = cat(1, w, Mim(begin(i):finish(i)));\n",
    "end;\n",
    "whos w\n",
    "# Now measure the base frequency of each piece\n",
    "d = zeros(size(begin,2),4);\n",
    "for i = 1:size(begin,2);\n",
    "d(i,:) = transpose(cell2mat(struct2cell(basecount(w(i,:)))));\n",
    "end;\n",
    "whos d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-12-adffd8fe97d2>, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-12-adffd8fe97d2>\"\u001b[0;36m, line \u001b[0;32m2\u001b[0m\n\u001b[0;31m    > d\u001b[0m\n\u001b[0m    ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "# The matrix should now contain frequencies of A, C, G and T \n",
    "# for each window\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# the frequency of cytosines alone is\n",
    "d(:,2)\n",
    "# And GC content is the sum of the middle two columns\n",
    "> GCfreq = d(:,2) + d(:,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# In MATLAB, for loops can be very slow\n",
    "# but matrix calculations are fast\n",
    "# so to speed things up, we can use `arrayfun' ,\n",
    "# which applies a function over a full array\n",
    "# then apply `basecount' to each window using `cellfun', like so\n",
    "w = arrayfun(@(x) Mim(x:(x+1000)), begin, `UniformOutput', false);\n",
    "d = arrayfun(@(x) basecount(x,:), w, `UniformOutput', false));\n",
    "gc = arrayfun(@(x) x.G + x.C, w, `UniformOutput', false);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# make a plot of GC content as a function of position\n",
    "plot(begin, gc); xlabel(`Position'); ylabel(`GC');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If all that worked for you, then how would you describe the scale and magnitude of the variability in GC content along this sequence?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}