<jupyter_start><jupyter_text>Exercise 1  (a)<jupyter_code>def single_read_locations(single_read,sequence):
    idx=0
    read_locations=[]
    idx=sequence.find(single_read)
    while idx>=0:
        read_locations.append(idx)
        idx=sequence.find(single_read,idx+1)
    return read_locations
    
def map_reads_naive(reads,sequence):
    all_locations={}
    for single_read in reads:
        locations = single_read_locations(single_read,sequence)
        if locations: all_locations[single_read] = locations
    return all_locations

import re

def map_reads_re(reads,sequence):
    all_locations={}
    for single_read in reads:
        locations = [m.start() for m in re.finditer(single_read,sequence)]
        if locations: all_locations[single_read] = locations
    return all_locations

seq="abracadabra"
reads = "a bra cada arba".split()
map_reads_naive(reads,seq)
map_reads_re(reads,seq)<jupyter_output><empty_output><jupyter_text>(b)<jupyter_code>import os
from glob import glob

filenames = glob(os.path.join(folder,"*.faa"))+glob(os.path.join(folder,r"drosophila\*.faa"))
word_list = [word.strip().upper() for word in open(os.path.join(folder,"words.txt")) if len(word)>=5]
from operator import itemgetter

def find_longest_most_word_naive(seq):
    locations = map_reads_naive(word_list,seq)
    return max(locations,key=len,default=""), locations

def find_longest_most_in_proteome_naive(proteome_file):
    result = [(find_longest_most_word_naive(seq),hd) for hd,seq in fasta_sequences(proteome_file)]
    return (max(map(lambda x: (x[0][0],x[1]),result),key=lambda x: len(x[0]),default=("","")),
            max(map(lambda x: (x[0][1],x[1]),result),key=lambda x: len(x[0]),default=(0,"")))

def find_longest_most_in_proteome_files_naive(file_list):
    for proteome_file in file_list:
        longest,most = find_longest_most_in_proteome_naive(proteome_file)
        print("The longest word in '{}' is:".format(os.path.basename(proteome_file)))
        print("'{}' in protein '{}'".format(*longest))
        
        print("The most unique words, namely {}, are found in '{}'".format(len(most[0]),most[1]))
        print(" ".join(most[0].keys()))
%time find_longest_most_in_proteome_files_naive(filenames)<jupyter_output>The longest word in 'ecoli-proteome.faa' is:
'RALLIERS' in protein 'gi|90111683|ref|NP_418505.2| membrane translocase (MDR) of MdtNOP efflux pump, PET family [Escherichia coli str. K-12 substr. MG1655]'
The most unique words, namely 46, are found in 'gi|16128891|ref|NP_415444.1| chromosome condensin MukBEF, ATPase and DNA-binding subunit [Escherichia coli str. K-12 substr. MG1655]'
ADDS AIVER ALAE ALEF ALGA ALGAL AREA AREAL ARES ARIA ARIAS DELE EFTS ELKS FART FIGS GALS GATS KALE KIER LADY LAVA LEAD LEADY LEAL LEAR MARE NEAR PENIS PENS PERK RADS RAKE RARE REAL RELAY SEAL SEAT SERF SLED SNARE SWEDE SWEDES TALI TINS VAIN
The longest word in 'chr2L.faa' is:
'LACTEALS' in protein 'gi|24585071|ref|NP_609917.2| short spindle 3 [Drosophila melanogaster]'
The most unique words, namely 253, are found in 'gi|320544542|ref|NP_001188694.1| Muscle-specific protein 300, isoform D [Drosophila melanogaster]'
AGEE AGER AINE AIRS ALAE ALAS ALEE ALES ALGA ALINE ALTS ANTA ARIL ARKS ARSE AVES CELL CLAW CRAM[...]<jupyter_text>Exercise 2  (a)<jupyter_code>import random

def generate_random_data(genome_size,nr_of_reads,read_length=None,seed=4711):
    genome_size=int(genome_size)
    nr_of_reads = int(nr_of_reads)
    alphabet= "ACGT"
    rng = random.Random(seed)
    rng.seed(seed)
    sequence = ''.join(rng.choices(alphabet,k=genome_size))
    if read_length==None: # If no length is given choose a reasonable random length
        read_length = int(math.log(sequence_size)/math.log(len(alphabet))+1)
        reads = ["".join(rng.choices(alphabet,k=read_length)) for _ in range(nr_of_reads)]
    else:
        reads=[]
        reads = [sequence[idx:idx+read_length] 
                    for idx in (rng.randint(0,genome_size-read_length) for _ in range(nr_of_reads))]
    return sequence,reads

