lengths.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lytir: Plotting lengths of scaffolds\n",
    "\n",
    "This Notebook is part of the [Lytir](https://github.com/karinlag/Lytir) genome analysis package.\n",
    "This notebook is covered by the BSD license.\n",
    "Enjoy!\n",
    "\n",
    "## How to use a jupyter notebook\n",
    "\n",
    "Jupyter Notebooks consist of annotated and commented runnable code, which analyse data and \n",
    "produce tables and figures. [You can find out more about notebooks here](http://jupyter-notebook-beginner-guide.readthedocs.io/en/latest/what_is_jupyter.html)\n",
    "\n",
    "Jupyter notebooks consists of cells (squares, as you will see below) with code that you can run.\n",
    "To run what is in a cell, click on a cell so that there is a blue bar on the left side, and \n",
    "then press either Shift-Enter or Control-Enter.\n",
    "\n",
    "## What does this notebook do\n",
    "This notebook has been developed to plot the lengths of the scaffolds of two different assemblies\n",
    "as histograms. To use this notebook, you need to have two assemblies available. You need to\n",
    "replace the file names that you see below. Then, start at the top of the notebook and execute\n",
    "each cell in order.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "!conda install --yes --prefix {sys.prefix} Biopython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import random\n",
    "import numpy\n",
    "from matplotlib import pyplot\n",
    "from Bio import SeqIO\n",
    "% matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Replace filenames here!\n",
    "Below, you will find the filenames that have to be replaced. Make sure\n",
    "to use the full path to where your files are!\n",
    "\n",
    "The first file will be \"A\" in the plots below, while the second will be the \"B\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fastaA = \"../testdata/95-Q03-48_ACTCGCTA-GTAAGGAG_spades_scaffolds.fasta\"\n",
    "fastaB = \"../testdata/97-Q04-03_ACTCGCTA-AAGGAGTA_spades_scaffolds.fasta\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lengths(fastafile):\n",
    "    lengths = []\n",
    "    with open(fastafile, \"r\") as handle:\n",
    "        for record in SeqIO.parse(handle, \"fasta\"):\n",
    "            lengths.append(len(record))  \n",
    "    lengths.sort()\n",
    "    return(lengths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lengthsA=lengths(fastaA)\n",
    "lengthsB=lengths(fastaB)\n",
    "maxval = max(lengthsA)\n",
    "if max(lengthsB) >= max(lengthsA):\n",
    "    maxval = max(lengthsB) \n",
    "bins = numpy.linspace(0, maxval, 50)\n",
    "pyplot.hist(lengthsA, bins, alpha=0.5, label='A')\n",
    "pyplot.hist(lengthsB, bins, alpha=0.5, label='B')\n",
    "pyplot.legend(loc='upper right')\n",
    "pyplot.xlabel(\"Lengths of scaffolds\")\n",
    "pyplot.title(\"Scaffold lengths for both assemblies, distributed over 50 buckets\")\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxval = 100000\n",
    "bins = numpy.linspace(0, maxval, 50)\n",
    "pyplot.hist(lengthsA, bins, alpha=0.5, label='A')\n",
    "pyplot.hist(lengthsB, bins, alpha=0.5, label='B')\n",
    "pyplot.legend(loc='upper right')\n",
    "pyplot.xlabel(\"Lengths of scaffolds\")\n",
    "pyplot.title(\"Scaffold lengths for both assemblies for scaffolds below 100 000 bps\")\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxval = 10000\n",
    "bins = numpy.linspace(0, maxval, 50)\n",
    "pyplot.hist(lengthsA, bins, alpha=0.5, label='A')\n",
    "pyplot.hist(lengthsB, bins, alpha=0.5, label='B')\n",
    "pyplot.legend(loc='upper right')\n",
    "pyplot.xlabel(\"Lengths of scaffolds\")\n",
    "pyplot.title(\"Scaffold lengths for both assemblies for scaffolds below 10 000 bps\")\n",
    "pyplot.show()"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}