Month: April 2020

ONA – Zadanie 2 – podstawowa liczba odtwarzania wirusa

Naszym zadaniem będzie implementacja funkcji, która będzie wykorzystywac metodę liniowych najmniejszych kwadratów, aby znaleźć podstawową liczbę odtwarzania dla przypadków COVID-19 w Polsce.

W naszym (niezmiernie uproszczonym i nie nadającym się do zastosowań praktycznych, ale ilustrującym kluczowe idee) modelu, przyjmiemy, że początkowy rozwój wirusa w populacji rośnie wg wykładniczego wzrostu:

eq1

gdzie t oznacza czas (w dniach), n_E(t) oznacza liczbę osób chorych w momencie t, R_0 jest szukaną podstawową liczbą odtwarzania, zaś \tau przyjmiemy za równą 1.

Nasz program będzie wykorzystywał fakt, że to równanie po obustronnym zlogarytmowaniu staje się liniowe:

eq2

Jeśli mamy dane o liczbie potwierdzonych chorych w pewnym okresie (dla każdego kolejnego dnia), znamy \tau, to możemy wykorzystać liniową ważoną metodę najmniejszych kwadratów do dopasowania rozwiązania, aby znaleźć log(R_0).

Warto  pamiętać, że w tym przypadku, jeśli nie zastosujemy wag w metodzie najmniejszych kwadratów, nasza metoda będzie minimalizować |log(x)-log(b)|^2, co nie będzie prawidłowe. Jak widać na poniższym rysunku, dopasowanie do danych jest lepsze (zwłaszcza z prawej strony) dla problemu ważonego.

Figure_1Dlatego warto przypomnieć sobie zadanie o ważonym zagadnieniu liniowych największych kwadratów, które ma podane rozwiązanie w materiałach do metod numerycznych. W naszym przypadku wagi w_i są równe n_E(i).

Punktacja:

  1. Napisanie funkcji wykorzystującej funkcję scipy.linalg.lstsq (proszę nie używać funkcji typu polyfit) i obliczającej ważony liniowy problem najmniejszych kwadratów (5 pktów)
  2. Napisanie funkcji wykorzystującej funkcję urlopen() i dane dostępne publicznie aby pobrać pierwsze 28 dni, kiedy w wybranym kraju (np. Polsce) występowały zachorowania (tj. pomijając pierwsze zera w linii odpowiadającej danemu krajowi) (5 pktów)
  3. Napisanie programu, który wykorzystuje te dane, aby obliczyć R_0 zgodnie z wzorami podanymi powyżej. Program powinien też generować wykres podobny do powyższego. (5 pktów)

Rozwiązania należy nadsyłać w ciągu dwóch tygodni (czyli do 13. maja włącznie) na mój adres e-mail z dopiskiem [ONA-2] w tytule. Jak zwykle, proszę o skomentowane pliki .py jako normalne załączniki.

ONA 9 – Interpolacja Funkcji

Tym razem zajmiemy się Interpolacją funkcji przy pomocy wielomianów i funkcji sklejanych.

Teoria z wykładu jeśli chodzi o wielomiany znajduje się tu a jeśli chodzi o funkcje sklejane tu.

Większość interesujących nas dzisiaj funkcji znajdziemy w module scipy.interpolate, ale najprostsze funkcje polyfit  i poly1d znajdują się w module numpy.  W module interpolate interesują nas funkcje:

Na laboratorium będziemy rozwiązywać następujące problemy:

1. Spróbuj zinterpolować funkcję sinus(x) na przedziale [0,math.pi] korzystając z równoodległych N węzłów (np. wygenerowanych używając np.linspace(0,math.pi,N)). Jak zachowuje się błąd średniokwadratowy tej interpolacji dla punktów np.linspace(0.math.pi,1000), gdy N rośnie?