def generate_random_data_alt(genome_size,nr_of_reads,read_length,seed=4711):
    genome_size=int(genome_size)
    nr_of_reads = int(nr_of_reads)
    alphabet= "acgt"
    random.seed(seed)
    sequence = ''.join(random.choices(alphabet,k=genome_size))
    reads=[]
    reads = [sequence[idx:idx+read_length] 
                for idx in (random.randint(0,genome_size-read_length) for _ in range(nr_of_reads))]
    return sequence,reads<jupyter_output><empty_output><jupyter_text>(b)<jupyter_code>import time
def time_fct(fct):
    start=time.perf_counter()
    fct()
    elapsed=time.perf_counter()-start
    return elapsed
#seed = "Holy Hand Grenade of Antioch"
titles=["set i.","set ii.", "set iii."]
genome_sizes = [10**6,10**6,10**7]
number_of_reads = [10**4,10**5,10**4]
read_lengths = [10,11,12]

for tt,gs,nr,rl in zip(titles,genome_sizes,number_of_reads,read_lengths):
    seq,reads = generate_random_data(gs,nr,rl,tt)
    elapsed = time_fct(lambda : map_reads_naive(reads,seq))
    print("For '{}' map reads took {:.1f} seconds".format(tt,elapsed))<jupyter_output>For 'set i.' map reads took 22.4 seconds
For 'set ii.' map reads took 227.8 seconds
For 'set iii.' map reads took 236.5 seconds<jupyter_text>Exercise 3<jupyter_code>from collections import defaultdict

def map_reads_lookup(reads,sequence,lookup_length):
    lookup=defaultdict(list)
    for i in range(len(sequence)-lookup_length+1):
        lookup[sequence[i:i+lookup_length]].append(i)
    all_locations=defaultdict(list)
    for single_read in reads:
        read_prefix = single_read[:lookup_length]
        read_suffix = single_read[lookup_length:]
        srl=len(single_read)
        if read_prefix in lookup:
            for idx in lookup[read_prefix]:
                if read_suffix==sequence[idx+lookup_length:idx+srl]:
                    all_locations[single_read].append(idx)
    return all_locations
for lookup_length in [4,8]:
    print("\nLookup length {}".format(lookup_length))
    for tt,gs,nr,rl in zip(titles,genome_sizes,number_of_reads,read_lengths):
        seq,reads = generate_random_data(gs,nr,rl,tt)
        elapsed = time_fct(lambda : map_reads_lookup(reads,seq,lookup_length))
        print("For '{}' map reads took {:.1f} seconds".format(tt,elapsed))<jupyter_output>Lookup length 4
For 'set i.' map reads took 8.4 seconds
For 'set ii.' map reads took 81.4 seconds
For 'set iii.' map reads took 90.1 seconds

Lookup length 8
For 'set i.' map reads took 0.4 seconds
For 'set ii.' map reads took 0.8 seconds
For 'set iii.' map reads took 4.1 seconds<jupyter_text>Exercise 4  (a)<jupyter_code>from itertools import islice

def lo_up(s,flag):
    return s.upper() if flag else s.lower()

class KeywordTree:
    def __init__(self):
        # The edges are represented by the dictionary
        # the key is the edge label and the value is the child node
        self.children={}
        # Will be set to the word if the node represents a valid word
        self.word = None 
       
    
    def add_word_rec(self,word,original):
        if word=="":
            self.word = original
        else:
            head = word[0]
            tail = word[1:]
            if head not in self.children:
                self.children[head] = KeywordTree()
            self.children[head].add_word_rec(tail,original)

    def add_word(self,word):
        self.add_word_rec(word,word)

    def add_word_list(self,words):
        for w in words:
            self.add_word(w)
            
    def lookup(self,word):
        if not word or word[0] not in self.children:
            return self.word is not None
        else:
            return self.children[word[0]].lookup(word[1:])
        
    def ascii_tree(self,max_degree,max_depth,first_prefix="*",prefix=" "):
        if max_depth==0 or not self.children:
            print(first_prefix)
            return
        items = [(i,label,child) for i,(label,child) in enumerate(islice(self.children.items(),max_degree))]
        last_idx = len(items)-1
        for idx,label,child in items:
            if idx == 0:
                new_first_prefix = first_prefix+"-+-"+lo_up(label,child.word)
            else:
                print(prefix+" |")
                new_first_prefix = prefix+" +-"+lo_up(label,child.word)
            new_prefix=prefix+("    " if idx == last_idx else " |  ")
            child.ascii_tree(max_degree,max_depth-1,new_first_prefix,new_prefix)
            joiner = True
            
