python-bioinformatics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Python for bioinformatics, and with bioinformatics packages\n",
    "\n",
    "## Python versus Bash for building toolchains/pipelines\n",
    "\n",
    "### The Bash approach\n",
    "\n",
    "- In Bash/Shell it's straightforward to build toolchains using pipes, and software found in your $PATH\n",
    "- This is useful to sequentially process outputs with specialist programs or operations for data manipulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hi\n",
      "4\n",
      "<class 'IPython.utils.text.SList'>\n",
      "<class 'str'>\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "# This environment runs on a Python interpreter, so it doesn't directly run shell commands.\n",
    "# But since it uses an IPython interpreter, prefixing with \"!\" allows shell commands to be used, e.g.:\n",
    "!echo \"Hi\"\n",
    "\n",
    "# Toolchains in Bash\n",
    "# Bash/Shell has the pipe character \"|\", which allows quite complex toolchains to be built up. Programs in $PATH can be directly called, e.g.:\n",
    "!echo -n \"1234\" | wc -m   # Count characters from echo command. -n suppresses the behaviour to also add a newline, ensuring only string characters counted.\n",
    "\n",
    "# Side note, we can integrate this with python, though it doesn't create a Python data type directly\n",
    "var = !echo \"1234\" # Python mixed with Shell\n",
    "print (type(var)) # Not a regular python data type, an Slist\n",
    "str_var = str(var[0]) # Converts first element of Slist (\"1234\") into a python string\n",
    "print(type(str_var)) # Evaluates as a python string, good!\n",
    "print(str_var[0]) # Works as expected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Python equivalent of pipes\n",
    "\n",
    "- In Python, instead of pipes we might use the \"subprocess\" module\n",
    "- This allows external commands to be run, but can be pretty wordy compared to a simple: \"echo -n \"1234\" | wc -m\" (see below example)\n",
    "- So sometimes in pipeline development, it can make sense to use Shell over Python for ease of both writing & reading, while bringing in Python for more advanced data manipulations & operations\n",
    "- Keep in mind that in the workflow paradigm Nextflow, you can integrate the two. Each process code block can have a shebang for setting the code interpreter (if none is added, default is Shell)\n",
    "- Alternatively, you can have a Nextflow Shell code block, and call Python as an external programme from within it, pointing at a script file, e.g.,: `python scripts/function_x.py`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b'4\\n' # This is a bytes object\n",
      "4\n"
     ]
    }
   ],
   "source": [
    "# The subprocess module\n",
    "import subprocess\n",
    "\n",
    "# Defines command 1, which would direct stdout into a pipe rather than printing to terminal (can now be accessed programatically)\n",
    "process1 = subprocess.Popen(['echo', '-n', '1234'], stdout=subprocess.PIPE) # process1 is an instance of .Popen class.\n",
    "\n",
    "# Defines second command, defines input from pipe, and output back to pipe\n",
    "process2 = subprocess.Popen(['wc', '-m'], stdin=subprocess.PIPE, stdout=subprocess.PIPE)\n",
    "\n",
    "# Actually runs the commands as a pipeline\n",
    "output1, _ = process1.communicate() # .communicate is a method, waits for the process to complete, capturing stdout + stderr in a tuple. stdout as first element, stderr as second. We use the underscore to throw away that value.\n",
    "output2, _ = process2.communicate(input=output1) # Sends captured stdout as input, and waits for completion\n",
    "\n",
    "# Final output\n",
    "out = output2.decode().strip() # Decode converts back from \"bytes object\" to string, strip clears leading/trailing whitespace\n",
    "print(output2, \"# This is a bytes object\")\n",
    "print(out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Python alternative to subprocess\n",
    "\n",
    "- The \"sh\" package is a subprocess replacement\n",
    "\n",
    "- Programs can be called as if they were functions\n",
    "\n",
    "- This can be used for any binary command on your system (i.e., in your Conda env)\n",
    "\n",
    "- Only works on Unix-like operating systems, as it uses underlying system calls (not python reimplementations)\n",
    "\n",
    "- It's in the repo, installed with: `conda install -c conda-forge sh`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "./python-bioinformatics.ipynb\n",
      "./python-machine-learning.ipynb\n",
      "./python-dashboards.ipynb\n",
      "./python-data-manipulation-and-graphics.ipynb\n",
      "./python-fundamentals.ipynb\n",
      "\n",
      "bedtools v2.30.0\n",
      "\n",
      "This is not the file you're looking for\n"
     ]
    }
   ],
   "source": [
    "import sh\n",
    "\n",
    "# Find as a \"sh\" argument\n",
    "print(sh.find(\".\",\"-name\",\"*ipynb\"))\n",
    "\n",
    "# Bedtools, an installed package in the Conda env\n",
    "print(sh.bedtools(\"-version\"))\n",
    "\n",
    "# Can build in python exception handling, for example if a file isn't real\n",
    "try:\n",
    "    sh.ls(\"./non-existant_file\")\n",
    "except sh.ErrorReturnCode_2:\n",
    "    print (\"This is not the file you're looking for\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bioinformatics packages built for Python\n",
    "\n",
    "### Some popular bioinformatics packages have reimplementions for direct python integration (e.g., Pybedtools)\n",
    "\n",
    "- To use Linux programmes, we need to install them and call via Sh or Subprocess\n",
    "- If there is a Python reimplementation available, these can be used directly; we will use Pybedtools as an example\n",
    "- Pybedtools is in the conda environment already, installed using: `conda install -c bioconda pybedtools`\n",
    "- Again, if already familiar with Bash, learning the python reimplementations can feel clunky"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A\n",
      "chr1\t1\t100\tfeature1\t0\t+\n",
      " chr1\t100\t200\tfeature2\t0\t+\n",
      " chr1\t150\t500\tfeature3\t0\t-\n",
      " chr1\t900\t950\tfeature4\t0\t+\n",
      " \n",
      "B\n",
      "chr1\t155\t200\tfeature5\t0\t-\n",
      " chr1\t800\t901\tfeature6\t0\t+\n",
      " \n",
      "A and B\n",
      "chr1\t155\t200\tfeature2\t0\t+\n",
      " chr1\t155\t200\tfeature3\t0\t-\n",
      " chr1\t900\t901\tfeature4\t0\t+\n",
      " "
     ]
    }
   ],
   "source": [
    "import pybedtools\n",
    "\n",
    "a = pybedtools.example_bedtool('a.bed')\n",
    "b = pybedtools.example_bedtool('b.bed')\n",
    "\n",
    "print(\"A\")\n",
    "a.head()\n",
    "print()\n",
    "print(\"B\")\n",
    "b.head()\n",
    "\n",
    "# Intersect a with b\n",
    "a_and_b = a.intersect(b)\n",
    "print()\n",
    "print (\"A and B\")\n",
    "a_and_b.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Biopython\n",
    "\n",
    "- [Biopython](https://biopython.org) is a collection of Python modules and functions for bioinformatics\n",
    "- Installed here via Conda (the command used was: `conda install -c conda-forge biopython`)\n",
    "- I'm not going to do a deep dive on this, just a small example\n",
    "- I am not yet sure of a clear advantage to using Biopython versus a Linux based install of certain original packages\n",
    "- One concern would be package maintenance & upgrade - these are likely less commonly updated than the original software\n",
    "- Additionally, many packages likely never make it to Biopython, which I suspect has a relatively small community\n",
    "- Where it probably does have advantages is in advanced programmatic manipulation or parsing of files and formats\n",
    "- Another possibility is for NCBI database accessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATCG Opposite\n",
      "GCTA Reverse\n"
     ]
    }
   ],
   "source": [
    "# Many modules, e.g. AlignIO is for multiple sequence alignment\n",
    "\n",
    "# We'll do a simple reverse complement with Seq\n",
    "from Bio.Seq import Seq\n",
    "dna_seq = Seq(\"TAGC\") # Returns a seq object\n",
    "print(dna_seq.complement(), \"Opposite\") # Returns the opposite strand\n",
    "print(dna_seq.reverse_complement(), \"Reverse\") # Returns the reverse complement"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python-data-science",
   "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.7.15"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}