Lesson3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison, looping, and plotting output\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## loops\n",
    "\n",
    "A lot of times we use computer programs to simplify, or speed up a task that is repetitive, or tedious to do by hand (remember translating RNA to protein, or finding START codons?)\n",
    "\n",
    "In order to do these things over and over and OVER, we use something called a 'loop'\n",
    "\n",
    "## for-loops"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# start with a list\n",
    "x = [1,2,3,4,5,6,7]\n",
    "\n",
    "# define a looping function:\n",
    "def times2(listofnumbers):\n",
    "    for number in listofnumbers:\n",
    "        print(number*2)\n",
    "        \n",
    "# call that function on our list\n",
    "times2(x)\n",
    "\n",
    "# Can you modify our times2 function to multiply by a different number?\n",
    "## can you modify it to take an argument, and multiply everything in the list by that argument?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ranges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = range(10)\n",
    "print(type(x))\n",
    "print(x)\n",
    "\n",
    "for thing in x:\n",
    "    print(thing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reuse, reuse, and reuse!!\n",
    "\n",
    "\n",
    "Crazy thing about changes to the genetic code: They add up over time and generations. So now what we need is a way to allow for multiple mutations, over a number of generations.\n",
    "\n",
    "We already wrote a function that would take a DNA sequence, and add a mutation at a random location.. But we would like to simulate *multiple* generations, and multiple DNA replication events.\n",
    "\n",
    "We can do this using a **loop!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## start by importing our mutate function from Lesson2.1\n",
    "from mutate import mutate\n",
    "\n",
    "## define a function that takes 2 arguments: an original sequence, and the number of generations\n",
    "def genetic_drift(seq, generations):\n",
    "    for generation in range(generations):\n",
    "        seq = mutate(seq)\n",
    "        # can you make this function display changes in each generation?\n",
    "    return seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## now we can test our new function!\n",
    "\n",
    "## start with a genetic sequence: \n",
    "dna1 = 'ATGAAAAAAAAAAA'\n",
    "from dna_analyze import dna_analyze\n",
    "#check our input dna sequence:\n",
    "dna_analyze(dna1)\n",
    "\n",
    "# now let's MUTATE that sequence!\n",
    "dna2 = genetic_drift(dna1,1)\n",
    "# check our new mutated sequence:\n",
    "dna_analyze(dna2)\n",
    "\n",
    "# Did you lose your start codon?\n",
    "# how many generations does it take?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Biopython!\n",
    "\n",
    "Biopython is a library that contains lots of functions relevant to biology and genetic computations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Alphabet import generic_dna\n",
    "from Bio.Seq import Seq\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio import SeqIO\n",
    "\n",
    "testseq = Seq('AAATGAAAA')\n",
    "record = SeqRecord(testseq,id='test1')\n",
    "\n",
    "print(record.seq)\n",
    "\n",
    "sequences = record  \n",
    "SeqIO.write(sequences, \"example.fasta\", \"fasta\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## rewrite your genetic drift function using conversions to biopython's bio.seq class\n",
    "\n",
    "from mutate import mutate\n",
    "\n",
    "## Need to add another input parameter:  name\n",
    "def genetic_drift(seq_record, generations, name):\n",
    "    raw_seq = seq_record.seq\n",
    "    for generation in range(generations):\n",
    "        new_seq = mutate(raw_seq)\n",
    "        new_record = SeqRecord(Seq(new_seq),id=name)\n",
    "    return new_record\n",
    "\n",
    "a = genetic_drift(record, 100, 'testmutation')\n",
    "print(a.seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import Entrez\n",
    "from Bio import SeqIO\n",
    "\n",
    "Entrez.email = \"A.N.Other@example.com\"\n",
    "with Entrez.efetch(db=\"nucleotide\", rettype=\"fasta\", retmode=\"text\", id=\"841163\") as handle:  #working: 6273291\n",
    "    seq_record = SeqIO.read(handle, \"fasta\")\n",
    "print(\"%s with %i features\" % (seq_record.id, len(seq_record.features)))\n",
    "print(seq_record.seq)\n",
    "entrez = SeqRecord(seq_record.seq,id=str(seq_record.id))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Phylogenetic trees!\n",
    "\n",
    "So now that we can mutate a DNA sequence, and store that sequence as a SeqRecord object.  Now let's see how those genetic changes add up over generations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## start with an original sequence (a common ancestor.)\n",
    "# OG_seq = SeqRecord(Seq('AATGACGAGTGAGCGTAACAACGTTTA'),id='OG')\n",
    "OG_seq = entrez\n",
    "## feed in, and change the number of generations\n",
    "# newDNA = function(originalDNA, number_of_gens, 'new name')\n",
    "a1 = genetic_drift(OG_seq, 10, 'a1')\n",
    "a2 = genetic_drift(OG_seq, 100, 'a2')\n",
    "b1 = genetic_drift(a1, 1000, 'b1')\n",
    "b2 = genetic_drift(a1,100, 'b2')\n",
    "c1 = genetic_drift(a2, 20, 'c1')\n",
    "c2 = genetic_drift(a2, 10, 'c2')\n",
    "\n",
    "records = OG_seq, a1,a2, b1, b2, c1, c2\n",
    "\n",
    "SeqIO.write(records, \"OG.fasta\", \"fasta\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering!\n",
    "\n",
    "To compare different sequences, we use statistics to line up our sequences, and score them based on how similar they are to each other\n",
    "Sometimes the program you want to run is not coded in Python directly.  Luckily, people have figured out ways to run non-python programs from the command line using python code and without ever leaving the Python coding space!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## The BioPython package contains commands that allow us to execute external programs\n",
    "\n",
    "import Bio.Align.Applications\n",
    "\n",
    "from Bio.Align.Applications import ClustalwCommandline\n",
    "# dir(Bio.Align.Applications)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cline = ClustalwCommandline(\"/data/clustalw-2.1/src/clustalw2\", infile=\"OG.fasta\")\n",
    "stdout, stderr = cline()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from Bio import AlignIO\n",
    "alignment = AlignIO.read(\"OG.aln\", \"clustal\")\n",
    "print(alignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trees!\n",
    "\n",
    "Aligning sequences is great, but it's difficult to see the actual differences between raw sequence alignments, so we use the differences in raw sequences to draw \"Phylogenetic trees\"  These trees help us visualize the differences between our sequences!\n",
    "\n",
    "We will be using the alignment file (OG.dnd) that our ClustalW program generated, and displaying it using `ascii` characters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from Bio import Phylo\n",
    "tree = Phylo.read(\"OG.dnd\", \"newick\")\n",
    "Phylo.draw_ascii(tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}