sample_tree=KeywordTree()
sample_tree.add_word_list("to tea ted ten a i in inn team and inning tennant anytime any intern toad torque".split())
sample_tree.ascii_tree(8,8)
word_tree=KeywordTree()
#word_list = [word.strip().upper() for word in open(os.path.join(folder,"words.txt"))]
word_tree.add_word_list(word_list)
word_tree.ascii_tree(2,16)<jupyter_output>*-+-a-+-a-+-h-+-e-+-D
  |   |   |   |
  |   |   |   +-i-+-n-+-G
  |   |   |
  |   |   +-l-+-i-+-I-+-S
  |   |       |
  |   |       +-S
  |   |
  |   +-b-+-a-+-c-+-A-+-S
  |       |   |   |
  |       |   |   +-I
  |       |   |
  |       |   +-f-+-T
  |       |
  |       +-b-+-a-+-c-+-i-+-e-+-S
  |           |   |   |
  |           |   |   +-Y
  |           |   |
  |           |   +-t-+-i-+-a-+-L
  |           |
  |           +-E-+-S-+-S-+-e-+-S
  |               |
  |               +-Y-+-S
  |
  +-b-+-a-+-a-+-e-+-D
      |   |   |
      |   |   +-i-+-n-+-G
      |   |
      |   +-b-+-A-+-S-+-s-+-U-+-S
      |       |
      |       +-b-+-i-+-t-+-T-+-e-+-D
      |           |           |
      |           |           +-i-+-n-+-G
      |           |
      |           +-l-+-E-+-D
      |               |   |
      |               |   +-R-+-S
      |               |
      |               +-i-+-n-+-G-+-S
      |
      +-d-+-e-+-l-+-l-+-i-+-u-+-M-+-S<jupyter_text>(b)<jupyter_code>def map_reads_tree(root,sequence):
    if not isinstance(root,KeywordTree):
        reads = root
        root = KeywordTree()
        root.add_word_list(reads)
    locations=defaultdict(list)
    n=len(sequence)
    for i in range(len(sequence)):
        node=root
        j=i
        while j<n and sequence[j] in node.children:
            node=node.children[sequence[j]]
            j=j+1
            if node.word:
                locations[node.word].append(i)
    return locations<jupyter_output><empty_output><jupyter_text>(c)<jupyter_code>from operator import itemgetter

def find_longest_most_word_tree(word_tree,seq):
    if isinstance(word_tree,KeywordTree):
        locations = map_reads_tree(word_tree,seq)
    else:
        locations = another_map_reads_tree(word_tree,seq)
    return max(locations,key=len,default=""), locations

def find_longest_most_in_proteome_tree(word_tree,proteome_file):
    result = [(find_longest_most_word_tree(word_tree,seq),hd) for hd,seq in fasta_sequences(proteome_file)]
    return (max(map(lambda x: (x[0][0],x[1]),result),key=lambda x: len(x[0]),default=("","")),
            max(map(lambda x: (x[0][1],x[1]),result),key=lambda x: len(x[0]),default=(0,"")))

def find_longest_most_in_proteome_files_tree(word_tree,file_list):
    for proteome_file in file_list:
        longest,most = find_longest_most_in_proteome_tree(word_tree,proteome_file)
        print("The longest word in '{}' is:".format(os.path.basename(proteome_file)))
        print("'{}' in protein '{}'".format(*longest))
        
        print("The most unique words, namely {}, are found in '{}'".format(len(most[0]),most[1]))
        print(" ".join(most[0].keys()))
