Čas 2: Pronalaženje početka replikacije

Gde počinje krug?

Bakterijska DNK je ciklična. Ipak, replikacija bakterijske DNK počinje od istog (istih regiona) - oriC regiona. Unutar oriC regiona nalaze se kratke DNK sekvence koje predstavljaju mesta za koja se protein DnaA vezuje i počinje proces replikacije. Te kratke sekvence, mesta vezivanja DnaA proteina, nazivaju DnaA boksovi. Evolucija je obezbedila višestruka ponavljanja DnaA boksova unutar oriC regiona, kako bi se sa većom verovatnoćom obezbedilo pronalaženje mesta i otpočela replikacija.

Problem: Pronalaženje DnaA boksova unutar oriC regiona

Problem pronalaženja DnaA boksova svodimo na problem pronalaženja višestruko ponovljenih kratkih sekvenci unutar teksta.

Rešenje 1.

Naivni pristup pronalaženja višestruko ponovljenih sekvenci zasniva se na izdvajanju svih podniski dužine k (k-grami) i proveri koliko puta su izdvojene podniske pronađene u tekstu.

Prednosti: Implementacija je jednostavna, predstavlja intuitivno rešenje navedenog problema
Mane: Loša efikasnost algoritma, vremenska složenost - k · |Tekst| · (|Tekst| - k) = O(k · |Tekst|2)

Funkcija pattern_count prebrojava pojavljivanje niske pattern u okviru niske text

In [21]:
def pattern_count(text, pattern):
    n = len(text)
    k = len(pattern)
    count = 0
    for i in range(0, n - k + 1):
        current_pattern = text[i:i + k]
        if current_pattern == pattern:
            count += 1 
    return count
In [40]:
pattern_count('XABABABX', 'AB')
Out[40]:
3

Funkcija frequent_words pronalazi najzastupljenije nukleotidne podsekvence dužine k u zadatom tekstu

In [43]:
def frequent_words(text, k):
    frequent_patterns = set([]) # Skup koji ce cuvati sve najzastupljenije sekvence, bez duplikata

    n = len(text)
    
    # Niz koja ce sluziti za brojanje pojavljivanja podniski. 
    # Indeksi su poceci podniski a vrednosti brojevi pojavljivanja podniske
    counts = [1 for i in range(0, n - k + 1)]

    # Izdvajanje svih podniski i njihovo prebrojavanje
    for i in range(0, n - k + 1):
        pattern = text[i:i + k]
        counts[i] = pattern_count(text, pattern)
    
    max_count = max(counts)
    
    for i in range(0, n - k + 1):
        if counts[i] == max_count:
            frequent_patterns.add(text[i:i + k])
            
    return list(frequent_patterns)
In [77]:
frequent_words('ATATGCTAGTGTCGATGTGCTA', 4)
Out[77]:
['TGCT', 'GCTA']
In [160]:
dna = 'ATCAATGATCAACGTAAGCTTCTAAGCATGATCAAGGTGCTCACACAGTTTATCCACAACCTGAGTGGATGACATCAAGATAGGTCGTTGTATCTCCTTCCTCTCGTACTCTCATGACCACGGAAAGATGATCAAGAGAGGATGATTTCTTGGCCATATCGCAATGAATACTTGTGACTTGTGCTTCCAATTGACATCTTCAGCGCCATATTGCGCTGGCCAAGGTGACGGAGCGGGATTACGAAAGCATGATCATGGCTGTTGTTCTGTTTATCTTGTTTTGACTGAGACTTGTTAGGATAGACGGTTTTTCATCACTGACTAGCCAAAGCCTTACTCTGCCTGACATCGACCGTAAATTGATAATGAATTTACATGCTTCCGCGACGATTTACCTCTTGATCATCGATCCGATTGAAGATCTTCAATTGTTAATTCTCTTGCCTCGACTCATAGCCATGATGAGCTCTTGATCA'
k = 9
print(faster_frequent_words(dna, k))

['TCTTGATCA', 'AAGCATGAT', 'AGCATGATC', 'CTCTTGATC', 'ATGATCAAG', 'GCATGATCA']

Rešenje 2.