2. korzystając z przykładu pokazanego na wykładzie, gdzie wartości y_i=[1,1,1,2], czy potrafisz dobrać pozycje x_i=[1,x2,x3,4] tak aby uzyskać dowolnie dużą amplitudę interpolowanego (polyfit(x_i,y_i,s=0) wielomianu na przedziale [1,4]? Czy błąd aproksymacji (polyfit, k=2) wielomianem stopnia 2 jest tak samo duży?

3. Spróbuj interpolować te same dane przy pomocy krzywych sklejanych. Czym różnią się splajny dla różnych k (k=1,2,3)?

4. Wygeneruj zaburzone obserwacje wg funkcji y=log_2(x) + scipy.stats.normal() w wielu (~1000) punktach na przedziale 1..100. Czy lepiej będzie interpolować, czy aproksymować aby znaleźć kształt funkcji?

 

WBO 9 – Tree reconciliation

English version below. Today’s lecture will be online, the link can be found in the chat channel.

Dzisiaj na wykładzie zajmowaliśmy się uzgadnianiem drzew.

Slajdy do wykładu:

Zadania na lab:

Dzisiaj zajmiemy się uzgadnianiem drzew. Rozważmy sytuację, gdy mamy 1 drzewo gatunków S i drzewo genów G zgodne z tym drzewem gatunków (dla każdego genu g\in G , s(g)\in S). Na początek możemy:

0. Drzewa możemy reprezentować jako słowniki, gdzie każdemu identyfikatorowi wierzchołka przypisujemy listę identyfikatorów jego potomków, w takiej reprezentacji, w pythonie, drzewa ze slajdu 5 z wykładu P. Góreckiego można zapisać tak:

G={"abcde" : ["abcd","e"], "abcd" : ["ab", "cd"], "e" : [], "ab" : ["a","b"], "cd" : ["c","d"], "a" : [], "b" : [],"c" : [],"d" : [] }
S={"abcde": ["abcd","e"], "abcd": ["abd","c"], "abd":["a","bd"], "bd":["b","d"],"a" : [], "b" : [],"c" : [],"d" : [] ,"e":[]}

1. Napisać algorytm LCA(G,S), który znajduje mapowanie LCA dla drzew G,S i zwraca je w postaci słownika

2. Napisać funkcje liczące koszty DC(G,S), D(G,S), a następnie wykorzystać wzór ze slajdu 10, aby wyliczyć  koszt strat L(G,S)

3. Napisać funkcję fat_tree(G,S), odtwarzającą scenariusz uzgadniający drzewaG i S (nie zawsze optymalny), w którym w korzeniu występują wszystkie duplikacje, zaś wszystkie straty występują możliwie nisko w drzewie gatunków. Scenariusz, to oczywiście etykietowanie drzewa gatunków zdarzeniami ewolucyjnymi, D i L. Zauważmy, że zarówno straty jak i duplikacje występują na krawędziach drzewa gatunków, co oznacza, że możemy zapisać wynik naszej funkcji, jako słownik, który każdej etykiecie węzła w drzewie S, przypisuje listę zdarzeń ewolucyjnych, które w wynikowym scenariuszu następują na krawędzi prowadzącej do tego węzła od jego rodzica.

4. Napisać funkcję optimal_tree(G,S), która znajduje jeden z optymalnych scenariuszy uzgadniających G i S. W tym celu musimy przenieść w drzewie “grubym”, które zrobiliśmy w zadaniu 3,  duplikacje możliwie daleko w dół drzewa, a straty w górę, jeśli to możliwe.

English version:

Today’s lecture was about gene and species tree reconciliation.

We have two sets of slides:

  • More theoretical approach by  Paweł Górecki (slides 1-20) : dlt
  • More applied approach by Terry Speed (slides 14-20): Terry-Speed-Reconciliation

For today’s practicals:

We will implement some algorithms for tree reconiliation. We will consider the scenario, wherewe have 1 species tree S  and one gene tree G, such that \forall g \in G  s(g) \in S).

0. We will represent trees in python as dictionaries, with tree node identifies as keys and lists of descendant node identifiers as values. In this representation, the sample trees from P. Górecki’s slide 5 can be written as below:

G={"abcde" : ["abcd","e"], "abcd" : ["ab", "cd"], "e" : [], "ab" : ["a","b"], "cd" : ["c","d"], "a" : [], "b" : [],"c" : [],"d" : [] }
S={"abcde": ["abcd","e"], "abcd": ["abd","c"], "abd":["a","bd"], "bd":["b","d"],"a" : [], "b" : [],"c" : [],"d" : [] ,"e":[]}

You can use these trees to test your implementations of the following assignments.

1. Implement the  LCA(G,S) function, that will return the LCA mapping for treees G and S as a dictionary, with keys being identifiers from the G tree

2. Implement the deep coalescence  [DC(G,S)] and duplication [D(G,S)] cost functions for trees G and S,  then use the formula from slide 10, to write a function calculating the loss cost L(G,S)

3. Write the function fat_tree(G,S), that returns the (not always optimal) reconciliation scenario for trees G and  S, that has all the needed duplications at the root node and all the losses, as low as possible. A scenario can be represented as a dictionary, assigning lists of events (of type Duplication  or  Loss) to edges of the species tree. In practice, it’s easiest to use identifiers of nodes from the S tree as keys and lists of events on the edge leading to this node from its parent as values.

