Category: teaching

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)

APB 6 – Workflow management software in bioinformatics – from Kepler to Galaxy

Dziś na wykładzie (online – link podam na rocket chat i usos-mail) będziemy mówić o rozwiązaniach do zarządzania przepływami pracy (ang. workflows). Jest to niejako kontynuacja poprzedniego wykładu o LIMS i reproducible research. slajdy są tu: wyk6-workflows

Na ostatnim laboratorium otrzymałem od wszystkich zespołów opisy projektów. Na dzisiejszym laboratorium będziemy znowu robić konsultacje nt. harmonogramów prac zespołów i podziału zadań. Wg poprzedniego harmonogramu, mieliśmy mieć teraz 3 tygodnie do ustalenia już konkretnego planu działań, ale biorąc pod uwagę, że mamy nieco inne warunki i kontakt zarówno wewnatrz zespołów, jak i między mną a Państwem jest nieco utrudniony to proponuję taki plan działań (jak zwykle – to jest program minimum – zachęcam do działania z wyprzedzeniem):

  • 15 IV – zgrubny podział zadań i ramowy harmonogram (wyznaczenie modułów do wykonania przez poszczególne osoby i ram czasowych głównych zadań tzw. kroków milowych – działające moduły, integracja modułów, dokumentacja i tutoriale)
  • 29 IV – szczegółowy opis funkcjonalności do stworzenia przez poszczególne osoby (tu już nie tylko moduły, ale i funkcje i ich wejścia wyjścia, zdefiniowanie punktów integracji tj. gdzie muszą się Państwo w swoich modułach trzymać ściśle specyfikacji, a gdzie mogą Państwo jeszcze zmieniać implementację swobodnie)
  • 13 V – prezentacja początkowych implementacji, dyskusja decyzji projektowych (tutaj oczekuję, że wszystkie zespoły będą miały już założone repozytoria i jakieś zręby początkowych działających implementacji)
  • 10 VI – prezentacje projektów

Tymczasem mamy jeszcze program wykładów i prezentacji:

  • 8 IV – piątek, nie ma zajęć
  • 15 IV – Wykład
  • 22 IV –  Wykład
  • 29 IV – referaty (przeniesione z 22 IV): Kiełek (Cytoscape) i Niedźwiecki (dokowanie)
  • 6 V – referaty (przeniesione z 29 IV): Domżał (?) i Goździewska (SNP)
  • 13 V – referaty (przeniesione z 6 V): Różycka (RNA Seq) i Grynkiewicz (Synthetic lethality testing)
  • 20 V – referaty: Kuśnierz (Dynamika Molekularna?) i Kiljan (metagenomika)
  • 27 V – referaty: Solak (MSA) i Szymczak (?)
  • 3 VI – referaty: Antonowicz (?) i Poziemski (Asemblacja DNA)
  • 10 VI – prezentacje zespołowe

 

WBO 6 – Searching for distant relatives- BLAST

The lecture is again online. I will send the link before. Labs are again after the lecture. English version below:

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

Na ćwiczeniach spróbujemy zająć się wykorzystaniem programu blast: 0. Wczytaj plik w formacie FastQ microbial_reads.fastq (jest już na naszym serwerze jupyter w katalogu WBO, nie trzeba go tam ładować), przy pomocy SeqIO.parse(). Są to odczyty z mikrobiomu jelitowego myszy.

Continue reading “WBO 6 – Searching for distant relatives- BLAST”

ONA 5 – analiza obrazu

Dzisiejsze zajęcia poświęcimy na analizę obrazów 2D i 3D. Slajdy są tu ONA5-Obrazy Wykład jak zwykle w trybie zdalnym, link podam na czacie wydziałowym.

Przydadzą się pakiety: scipy.ndimage i matplotlib.

Zadania na dziś:

0. Pozyskaj przy pomocy kinect swoje zdjęcie w pracowni komputerowej zarówno w postaci macierzy głebokości, jak i obrazów rgb Niestety, z powodu kwarantanny, musimy polegać na przykladowych obrazach z kamery głębokości tu, a dalmierza tu. Instrukcje jak je wczytać  i wyświetlić mogą Państwo znaleźć na naszym serwerze jupyter tutaj

1. Przekształć obraz RGB do skali szarości poprzez uśrednienie składowych.

2. Wyświetl te obrazy na swoim komputerze przy pomocy biblioteki matplotlib

3. Wyświetl histogram obrazu i dokonaj jego “wyrównania” pisząc program w  języku python,  gdzie nowa wartość pixeli o jasności x wynosi n(x) = (x-min)/(max-min)*255 dla parametrów min i max będących maksymalną i minimalną jasnością punktu. To zadanie lepiej wykonywać na pliku o zaburzonym histogramie, np. tym prześwietlonym plikuobrazek.png

4. Zastosuj filtr Gaussowski z biblioteki ndimage rozmiaru k, aby “wygładzić” obraz i obejrzyj wyniki dla różnych k (3,5,7,9,11,…)

