{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Pathogenicity islands\n",
        "\n",
        "This exercise was inspired by [Libeskind-Hadas and Bush, *Computing for Biologists*, Cambridge University Press, 2014](https://www.cs.hmc.edu/CFB).\n",
        "\n",
        "[Data set download](https://raw.githubusercontent.com/justinbois/bootcamp/refs/heads/main/data/salmonella_spi1_region.fna)\n",
        "\n",
        "<hr>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For this and the next problem, we will work with real data from the *Salmonella enterica* genome.  The section of the genome we will work with is downloadable [here](https://raw.githubusercontent.com/justinbois/bootcamp/refs/heads/main/data/salmonella_spi1_region.fna). I cut it out of the [full genome](http://www.ncbi.nlm.nih.gov/nucleotide/821161554). It contains *Salmonella* pathogenicity island I (SPI1), which contains genes for surface receptors for host-pathogen interactions.\n",
        "\n",
        "Pathogenicity islands are often marked by different GC content than the rest of the genome. We will try to locate the pathogenicity island(s) in our section of the *Salmonella* genome by computing GC content.\n",
        "\n",
        "**a)** Write a function that divides a sequence into blocks and computes the GC content for each block, returning a tuple. The function signature should look like\n",
        "\n",
        "    gc_blocks(seq, block_size)\n",
        "    \n",
        "To be clear, if `seq = 'ATGACTACGT'` and `block_size = 4`, the blocks to be considered are\n",
        "\n",
        "    ATGA\n",
        "    CTAC\n",
        "    \n",
        "and the function should return `(0.25, 0.5)`. Note that the blocks are non-overlapping and that we don't bother with the fact that end of the sequence that does not fit completely in a block."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**b)** Write a function that takes as input a sequence, block size, and a threshold GC content, and returns the original sequence where every base in a block with GC content above threshold is capitalized and every base below the threshold is lowercase. You would call the function like this:\n",
        "\n",
        "    mapped_seq = gc_map(seq, block_size, gc_thresh)\n",
        "\n",
        "For example, \n",
        "\n",
        "    gc_map('ATGACTACGT', 4, 0.4)\n",
        "\n",
        "returns `'atgaCTAC'`. Note that bases not included in GC blocks (in this case the `GT` at the end of the sequence) are not included in the output sequence."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**c)** Use the `gc_map()` function to generate a GC content map for the *Salmonella* sequence with `block_size = 1000` and `gc_thresh = 0.45`. Where do you think the pathogenicity island is?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**d)** Write the GC-mapped sequence (with upper and lower characters) to a new FASTA file. Use the same description line (which began with a `>` in the original FASTA file), and have line breaks every 60 characters in the sequence."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<br />"
      ]
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "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.13.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}