ONA 10 – Wartości i wektory własne macierzy

Dziś poznaliśmy numeryczne rozwiązania problemu znajdowania dominujacych wartości własnych metodą potęgową i wartości własnych bliskich zadanej wartości metodą odwrotną potęgową. Dużo więcej materiału można znaleźć na tej stronie

Jeśli chodzi o funkcje przydatne w pythonie, to przede wszystkim interesują nas funkcje:

scipy.linalg.eig(A)

oraz

scipy.linalg.eigvals(A)

Zadania na dziś:

1. Zaimplementuj metodę potęgową z wykładu i zastosuj ją dla rozważanej na wykładzie macierzy: M=array([[1,2.],[4,3]]). Czy Twoje rozwiązanie różni się od rozwiązania z funkcji eig? Jakie równanie spełniają wektory własne zwrócone przez tę funkcję i Twoją metodę?

2. Zaimplementuj metodę Rayleigha RQI (z wykładu ). zastosuje ją do znalezienia najmniejszej wartości własnej macierzy M z zad. 1.

3. Rozważ macierz A=array([[-1,2,2],[2,2,-1.],[2,-1,2]]). Spróbuj znaleźć jej wartości własne zaimplementowanymi przez siebie metodami. Skąd biorą się problemy? Czy podobne problemy spotkają nas dla macierzy rand(3,3)?

4. Rozważ macierz M=array([[a,1.],[0,b]]). Gdzie a i b są dodatnimi liczbami całkowitymi. Jak zmienią się wartości własne, gdy zamiast zera wstawimy do M[1,0] wartośći 1/k dla k rosnącego wykładniczo?

WBO 10 – gene families and functions

Slides for today’s lecture are here: WBO-10-gene-families

English version below

Zadania na dziś:

1. Zapoznaj się z testem Fishera  (np w pakiecie scipy.stats.fisher_exact) i metodą GSEA zaimplementowaną w Istytucie Broad’a 

2. Zapoznaj się z formatami plików OBO i GAF z projektu gene ontology: opis formatówprzykładowy plik GAFprzykładowy plik OBO

3. Napisz prosty program, który testuje testem Fishera, wraz z poprawką Bonferroniego wzbogacenie funkcji GO w zadanym zbiorze genów.

4. Przetestuj to na połączonym zbiorze histonów i homologów PAH

5 (*). Zaimplementuj uproszczoną  (wyliczanie statystyki ES dla kolejnych fragmentów rankingu, wykreślanie krzywej wzbogacenia i istotności, bez poprawki FDR) metodę GSEA na podstawie opisu  z suplementu do pracy o GSEA 06580SuppText

 

English version of today’s assignments:

1. Take a look at the  Fisher’s exact test  as implemented in the  scipy.stats.fisher_exact and with the  GSEA   method implemented at the Broad Institute

2. Look at the ontology ( OBO ) i and gene annotation ( GAF) file formats from the gene ontology project: format descriptionsample GAF filesample OBO file

3. Write a simple program that implements gene set enrichment testing on a given set of gene identifiers for a given set of annotations (assume a GAF file is given, and test all the annotated terms in this file against the given set of gene identifiers, score using FIsher’s exact test with bonferroni correction).

4. Test it on the joined set of PAH and histone genes

5 (*). Implement a simplified GSEA approach (computing the ES score and the enrichment curve  for a single gene set and gene ranking, without the FDR)based on supplementary methods to the  GSEA paper 06580SuppText

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.