MiASD/L/z3

From WikiZMSI

Na zajęciach

Problem globalnego porównywania sekwencji:

Dane są dwie sekwencje nukleotydowe v = (v_1, ..., v_n) i w = (w_1, ..., w_n) oraz macierz oceniajaca delta. Należy znaleźć takie zestawienie (alignment) podanych sekwencji, dla którego suma wypłat wg podanej macierzy ocen jest największa.

Wskazówki:

• Wykorzystać (dodać do nowego skryptu) z poprzedniego laboratorium zmienną pomocniczą bases = {"a", "c", "g", "t"} oraz funkcje pomocnicze BToI oraz IToB.

• Wprowadzić w Mathematice macierz ocen delta o wymiarach 5 x 5, gdzie 5-ty wiersz i 5-ta kolumna będą związane z indelem (insercja lub delecja) a pierwsze cztery wiersze / kolumny z nukleotydami a, c, g, t. Początkowo macierz może być wprowadzona w sposób ogólny tj. z parametrami mi (kara z mismatch) i sigma (kara za indel). A następnie można ustalić te parametry na wartości mi = 0.5 i sigma = 0.7. Nagroda za match niech wynosi 1.

• Napisać funkcję GlobalAlignmtent (jako Module) przyjmującą na wejście sekwencje v, w oraz macierz delta i realizującą globalne porównanie tych sekwencji za pomocą programowania dynamicznego. Na czas laboratoriów wystarczy, aby funkcja ta zwracała przeliczony grid liczb s_ij oraz ewentualnie dodatkowo wartość globalnego dopasowania tj. wartość ostatniej komórki s_nm.

Do domu

• Dane są trzy sekwencje DNA dotyczące immunoglobulin dla człowieka, myszy domowej oraz szczura wędrownego pobrane z repozytorium NCBI (źródło - [1]):

Homo sapiens (człowiek) mRNA immunoglobulin heavy chain variable region to JH segment, BA-1F (jakiś fragment) - dł. 418 Link do formatu FASTA - [2]

1 ctcaggatgg actggacctg gaggatcctc ttcttggtgg cagcagcaac aggtgcccac

61 tcccaggtgc acctggtgca atctgggtct gagttgaaga agcctggggc ctcggtgaaa

121 atttcctgca aggcttctgg atacaccttc actagatttg ctgtgaattg ggtccgccag

181 gcccctggac aagggcttga gtggatggga tggatcaaca ctaacactgg aaagccaacc

241 tatgcccagg gcttcacagg acggtttgtc ttttccttgg acacctctgt caccacgtca

301 tatctgcaga tcagcacgct aaaggctgag gacactgccg tttatttctg tgcgagagat

361 ttctttcgag actactttga ctactggggc cagggaaccc tggtcaccgt ctcctcag

Mus musculus (mysz domowa) mRNA immunoglobulin heavy chain - fragment o dł. 492 [3]

1 caggtccagc tgcagcagtc tggagctgag ttggtaaggc ctgggacttc agtgaagatg

61 tcctgcaagg ctgctggata caccttcact aaccactgga taggttgggt aaagcagagg

121 cctggacatg gccttgagtg gattggagat atttaccctg gaagtggtta tactaactac

181 aatgagaagt tcaagggcaa ggccacactg actgcagaca catcctccag cacagcctac

241 atgcagctca gcagcctgac atctgaggac tctgccatct attactgtgc aagatccgat

301 tactacggct cctggtttgc ttactggggc caagggactc tggtcactgt ctctgcagcc

361 aaaacgacac ccccatctgt ctatccactg gcccctggat ctgctgccca aactaactcc

421 atggtgaccc tgggatgcct ggtcaagggc tatttccctg agccagtgac agtgacctgg

481 aactctggat cc

Rattus norvegicus (szczur wędrowny) immunoglobulin heavy chain variable region (VHHAR.2) gene, partial cds - dł. 468 [4]

