running_software_solution.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:b5d82e425367c89c9a7f5bac16c493b43f8793653a1d415416fc017e7e755b39"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Running bioinformatics software\n",
      "\n",
      "## You try\n",
      "\n",
      "1. Use the help of cd-hit to learn how to cluster proteins at different sequence identity thresholds.  Run cd-hit to cluster the proteins at 100%, 80%, and 70%. \n",
      "\n",
      "2. The output fasta file only includes the consensus sequence for each cluster. How many clusters are there for each protein identity threshold?\n",
      "\n",
      "**Hints:**\n",
      "  1. Each protein has one header line that starts with '>'\n",
      "  2. You can use the unix command '`wc -l`' and '`grep`' to count the protein.\n",
      "  \n",
      "## Solution\n",
      "\n",
      "1. Use the help of cd-hit to learn how to cluster proteins at different sequence identity thresholds.  To print the help, type:\n",
      "\n",
      "        $ cd-hit -help\n",
      "        ====== CD-HIT version 4.6 (built on Jan  5 2015) ======\n",
      "        \n",
      "        Usage: cd-hit [Options]\n",
      "        \n",
      "        Options\n",
      "        \n",
      "           -i   input filename in fasta format, required\n",
      "           -o   output filename, required\n",
      "           -c   sequence identity threshold, default 0.9\n",
      "                this is the default cd-hit's \"global sequence identity\" calculated as:\n",
      "                number of identical amino acids in alignment\n",
      "                divided by the full length of the shorter sequence\n",
      "\n",
      "  Thus using the option '`-c #`' will cluster the sequences at `#` percent identity.\n",
      "  \n",
      "2. Run cd-hit with each percent identity.  Use unique output filenames so the commands don't overwrite your files.\n",
      "    1. 100%\n",
      "            $ cd-hit -i proteins.faa -o clustered_proteins_100.faa -c 1.0\n",
      "\n",
      "    2. 80%\n",
      "            $ cd-hit -i proteins.faa -o clustered_proteins_80.faa -c 0.8\n",
      "            \n",
      "    3. 70%\n",
      "            $ cd-hit -i proteins.faa -o clustered_proteins_70.faa -c 0.7\n",
      "\n",
      "3. Since each fasta output file contains one sequence per cluster, to count the number of clusters, we need only count the number of sequences.  Since each sequence starts with '`>`', use `grep` to find the first line of each sequence and then count the number of lines using `wc`.  Use `wc`'s '`-l`' option to count the number of lines.  To send the output of `grep` to the input of `wc`, use a pipe ('`|`').  The result of the command will be the number of clusters.\n",
      "    1. 100%\n",
      "            $ grep \">\" clustered_proteins_100.faa | wc -l\n",
      "            199\n",
      "            \n",
      "    2. 80%\n",
      "            $ grep \">\" clustered_proteins_80.faa | wc -l\n",
      "            125\n",
      "            \n",
      "    3. 70%\n",
      "            $ grep \">\" clustered_proteins_70.faa | wc -l\n",
      "            118   \n",
      "            \n",
      "4. Thus there are 199 clusters if a 100% protein identity threshold is used, 125 clusters if a 80% protein identity threshold is used, and 118 clusters if a 70% protein identity threshold is used."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}