day_3.ipynb

{
 "metadata": {
  "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.5"
  },
  "orig_nbformat": 4,
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.8.5 64-bit ('myenv': conda)"
  },
  "interpreter": {
   "hash": "f5c9dfabb21bd2a91b63810df81acbdc6b5e617e45414f0ef050ca96090c868f"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "source": [
    "# Welcome to Day 3! \n",
    "\n",
    "## More sequence modification and data extraction\n",
    "\n",
    "### Section 1: Filtering for sequences matching certain conditions\n",
    "\n",
    "### Section 2: Storing filtered sequences to files\n",
    "\n",
    "### Section 3: Reformating and storing sequences\n",
    "\n",
    "---\n",
    "\n",
    "## Session summary\n",
    "\n",
    "\n",
    "Today things start to get a little serious. Building on concepts we learned in Day 1 and Day 2, we can do things like filter for sequences of a certain length or containing a motif of interest, remove undesired characters, or rename headers. We then go over how to store all of our modifications and reformats.\n",
    "\n",
    "---\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "---\n",
    "## For Google colab users only\n",
    "\n",
    "Run the following commands"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install Biopython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/to_reformat.fasta\n",
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/mixed_args.fasta\n",
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/genes_from_ocean.fasta\n",
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/genes_from_desert.fasta"
   ]
  },
  {
   "source": [
    "---\n",
    "\n",
    "# Section 1: Filtering for sequences matching certain conditions\n",
    "\n",
    "Let's say a collaborator gives you fasta file containing 13 sequences of two types of antibiotic resistance genes: one gene type, ermA, confers erythromycin resistance and the other gene type, blaC, imparts beta lactam antibiotic resistance. \n",
    "\n",
    "You know that there are 8 erm sequences and 5 bla sequences.\n",
    "\n",
    "How would you separate the two sets of antbiotic resistance gene classes into two files? "
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import our two Biopython modules that we will be using\n",
    "from Bio import SeqIO\n",
    "from Bio.Seq import Seq"
   ]
  },
  {
   "source": [
    "One way to distinguish between the different gene types is their nucleotide lengths. A useful way to examine their lengths is to parse a fasta file and print the length (using the `len()` function), and the gene name in a single line. "
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    print(len(record.seq), ' ',record.id)\n",
    "    print('\\n')"
   ]
  },
  {
   "source": [
    "\n",
    "\n",
    "From the above output, we can see that ermA genes are 729 nucleotdies long while Bla genes are at least 921 nucleotides long. \n",
    "\n",
    "However, there are some gene names that are not especially clear. Are they erm or bla? \n",
    "\n",
    "We can access additional information stored in each gene's header using record.description \n",
    "\n",
    "---\n",
    "### Exercise 1\n",
    "\n",
    "Using the above code as an example, use `SeqIO.parse()` to parse `mixed_args.fasta` and `print()` the `record.description` of each record followed by the length of the `record.seq`"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in ___.parse('mixed_args.fasta', '___'):\n",
    "    print(record.___, ' ', len(___))\n",
    "    print('\\n')"
   ]
  },
  {
   "source": [
    "---\n",
    "\n",
    "The output now clearly shows that there is a nucleotide length difference between ermA and blaC that is sufficiently different to separate the two gene types.\n",
    "\n",
    "Now we will filter for sequences that match the ermA length and print them."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    if len(record.seq) == 729:\n",
    "        print(record.description)\n",
    "        print('\\n')"
   ]
  },
  {
   "source": [
    "The blaC genes have lengths that range from 921 to 1143. What if we want to filter for only blaC genes? We can do this in a couple of different ways.\n",
    "\n",
    "The simplest is to search for anything larger than 900\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    if len(record.seq) > 900:\n",
    "        print(record.description)\n",
    "        print('\\n')"
   ]
  },
  {
   "source": [
    "Another is to set length boundaries between 921 and 1143"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    if len(record.seq) >= 921:\n",
    "        if len(record.seq) <= 1143:\n",
    "            print(record.description)\n",
    "            print('\\n')"
   ]
  },
  {
   "source": [
    "---\n",
    "### Exercise 2\n",
    "\n",
    "Show only Bla genes less than 1143 nucleotides.\n",
    "\n",
    "`print()` the `record.description` only if the `record.seq` has a `len()` **less** than `1143` nucleotides."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ___ in SeqIO.___('mixed_args.fasta','___'):\n",
    "    if len(___) < ___:\n",
    "        print(record.___)\n",
    "        print('\\n')"
   ]
  },
  {
   "source": [
    "---\n",
    "\n",
    "Another way to filter for genes of interest is to select only records that contain key words in their headers. Such as all sequences that have the word Beta-lactamase in `record.description`. We can use the base Python function `find()` for filtering.\n",
    "\n",
    "`find()` returns the position that matches the start of a string of interest. If the string is not found, it retuns -1.\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# where starts at position 0 of the string 'where am I?'\n",
    "'where am I?'.find('where')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ? starts at position 10 of the string 'where am I?'\n",
    "'where am I?'.find('?')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the first occurrence of e is at position 2 of the string 'where am I?'\n",
    "'where am I?'.find('e')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# q is not in the string 'where am I?' and so it returns a -1 position\n",
    "'where am I?'.find('q')"
   ]
  },
  {
   "source": [
    "We will filter for sequences that have the word Beta-lactamase in `record.description`. "
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    # if the record.description does have Beta-lactamase in it, then the character position will be greater than -1.\n",
    "    if record.description.find('Beta-lactamase') > -1:\n",
    "        print(record.description)"
   ]
  },
  {
   "source": [
    "---\n",
    "\n",
    "### Exercise 3a\n",
    "\n",
    "Filter for sequences with the term `'Erythromycin resistance'` in their `record.description` using `find()`\n",
    "\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ___ in SeqIO.___('mixed_args.fasta','___'):\n",
    "    if record.description.___('___') > __:\n",
    "        print(record.description)\n",
    "        print('\\n')"
   ]
  },
  {
   "source": [
    "### Exercise 3b\n",
    "\n",
    "Filter for sequences that do **not** contain `'Erythromycin resistance'` in their `record.description` using `find()`\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ___ in SeqIO.___('mixed_args.fasta','___'):\n",
    "    if record.description.___('___') ___ -1:\n",
    "        print(record.description)\n",
    "        print('\\n')"
   ]
  },
  {
   "source": [
    "___\n",
    "\n",
    "From exercise 3, we can tell that sometimes gene names and descriptions are not always useful for filtering. \n",
    "\n",
    "We can also try filtering by searching for sequence characters (nucleotide or amino acid) of interest.\n",
    "\n",
    "For example, we can look for all sequences that contain the DNA motif `'GCTGTAAATC'`\n",
    "\n",
    "\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    if record.seq.find('GCTGTAAATC') > -1:\n",
    "        print(record.description)"
   ]
  },
  {
   "source": [
    "Now let's try filtering for sequences using protein motifs. \n",
    "\n",
    "This is a two step-challenge: first we must translate the DNA sequence and afterwards we need to search for the motif. \n",
    "\n",
    "We will search by filtering for sequences that carry the blaC motif, GAKR"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    translated_record = record.translate()\n",
    "    if translated_record.seq.find('GAKR') > -1:\n",
    "        print(record.description)"
   ]
  },
  {
   "source": [
    "What if we want to use more than one motif to more clearly differentiate between sequences?\n",
    "\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    translated_record = record.translate()\n",
    "    if translated_record.seq.find('GAKR') > -1:\n",
    "        if translated_record.seq.find('DGTA') > -1:\n",
    "            print(record.description)\n",
    "# compare result to the search with only GAKR motif"
   ]
  },
  {
   "source": [
    "With the double motif filter we end up losing the blaC gene belonging to citroboacter freundii, which seems to be more phylogenetically distant than the other blaC (and 1143 nucleotides long compared to 921)\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "Many protein motifs actually have stretches of unconserved amino acids spaced between conserved amino acids. `find()` will not be sufficient here. \n",
    "\n",
    "Instead we can use [regex](https://docs.python.org/3/howto/regex.html), which is designed to search for character patterns in strings.\n",
    "\n",
    "Regex is a very powerful tool and only its most basic usage will be employed here to demonstrate its utility in identifying motifs\n",
    "\n",
    "Let's search for the protein motif `'T-X-GAKR'` where X is a single occurrence of any amino acid, in the string `'WWWWTLGAKRYYYYY'`"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regex_match_list = re.findall(r\"T.{1}GAKR\", 'WWWWTLGAKRYYYYY')\n",
    "regex_match_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regex_match_list = re.findall(r\"T.{1}GAKR\", 'WWWWTLGAKRYYYYYTVGAKRREHE')\n",
    "regex_match_list"
   ]
  },
  {
   "source": [
    "Now let's let's use the same motif to find only bla genes in `mixed_args.fasta`"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    translated_record = record.translate()\n",
    "\n",
    "    regex_match_list = re.findall(r\"T.{1}GAKR\", str(translated_record.seq))\n",
    "\n",
    "    if len(regex_match_list) > 0:\n",
    "        print(record.description)\n",
    "        print('\\n')"
   ]
  },
  {
   "source": [
    "Success! We correctly printed only and all five bla genes. \n",
    "\n",
    "One noticeable thing we did here was we had to convert translated_record.seq to a string using `str()`. This is because `regex` works on strings and not on `Seq()` objects"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "---\n",
    "\n",
    "### Exercise 4a\n",
    "\n",
    "`translate()` the `record` object and assign it to the variable `translated_record`. Afterwards, filter for sequences that contain the ermA sequence motif `'QRALGL'` using `find()`"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ___ in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "    translated_record = ___.___\n",
    "    if translated_record.___.___('___') > ___:\n",
    "        print(record.description)\n",
    "        print('\\n')"
   ]
  },
  {
   "source": [
    "### Exercise 4b\n",
    "\n",
    "Use `re.findall()` to search for DNA motif `'C.{3}TTGC'`. \n",
    "\n",
    "Make sure to convert the `Seq()` object, `record.seq`, to a string with `str()` prior to using regex's `re.findall()` function.\n",
    "\n",
    "if the motif is detected, `print()` the `record.description`\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ___ in SeqIO.parse('mixed_args.fasta','fasta'):\n",
    "\n",
    "    regex_match_list = ___.___(r\"___\", ___(record.seq))\n",
    "    if len(___) > 0:\n",
    "        print(___.___)\n",
    "        print('\\n')"
   ]
  },
  {
   "source": [
    "### Exercise 4c\n",
    "\n",
    "You wonder if your DNA motif `'C.{3}TTGC'` is specific to only the forward strand of blaC genes.\n",
    "\n",
    "Test for its detection in the forward, `complement()`, reverse `[::-1]`, or `reverse_complement()` oritentation.\n",
    "\n",
    "`print()` the `record.description` each time the DNA motif is detected\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ___ in ___.___('mixed_args.fasta','fasta'):\n",
    "\n",
    "    # search record.seq\n",
    "    ___ = re.findall(___\"C.{3}TTGC\", ___(record.___))\n",
    "    if len(regex_match_list) > 0:\n",
    "        print('found in forward sequence')\n",
    "        print(record.description)\n",
    "        print('\\n')\n",
    "    \n",
    "    # search the reverse of record.seq\n",
    "    ___ = re.findall(___\"C.{3}TTGC\", ___(record.___[::-1]))\n",
    "    if len(regex_match_list) > 0:\n",
    "        print('found in reverse sequence')\n",
    "        print(record.description)\n",
    "        print('\\n')\n",
    "\n",
    "    # search the complement of record.seq\n",
    "    ___ = re.findall(___\"C.{3}TTGC\", ___(record.___()))\n",
    "    if len(regex_match_list) > 0:\n",
    "        print('found in complement sequence')\n",
    "        print(record.description)\n",
    "        print('\\n')\n",
    "\n",
    "    # search the reverse complement of record.seq\n",
    "    ___ = re.findall(r___\"C.{3}TTGC\", ___(record.___()))\n",
    "    if len(regex_match_list) > 0:\n",
    "        print('found in reverse complement sequence')\n",
    "        print(record.description)\n",
    "        print('\\n')\n",
    "\n"
   ]
  },
  {
   "source": [
    "---\n",
    "\n",
    "# Section 2: Storing filtered sequences to files\n",
    "\n",
    "After we parse files and look for features of interest, we typically want to store them in a new location, such as a new fasta file. \n",
    "\n",
    "We will first store the sequences we have filtered for in a list (like we learned to do in Day 2). Afterwards, we will write that list to a file using `SeqIO.write()`, which we also learned to do in Day 2.\n",
    "\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "# empty list to store modified sequence records\n",
    "filtered_records_list = []\n",
    "\n",
    "for record in SeqIO.parse('mixed_args.fasta', 'fasta'):\n",
    "    # translate the sequence\n",
    "    translated_seq = record.seq.translate()\n",
    "    # search for the ermA motif\n",
    "    if translated_seq.find('QRALGL') > -1:\n",
    "        # if the sequence contains the motif, append it to our list\n",
    "        filtered_records_list.append(record)\n",
    "\n",
    "# examine the contents (will look a little messy)\n",
    "print(filtered_records_list)\n"
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": []
  },
  {
   "source": [
    "Now that we have the list ready, we can write it a new fasta file using `SeqIO.write()`\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SeqIO.write(filtered_records_list, 'erma.fasta','fasta')\n",
    "# open and inspect!"
   ]
  },
  {
   "source": [
    "---\n",
    "\n",
    "### Exercise 5a\n",
    "\n",
    "For each sequence, keep only those with with a blaC motif. `append()` each record that passes these conditions to the list `extracted_seqs`. Use `SeqIO.write()` to write them to a new file called `'day3_5a.fasta'`\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# empty list to store selected sequence records\n",
    "extracted_seqs = ___\n",
    "\n",
    "for ___ in ___('mixed_args.fasta', 'fasta'):\n",
    "    # translate the sequence\n",
    "    translated_seq = record.seq.translate()\n",
    "    # search for the blaC motif\n",
    "    if translated_seq.find('GAKR') > -1:\n",
    "        # if the sequence contains the motif, append it to our list\n",
    "        extracted_seqs.___(record)\n",
    "\n",
    "# examine the contents\n",
    "print(extracted_seqs)\n",
    "\n",
    "#write to file\n",
    "SeqIO.___(___, 'day3_5a_blac.fasta','___')"
   ]
  },
  {
   "source": [
    "### Exercise 5b\n",
    "\n",
    "For each sequence with a blaC motif within `record.seq`, select only those with a `record.seq` `len()` shorter than 1143 nucleotides. `append()` each record that passes these conditions to the list `extracted_seqs`. Use `SeqIO.write()` to write them to a new file called `'day3_5b.fasta'`"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# empty list to store selected sequence records\n",
    "extracted_seqs = ___\n",
    "\n",
    "for ___ in ___('mixed_args.fasta', 'fasta'):\n",
    "    # translate the sequence\n",
    "    translated_seq = record.seq.translate()\n",
    "    # search for the blaC motif\n",
    "    if translated_seq.find('GAKR') > -1:\n",
    "        # check if the nucleotide sequence is less than 1143 nucleotides\n",
    "        if ___(___.seq) < 1143:\n",
    "            # append to sequence\n",
    "            extracted_seqs.___(record)\n",
    "\n",
    "# write to file\n",
    "SeqIO.___(___, '___.fasta','fasta')"
   ]
  },
  {
   "source": [
    "___\n",
    "\n",
    "# Section 3: Reformating and storing sequences\n",
    "\n",
    "It is sometimes useful to rename or reformat the headers of sequences.\n",
    "\n",
    "You can use string manipulation to modify record.id or record.description.\n",
    "\n",
    "Let's look at the unsightly headers in a file called `'to_reformat.fasta'`"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('to_reformat.fasta', 'fasta'):\n",
    "    print('record id:',' ',record.id)\n",
    "    print('record description:',' ',record.description)\n",
    "    print('\\n')"
   ]
  },
  {
   "source": [
    "Let's say we don't want our record id to say UgLyGenE_ because we already have a nice gene name next to it. We can use the `replace()` function we learned in Day 1"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('to_reformat.fasta', 'fasta'):\n",
    "    record.id = record.id.replace('UgLyGenE_', '')\n",
    "    \n",
    "    print('record id:',' ',record.id)\n",
    "    print('record description:',' ',record.description)\n",
    "    print('\\n')"
   ]
  },
  {
   "source": [
    "Now we also want to change the description. Maybe all these genes are from *E. coli* and so we want to include that in the description"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('to_reformat.fasta', 'fasta'):\n",
    "    record.id = record.id.replace('UgLyGenE_', '')\n",
    "    record.description = 'this gene is from Escherichia coli'\n",
    "\n",
    "    print('record id:',' ',record.id)\n",
    "    print('record description:',' ',record.description)\n",
    "    print('\\n')"
   ]
  },
  {
   "source": [
    "We also notice one more thing: All of our sequences have an asterik at the end. Let's get rid of it!"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('to_reformat.fasta', 'fasta'):\n",
    "    print(record.seq)\n",
    "    print('\\n')"
   ]
  },
  {
   "source": [
    "Let's get rid of the asterik in the sequences."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for record in SeqIO.parse('to_reformat.fasta', 'fasta'):\n",
    "    # create a string version of record.seq. This is because replace only works on strings\n",
    "    sequence_string = str(record.seq)\n",
    "    sequence_string = sequence_string.replace('*','')\n",
    "    # afterwards, convert the modified string into a Seq object\n",
    "    record.seq = Seq(sequence_string)\n",
    "    print(record.seq)\n",
    "    print('\\n')"
   ]
  },
  {
   "source": [
    "Putting it all together:"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "store_modifed_record = []\n",
    "\n",
    "for record in SeqIO.parse('to_reformat.fasta', 'fasta'):\n",
    "    # change record.id\n",
    "    record.id = record.id.replace('UgLyGenE_', '')\n",
    "    # change record.description\n",
    "    record.description = 'this gene is from Escherichia coli'\n",
    "    # change sequence\n",
    "    sequence_string = str(record.seq)\n",
    "    sequence_string = sequence_string.replace('*','')\n",
    "    record.seq = Seq(sequence_string)\n",
    "    # store \n",
    "    store_modifed_record.append(record)\n",
    "\n",
    "store_modifed_record\n",
    "# write to file\n",
    "SeqIO.write(store_modifed_record,'reformatted_sequences.fasta','fasta')\n",
    "# open and view file\n"
   ]
  },
  {
   "source": [
    "---\n",
    "\n",
    "### Exercise 6a\n",
    "\n",
    "Modify `record.id` so that you `replace()` the string `'UgLyGenE_'` found in `record.id` with the emtpy string `''`.\n",
    "\n",
    "Assign `record.description` to equal an uppercased `record.id` using `upper()`. \n",
    "\n",
    "Afterward, add the string `'_Escherichia coli'` to `record.description."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "for record in SeqIO.parse('to_reformat.fasta', 'fasta'):\n",
    "    ___.id = ___.id.___('___', '')\n",
    "    ___.description = record.___.upper()\n",
    "    ___.description = ___.description + '___'\n",
    "\n",
    "    print(record.description)\n",
    "    print('\\n')"
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": []
  },
  {
   "source": [
    "### Exercise 6b\n",
    "\n",
    "Make an empty list called `store_records`\n",
    "\n",
    "Modify `record.id` so that you `replace()` the string `'UgLyGenE_'` found in `record.id` with the emtpy string `''`.\n",
    "\n",
    "Assign `record.description` to equal an uppercased `record.id` using `upper()`. \n",
    "\n",
    "Afterward, add the string `'_Escherichia coli'` to `record.description.\n",
    "\n",
    "`append()` the modified record to `store_records`\n",
    "\n",
    "Use `SeqIO.write()` to write the contents to the fasta file `'day3_6b.fasta'`\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "___ = []\n",
    "\n",
    "for record in SeqIO.parse('to_reformat.fasta', 'fasta'):\n",
    "    ___.id = ___.id.___('___', '')\n",
    "    ___.description = record.___.upper()\n",
    "    ___.description = ___.description + '___'\n",
    "\n",
    "    store_records.___(___)\n",
    "\n",
    "    print(record.description)\n",
    "    print('\\n')\n",
    "\n",
    "SeqIO.write(store_records, '___.fasta', 'fasta')"
   ]
  },
  {
   "source": [
    "Sometimes the fasta files themselves have important information that can be useful to place in the `record.description` of a sequence when combining multiple fasta files into one.\n",
    "\n",
    "We can use `os.listdir()` to help do this"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "current_directory = os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for file in os.listdir(current_directory):\n",
    "    if file.endswith('.fasta'):\n",
    "        print(file)"
   ]
  },
  {
   "source": [
    "### Exercise 7a\n",
    "\n",
    "loop through each `file` in `current_directory` using `os.listdir()`\n",
    "\n",
    "for each `file`, make sure it `startswith()` `'genes_from_'` and that it `endswith()` `'.fasta'` \n",
    "\n",
    "if the `file` meets these criteria, `replace()` the string `'.fasta'` in the `file` string with an empty string `''`, store it to the variable `file_description`\n",
    "\n",
    "afterward, `print()`the `file_description`"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ___ in ___.___(current_directory):\n",
    "    if file.___('genes_from_'):\n",
    "        if file.endswith('___'):\n",
    "            ___ = file.replace('___', '')\n",
    "            print(file_description)"
   ]
  },
  {
   "source": [
    "### Exercise 7b\n",
    "\n",
    "build on the code in 7a\n",
    "\n",
    "start by making an empty `[]` list called `store_records`\n",
    "\n",
    "afterwards, for each file that passess the above criteria, modify the `record.description` by adding the `file_description` to it.\n",
    "\n",
    "`append()` each modified record to `store_records`\n",
    "\n",
    "afterewards, use `SeqIO.write()` to write `store_records` to a new `'fasta'` file called `'day3_7b.fasta'`"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "___ = ___\n",
    "\n",
    "for ___ in ___.___(current_directory):\n",
    "    if file.___('genes_from_'):\n",
    "        if file.endswith('___'):\n",
    "            for ___ in SeqIO.parse(file, 'fasta'):\n",
    "                ___ = file.replace('___', '')\n",
    "                record.description = record.___ + ' location: '+ ___\n",
    "                ___.append(record)\n",
    "\n",
    "___.___(store_records, 'day3_7b.fasta','___')\n"
   ]
  },
  {
   "source": [
    "---\n",
    "\n",
    "That's it for Day 3! "
   ],
   "cell_type": "markdown",
   "metadata": {}
  }
 ]
}