4. Finally, using all the previous functions, write the optimal_tree(G,S) function, returning one of the optimal reconciliation scenarios for the trees G and S. This is done by taking the “fat” tree, from task 3, and moving all duplication events as low as possible, and the lossess as high as possible.

 

ONA 8 – Aproksymacja funkcji

Dzisiaj na wykładzie omówiliśmy zasadniczo tematy zawarte w wykładzie 12. z metod numerycznych. Jeśli ktoś chciałby doczytać to znajdzie materiały tutaj

Zadania na lab:

0. Rozważmy cztery punkty na płaszczyźnie: (0,4), (1,0), (2,0) i (3,0). Jaka prosta przechodzi najbliżej nich? ułóż układ równań liniowych, który można rozwiązać metodą najmniejszych kwadratów. Wykorzystaj funkcje scipy.linalg.qr aby zobaczyć rozkład QR macierzy A. Użyj funkcji scipy. linalg.lstsq aby znaleźć rozwiązanie. Jakie jest znaczenie wartośći zwróconych przez tę funkcje?

1. Rozważmy dane:

x f(x)
0.00 4.00000000000000e+00
1.25 3.28650479686019e+00
2.50 3.08208499862390e+00
3.75 3.02351774585601e+00
5.00 3.00673794699909e+00
6.25 3.00193045413623e+00
7.50 0.00055308437015e+00
8.75 3.00015846132512e+00
10.00 3.00004539992976e+00

Jak będzie wyglądało dopasowanie met. najlepszych kwadratów funkcji f(x)=a*1+b*exp(-x) do tych danych?

2. Rozważmy dane o obwodzie pnia, trees-stripped ( do wczytywania przyda się funkcja numpy.loadtext). Kolejne kolumny oznaczają tu:

  • obwód pnia
  • wysokość drzewa
  • objętość pozyskanego drewna.

Spróbuj dopasować (metodą najmniejszych kwadratów) objętość drzewa jako funkcję:

  • kombinację liniową obwodu pnia i wysokości drzewa
  • iloczynu wysokości przez obwód
  • kombinację liniową powyższych

Gdzie uzyskujemy najmniejszy błąd przybliżenia?

3(*). Przedstaw na wykresie wynik zadania 1 (wykres punktowy obserwacji i dużo gęściej próbkowany wykres liniowy znalezionej funkcji). Jaki jest problem naszego rozwiązania? czy można jakoś pomóc sobie używając ważonego problemu średnich kwadratów? jak to zrobić dla wagi jednej z obserwacji=0? a jak dla wagi 0.1?

 

WBO 8 – profile HMMs and protein domains

English version below. Today’s lecture is online again, the link was sent by USOS-mail.

Dzisiaj mówiliśmy o zastosowaniach HMM do reprezentacji profilów uliniowień oraz wyszukiwania genów.

Slajdy są tu: WBO-8-profile_hmm

W ramach ćwiczeń, warto jest zapoznać się z bazą PFAM  i możliwością wyszukiwania profilów domen przy pomocy programu HMMER. Mieliśmy też wykorzystać moduł biopythona HmmerIO, ale niestety nie działa on po niedawnych zmianach formatu HMMER.

Możemy jednak wykonać ćwiczenia nie korzystając z python’a ale  wykorzystując narzędzia pfam+hmmer do wyszukania domen i znalezienia ich modeli. Proponuję zacząć od białka tramtrack, które posiada 2 istotnie różne wersje TTK69 i TTK88. Proszę znaleźć te dwie sekwencje aminokwasowe w bazie flybase, ściągnąć je na swój komputer i wrzucić do wyszukiwarki domen pfam/hmmer. Następnie porównać struktury domen tych białek i zobaczyć jak to się ma do struktury egzonów genu tramtrack. Warto zobaczyć jak wyglądają opisy domen w postaci profilów hmm.

English version:

Today we discussed HMM application to protein domain modeling and gene finding.

The slides are available here: WBO-8-profile_hmm

As the practical assignments for today, we will look at the PFAM database and the HMMER program functionality. In theory, we could also use the biopython’s HmmerIO, however it is buggy after recent changes in HMMER, so we will not use python this week.

We will just use the web interfaces to the pfam db and hmmer to search for protein domains. I would suggest starting from the Tramtrack protein,  that has 2 major protein isoforms: TTK69 and TTK88. Your assignment is to find these protein sequences in the flybase database , download them to your computer and submit it for domain search in PFAM hmmer search. Then you can compare the pfam domain search results with the dna gene exon structure and identify which exons correspond to the differences in the identified domains.