5. Napisz program wykrywający krawędzie przy pomocy filtru Sobel’a, porównaj wyniki swojej implementacji (na podstawie wykładu) z wynikam funkcji sobel 

6(*). Napisz program wykrywający cienie na obrazie z miernika odległości i wyświetlający tylko tę część obrazu z kamery, która jest w określonym zakresie odległości od kamery.

7(*). Użyj biblioteki pydicom do wczytania przekrojów przez głowę z projektu visible human. Napisz program, który zamieni te przekroje poziome na prekroje pionowe. Użyj wygładzania krawędzi.

APB – zajęcia 5 – Pakiety LIMS

Na dzisiejszym wykładzie (także online, link podam przez chat wydziałowy) będziemy rozmawiać o problemach z powtarzalnymi badaniami naukowymi (w bioinformatyce, ale nie tylko) i powiemy co nieco o pakietach takich jak BASE czy Galaxy w kontekście ich roli jako rozszerzenie pakietów LIMS. Slajdy są dostępne tu: wyk5-lims

W czasie labów, znów zrobimy sobie konsultacje z zespołami (proszę o rejestrację w tym doodle ). Zakładam, że podczas dzisiejszego spotkania będziemy mogli już omawiać Państwa projekty (przesłane na mój e-mail w formie załącznika – a nie google docs, czy innej formy edytowalnej zdalnie). Oprócz komentarzy z mojej strony nt. samych projektów, będziemy już też rozmawiać o planowanym podziale zadań i harmonogramie.

WBO 5 – multiple sequence alignment methods

English version below:

Dzisiaj mówiliśmy o metodach uliniowienia wielu sekwencji. Slajdy są tu: wyk5

Na zajęciach naszym celem jest zapoznanie się z metodami uliniowienia wielu sekwencji w pakiecie Biopython oraz metodami uliniowienia progresywnego. Warto przeczytać odpowiedni rozdział tutorialu biopythona

0. Tutaj warto, aby każdy z Państwa korzystających z serwera jupyter, założył sobie na serwerze Jupyter własny folder, żeby pliki nam się nie mieszały. Proponuję numer indeksu, albo inny identyfikator jednoznaczny.

1.Weźmy na początek naszą rodzinę sekwencji białek histonowych. Jest dostępna na naszym serwerze jupyter w pliku histones.fa Na początek potrzebujemy wczytać je z pliku (SeqIO.parse) przetłumaczyć je na sekwencje białkowe (metoda translate()) i zapisać do pliku fasta (SeqIO.write) wraz z odpowiednimi identyfikatorami i usuniętymi kodonami stopu (*, która powstaje w wyniku translate)

2. Wykorzystajmy program clustalw (Bio.Align.Applications.ClustalwCommandline) do wykonania mulituliniowienia tych sekwencji na naszym serwerze plik wykonywalny jest zainstalowany w katalogu /usr/bin. Po wykonaniu zadania wczytaj uliniowienie z pliku .aln (AlignIO.read format “clustal”) oraz drzewo filogenetyczne z pliku .dnd (Phylo.read format=”newick”) (osoby pracujące na mac’ach mogą zainstalować sobie clustal Omega)

3. Wykorzystaj teraz program muscle (Bio.Align.Applications.MuscleCommandline) do wykonania multiuliniowienia tych samych sekwencji. Wczytaj je z pliku w formacie fasta

4. Stwórz drzewo filogenetyczne na podstawie uzyskanego uliniowienia metodą neighbor joining jak na poprzednich zajęciach. Porównaj to drzewo do drzewa otrzymanego z programu clustalw

5. Powtórz ćwiczenia 2-4 dla większego zbioru sekwencji (np. Human_PAH_homologs). Czy tutaj widoczne są większe różnice między metodami?

6(*). Praca domowa za 2 punkty. Napisz prostą implementację uliniowienia progresywnego, która wykorzystuje globalne uliniowienie z modułu Bio.pairwise2. To będzie wymagać implementacji funkcji call-back dla substytucji profili i niewielkich zmian w samym module pairwise2

English version

0. If you are using our Jupyter server, please make sure that you work in a separate folder that identifies you by name.

1.Firstly, let us take the histones.fa file. We need to parse it (using Bio.SeqIO.parse method), then translate each of them into protein alphabet (using the  .translate() method) and save to another file using SeqIO.write( remember about removing the * symbols that represent translations of a stop codon)

2. Let us then use the clustalw (Bio.Align.Applications.ClustalwCommandline) to create a multiple alignment of these sequences (the clustalw executable should be in  /usr/bin on the server). After you are able to run it, you need to read the  .aln file with the alignment (AlignIO.read(file, format=”clustal”)) and the guide tree from the  .dnd file (Phylo.read format=”newick”) (students using apple macs can install clustal Omega)

