<jupyter_start><jupyter_text>2. Needleman-Wunschov algoritmus s afinným skóre za diery> Implementujte Needleman-Wunshov algoritmus na globálny alignment s afinnými pokutami za diery. Vaša funkcia bude volaná nasledovne    NWAffine(seqA, seqB, match, mismatch, gapopen, gapext)kde - `seqA` je string obsahujúci prvú sekvenciu, - `seqB` je string obsahujúci druhú sekvenciu, - `match` je kladné skóre za zhodu, - `mismatch` je záporná pokuta za nezhodu, - `gapopen` je záporná pokuta za otvorenie diery, - `gapext` je záporná pokuta za predĺženie diery o jeden znak.Funkcia bude vracať skóre globálneho alignmentu sekvencií `seqA` a `seqB`. Volaním tejto funkcie s parametrami: - `match = 1` - `mismatch = -1` - `gapopen = -2` - `gapext = -1`spočítajte skóre pre všetky dvojice sekvencií: - Homo sapiens insulin (INS), transcript variant 1, mRNA - Rattus norvegicus insulin 1 (Ins1), mRNA - Mus musculus insulin II (Ins2), transcript variant 1, mRNAUvedené sekvencie môžete získat' z "NCBI Entrez sequence retrieval system" na adrese http://www.ncbi.nlm.nih.gov/nucleotide?cdb=nucleotide. Hľadajte iba v sekcii nukleových sekvencií. Sekvencie by mali mať dĺžku medzi 400 a 500 báz.<jupyter_code>import numpy as np
import gzip

def loadfasta(filename,verbose=False):
    """ Parses a classically formatted and possibly 
        compressed FASTA file into a dictionary where the key
        for a sequence is the first part of its header without 
        any white space; if verbose is nonzero then the identifiers 
        together with lengths of the read sequences are printed"""
    if (filename.endswith(".gz")):
        fp = gzip.open(filename, 'rt')
    else:
        fp = open(filename, 'r')
    # split at headers
    # data = fp.read().split('>')
    data = fp.read()
    data = data.split('>')
    fp.close()
    # ignore whatever appears before the 1st header
    data.pop(0)     
    # prepare the dictionary
    D = {}
    for sequence in data:
        lines = sequence.split('\n')
        header = lines.pop(0).split()
        key = header[0]
        D[key] = ''.join(lines)
        if verbose:
            print("Sequence %s of length %d read" % (key,len(D[key])))
    return D
def score(a, b, match, mismatch):
    return match if a == b else mismatch
def NWAffine(seqA, seqB, match, mismatch, gapopen, gapext, verbose=False):
    """Needleman-Wunsh algorithm for global alignment with affine penalties
    """
    
    # Init
    A = len(seqA)
    B = len(seqB)
    M = np.zeros((A+1, B+1))
    I = np.zeros((A+1, B+1))
    
    for i in range(A+1):
        M[i][0] = np.NINF
        I[i][0] = gapopen+(i-1)*gapext
        
    for i in range(B+1):
        M[0][i] = np.NINF
        I[0][i] = gapopen+(i-1)*gapext
    
    I[0][0] = M[0][0] = 0 # This needs to be there for the first step
    
    for x in range(1,A+1):
        for y in range(1,B+1):
            I[x][y] = max([
                M[x-1][y]+gapopen,
                I[x-1][y]+gapext,
                M[x][y-1]+gapopen,
                I[x][y-1]+gapext
            ])
            M[x][y] = max([
                M[x-1][y-1]+score(seqA[x-1], seqB[y-1], match, mismatch),
                I[x-1][y-1]+score(seqA[x-1], seqB[y-1], match, mismatch)
            ])
    
    if verbose:
        print(M)
        print(I)
            
    return max(M[A][B], I[A][B])
sequences = [
    loadfasta('NC_000011.10.fasta')['NC_000011.10:c2161209-2159779'],
    loadfasta('NC_000073.6.fasta')['NC_000073.6:c142679726-142678656'],
    loadfasta('NC_005100.4.fasta')['NC_005100.4:272799784-272800351']
]

S = len(sequences)
results = np.zeros((S, S))
for i in range(S):
    for j in range(i+1, S):
        results[i][j] = NWAffine(sequences[i], sequences[j], 1, -1, -2, -1)