ONA 7 – Układy równań liniowych

Dziś zajmujemy się układami równań liniowych w reprezentacji macierzowej.

Wykład będzie prowadzony prz tablicy, więc slajdów nie ma, ale potrzebne materiały ( i dużo więcej niż nam potrzeba) są dostępne w materiałach do wykładu z Metod numerycznych . Nas w szczególności interesują wykłady 5 (eliminacja Gaussa) i 7 (uwarunkowanie problemu). Można też zajrzeć do wykładu 8 (macierze rzadkie), ale jest on o dla nas zdecydowanie zbyt obszerny, a  dla nas to tylko temat poboczny.

Jeśli chodzi o operacje algebry liniowej w pythonie, podczas laboratorium przydadzą nam się następujące biblioteki:

numpy.linalg – podstawowe operacje na macierzach, m.in. cond (A), det(A)

scipy.linalg – nieco więcej operacji na macierzach, m.in. hilbert(N), solve(A,B)

scipy.sparse.linalg – operacje na macierach rzadkich m.in. spsolve(A,B)

scipy.io – wczytywanie różnych ciekawych formatów – np. macierzy rzadkich w formacie matlaba loadmat(f)

Zadania na laboratorium :

  1. Stwórz macierz Hilberta H[i,j]=1/i+j-1 dla sensownie dużego rozmiaru macierzy (100×100, 200×200) i wyświetl ją przy pomocy imshow() z biblioteki matplotlib
  2. Przypomnij sobie rozwiązywanie ukł. równań poprzez rozkład LU. Zaimplementuj ręcznię tę metodę (w wersji bez wyboru el. głównego) i spróbuj rozwiązać tym sposobem układ postaci A= np.array([[1e-18,1.0],[1.0,2.0]]) B=np.array([1.0,4.0]). Porównaj wynik z wynikiem scipy.linalg.solve(A,B). Dlaczego w przypadku tej macierzy wybór elementu głównego ma takie znaczenie?
  3. Korzystając z metody scipy.rand(1000,1000) stwórz macierz losową. Użyj jej do przetestowania rozkładu macierzy  LU (scipy.linalg.lu) w wersji z permutacją i bez. Rozwiąż taki układ (dla b=scipy.rand(1000)) przy pomocy scipy.linalg.solve()
  4. Załóżmy, że mamy teraz jedną macierz A (my użyjemy losowej) i bardzo wiele różnych warunków brzegowych B (my wylosujemy 1000). Spróbuj rozwiązać wszystkie układy równań powstałe z przyrównania tej samej macierzy A do wielu różnych wektorów B. Czy można zamiast korzystać wielokrotnie z funkcji solve() coś przyspieszyć? Np. korzystając z funkcji lu()?
  5. Duże macierze równań liniowych często powstają w problemach inżynierskich. Pobierz jedną z macierzy z kolekcji układów powstałych przy projektowaniu elementów samolotów Boeing (np. nr 38 ) i spróbuj wczytać ją jako macierz rzadką do ipython’a. Następnie rozwiąż ją dla losowych warunków brzegowych b.

WBO 7 – Hidden Markov Models

Today, we will be talking about Hidden Markov Models. The slides are here: wyk7-hmm. Polish version of the notes can be found below. The lecture will be online again, check the chat.mimuw.edu.pl for the link. 

Today’s practicals are here:

0. Take a look at the  hmmlearn package. We will be using it to train hidden Markov Models

1.CpG islands are found frequently in diverse genomes, and human genome in particular. Try to define a 2-state HMM (CpG and background) and train it on the  cpg.fa input file. You can try and do it with four possible emissions (ACGT), and for 16 di-nucleotide emissions (AA,AC,AG,AT,  etc.) Can you interpret the emission matrix and assign the learned states to the CpG and background sequences? Perform the learning process  (Baum-Welch) several times and see if the results are similar. Can you interpret the transition matrix probabilities? What does it tell us about the average length of CpG islands?

It might be worth to take a look at the  MultinomialHMM usage example

2. We can also use HMMs for quantitative signal segmentation. Let us take a sample file – histone methylation data from A. thaliana plant (H3K9me2-crh6-3-chr1.sgr) giving us methylation quantification for different positions in the genome (this is a text file with a position and a value in each line).Use GaussianHMM class to identify 2-state and 3-state models from this data. What are the mean emission values for these data? You can find an example of GaussianHMM here.

Homework (2 pts):

Use the trained CpG model on the following 30 sequences from a file cpg_test.fa. As a result, send me the code (in .py) and the results (in a text file). This time you have 2 weeks for this task.

