{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Homework Problem \\#1 - DNA sequencing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Hand-in format:** IPython Notebook or python program. Submit via email.\n", "\n", "As a reminder: please make sure your code is clean, documentated, and understandable. Make sure it runs without errors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Background" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[DNA](http://en.wikipedia.org/wiki/DNA) is a molecule that encodes genetic instructions for living organisms. DNA typically takes the form of double-stranded helix: each strand corresponds to a sequence of *nucleotides*.\n", "\n", "Each nucleotide consists of a nucleobase (guanine, adenine, thymine, or cytosine) attached to sugars, which are in turn separated from each other by phosphate groups:\n", "\n", "![DNA](http://upload.wikimedia.org/wikipedia/commons/thumb/4/4c/DNA_Structure%2BKey%2BLabelled.pn_NoBB.png/492px-DNA_Structure%2BKey%2BLabelled.pn_NoBB.png)\n", "\n", "The nucleobases are commonly referred to with the letters:\n", "* `G` (guanine)\n", "* `A` (adenine)\n", "* `T` (thymine) and \n", "* `C` (cytosine).\n", "\n", "The nucleobases form pairs (or \"complements\"): `G` pairs with `C`, and `T` pairs with `A`. \n", "\n", "The information content of a piece of DNA is therefore commonly represented as a string, giving a sequence of nucleobases, such as `GATTACACCTCATTATAAA`.\n", "\n", "In this problem set, we will be working with fake genetic data, since real-life data often has added complications, but the functions and techniques uses here are the same or very similar to techniques used on real data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The **reverse complement** of DNA is found by reversing the DNA sequence, then replacing each base by its complement. For example, the reverse complement of `ATGCGGC` is `GCCGCAT`.\n", "\n", "Write a Python function `reverse_complement()` that takes a DNA sequence as a string, and returns the reverse complement. Test your function, then find and print the reverse complement of the following sequence:\n", "\n", " ATGCGCGGATCGTACCTAATCGATGGCATTAGCCGAGCCCGATTACGC\n", "\n", "Note that this function is not needed for the remaining questions below." ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CGCATTAGCCCGAGCCGATTACGGTAGCTAATCCATGCTAGGCGCGTA\n" ] } ], "source": [ "# your solution here\n", "#DNA = \"ATGCGCGGATCGTACCTAATCGATGGCATTAGCCGAGCCCGATTACGC\"[::-1]\n", "#print(DNA)" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['C', 'G', 'C', 'A', 'T', 'T', 'A', 'G', 'C', 'C', 'C', 'G', 'A', 'G', 'C', 'C', 'G', 'A', 'T', 'T', 'A', 'C', 'G', 'G', 'T', 'A', 'G', 'C', 'T', 'A', 'A', 'T', 'C', 'C', 'A', 'T', 'G', 'C', 'T', 'A', 'G', 'G', 'C', 'G', 'C', 'G', 'T', 'A']\n" ] } ], "source": [ "#DNA_1 = \"CGCATTAGCCCGAGCCGATTACGGTAGCTAATCCATGCTAGGCGCGTA\"\n", "#DNA_1 = DNA_1.replace('C', 'G')\n", "#DNA_1 = DNA_1.replace('T', 'A')\n", "#print(DNA1)" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CGCATTAGCCCGAGCCGATTACGGTAGCTAATCCATGCTAGGCGCGTA\n", "GCGTAATCGGGCTCGGCTAATGCCATCGATTAGGTACGATCCGCGCAT\n" ] } ], "source": [ "DNA='ATGCGCGGATCGTACCTAATCGATGGCATTAGCCGAGCCCGATTACGC'\n", "#def reverse_sequence(sequence):\n", " #reverse = sequence[::-1]\n", " #return reverse\n", "\n", "#reverse=reverse_sequence(DNA)\n", "#print(reverse)\n", "\n", "#def complement_sequence(sequence):\n", " #base_complements={'A':'T','T':'A','G':'C','C':'G'}\n", " #complement_of_sequence=[base_complements[base] for base in sequence]\n", " #complement_of_sequence=''.join(complement_of_sequence)\n", " #return complement_of_sequence\n", "\n", "#reverse_complement=complement_sequence(reverse)\n", "#print(reverse_complement)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Ribonucleic acid (RNA)](http://en.wikipedia.org/wiki/RNA) is a family of large biological molecules that are *transcribed* from DNA by the \"RNA polymerase\" enzyme. \n", "\n", "[Messenger RNA](http://en.wikipedia.org/wiki/Messenger_RNA) molecules (or mRNA) are a subset of RNA molecules that are used to pass information from DNA to ribosomes, which then translates the mRNA to protein sequences.\n", "\n", "mRNA consists of a single strand of nucleotides that are identical to the ones found in DNA, with the exception of uracil (`U`), which replaces thymine (`T`).\n", "\n", "Write a Python function `dna_to_mrna()` that takes a DNA sequence and returns the corresponding mRNA sequence. For example, the DNA sequence `ATCGCGAT` should produce the mRNA sequence `AUCGCGAU`." ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "None\n" ] } ], "source": [ "# your solution here\n", "DNA_1='ATCGCGAT'\n", "def complement_sequence_1(sequence):\n", " base_complements_1={'T':'U'}\n", " complement=[base_complements_1[base] for base in sequence]\n", " complement=''.join(complement)\n", " return complement\n", "\n", "dna_to_mrna = complement(DNA_1)\n", "print(dna_to_mrna)" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "dna='ATCGCGAT'\n", "def DNA_strand(dna):\n", " my_dict = {'T':'U'}\n", " my_str = ''\n", " for i in dna:\n", " for key, value in my_dict.items():\n", " if i == key:\n", " my_str +=value\n", " elif i == value:\n", " my_str+=key\n", " return my_str\n", "print(DNA_strand(dna))\n" ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "scrolled": true }, "outputs": [ { "ename": "KeyError", "evalue": "'A'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", "Input \u001b[0;32mIn [62]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 3\u001b[0m sequence_complement\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m base \u001b[38;5;129;01min\u001b[39;00m dna:\n\u001b[0;32m----> 5\u001b[0m base_complement\u001b[38;5;241m=\u001b[39m\u001b[43mbase_complement\u001b[49m\u001b[43m[\u001b[49m\u001b[43mbase\u001b[49m\u001b[43m]\u001b[49m\n\u001b[1;32m 6\u001b[0m sequence_complement\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39mbase_complement\n\u001b[1;32m 7\u001b[0m \u001b[38;5;28mprint\u001b[39m(sequence_complement)\n", "\u001b[0;31mKeyError\u001b[0m: 'A'" ] } ], "source": [ "dna='ATCGCGAT'\n", "base_complement={'T':'U'}\n", "sequence_complement=''\n", "for base in dna:\n", " base_complement=base_complement[base]\n", " sequence_complement+=base_complement\n", "print(sequence_complement)" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "ename": "KeyError", "evalue": "'A'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", "Input \u001b[0;32mIn [63]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 4\u001b[0m dna_to_mrna[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mT\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m nucleotide \u001b[38;5;129;01min\u001b[39;00m dna:\n\u001b[0;32m----> 6\u001b[0m rna_seq\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m\u001b[43mdna_to_mrna\u001b[49m\u001b[43m[\u001b[49m\u001b[43mnucleotide\u001b[49m\u001b[43m]\u001b[49m\n\u001b[1;32m 7\u001b[0m \u001b[38;5;28mprint\u001b[39m(rna_seq)\n", "\u001b[0;31mKeyError\u001b[0m: 'A'" ] } ], "source": [ "dna='ATCGCGAT'\n", "rna_seq = str()\n", "dna_to_mrna={'T':'U'}\n", "dna_to_mrna['T']\n", "for nucleotide in dna:\n", " rna_seq+=dna_to_mrna[nucleotide]\n", "print(rna_seq)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dna='ATCGCGAT'\n", "dna_to_mrna=dna.maketrans('T','U')\n", "print(dna.translate(dna_to_mrna))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "i,j='T','U'\n", "dna_mrna=str.maketrans(i,j)\n", "my_str='ATCGCGAT'\n", "print(my_str.translate(dna_mrna))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When the mRNA is translated to a protein sequence, each set of three nucleotides, called a **codon**, is translated into a single amino acid. For example, the codon ``UUC`` translates to the amino acid *Phenylalanine*. Each amino acid can be represented by a single letter - for example Phenylalanine is represented by the letter ``F``.\n", "\n", "> A protein, which is formed from a sequence of amino acids, can therefore be written as a sequence of letters in the same way as DNA or mRNA, but using more of the letters of the alphabet since there are more than four amino acids.\n", "\n", "The [data/p1_codons.txt](data/p1_codons.txt) file contains a list of codon-amino acid pairs. There are two columns:\n", "* first column: codon (represented by three letters).\n", "* second column: corresponding amino acid (represented by a single letter).\n", "\n", "Certain codons do not correspond to an amino acid, but instead indicate that the amino acid sequence is finished. These are indicated by `Stop`.\n", "\n", "Write a function `mrna_to_protein()` that takes an mRNA sequence (as a string) and returns the sequence of amino acids (as a string), stopping the first time a `Stop` codon is encountered.\n", "\n", "> Make sure that the codon map file is only read once when running the script (and not every time you want to translate a codon).\n", "\n", "Then, write a function `dna_to_protein()` that takes a DNA sequence (as a string) and returns the sequence of amino acids (as a string), making use of the functions that you wrote previously.\n", "\n", "Print out the amino acid sequence for the following DNA sequence:\n", "\n", " AATCTCTACGGAAGTAGGTCAGTACTGATCGATCAGTCGATCGGGCGGCGATTTCGATCTGATTGTACGGCGGGCTAG\n", " " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# your solution here\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the previous questions, we have been specifying the DNA sequence by hand, but DNA sequences are usually long and are stored in files. A common file format is the [FASTA format](http://en.wikipedia.org/wiki/FASTA_format) which is a text-file of the form:\n", "\n", " >label1\n", " ACTGTATCGATGCTAGCTACGTAGCTAGCTAGCTAGCTGACGTA\n", " ACGATGTGCGAGGGTCATGGGACGCGAGCGAGTCTAGCACGATC\n", " >label2\n", " ACTGGGCTTGACTACGGCGGTATCTGACGGGCGAGCTGTACGAG\n", " ACGGACTAGGGCGCGGCGGGGCGGATTTTCGAGTCGAGCGTTAT\n", " \n", "The first line starts with a ``>`` which is immediately followed by a label (an arbitrary string, e.g. the name of the gene). The sequence then starts on the second line, and may continue on several lines. It is common to limit the length of each line to 80, but this may vary from file to file. The sequence stops once either the file ends, or a line starts with ``>``, which indicates that a new sequence is being given. There may be any number of sequences in a file.\n", "\n", "Write a function `read_fasta()`, that takes the name of a file (as a string) and returns a Python dictionary containing all the sequences from the file, with the keys in the dictionary corresponding to the label. \n", "\n", "Use this function and the functions you have written above to read in the [data/p1_fasta_q4.txt](data/p1_fasta_q4.txt) file and print out, for each sequence, the label, followed by the **amino acid** sequence (not the DNA sequence!)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# your solution here\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given several sequences with the same length, some of which may include **point mutations** (i.e. individual nucleotides are changed), our goal is to try and find the most likely original sequence. For example, if we have the following sequences:\n", "\n", " sequence 1: A C T C T\n", " sequence 2: A C T C G\n", " sequence 3: G C C C T\n", " sequence 2: A C T C T\n", " sequence 4: A T G C T\n", " \n", "we can go through each position (of 5, in the both example) and find the most common nucleotide. In this case:\n", " \n", " A: 4 0 0 0 0\n", " C: 0 4 1 5 0\n", " G: 1 0 1 0 1\n", " T: 0 1 3 0 4\n", "\n", "where the five columns represent the five positions in the sequence, and the value in each gives the count of that nucleotide, summing across all the sequences.\n", "\n", "In this example, the most common \"first base\" is A, the most common second base is C, and so on. The most common sequence is then `ACTCT`. This is called the [consensus sequence](http://en.wikipedia.org/wiki/Consensus_sequence).\n", "\n", "Write a function `consensus_sequence()` that takes a dictionary of sequences returned by `read_fasta()` and computes the consensus sequence. \n", "\n", "Calculate the consensus sequence from the twenty sequences in [data/p1_fasta_q5.txt](data/p1_fasta_q5.txt)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# your solution here\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 6" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One common task is to identify the longest common sub-sequence between two sequences. For example, in the sequences `ACTGCT` and `TGCCCT`, the longest common sub-sequence is `TGC` (AC**TGC**T and **TGC**CCT). Note that these do not have to be at the same positon in each sequence.\n", "\n", "Write a function `longest_common_sequence()` that takes a dictionary of sequences and returns the longest common sub-sequence found in all the sequences.\n", "\n", "Calculate the longest common sub-sequence from the twenty sequences in [data/p1_fasta_q6.txt](data/p1_fasta_q6.txt)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# your solution here\n" ] } ], "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.9.12" } }, "nbformat": 4, "nbformat_minor": 4 }