homework04_test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 4: Practical analysis with BioPython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the homework, you are going to extend the code from real biological analysis of our FASTQ files in lectures 8 and 9.\n",
    "\n",
    "As described in the Jupyter notebook for that lecture, the FASTQ reads can originate from **either** HA or NA, and that will be distinguished by the most 3' end of the read.\n",
    "But in our example, we did not distinguish among reads matching to HA and NA, as we didn't even look far enough into the read to tell the identity.\n",
    "\n",
    "For the homework, your goal is to write code that extends the material from lectures 8 and 9 to also distinguish between HA and NA.\n",
    "\n",
    "Please include code to address each of the following questions. Please include code comments to explain what your code is attempting to accomplish. Don't forget to include references to the sources you used to obtain your answer, including your classmates (if you are working in groups).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A real biological analysis: parsing barcodes\n",
    "The reads that we just read as `seqreads_str` come from a real sequencing run of influenza virus HA and NA genes.\n",
    "\n",
    "The sequences are as follows:\n",
    "\n",
    "    5'-[end of HA]-AGGCGGCCGC-[16 X N]-3'\n",
    "    \n",
    "or \n",
    "\n",
    "    5'-[end of NA]-AGGCGGCCGC-[16 X N]-3'\n",
    "    \n",
    "The end of NA is:\n",
    "\n",
    "    ...CACGATAGATAAATAATAGTGCACCAT\n",
    "    \n",
    "The end of HA is:\n",
    "\n",
    "    ...CCGGATTTGCATATAATGATGCACCAT\n",
    "    \n",
    "The sequencing run reads from the reverse end of the molecules, so the first thing in the sequencing reads is the barcode followed by the constant sequence and the end of HA or NA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. How many reads map to HA, and how many reads map to NA?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "#taken from lecture09\n",
    "\n",
    "import re #imports regular expression package onto jpynb\n",
    "\n",
    "import Bio.SeqIO #imports Bio python package onto jpynb\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 10000 sequencing reads in barcpdes_R1 file\n",
      "Found 10000 sequences in seqreads\n"
     ]
    }
   ],
   "source": [
    "#taken from lecture09\n",
    "\n",
    "seqreads = list(Bio.SeqIO.parse('barcodes_R1.fastq', format='fastq')) #parses through barcodes_R1.fastq file and stores every entry in a list called seqreads\n",
    "seqreads_str = [str(s.seq) for s in seqreads] #iterates through all entries from seqreads and stores every sequence into a list of strings called seqreads_str\n",
    "\n",
    "print(f'Found {len(seqreads)} sequencing reads in barcpdes_R1 file') #prints total number of entries in seqreads\n",
    "print (f'Found {len(seqreads_str)} sequences in seqreads') #prints total number of sequences in seqreads_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-52-51e189454aef>, line 120)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-52-51e189454aef>\"\u001b[0;36m, line \u001b[0;32m120\u001b[0m\n\u001b[0;31m    invalid_matcher = re.compile(f'^(?P<barcode>[ACGT]{{{bclen}}})' + reverse_complement ((?!^ABC$(upstream)))\u001b[0m\n\u001b[0m                                                                                           ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "#defining reverse_complement function to get the reverse complement, append to the revcomplement list and join the nucleotides\n",
    "\n",
    "def reverse_complement(seq): #\"seq\" is what the user will input, which will be fed through the reverse_complement function\n",
    "    \"\"\"Get reverse complement of a DNA sequence.\n",
    "    \n",
    "    Parameters\n",
    "    -----------\n",
    "    seq : str\n",
    "        Uppercase DNA sequence.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    str\n",
    "        Reverse complement of the sequence in upper case.\n",
    "        \n",
    "    Example\n",
    "    --------\n",
    "    >>> reverse_complement('ATGCAC')\n",
    "    'GTGCAT'\n",
    "    \n",
    "    \"\"\"\n",
    "    rc_dict = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G', 'N': 'N'} #rc_dict is the dictionary that will be iterated through to get the reverse_complement of each sequence\n",
    "    revcomplement = [] #revcomplement = [] is the empty list that each reverse complement sequence will be added to\n",
    "    for nt in reversed(seq.upper()): #seq.upper makes all the nucleotides in each sequence into an uppercase letter\n",
    "        revcomplement.append(rc_dict[nt]) #adds each reverse complemented sequence to the revcomplement list\n",
    "    return ''.join(revcomplement) #joins each nucleotide in each specific sequence together\n",
    "\n",
    "\n",
    "#####################################################################################################################\n",
    "\n",
    "\n",
    "\n",
    "#defining HA_barcode function to find HA specific pattern (through regex) and search through seqread\n",
    "def HA_barcode(seqread, bclen=16, upstream='AGGCGGCCGC', HA='CCGGATTTGCATATAATGATGCACCAT'):\n",
    "    \"\"\"Identify barcode with known upstream sequence.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    seqread : str\n",
    "        Nucleotide sequence matching UPSTREAM-BARCODE read in reverse orientation.\n",
    "    bclen : int\n",
    "        Length of barcode\n",
    "    upstream: str\n",
    "        Sequence upstream of the barcode.\n",
    "    HA : str\n",
    "        Distinguishing HA genes within seqreads_string\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    str or None\n",
    "        Sequence of the barcode in the forward orientation, or `None` if no match to expected barcoded sequence.\n",
    "        \n",
    "\n",
    "        \n",
    "    \"\"\"\n",
    "    HA_matcher = re.compile(f'^(?P<barcode>[ACGT]{{{bclen}}})' + reverse_complement(upstream) + reverse_complement(HA)) #finding HA pattern through regex\n",
    "    HA_search = HA_matcher.search(seqread) #searching seqread using the regex pattern stated above\n",
    "    if not HA_search: #if it's not HA pattern, skip it\n",
    "        return None\n",
    "    else:\n",
    "        return reverse_complement(seqread[: bclen]) #else, print the barcode\n",
    "\n",
    "\n",
    "######################################################################################################################\n",
    "\n",
    "\n",
    "# defining NA_barcode function to find NA specific pattern (through regex) and search through seqread\n",
    "def NA_barcode(seqread, bclen=16, upstream='AGGCGGCCGC', NA='CACGATAGATAAATAATAGTGCACCAT'):\n",
    "    \"\"\"Identify barcode with known upstream sequence.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    seqread : str\n",
    "        Nucleotide sequence matching UPSTREAM-BARCODE read in reverse orientation.\n",
    "    bclen : int\n",
    "        Length of barcode\n",
    "    upstream: str\n",
    "        Sequence upstream of the barcode.\n",
    "    HA : str\n",
    "        Distinguishing HA genes within seqreads_string\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    str or None\n",
    "        Sequence of the barcode in the forward orientation, or `None` if no match to expected barcoded sequence.\n",
    "        \n",
    "        \n",
    "    \"\"\"\n",
    "    NA_matcher = re.compile(f'^(?P<barcode>[ACGT]{{{bclen}}})' + reverse_complement(upstream) + reverse_complement(NA))\n",
    "    NA_search = NA_matcher.search(seqread)\n",
    "    if not NA_search:\n",
    "        return None\n",
    "    else:\n",
    "        return reverse_complement(seqread[: bclen])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the sum of HA is 5245\n",
      "the sum of NA is 3907\n"
     ]
    }
   ],
   "source": [
    "#got help from Tim Yu\n",
    "HA_dict = {}\n",
    "NA_dict = {}\n",
    "\n",
    "\n",
    "HA_count = 0\n",
    "NA_count = 0\n",
    "\n",
    "for seq in seqreads_str:\n",
    "    HA = HA_barcode(seq)\n",
    "    if HA is not None:\n",
    "        if HA not in HA_dict:\n",
    "            HA_dict[HA] = 1 #add the barcode as the key and give it a value of one\n",
    "        else:\n",
    "            HA_dict[HA] += 1 #updates the value of each key\n",
    "    NA = NA_barcode (seq)\n",
    "    if NA is not None:\n",
    "        if NA not in NA_dict:\n",
    "            NA_dict[NA] = 1 #add the barcode as the key and give it a value of one\n",
    "        else:\n",
    "            NA_dict[NA] += 1 #updates the value of each key\n",
    "\n",
    "            \n",
    "\n",
    "HA_dict\n",
    "print (f'the sum of HA is {sum(HA_dict.values())}')\n",
    "NA_dict\n",
    "print (f'the sum of NA is {sum(NA_dict.values())}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. What is the HA barcode with the most counts (and how many counts)? Also answer the same question for NA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum HA value is 155 and the barcode is CCCGACCCGACATTAA\n",
      "The maximum NA value is 152 and the barcode is ACCAGTTCTCCCCGGG\n"
     ]
    }
   ],
   "source": [
    "all_HA_values = HA_dict.values()\n",
    "max_HA_value = max(all_HA_values)\n",
    "\n",
    "all_HA_keys = HA_dict.keys()\n",
    "max_HA_keys = max(HA_dict, key=lambda key: HA_dict[key])\n",
    "\n",
    "print(f'The maximum HA value is {max_HA_value} and the barcode is {max_HA_keys}')\n",
    "\n",
    "\n",
    "all_NA_values = NA_dict.values()\n",
    "max_NA_value = max(all_NA_values)\n",
    "\n",
    "all_NA_keys = NA_dict.keys()\n",
    "max_NA_keys = max(NA_dict, key=lambda key: NA_dict[key])\n",
    "\n",
    "print (f'The maximum NA value is {max_NA_value} and the barcode is {max_NA_keys}')\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "40d3a090f54c6569ab1632332b64b2c03c39dcf918b08424e98f38b5ae0af88f"
  },
  "kernelspec": {
   "display_name": "Python 3.8.8 64-bit ('base': conda)",
   "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.8.8"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {
    "height": "623.8px",
    "left": "0px",
    "right": "829.6px",
    "top": "33.2px",
    "width": "322.4px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}