Kako se DNK sastoji od 4 različitih tipova nukleotida (Adenin, Timin, Citozin, Guanin) to se i niska koja predstavlja DNK sastoji samo od karaktera 'A', 'T', 'C' i 'G', pa je broj različitih DNK sekvenci dužine k jednak 4k. Kako su, u konkretnom problemu pronalaženja početka replikacije, DnaA boksovi kratke sekvence (do 10-tak nukleotida) sledi da je k mali broj. Iz tog razloga, poboljšanje efikasnosti se može dobiti mapiranjem svih nukleotidnih sekvenci u broj 0 - 4k i čuvanjem njihovog broja pojavljivanja.

Prednosti: Bolja vremenska složenost - O(|Tekst|)
Mane: Lošija prostorna složenost - O(4k), za male vrednosti k ova složenost je podnošljiva

Transformacija DNK sekvence u broj i broja u DNK sekvencu

Ako se za karakterima koji označavaju nukleotide sekvence uvede preslikavanje koje svaki karakter preslikava u cifru od 0 do 3, sekvenca zapisana ovim karakterima se može tumačiti kao broj u brojčanom sistemu sa osnovom 4. Ukupan broj takvih brojeva koji se mogu zapisati sa k cifara je 4k.

Neka je preslikavanje definisano kao:

A -> 0
T -> 1
C -> 2
G -> 3

Sekvenca AAA (dužine 3) bi se preslikavala u broj 0 (000) dok bi se sekvenca GGG preslikavala u broj 3 · 4 2 + 3 · 4 1 + 3 · 4 0 = 4 3 - 1 = 4 k - 1 = 64

Primer

ATCG -> (0123)4 = 0 · 43 + 1 · 42 + 2 · 4 + 3 = 27

Funkcija symbol_to_number vrši preslikavanje simbola nukleotida u odgovarajuće numeričke vrednosti, cifre iz brojčanog sistema sa osnovom 4.

In [37]:
def symbol_to_number(symbol):
    mapping = {
        'A': 0,
        'T': 1,
        'C': 2,
        'G': 3
    }
    return mapping[symbol]
In [38]:
symbol_to_number('C')
Out[38]:
2

Funkcija pattern_to_symbol vrši preslikavanje cele nukleotidne sekvence u odgovarajući dekadni broj.
Prevođenje se vrši rekurzivno primenom Hornerove sheme:

pattern_to_number('ATCG') = pattern_to_number('ATC') * 4 + symbol_to_number('G')

In [30]:
def pattern_to_number(pattern):
    if (len(pattern) == 1):
        return symbol_to_number(pattern)
    else:
        last_symbol = pattern[-1]
        return pattern_to_number(pattern[:-1]) * 4 + symbol_to_number(last_symbol)
In [33]:
pattern_to_number('GGG')
Out[33]:
63

Slično, kao u slučaju prevođenja sekvence u broj, prevođenje broja u sekvencu vrši se prevođenjem vrednosti iz dekadnog u vrednost u sistemu sa osnovom 4. Ukoliko je dužina dobijene sekvence kraća od tražene sekvence dužine k, sekvenca se može dopuniti ciframa 0 (karakterom 'A') sa leve strane

Primer

ATC = AATC = AAATC

Funkcija number_to_symbol za unetu cifru iz sistema sa osnovom 4 vraća odgovarajući simbol nukleotida.

In [49]:
def number_to_symbol(number):
    mapping = {
        0: 'A',
        1: 'T',
        2: 'C',
        3: 'G'
    }
    return mapping[number]

Funkcija number_to_pattern prevodi broj u nukleotidnu sekvencu dužine k.
Prevođenje se vrši rekurzivno:

