LCSM-Finding_a_shared_motif.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finding a shared motif, from [rosalind.info](http://rosalind.info/)\n",
    "\n",
    "(text copied from http://rosalind.info/problems/lcsm/)\n",
    "\n",
    "> **Problem**\n",
    ">\n",
    "> A common substring of a collection of strings is a substring of every member of the collection. We say that a common substring is a longest common substring if there does not exist a longer common substring. For example, \"CG\" is a common substring of \"A**CG**TACGT\" and \"AAC**CG**TATA\", but it is not as long as possible; in this case, \"CGTA\" is a longest common substring of \"A**CGTA**CGT\" and \"AAC**CGTA**TA\".\n",
    ">\n",
    "> Note that the longest common substring is not necessarily unique; for a simple example, \"AA\" and \"CC\" are both longest common substrings of \"AACC\" and \"CCAA\".\n",
    ">\n",
    "> **Given**: A collection of k (k≤100) DNA strings of length at most 1 kbp each in FASTA format.\n",
    ">\n",
    "> **Return**: A longest common substring of the collection. (If multiple solutions exist, you may return any single solution.)\n",
    ">\n",
    "> **Sample Dataset**\n",
    ">\n",
    "> ```\n",
    "> >Rosalind_1\n",
    "> GATTACA\n",
    "> >Rosalind_2\n",
    "> TAGACCA\n",
    "> >Rosalind_3\n",
    "> ATACA\n",
    "> ```\n",
    ">\n",
    "> **Sample Output**\n",
    ">\n",
    "> ```\n",
    "> AC\n",
    "> ```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## My interpretation/reasoning\n",
    "\n",
    "1. I am going to look for common substrings of DNA sequences in fasta format.  \n",
    "    From the example, it seems that 'common' means 'shared by _all_ sequences'!\n",
    "\n",
    "2. Of all the common substrings that exist, I only want to find the longest.\n",
    "\n",
    "3. If multiple substrings are the longest, I may return a random one. (Any correct answer is good.)\n",
    "\n",
    "4. There will be no more than 100 DNA sequences, with a length up to 1,000 bp.\n",
    "\n",
    "Now practically, this means the code is to:\n",
    "\n",
    "  - Open and read a fasta text file\n",
    "  - Take the first sequence\n",
    "  - In the second sequence, look for overlaps (>1) with the first and save all of them\n",
    "  - For each following sequence, compare to all overlapping sequences and overwrite when necessary (when the new overlap is shorter)\n",
    "  - In the end, find the longest remaining overlapping sequence. If multiple exist, pick the first.\n",
    "  \n",
    "The most difficult part will be in finding overaps/common substrings. \n",
    "It looks like the `regex` library has nice functionality for that.\n",
    "(Found thanks to David C at https://stackoverflow.com/a/18966891)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "import regex\n",
    "\n",
    "def read_fasta_file(fasta_file):\n",
    "    \"\"\"\n",
    "Takes a fasta file as input to return a dictionary with sequence IDs as keys and sequences as values.\n",
    "    \"\"\"\n",
    "    sequence_dict = {}\n",
    "\n",
    "    for record in SeqIO.parse(fasta_file, \"fasta\"):\n",
    "        sequence_dict[record.id] = record.seq\n",
    "        \n",
    "    return(sequence_dict)\n",
    "\n",
    "def find_initial_overlaps(seq1, seq2):\n",
    "    \"\"\"\n",
    "List all overlaps between the first two sequences.\n",
    "    \"\"\"\n",
    "    substring_lengths = range(2, len(seq1))[::-1]\n",
    "    #List all possible substring lengths, \n",
    "    # from the length of sequence 1 minus 1 down to 2.\n",
    "    \n",
    "    common_substrings = []\n",
    "    #Save all common substrings in a list\n",
    "\n",
    "    for length in substring_lengths:\n",
    "        #For each possible length of substrings...\n",
    "        substring_list = regex.findall(r\"\\w{%i}\" % length, seq1, overlapped=True)\n",
    "        #... list all possible substrings.\n",
    "        \n",
    "        for substring in substring_list:\n",
    "            #For each of the substrings\n",
    "            if substring in seq2 and substring not in common_substrings:\n",
    "                #if it is in sequence 2 and not yet in the list\n",
    "                common_substrings.append(substring)\n",
    "                #add it to the list.\n",
    "            else:\n",
    "                #Either the substring is not shared between the sequences,\n",
    "                # or it has been added to the list already.\n",
    "                pass\n",
    "    \n",
    "    #Return all substrings common to sequences 1 and 2:\n",
    "    return(common_substrings)\n",
    "\n",
    "def compare_overlaps(motif_list, seq):\n",
    "    \"\"\"\n",
    "Compare the current sequence (from 3 to k) to the list of overlaps (motifs) that have been found to far.\n",
    "    \"\"\"\n",
    "    common_substrings = []\n",
    "    #Keep track of the substrings that are common between the\n",
    "    # motif list and the current sequence.\n",
    "    \n",
    "    for motif in motif_list:\n",
    "        #For each motif previously identified\n",
    "        if motif in seq:\n",
    "            #Check if it exists in the current sequence\n",
    "            common_substrings.append(motif)\n",
    "            #and save it if it is.\n",
    "        else:\n",
    "            #Otherwise, it is no longer common between all sequences,\n",
    "            # so pass. Do not track this motif any longer.\n",
    "            pass\n",
    "    \n",
    "    #Return only substrings/motifs that are shared with the current sequence.\n",
    "    return(common_substrings)\n",
    "\n",
    "def pick_longest(motif_list):\n",
    "    \"\"\"\n",
    "Given a list of motifs take the longest motif. \n",
    "Takes the first occurrence if multiple tie for longest, i.e. in alphabetical order.\n",
    "    \"\"\"\n",
    "    motif_list = sorted(motif_list) #sort alphabetically first\n",
    "    longest = sorted(motif_list, key = len, reverse = True)[0]\n",
    "    return(longest)\n",
    "\n",
    "def find_a_shared_motif(input_file):\n",
    "    \"\"\"\n",
    "The complete program as one function:\n",
    " 1. read the fasta file\n",
    " 2. find overlaps between sequences 1 and 2 (create initial motif list)\n",
    " 3. compare motif list with all other sequences, keeping only those they\n",
    "    have in common\n",
    " 4. return the longest motif common to all sequences\n",
    "    \"\"\"\n",
    "    sequence_dict = read_fasta_file(input_file)\n",
    "    \n",
    "    initial_motifs = find_initial_overlaps(\n",
    "                       str(list(sequence_dict.values())[0]),\n",
    "                       str(list(sequence_dict.values())[1]))\n",
    "    \n",
    "    motif_list = initial_motifs\n",
    "    \n",
    "    for other_sequence in list(sequence_dict.values())[2:]:\n",
    "        motif_list = compare_overlaps(motif_list, str(other_sequence))\n",
    "\n",
    "    return(pick_longest(motif_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['DEFG', 'DEF', 'EFG', 'DE', 'EF', 'FG']\n"
     ]
    }
   ],
   "source": [
    "print(find_initial_overlaps(\"ABCDEFG\", \"DEFGHI\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['ABC', 'BCD']"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regex.findall(r\"\\w{3}\", \"ABCD\", overlapped=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'AC'"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "find_a_shared_motif(\"data/Example_finding_a_shared_motif.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I think I have a pretty decent algorithm now. Let's try the 'real' dataset and see how that goes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'TTGGACTCGAAGCCGAGTCCACGCGTAGCGCCTTACTGAGATGATAGTAAGTGCCCCGCTGTCTGGAAGGCTGA'"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "find_a_shared_motif(\"data/rosalind_lcsm.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Success!!\n",
    "\n",
    "That worked! Now let's save this notebook and commit to git."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}