word_tree = KeywordTree()
word_tree.add_word_list(word_list)
%time find_longest_most_in_proteome_files_tree(word_tree,filenames)<jupyter_output>The longest word in 'ecoli-proteome.faa' is:
'RALLIERS' in protein 'gi|90111683|ref|NP_418505.2| membrane translocase (MDR) of MdtNOP efflux pump, PET family [Escherichia coli str. K-12 substr. MG1655]'
The most unique words, namely 46, are found in 'gi|16128891|ref|NP_415444.1| chromosome condensin MukBEF, ATPase and DNA-binding subunit [Escherichia coli str. K-12 substr. MG1655]'
FART TALI GATS TINS PENS SEAT KALE ALEF MARE LEAD LEADY KIER NEAR ELKS LADY RAKE VAIN DELE LEAL LEAR ARIA ARIAS SNARE EFTS AIVER ALAE SLED ARES SERF FIGS LAVA ADDS RARE AREA AREAL REAL RADS RELAY ALGA ALGAL PERK SEAL GALS SWEDE SWEDES PENIS
The longest word in 'chr2L.faa' is:
'LACTEALS' in protein 'gi|24585071|ref|NP_609917.2| short spindle 3 [Drosophila melanogaster]'
The most unique words, namely 253, are found in 'gi|320544542|ref|NP_001188694.1| Muscle-specific protein 300, isoform D [Drosophila melanogaster]'
MADS GAGA LETS LEGER EGER SERS TEEN NIDI STEEP TITS PASS HIKE WINS ANTA FLASK LARS ARSE LIES GH[...]<jupyter_text>(d)<jupyter_code>for tt,gs,nr,rl in zip(titles,genome_sizes,number_of_reads,read_lengths):
        seq,reads = generate_random_data(gs,nr,rl,tt)
        elapsed = time_fct(lambda : map_reads_tree(reads,seq))
        print("For '{}' map reads took {:.1f} seconds".format(tt,elapsed))<jupyter_output>For 'set i.' map reads took 1.8 seconds
For 'set ii.' map reads took 4.2 seconds
For 'set iii.' map reads took 20.8 seconds<jupyter_text>4 Alternative solution using nested dictionaries  (a)<jupyter_code>class AnotherKeywordTree:
    def __init__(self,words=[]):
        # The edges are represented by the dictionary
        # the key is the edge label and the value is the child node
        self.root = {} 
        self.add_word_list(words)
    
    def add_word(self,word):
        node = self.root
        for c in word:
            node = node.setdefault(c,{})
        node['_end'] = None

    def add_word_list(self,words):
        for w in words:
            self.add_word(w)
            
    def lookup(self,word):
        node = self.root
        for c in word:
            if c in node:
                node = node[c]
            else:
                return False
        return '_end' in node
        
    def ascii_tree(self,max_degree,max_depth,node=None,first_prefix="*",prefix=" "):
        if node is None:
            node = self.root
        if max_depth==0 or len(node)==('_end' in node): # A leaf should only have the '_end' key
            print(first_prefix)
            return
        items = [(i,label,child) for i,(label,child) in 
                 enumerate([item for item in islice(node.items(),max_degree) if item[0]!='_end'])]
        last_idx = len(items)-1
        for idx,label,child in items:
            if idx == 0:
                new_first_prefix = first_prefix+"-+-"+lo_up(label,'_end' in child)
            else:
                print(prefix+" |")
                new_first_prefix = prefix+" +-"+lo_up(label,'_end' in child)
            new_prefix=prefix+("    " if idx == last_idx else " |  ")
            self.ascii_tree(max_degree,max_depth-1,child, new_first_prefix,new_prefix)
            joiner = True
            