print(results)<jupyter_output>[[   0.   43. -545.]
 [   0.    0.  -75.]
 [   0.    0.    0.]]<jupyter_text>Results|      | INS  | Ins1 | Ins2 || :--: | :--: | :--: | :--: || INS  | ∞    | 43   | -545 || Ins1 | 43   | ∞    |  -75 || Ins2 | -545 |  -75 | ∞    |  1. Globálny a lokálny alignment> Implementujte Needleman-Wunshov algoritmus na globálny alignment s lineárnymi pokutami za diery. Vaša funkcia bude volaná nasledovne    NeedlemanWunsh(seqA, seqB, match, mismatch, gap)kde - `seqA` je string obsahujúci prvú sekvenciu, - `seqB` je string obsahujúci druhú sekvenciu, - `match` je kladné skóre za zhodu, - `mismatch` je záporná pokuta za nezhodu, - `gap` je záporná pokuta za dieru dĺžky 1.Funkcia bude vracať skóre globálneho alignmentu sekvencií `seqA` a `seqB`.<jupyter_code>def NeedlemanWunsh(seqA, seqB, match, mismatch, gap, verbose=False):
    """ NW is the same as NWAffine but it penalizes gap extensions same as gap openings
    """
    return NWAffine(seqA, seqB, match, mismatch, gap, gap, verbose)<jupyter_output><empty_output><jupyter_text>> Implementujte Smith-Watermanov algoritmus na lokálny alignment s lineárnymi pokutami za diery. Vaša funkcia bude volaná nasledovne    SmithWatermann(seqA, seqB, match, mismatch, gap)kde - `seqA` je string obsahujúci sekvenciu, - `seqB` je string obsahujúci sekvenciu, - `match` je kladné skóre za zhodu, - `mismatch` je záporná pokuta za nezhodu, - `gap` je záporná pokuta za dieru dĺžky 1.Funkcia bude vracať skóre lokálneho alignmentu sekvencií `seqA` a `seqB`.<jupyter_code>def SmithWatermann(seqA, seqB, match, mismatch, gap, verbose=False):
    """Smith-Watermann algorithm for local alignment with linear penalties
    """
    
    A = len(seqA)
    B = len(seqB)
    F = np.zeros((A+1, B+1))
    
    high_score = 0
    high_x = 0
    high_y = 0
    
    for i in range(1, A+1):
        for j in range(1, B+1):
            F[i][j] = max([
                0,
                F[i-1][j-1]+score(seqA[i-1], seqB[j-1], match, mismatch),
                F[i-1][j]+gap,
                F[i][j-1]+gap
            ])
            if F[i][j] > high_score:
                high_score = F[i][j]
                high_x = i
                high_y = j
    
    if verbose:
        print(F)
    
    return high_score<jupyter_output><empty_output><jupyter_text>Pomocou oboch funkcií s parametrami match = 1, mismatch = -1 a gap = -2 spočítajte skóre prevšetky dvojice sekvencií: - H.fasta - I.fasta - J.fastaVýsledné skóre porovnajte.<jupyter_code>sequences = [
    loadfasta('H.fasta')['seq'],
    loadfasta('I.fasta')['seq'],
    loadfasta('J.fasta')['seq']
]

S = len(sequences)
results = np.zeros((2, S, S))
for i in range(S):
    for j in range(i+1, S):
        results[0][i][j] = NeedlemanWunsh(sequences[i], sequences[j], 1, -1, -2)
        results[1][i][j] = SmithWatermann(sequences[i], sequences[j], 1, -1, -2)

print(results)<jupyter_output>[[[   0. -214. -206.]
  [   0.    0. 1020.]
  [   0.    0.    0.]]

 [[   0. 1292.  727.]
  [   0.    0. 1042.]
  [   0.    0.    0.]]]
{'fraction_non_alphanumeric': 0.13693374649549947, 'fraction_numerical': 0.03659436328759038, 'mean_word_length': 2.6618044300378174, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '272240', 'license': 'MIT', 'n_tokens_mistral': 2841, 'n_tokens_neox': 2660, 'n_words': 690, 'path': 'Bioinformatics Algorithms/Homework 4/Task 4.ipynb', 'repo_name': 'JakubMifek/mff-docs', 'size': '6777'}