number_to_pattern (27, 4) = number_to_pattern(27 // 4, k - 1) + number_to_symbol(27 % 4)

In [55]:
def number_to_pattern(number, k):
    if k == 1:
        return number_to_symbol(number)
    else:
        prefix = number // 4
        remainder = number % 4
        return number_to_pattern(prefix, k - 1) + number_to_symbol(remainder)
In [57]:
number_to_pattern(27, 4)
Out[57]:
'ATCG'

Funkcija computing_frequencies formira niz brojaca pojava podniski duzine k u tekstu. Podniske se prevode u numeričku formu i pri svakoj pojavi se odgovarajući brojač uvećava za 1.

In [177]:
def computing_frequencies(text, k):
    frequency_array = [0 for x in range(4**k)]
    for i in range(0, len(text) - k + 1):
        pattern = text[i:i + k]
        j = pattern_to_number(pattern)
        frequency_array[j] += 1
    return frequency_array

Funkcija faster_frequent_words pronalazi najzastupljenije nukleotidne podsekvence u zadatom tekstu korišćenjem optimizacije preslikavanjem podsekvenci dužine k u numeričke vrednosti.

In [178]:
def faster_frequent_words(text, k):
    frequent_patterns = set([])
    frequency_array = computing_frequencies(text, k) # Prebrojavanje podsekvenci

    max_count = max(frequency_array)
    
    for i in range(0, 4**k):
        if frequency_array[i] == max_count:
            frequent_patterns.add(number_to_pattern(i, k))
            
    return list(frequent_patterns)
In [179]:
faster_frequent_words('ATATGCTAGTGTCGATGTGCTA', 4)
Out[179]:
['TGCT', 'GCTA']

Navedeni kod može imati složenost O(|Tekst|) samo u slučaju kada je tekst dužine 4k ili više, zbog petlje koja prolazi 4k zato je potrebno dodatno optimizovati kod. Pravo rešenje vremenske složenosti O(|Tekst|) dobija se korišćenjem sufiksnog stabla.

In [180]:
dna = 'ATCAATGATCAACGTAAGCTTCTAAGCATGATCAAGGTGCTCACACAGTTTATCCACAACCTGAGTGGATGACATCAAGATAGGTCGTTGTATCTCCTTCCTCTCGTACTCTCATGACCACGGAAAGATGATCAAGAGAGGATGATTTCTTGGCCATATCGCAATGAATACTTGTGACTTGTGCTTCCAATTGACATCTTCAGCGCCATATTGCGCTGGCCAAGGTGACGGAGCGGGATTACGAAAGCATGATCATGGCTGTTGTTCTGTTTATCTTGTTTTGACTGAGACTTGTTAGGATAGACGGTTTTTCATCACTGACTAGCCAAAGCCTTACTCTGCCTGACATCGACCGTAAATTGATAATGAATTTACATGCTTCCGCGACGATTTACCTCTTGATCATCGATCCGATTGAAGATCTTCAATTGTTAATTCTCTTGCCTCGACTCATAGCCATGATGAGCTCTTGATCA'
k = 9
print(faster_frequent_words(dna, k))

['TCTTGATCA', 'AAGCATGAT', 'AGCATGATC', 'CTCTTGATC', 'ATGATCAAG', 'GCATGATCA']

Kako DnaA boksovi mogu sadržati pojedine mutacije, pronalaženje najzastupljenijih podeskvenci mora uzeti u obzir i eventualna dozvoljena odstupanja između sadržaja ponovljenih podsekvenci.

Hamingovo rastojanje definise se kao broj pozicija na kojima se dve niske razlikuju.

In [85]:
def hamming_distance(string1, string2):
    n = len(string1)
    m = len(string2)
    distance = 0
    for i in range(min(n, m)):
        if string1[i] != string2[i]:
            distance += 1
    return distance + max(n, m) - min(n, m)
In [89]:
hamming_distance('ATC','AGG')
Out[89]:
2

Funkcija neighbors generiše listu niski koje su na najviše d Hamingovom rastojanju od zadate niske.

In [100]:
def neighbors(pattern, d):
    if d == 0:
        return {pattern}
    if len(pattern) == 1:
        return {'A', 'C', 'G', 'T'}
    neighborhood = set([])
    sufix_neighbors = neighbors(pattern[1:], d)
    for text in sufix_neighbors:
        if hamming_distance(pattern[1:], text) < d:
            for x in ['A', 'C', 'G', 'T']:
                neighborhood.add(x + text)
        else:
            neighborhood.add(pattern[0] + text)
    return list(neighborhood)
In [102]:
neighbors('AAAA',1)
Out[102]:
['ACAA',
 'AAAC',
 'AAGA',
 'AATA',
 'GAAA',
 'ATAA',
 'CAAA',
 'AAAG',
 'TAAA',
 'AAAT',
 'AGAA',
 'AACA',
 'AAAA']

Funkcija immediate_neighbors generiše listu niski koje su na Hamingovom rastojanju 1 od zadate niske.

In [3]:
def immediate_neighbors(pattern):
    neighborhood = {pattern}
    pattern_list = list(pattern) # Konverzija zbog imutabilnosti niski
    
    for i in range(len(pattern)):
        for nucleotide in ['A', 'T', 'C', 'G']:
            if nucleotide != pattern[i]:
                neighbor = pattern_list[:]
                neighbor[i] = nucleotide
                neighborhood.add(''.join(neighbor))
    return list(neighborhood)
In [4]:
immediate_neighbors('AAA')
Out[4]:
['GAA', 'AAG', 'ATA', 'AGA', 'ACA', 'AAT', 'TAA', 'AAC', 'CAA', 'AAA']

Funkcija iterative_neighbors predstavlja iterativnu implementaciju funkcije neighbors.

In [5]:
def iterative_neighbors(pattern, d):
    neighborhood = {pattern}
    for j in range(d):
        new_neighbors = []
        for current_pattern in neighborhood:
            new_neighbors += immediate_neighbors(current_pattern)
        for neighbor in new_neighbors:
            neighborhood.add(neighbor)
    return list(neighborhood)
In [6]:
iterative_neighbors('AAAA', 1)
Out[6]:
['AAAA',
 'AAAT',
 'AGAA',
 'AAGA',
 'TAAA',
 'AAAC',
 'ATAA',
 'AACA',
 'AAAG',
 'ACAA',
 'AATA',
 'CAAA',
 'GAAA']

Funkcija approximate_pattern_count prebrojava pojavljivanje podniski u tekstu koje su na d Hamingovom rastojanju od zadate niske.

In [141]:
def approximate_pattern_count(text, pattern, d):
    k = len(pattern)
    count = 0
    for i in range(len(text) - k + 1):
        if hamming_distance(pattern, text[i:i + k]) <= d:
            count += 1
    return count

Funkcija frequent_words_with_missmatches pronalazi listu svih najzastupljenijih podniski dužine k u zadatom tekstu, dopuštajući pri poređenju niski d Hamingovo rastojanje. Drugim rečima, funkcija pronalazi najzastupljenije podsekvence dužine k sa najviše d dozvoljenih mutacija.

In [144]:
def frequent_words_with_missmatches(text, k, d):
    frequent_patterns = set([])
    close = [0 for i in range(4**k)]
    frequency_array = [0 for i in range(4**k)]
    
    for i in range(len(text) - k + 1):
        neighborhood = neighbors(text[i:i + k], d)
        for pattern in neighborhood:
            index = pattern_to_number(pattern)
            close[index] = 1
    
    for i in range(4**k):
        if close[i]:
            pattern = number_to_pattern(i, k)
            frequency_array[i] = approximate_pattern_count(text, pattern, d)
    
    max_count = max(frequency_array)
    
    for i in range(4**k):
        if frequency_array[i] == max_count:
            pattern = number_to_pattern(i, k)
            frequent_patterns.add(pattern)
        
    return list(frequent_patterns)
In [155]:
frequent_words_with_missmatches('ATATGCTAGTGTCGATGTGCTA', 4, 2)
Out[155]:
['TTGG', 'GGTT', 'GTTG']

G-C Skew

In [182]:
def calculate_skew(dna_sequence):
    skew = [0 for c in dna_sequence]
    last = 0
    for i in range(0,len(dna_sequence)):
        if dna_sequence[i] == 'G':
            skew[i] = last + 1
        elif dna_sequence[i] == 'C':
            skew[i] = last - 1
        else:
            skew[i] = last
        last = skew[i]
    return skew
In [184]:
print(calculate_skew('ACTAGTAGTGCACTAGCTAGCTA'))

[0, -1, -1, -1, 0, 0, 0, 1, 1, 2, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0]
In [185]:
import matplotlib.pyplot as plt

def draw_skew(dna_sequence):
    skew = calculate_skew(dna_sequence)
    n = len(skew)
    positions = [i for i in range(n)]
    plt.plot(positions, skew)
In [186]:
draw_skew('ACTAGTAGTGCACTAGCTAGCTA')
In [189]:
input_file = open('ecoli.txt')
dna_sequence = input_file.read()
draw_skew(dna_sequence)