1 atggacatca ggctcagctt agttttcctt gtccttttca taaaaggtaa ttgataaaag

61 tgtgattatc tctgtggtgt tcacatgaga ataagagagt ttattttgtt ttgttgtgtt

121 agtgatggtt ttctaaccag tattctctgt ttgcaggtgt ccagtgtgag gtgcagctgg

181 tggagtctgg gggaggctta gtgcagcctg gaaggtccct gaaactctcc tgtgcagcct

241 caggattcac tttcagtgac tataacatgg cctgggtccg ccaggctcca aagaagggtc

301 tggagtgggt cgcaaccatt agttatgatg gtagtagcac ttactatcga gactccgtga

361 agggccgatt caccatctcc agagataatg caaaaagcac cctatacctg caaatggaca

421 gtctgaggtc tgaggacacg gccacttatt actgtacaag acacagtg

• Dla ułatwienia sekwencje te są dodatkowo umieszczone plikach tekstowych (bez spacji i łamań): Czlowiek_(homo_sapiens).txt, Mysz_domowa_(mus_musculus).txt,Szczur_wedrowny_(ratus_norvegicus).txt, i można je wygodnie wczytać w Mathematice do zmiennych typu string np. poprzez polecenie typu v = Import["c:/czlowiek_(homo_sapiens).txt"] po ustawieniu odpowiedniej ścieżki.

• Rozrzeszyć funkcję GlobalAlignment o zapamiętywanie dla każdego węzła grida odpowiedniego wskazania (lub kilku wskazań) wstecznego - tzw. backpointery. Innymi słowy oprócz dotychczasowego grida (tablicy) s, funkcja powinna przechowywać i napełniać dodatkową drugą tablicę (np. o nazwie b) przechowującą wskazania wsteczne. Sugestia: zamiast przechowywać w tej tablicy jakieś odpowiedniki strzałek, prościej jest po prostu przechowywać współrzędne węzła grida, do którego należy wrócić. Np. jeżeli b[[10, 20]] = {{9, 20}, {9, 19}}, to oznacza to, że dla węzła o współrzędnych (10, 20) są dwa tak samo dobre (o tej samej wypłacie) wskazania wsteczne - jedno w lewo (indel) do węzła {9, 20} oraz jedno po skosie (match lub mismatch) do węzła (19, 20). Inne własne pomysły na sposób przechowywania backpointerów są oczywiście dopuszczalne.

• Napisać dodatkową funkcję Backtrack (także jako Module) przyjmującą na wejście sekwencje v, w oraz tablicę backpointerów b. Funkcja powinna w pętli przejść wstecz po kolejnych backpointerach i zwrócić na wyjściu napis (lub macierz 3 wierszową) reprezentujący najlepsze globalne zestawienie v i w, np. w następującej lub zbliżonej formie (znakując mismatch jako 'x' oraz indel jako '-'):

{"ctcaggatgga-ctgga-c---ctggag--ga-tcctc-"},
{"        xx     x                     x "},
{"--cagg-tccagctgcagcagtctggagctgagt--tgg"}

Uwaga 1: w trakcie przechodzenia wstecz po backpointerach, w przypadku natrafienia na sytuację, gdzie możliwych jest więcej niż jedno przejście wstecz, dla uproszczenia można przyjąć, że wykonujemy tylko jedno z nich np. dowolne lub pierwsze z brzegu (śledzenie rekurencyjne i generowanie wszystkich możliwych wynikowych zestawień - alignmentów - o tej samej maksymalnej wypłacie poprzez podążanie za wszelkimi możliwymi backpointearami może potencjalnie prowadzić do złożoności wykładniczej). Uwaga 2: w razie podjęcia takiej próby (opcjonalne) może okazać się potrzebne ustawienie w Mathematice zmiennej $RecursionLimit na nieskoczoność.

• Porównać parami trzy podane wyżej sekwencje, oraz wskazać (podczas pokazywania programu prowadzącemu) najbardziej podobną i niepodobną parę.