3. Now, for comparison, you can use the  muscle  program (Bio.Align.Applications.MuscleCommandline) to create a multiple alignment of the same sequences. It should give you the alignment in a  fasta format.

4. You can use neighbor joining algorithm we discussed last week to create the tree based on the muscle alignment and compare it visually to the guide tree from clustalw. Do you see differences?

5. If you still have time, you can repeat assignments 2-4 for a larger set of sequences (e.g. Human_PAH_homologs). Are the differences between trees more visible here?

6(*). Homework assignment(2pts). Please implement a simple progressive alignment method that uses Bio.pairwise2 module. It is not as simple as it seems, as you need to implement profile substitution matrix using a callback function and possibly some changes in bio.pairwise2 module.

 

ONA – Zadanie 1- stochastyczny automat komórkowy

Podobnie jak w zeszłym roku, naszym pierwszym zadaniem zaliczeniowym będzie wykonanie programu w pythonie, który symuluje automat komórkowy. Tym razem będzie to tzw. stochastyczny automat komórkowy, gdzie zmiany stanów następują losowo, zgodnie z pewnymi regułami. Będziemy to robić na przykładzie bardzo uproszczonego modelu epidemii podczas kwarantanny. Będziemy używali adaptacji uproszczonego modelu Kermack’a-McKendricka. Osoby zainteresowane mogą poczytać wicej o podobnych modelach w tym artykule:schneckenreither2008

W naszym przypadku model automatu komórkowego to dwuwymiarowa macierz o wymiarach NxN (gdzie N jest parametrem wykonania programu). Każda komórka macierzy reprezentuje jedną osobę i może przyjmować jedną z 3 możliwych wartości:

  • 0 – podatny na zarażenie,
  • 1 – zarażony,
  • 2 – ozdrowieniec (niepodatny na zarażenie).

Macierz jest w większości inicjowana na 0 (cała populacja jest na razie zdrowa i podatna na zarażenie), oprócz k pierwszych zarażonych, rozłożonych losowo w różnych miejscach macierzy (k to kolejny parametr wykonania)

W każdym kroku symulacji, wartości w macierzy będą zmieniały się zgodnie z następującymi regułami:

  • ozdrowieńcy (komórki o wartości 2) nie zmieniają już swojego stanu
  • każda komórka w stanie zarażonym (1) ma szanse na wyzdrowienie (przejście do stanu 2) w każdym kroku równe p_w (parametr)
  • każda komórka w stanie podatnym (0) ma szanse na zarażenie się od każdego chorego sąsiada równe p_z (parametr)

Sąsiedztwo wyznaczamy wg. metody Moore’a, przy czym rozważamy sąsiedztwo o promieniu r, (parametr, domyślnie 1), czyli sąsiadami są pola o odległości Czebyszewa nie większej niż r.

W związku z tym, w każdym kroku możemy obliczyć dla każdego pola liczbę sąsiadów zarażonych Z (liczba pól w zadanym sąsiedztwie, które w obecnym stanie mają wartość 1). Wtedy prawdopodobieństwo zarażenia w danym kroku (zakładamy niezależne zarażenia od wszystkich sąsiadów) to 1- (1-p_z)**Z.

Interesuje nas ewolucja takich układów dla zadanych parametrów N,k,p_w,p_z,r.

Państwa zadaniem jest napisanie programu, który pozwala na:

  1. Przeprowadzenie symulacji (dla zadanych parametrów), która wykonuje maksymalnie M kroków, chyba, że wcześniej zostanie osiągnięty stan, w którym wszyscy są już niepodatni na zarażenie ozdrowieńcami (5 pkt)
  2. Wizualizacja poszczególnych kroków algorytmu w postaci:
    a) wykresu liniowego (liczba komórek o wartościach 0,1,2 w zależności od czasu – 1 pkt)
    b) wizualizacji stanów systemu w kolejnych krokach (przy pomocy funkcji imshow() – 1pkt)
    c) zapisu tych wizualizacji do formatu .png (1 pkt)
    c) animacji przy pomocy klasy animation (2 pkt).

Mają Państwo na to 2 tygodnie, czyli spodziewam się przesłania do mnie e-mailem, z tagiem [ONA] w tytule, rozwiązania w postaci pojedyńczego, skomentowanego pliku pythonowego, jako załącznik, do dnia 1. kwietnia 2020.

Pytania najlepiej przesyłać w postaci komentarzy do tego posta.

ONA 4 – przetwarzanie sygnałów

Dzisiejszy wykład także będzie przeprowadzony zdalnie na platformie google meet. Link wyślę Państwu przez e-mail i chat wydziałowy.

Dziś zajmujemy się przetwarzaniem sygnałów, slajdy są dostępne tu: ONA4-DSP

Animacja transformaty Fourier’a z wikipedii (ze slajdów) jest tutaj, zaś animacje splotu, można znaleźć na stronie wikipedii tu: Convolution

Continue reading “ONA 4 – przetwarzanie sygnałów”