Polish version below:

Dziś na wykładzie omówimy ukryte modele Markowa i algorytmy rekonstrukcji parametrów z danych. Slajdy są tu: wyk7-hmm

Zadania na dziś:

0. Zapoznaj się z modułem hmmlearn. Będziemy go wykorzystywać do uczenia modeli z emisjami

1. Wyspy CpG znajdują się często w genomach, w szczególności genomie ludzkim. Spróbuj zdefiniować ukryty model Markowa, który ma 2 stany i spróbuj nauczyć go na sekwncji zawartej w pliku cpg.fa. Zrób to zarówno dla sekwencji liter (ACGT), jak i dla sekwencji dwunukleotydów (AA,AC,AG,AT, itp…) Czy możesz zinterpretować macierze emisji i przypisać jeden ze stanów do wysp CpG? Wykonaj kilkakrotnie proces uczenia (Baum-Welch) i zobacz czy wyniki są podobne. Jak interpretujesz prawdopodobieństwa w macierzy przejść. Czy coś możesz powiedzieć o średniej długości wysp CpG?

Warto przyjrzeć się przykładowi użycia klasy MultinomialHMM

2. czasami potrzebujemy użyć łańcuchów Markowa do segmentacji sygnału. Weżmy przykładowy plik – dane o methylacji histonów w rzodkiewniku (H3K9me2-crh6-3-chr1.sgr) zawierający poziom metylacji w różnych pozycjach (plik jest tekstowy, każda linia zawiera informacje o pozycji i wartośći). Wykorzystaj HMM z emisjami Gaussowskimi do segmentacji genomu na 2 lub 3 stany. Jakie są wysetymowane wartości średnie dla różnych stanów? Przykład wykorzystania HMMów z emisjami Gaussowskimi można znaleźć tu 

praca domowa (2 pkt):

Wykorzystaj model nauczony na danych o CpG i przetestuj które z 30 sekwencji w pliku cpg_test.fa są naprawdę wyspami CpG. Jako wynik proszę przysłać program w pythonie i wynik w pliku tekstowym.

ONA 6 – kompresja

Dzisiaj wykład też online. Slajdy są tutaj ONA6-kompresja

Zadania na dziś:

1. Zapoznaj się z modułem gzip . Weź plik tekstowy z sekwencją wirusa COVID-19 (covid-19.fasta), napisz skrypty pythona, które używając modułu gzip zapiszą jego wersję skompresowaną, a następnie zdekompresują. Jaki współczynnik kompresji osiąga się przy takim pliku?

2. Weźmy takie pomocnicze funkcje w pythonie opisujące dyskretną transformatę kosinusową:

import urllib,io
from urllib import request
from PIL import Image
from scipy import fftpack
import numpy as np

URL='http://regulomics.mimuw.edu.pl/wp/wp-content/uploads/2019/04/palmy.png'

def get_image_from_url(image_url=URL, size=(256, 256)):
    file_descriptor = request.urlopen(image_url)
    image_file = io.BytesIO(file_descriptor.read())
    image = Image.open(image_file)
    img_color = image.resize(size, 1)
    img_grey = img_color.convert('L')
    img = np.array(img_grey, dtype=np.float)
    return img

def get_2D_dct(img):
    """ Get 2D Cosine Transform of Image
    """
    return fftpack.dct(fftpack.dct(img.T, norm='ortho').T, norm='ortho')

def get_2d_idct(coefficients):
    """ Get 2D Inverse Cosine Transform of Image
    """
    return fftpack.idct(fftpack.idct(coefficients.T, norm='ortho').T, norm='ortho')

def get_reconstructed_image(raw):
    img = raw.clip(0, 255)
    img = img.astype('uint8')
    img = Image.fromarray(img)
    return img

wykonaj program, który dokonuje Ntego przybliżenia obrazu, przy pomocy zerowania wartośći powyżej Ntego wiersza i Ntej kolumny macierzy get_2D_dct(img) i wyświetl kilkanaście pierwszych przybliżeń.

3. Wyrysuj wykres Entropii dla kanału binarnego w zależności od wartości P(1), na przedziale (0,1)

4. (Praca domowa za 1 pkt) Napisz program, który dla zadanego pliku tekstowego (covid-19) tworzy tablicę częstotliwości znaków, kody Huffmana dla wszystkich znaków oraz wylicza entropię tego pliku i długość kodu Huffmana, który opisywałby cały plik (żeby dostać dodatkowy punkt, oczywiście przesyłamy plik .py jako załącznik maila z dopiskiem ONA w ciągu tygodnia)