WBO 8 – Algorytm BLAST

Dziś rozmawialiśmy o algorytmach przybliżonego wyszukiwania sekwencji w bazach danych. Slajdy są tu: wyk8

Na ćwiczeniach spróbujemy zająć się wykorzystaniem programu blast:

0. Wczytaj plik w formacie FastQ microbial_reads.fastq (jest już na naszym serwerze jupyter, nie trzeba go tam ładować), przy pomocy SeqIO.parse(). Są to odczyty z mikrobiomu jelitowego myszy.

1. Wybierz kilka losowych, dość długich sekwencji DNA i uruchom dla nich program BLAST online, obejrzyj wyniki (jeśli nic się nie “trafiło”, możesz wybrać inne, dłuższe sekwencje)

2. Wykonaj wyszukiwanie dla tych samych sekwencji przy pomocy interfejsu API biopython’a do blasta online (NCBIWWW) i parsera xml (NCBIXML)

3. Znajdź, która z twoich wybranych sekwencji ma najistotniejsze trafienia do bazy NCBI

4. Przeanalizuj “trafioną” sekwencję. Czy to możliwe, aby ta sekwencja była dokładnie z tego organizmu, który był badany? Porównaj z wynikami algorytmu Smith’a-Waterman’a

 

 

ONA 8 – Metoda najmniejszych kwadratów

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 trzy punkty na płaszczyźnie: (0,6), (1,0) i (2,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+b*exp(-x) do tych danych?

2. Rozważmy dane o obwodzie pnia, trees-stripped ( do wczytywania przyda się funkcja scipy.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 7 – zastosowania HMM

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

Slajdy są tu: wyk7

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ł HmmerIO, ale niestety nie działa on po niedawnych zmianach formatu HMMER.

Zasadniczo powinniśmy móc wykorzystać to narzędzie (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 do siebie 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.

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łądy 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 jest to dla nas tylko temat niejako poboczny.

Jeśli chodzi o operacje algebry liniowej w pythonie, 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 (ew. nierozwiązane jako praca domowa):

  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()/
  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óbój 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.

APB 3 – testy oprogramowania

Slajdy do wykładu są tu: Wyk3-testy

Umówiliśmy się co  do terminów postępowania z projektem zaliczeniowym:

  • ustalenie składu zespołu(ów) +
  • ustalenie tematu projektu –
  • ustalenie zakresu projektu –
  • ustalenie ról w projekcie i podział zadań na osoby –
  • ustalenie harmonogramu uwzględniającego rozwój, testy i dokumentację –
  • raporty z postępów –
  • prezentacja –

WBO 6 – Ukryte modele Markowa

Dziś na wykładzie omówiliśmy ukryte modele Markowa i algorytmy rekonstrukcji parametrów z danych. Slajdy są tu: wyk6

Przy okazji – na naszym serwerze jupyter doinstalowalem mozliwosc uzywania python’a 3, ale niestety przestala dzialac funkcjonalnosc python’a 2. Jesli ktos chce korzystac ze stworzonych wczesniej notebooków, to trzeba wybrac opcję: kernel->change kernel-> python 3

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:

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 – Zadanie zaliczeniowe 2 – operacje na plikach sgr

Naszym zadaniem będzie napisanie programu, który będzie przetwarzał strumienie danych genomicznych zapisane w plikach sgr . Są to bardzo proste pliki, które pozwalają opisywać zmienność parametru wzdłuż chromosomów. Każda linijka ma tylko 3 wartości:

chromosom pozycja wartość

Wartości są oddzielane znakiem tabluacji (“\t”), pozycja jest liczbą całkowitą dodatnią, a wartość jest liczbą zmiennoprzecinkową.

Zakładamy, że wartości pomiędzy punktami w pliku są odcinkami prostoliniowymi, co pozwala nam na podstawie pliku sgr, wyliczyć wartości dla wszystkich pozycji chromosomu.

Przykładowe pliki sgr dla dwóch sygnałów i dwóch filtrów są w pliku Zad2

Pliki te mogą być duże – np. większe niż dostępna pamięć, co oznacza, że nie możemy ich wczytać w całości do wektora.

Będziemy chcieli, żeby nasz program oferował następujące funkcje:

– Suma, różnica, iloczyn, iloraz – zwróć wynik operacji arytmetycznej na  dwóch sygnałach dla tych samych chromosomów. Zakładamy,  że na wejściu mamy dwa pliki, które opisują sygnał wzdłuż tych samych chromosomów, choć niekoniecznie w tych samych pozycjach. OPeracje są wykonywane “po pozycjach”. Zwracamy plik sgr, który opisuje funkcję wynikową (4 pkt)

– wygładzanie – przy pomocy średniej kroczącej o zadanej długości (w parach zasad), zwróć sygnał wygładzony (2 pkt)

– splot – Tym razem na wejściu bierzemy jeden plik sgr z sygnałem (wzdłuż chromosomów) i drugi plik sgr z wartościami funkcji filtra (np. filtr prostokątny, lub Gaussowski), wzdłuż sztucznego chromosomu o nazwie “filtr”. Wynikiem powinien być splot funkcji, czyli przefiltrowany sygnał 6  pkt)

– upraszczanie – dla zadanego pliku sgr wykryj, które punkty są niepotrzebne (tzn po ich usunięciu sygnał nie zmieni się) i zwróć sygnał pozbawiony tych linijek (3 pkt).

Rozwiązanie polega na napisani 7 funkcji, każda z nich powinna działać na 2 plikach (+ parametr szerokość dla wygładzania), wszystko w jednym pliku .py

Rozwiązania podpisane, z dopiskiem [ONA-2018-2] w temacie proszę wysyáć mailem do mnie, do 22. kwietnia 2018.

ONA 6 – Kompresja

Slajdy są tutaj ONA6-kompresja

Zadania na dziś:

1. Zapoznaj się z modułem gzip  zapisz plik tekstowy w pliku skompresowanym ze skryptu pythona, zdekompresuj go przy pomocy programu gunzip i odwrotnie

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

import urllib2,io
import Image
from scipy import fftpack

image_url='http://i.imgur.com/8vuLtqi.png'

def get_image_from_url(image_url='http://i.imgur.com/8vuLtqi.png', size=(128, 128)):
    file_descriptor = urllib2.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 P(1)

4. Napisz program, który dla zadanego pliku tekstowego 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