sample_tree=AnotherKeywordTree()
sample_tree.add_word_list("to tea ted ten a i in inn team and inning tennant anytime any intern toad torque".split())
sample_tree.ascii_tree(8,8)
another_word_tree = AnotherKeywordTree(word_list)
another_word_tree.ascii_tree(2,16)<jupyter_output>*-+-a-+-a-+-h-+-e-+-D
  |   |   |   |
  |   |   |   +-i-+-n-+-G
  |   |   |
  |   |   +-l-+-i-+-I-+-S
  |   |       |
  |   |       +-S
  |   |
  |   +-b-+-a-+-c-+-A-+-S
  |       |   |   |
  |       |   |   +-I
  |       |   |
  |       |   +-f-+-T
  |       |
  |       +-b-+-a-+-c-+-i-+-e-+-S
  |           |   |   |
  |           |   |   +-Y
  |           |   |
  |           |   +-t-+-i-+-a-+-L
  |           |
  |           +-E-+-S-+-S-+-e-+-S
  |
  +-b-+-a-+-a-+-e-+-D
      |   |   |
      |   |   +-i-+-n-+-G
      |   |
      |   +-b-+-A-+-S-+-s-+-U-+-S
      |       |
      |       +-b-+-i-+-t-+-T-+-e-+-D
      |           |
      |           +-l-+-E-+-D
      |               |
      |               +-i-+-n-+-G-+-S
      |
      +-d-+-e-+-l-+-l-+-i-+-u-+-M-+-S<jupyter_text>(b)<jupyter_code>def another_map_reads_tree(tree,sequence):
    if not isinstance(tree,AnotherKeywordTree):
        reads = tree
        tree = AnotherKeywordTree()
        tree.add_word_list(reads)
    locations=defaultdict(list)
    n=len(sequence)
    for i in range(len(sequence)):
        node=tree.root
        j=i
        while j<n and sequence[j] in node:
            node=node[sequence[j]]
            j=j+1
            if '_end' in node:
                locations[sequence[i:j]].append(i)
    return locations<jupyter_output><empty_output><jupyter_text>(c)<jupyter_code>another_word_tree = AnotherKeywordTree(word_list)
%time find_longest_most_in_proteome_files_tree(another_word_tree,filenames)<jupyter_output>The longest word in 'ecoli-proteome.faa' is:
'RALLIERS' in protein 'gi|90111683|ref|NP_418505.2| membrane translocase (MDR) of MdtNOP efflux pump, PET family [Escherichia coli str. K-12 substr. MG1655]'
The most unique words, namely 46, are found in 'gi|16128891|ref|NP_415444.1| chromosome condensin MukBEF, ATPase and DNA-binding subunit [Escherichia coli str. K-12 substr. MG1655]'
FART TALI GATS TINS PENS SEAT KALE ALEF MARE LEAD LEADY KIER NEAR ELKS LADY RAKE VAIN DELE LEAL LEAR ARIA ARIAS SNARE EFTS AIVER ALAE SLED ARES SERF FIGS LAVA ADDS RARE AREA AREAL REAL RADS RELAY ALGA ALGAL PERK SEAL GALS SWEDE SWEDES PENIS
The longest word in 'chr2L.faa' is:
'LACTEALS' in protein 'gi|24585071|ref|NP_609917.2| short spindle 3 [Drosophila melanogaster]'
The most unique words, namely 253, are found in 'gi|320544542|ref|NP_001188694.1| Muscle-specific protein 300, isoform D [Drosophila melanogaster]'
MADS GAGA LETS LEGER EGER SERS TEEN NIDI STEEP TITS PASS HIKE WINS ANTA FLASK LARS ARSE LIES GH[...]<jupyter_text>(d)<jupyter_code>for tt,gs,nr,rl in zip(titles,genome_sizes,number_of_reads,read_lengths):
        seq,reads = generate_random_data(gs,nr,rl,tt)
        elapsed = time_fct(lambda : another_map_reads_tree(reads,seq))
        print("For '{}' map reads took {:.1f} seconds".format(tt,elapsed))<jupyter_output>For 'set i.' map reads took 1.2 seconds
For 'set ii.' map reads took 2.7 seconds
For 'set iii.' map reads took 17.8 seconds
{'fraction_non_alphanumeric': 0.13360637591321672, 'fraction_numerical': 0.022138587558113793, 'mean_word_length': 2.8346774193548385, 'pattern_counts': {'":': 1, '<': 46, '<?xml version=': 0, '>': 46, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '411071', 'license': 'MIT', 'n_tokens_mistral': 6970, 'n_tokens_neox': 6420, 'n_words': 1374, 'path': 'sample solutions/Handout 7.ipynb', 'repo_name': 'RSRamKumar/Programming_Lab_1', 'size': '18068'}