teaching – regulatory genomics

ONA Wykład 1 – Arytmetyka komputerów

Slajdy do wykładu można znaleźć tu: ONA1-Arytmetyka

Kiedy tylko LabKom uruchomi dla nas kursy w moodle’u, przeniesiemy się tam z materiałami do tego przedmiotu. Będę Państwa informował o tym e-mailem.

Pracujemy w narzędziu python3 (nie w octave) i wykonujemy zadania częściowo inspirowane zadaniami z labu 3 do metod numerycznych. Zachęcam do użycia interpretera ipython3, jeśli pracują Państwo na domowym komputerze. Jeśli ktoś ma kłopoty z interpreterem w domu, to może korzystać z serwera jupyter, używając hasła podanego na wykładzie. Proszę wtedy pracować w osobnych podpisanych notebookach. Nieco więcej technicznych/matematycznych szczegółów jest też dostępne w materiałach do wykładu z metod numerycznych.

Będziemy też korzystać z modułu numpy (pisząc “import numpy as np”) a w szczególności z jego typów danych float64 i float128

Uwaga nt. liczb całkowitych. Obecnie python3 ukrywa zupełnie przed nami fakt istnienia ograniczeń wielkości rejestrów przy operacjach na liczbach całkowitych. Aby zobaczyć jak to działa, warto uruchomić interpreter python2.7 i wykonać operacje na dużych liczbach.

W szczególności:

1. Policz epsilon maszynowy (zarówno jako wartość dziesiętną ale też liczbe bitów cechy binarnej przy pomocy sprawdzania dla jakich n 1.0+2.0**(-n) >1.0. Co by się stało, gdybyśmy wykonali tę próbę w okolicy 0.0 zamiast 1.0. Czy taki epsilon maszynowy zależy od liczby bitów cechy, czy mantysy?

2. Oblicz wartość funkcji (pochodzącej od Siegfrieda Rump’a): f(x,y)= 333.75*y**6+x**2*(11*x**2*y**2-y**6-121*y**4-2)+5.5*y**8+x/(y*2) w punkcie x=77617, y=33096. Wykonaj obliczenia w pythonie2, pythonie3, np.float64, np.float32, np.float128 i module Fractions. Jakie błędy obliczeń tam zachodzą? Jeśli ktoś napisze program obliczający tę wartość prawidłowo, i prześle mi przez moodle program i prawidłowy wynik, to dostanie dodatkowy (bonusowy) punkt do zaliczenia.

3(*). Zastanów się nad zadaniem z szeregami zbieżnymi (lub rozbieżnymi) numerycznie ( w materiałach do labu 3). Zaimplementuj liczenie sum S(n) w pythonie, dla dowolnego n i dla kilku szeregów podanych w rozwiązaniu. Które z nich można dokładnie policzyć pry użyciu modułu Fractions?

4(*). Policz wartości e(x) wg rozwinięcia wzoru Taylora dla zadanego maksymalnego stopnia. Porównaj różnicę wyniku jaki otrzymujesz dla x=-90 i x=2.

5(*). Policz sumy częściowe szeregu harmonicznego przy użyciu liczb zmiennoprzecinkowych i przy użyciu modułu Fractions. Porównaj wyniki.

WBO – final results

Sprawdziłem Zadanie 2. Wyniki są całkiem niezłe.

Jedynym powtarzającym się błędem jest to, że wielu spośród Państwa nie używało parametru alternative=”greater” w teście dwumianowym, co prowadziło np. do “wzbogacenia” motywu w obu grupach nawzajem. Oczywiście nie jest to dobry wynik i musiałem obciąć po kilka punktów osobom, które tak zrobiły.

W związku z tym, że kolokwium było dość niemiarodajne (można było zgadnąć trochę, albo głupio się pomylić w pośpiechu), to proponuję następującą skalę ocen minimalnych:

55 punktów lub więcej – ocena dobra+ i zaproszenie na egzamin ustny na 5
50 punktów lub więcej – ocena dobra i zaproszenie na egzamin ustny na 4+ lub 5
46 punktów lub więcej – ocena dostateczna+ i zaproszenie na egzamin ustny na 4, lub 4+
42 punkty lub więcej – ocena dostateczna i zaproszenie na egzamin ustny na 3+ lub 4
35 punktów lub więcej – zaproszenie na egzamin ustny na 3 lub 3+

Proponuję egzamin ustny w dniach 22-25 Czerwca (można zapisywać się przez doodle) ale jeśli komuś ten termin nie pasuje, to można się umawiać indywidualnie. Zakres materiału taki jak na kolokwium.

Osoby niedopuszczone do egzaminu (poniżej 35 punktów) mogą zdawać egzamin ustny we wrześniu (podobnie jak osoby, które nie będą zadowolone ze swojej oceny po egz. ustnym).

ONA zadanie 3 – Steganografia

W naszym ostatnim i zarazem największym zadaniu zaliczeniowym odejdziemy od tematyki epidemiologicznej i zajmiemy się przetwarzaniem obrazów i ich kompresją. Konkretnie, naszym zadaniem będzie napisanie programu pozwalającego na “ukrywanie” informacji w obrazach, czyli steganografię. Wykorzystamy najprostszą metodę steganografii, czyli wykorzystanie modyfikacji najmniej znaczącego bitu do zakodowania “ukrytej informacji”.

Zasada działania:

Załóżmy, że mamy do przekazania komuś wektor bitów S, długości N (może to być obraz, lub tekst – jest to dla nas w tej chwili bez znaczenia). Możemy wziąć jakiś obraz reprezentowany w standardowym formacie RGB, z 8 bitami na kanał i powiedzmy wymiarach 2000×1000 pikseli. Obraz ten, możemy wczytać do pamięci w postaci macierzy M o wymiarach 2000x1000x3 z typem danych ~~int8~~uint8. Zauważmy, że każda z liczb w tej macierzy zawiera intensywność jednego z kolorów podstawowych w postaci wartości od 0 do 255. Niewielkie zmiany tych wartości (np. o 1) nie są łatwo postrzegalne dla ludzkich oczu, ale program komputerowy może bardzo łatwo je wykryć.

Nasze zadanie będzie polegało na tym, że każdy kolejny bit z wektora S “wstawiać” jako ostatni bit kolejnego bajtu macierzy M. Jak go wstawiać? otóż wystarczy, że sprawdzimy czy reszta z dzielenia kolejnego bajtu macierzy M przez 2 jest równa wartości odpowiedniego bitu z wektora S.

M.reshape((2000*1000*3,))[i]%2 == S[i]

Jeśli tak, to nic nie musimy robić. Jeśli nie, to musimy odjąć lub dodać 1 do M.reshape((2000*1000*3,))[i], żeby nasze równanie było prawdziwe (uważając przy tym aby nie odjąć od zera ani nie dodawać do 255).

W ten sposób otrzymujemy nową macierz M, która ma tę własność, że ostatnie bity reprezentacji kolejnych bajtów M są równe kolejnym bitom S. Taki plik możemy zapisać jako obraz i będzie on bardzo subtelnie wizualnie zmieniony, natomiast odbiorca tego pliku, który wie czego szukać może bardzo łatwo odczytać wektor S

Aby rozwiązanie było praktyczne, musimy jeszcze być w stanie podać długość wektora S oraz typ danych do zapisu tego wektora. W naszym rozwiązaniu, pierwszy zakodowany bit (M[0]%2) będzie określał, czy mamy do czynienia z tekstem w alfabecie DNA (M[0]%2==0) czy z obrazem w skali szarości (M[0]%2==1). Kolejne 32 bity będą podawać nam długość ciągu DNA (w postaci jednej 32 bitowej całkowitej liczby dodatniej w kodzie dwójkowym) lub wymiary obrazka (w postaci dwóch 16bitowych całkowitych liczb dodatnich oznaczających odpowiednio szerokość x i wysokość y zakodowanego obrazka).

Zadanie:

Napisz moduł pythonowy (z komentarzami!) o nazwie stegano.py, który realizuje następujące funkcje:

read_image_greyscale (image_file) – zwracający wymiary obrazka oraz ciąg bitów wartości na podstawie pliku w formacie png (2 pkt)
write_image_greyscale (x,y,bit_vector,output_image_file) – zapisujący obrazek do formatu png w skali szarości (2pkt)
read_DNA(text_file) – zwracający ciąg bitów na podstawie pliku tekstowego zawierającego wyłącznie alfabet ACGT (kodowanie bitowe: A=00 ,C=01, G=10, T=11) (2 pkt)
write_DNA(bit_vector,output_text_file) – zapisuje ciąg liter DNA do pliku tekstowego na podstawie ciągu bitów wg powyższego kodowania (2 pkt)
encode_DNA(bit_vector,image_file, output_file) – wczytujący obraz RGB z pliku image file i zakodowujący wektor bit_vector reprezentujący DNA w ostatnich bitach obrazu (zgodnie z powyższym opisem, łącznie z 33 bitami nagłówka). Zmodyfikowany obraz powinien być zapisany do pliku output_file. Jeśli obraz jest za mały aby zakodować cały wektor, program powinien zwrócić błąd. (5 pkt)
encode_image(bit_vector,x,y,image_file,output_file) – wczytujący obraz RGB z pliku image file i zakodowujący wektor bit_vector reprezentujący obraz w ostatnich bitach obrazu (zgodnie z powyższym opisem, łącznie z 33 bitami nagłówka). Zmodyfikowany obraz powinien być zapisany do pliku output_file. Jeśli obraz jest za mały aby zakodować cały wektor, program powinien zwrócić błąd. (5 pkt)
decode(image_file, output_file) – powinien wczytywać obraz RGB z pliku, odczytywać z nagłówka (33 ostatnie bity z pierwszych 33 bajtów) informację o tym, czy zakodowane jest DNA, czy obraz, następnie odkodowywać informację i zapisać (korzystając z write_DNA lub write_image_grayscale) wynik do pliku output_file. (5 pkt)
main(), która powinna umożliwiać zakodowywanie i odkodowywanie informacji przy pomocy parametrów z linii komend obsługiwanych przy pomocy modułu Argparse. (2 pkt)

Terminy:

Rozwiązania (w postaci e-maila z normalnym załącznikami w postaci pliku stegano.py oraz zakodowanych przykładowych danych i dopiskiem [ONA-3] w temacie) w ciągu 3 tygodni od dzisiaj, czyli do 18. VI 2020

Przykładowe dane:

Można popróbować np. zakodować taki obrazek:

lub kod DNA koronawirusa covid-19.fasta

W jakimś ulubionym pliku ze zdjęciem przedstawiającym swoją fizjonomię.

Wersja bonusowa (+5 pkt)

Można wyposażyć nasz program w dodatkową funkcjonalność polegającą na kompresji ciągu bitów przed zakodowaniem. Jeśli ktoś zaimplementuje kompresję (i dekompresję) metodą LZ77 to może otrzymać dodatkowe 5 punktów).

ONA 13 – Równania nieliniowe i metody Monte Carlo

Dziś przede wszystkim zajmowaliśmy się zerami funkcji nieliniowych, ale zrobiliśmy też dygresję (w zasadzie uzupełnienie do wykładu o całkowaniu numerycznym) n.t. metod Monte Carlo. Notatki do równań nieliniowych znajdą Państwo tu. Opis całkowania metodą MC znajdą Państwo np. tu

Zadania na dziś:

1. Napisz program, który znajduje przybliżenie liczby Pi przy pomocy równania x**2+y**2<=1

2. Napisz program, który korzysta z metody bisekcji i oblicza wartość pierwiastka kwadratowego z 2. Jakiego równania zera szukamy? ile iteracji potrzeba?

3. Napisz program, który korzysta z metody stycznych, aby obliczyć tę samą wartość. O ile szybciej możemy znaleźć dokładną wartość?

4 (*). Zaimplementuj metodę iteracji prostej Banacha dla rozwiązania równania Keplera. Możesz posiłkować się tym artykułem

Materiały do powtórki przed kolokwium: ONA_powtorka

WBO 13 – Genome assembly and metagenomics

Today’s lecture is making a full circle to the topics from our first meeting. We’re back at DNA sequence identification from short reads, i.e. genome assembly. The slides are here: wyk13-assembly-meta

We do not have any practical assignments for today (we have done simple exercises on k-mers and DeBruijn graphs in the first lab, so you are welcome to look at them again here ), however we will spend some time reviewing the material.

I have prepared some slides with some keywords from this semester here: wyk14-review

Next week, during the lecture, we will have the test. It will be taken online

ONA 12 – obliczenia symboliczne i pakiety do obliczeń

Kwestie organizacyjne – Na ostatnim wykładzie za dwa tygodnie (4. czerwca) będziemy mieli kolokwium (przeprowadzone online, podczas wykładu). Będzie to test, podobny do wcześniejszych: kolokwium-2017, z tym, że przeprowadzony przy pomocy systemu online. Podobnie zrobimy z egzaminem, w czasie sesji.

Ostatnie, trzecie, zadanie zaliczeniowe opublikuję dopiero za tydzień, ale będą Państwo mieli nadal 3 tygodnie na wykonanie go, czyli do 18. czerwca.

Dzisiejszy wykład porównujący cechy różnych pakietów do obliczeń jest tu: ONA12-Obliczenia-symboliczne-i-pakiety

Zadania na dziś są dość proste i dotyczą przede wszystkim obliczeń symbolicznych i pakietu sympy (można to robić w platformie jupyter)

Przy użyciu pakietu sympy, zbadaj podstawowe możliwości obliczeń symbolicznych: stwórz zmienną sumboliczną x (Symbol(‘x’)), stwórz wyrażenie wielomianowe (np. x**2-2) oraz wyrażenie trygonometryczne (np. sin(x)) i policz dla nich pochodne, całki (oznaczone i nieoznaczone) oraz granice (w 0, z prawej i lewej strony oraz , granice 1/x w nieskonczonosci: sympy.oo,).
Zobacz jak wyniki tego typu zadań wyglądają w Wolfram Alpha
Zarejestruj się w systemie SageMath i tam spróbuj rozwiązać te same problemy
(*) Spróbuj rozwiązać symbolicznie układ równań różniczkowych zwyczajnych z poprzednich zajęć (model Lotki-Volterry albo Zombie Invasion) w jednym z wymienionych wcześniej pakietów

WBO 12 – Short read mapping

Today, we are concerned with mapping short reads to long genomes. the lecture will be, as usual, online and the slides are here: wyk12-ngs

Today’s assignments are mostly related to fast text indexes:

1. Write a function compute_BWT(txt), that returns the Burrows-Wheeler Transform of the given text (as a string) and the arrays C(x) and OCC(i,x) as dictionaries.

2. Apply the function to the following DNA sequence:
CGAGCCGCTTTCCATATCTATTAACGCATAAAAAACTCTGCTGGCATTCACAAATGCGCAGGGGTAAAACGTTTCCTGTAGCACCGTGAGTTATACTTTGT

3. Using the method described during the lecture, compute the LTF mapping, and write a method for identifying short strings. Find the occurrences of the string AAAC in the text above using your function.

4. Write another function (that does splitting the query into two and matching at least one half with no errors) to find occurrences of AAAACTCCGCTGGCATTCACAAAT in the text with at most one error

Homework (sixth out of 5 Write a code that uses simplified BWT matching, based on today’s assignments, to create a BWT of the E. coli genome (Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.dna.chromosome.Chromosome.fa) and find positions of occurrences of sequences from this file ( ecoli_proms.fa ) in it. Please send the code and the resulting positions to me by e-mail before the next lecture. Note – we should absolutely avoid any methods that use quadratic memory size (like generating all suffixes of the genome in memory).

ONA 11 – różniczkowanie i całkowanie numeryczne

Dziś będziemy zajmować się różniczkowaniem i całkowaniem numerycznym. Slajdy są tu: ONA11-calkowanie ONA11-calkowanie

Przydatny będzie nam przede wszystkim pakiet scipy.integrate ale będziemy też korzystać z nowych funkcji pakietu matplotlib, takich jak wykresy strzałkowe

0. Na początek proponuje obejrzeć sobie program opisujący model Lotki-Volterry ( lotka_volterra ) i spróbować wykonać go po kawałku ze zrozumieniem. Jeśli są pytania, to proszę pytać

1. Naszym głównym zadaniem zwykle było zadanie o zombie (tutaj wersja zeszłoroczna). W związku z okolicznościami przyrody, może lepiej zajmiemy się modelem SEIR ewolucji epidemii. Jest to model opisany równaniami:

Gdzie:

S to liczba osób niezarażonych,
E to liczba osób zarażonych, ale jeszcze nie zarażających,
I to liczba osób zarażonych i zarażających
R to liczba osób, które już ozdrowiały i nie są podatne na zarażenie

Model ma następujące parametry:

\beta to współczynnik “skuteczności zarażania”
\sigma to współczynnik “inkubacji” im większy, tym krótsza inkubacja
\gamma to współczynnik “zdrowienia” im większy, tym krótszy czas “zarażania”
\nu to współczynnik “zyskiwania/utraty odporności”, jeśli jest dodatni, część populacji zyskuje odporność bez chorowania, jeśli jest ujemny, część populacji traci odporność po ozdrowieniu
\mu to współczynnik normujący związany ze śmiertelnością populacji. Ten model zakłada stałe N w czasie, ale zawiera przepływy normujące \mu, które związane są z efektami śmiertelności na populacje

Interesują nas następujące pytania:

a) Jak wygląda ewolucja systemu w czasie (w zakresie 0..50), jeśli na początku jest tylko 1000 ludzi niezarażonych, i jeden zarażony i mamy następujące prametry:

\beta=0.9
\gamma=0.2
\sigma=0.5
\mu=0
\nu=0

Rozwiązanie przedstaw na wykresie 4 zmiennych (S, E, I, R) od czasu

b) Jak na sytuację wpłyną zmiany współczynnika \nu? Co jeśli \nu jest ujemne, jak w przypadku grypy? Co jeśli \nu jest dodatnie, co oznaczałoby spontaniczne zyskiwanie odporności, bez fazy zarażania?

c) Zaprezentuj pole wektorowe populacji S i I względem parametrów \beta i \gamma?

d*) Zaprezentuj pole wektorowe 3d dla zmiennych SEI

Praca domowa (za 1 pkt)

Trzeba zaimplementować model z niezerową śmiertelnością (powiedzmy \mu=0.01), wykonać symulacje i opisać jednym zdaniem co się zmienia jeśli chodzi o liczbę chorych infekujących (I) w stosunku modelu z \mu=0.

Żeby dostać ten punkt, trzeba wysłać do mnie jako załącznik plik .py generujący wykresy (SEIR od czasu) dla dwóch modeli (\mu=0 i \mu=0.01) a w treści wiadomości napisać co zaobserwowaliśmy. Czas, jak zwykle, do następnych zajęć.

WBO – Assignment 2 – Zadanie 2

In our second assignment, we will try to combine a few of the problems we have encountered: Searching for related sequences, identifying regulatory regions and motifs and enrichment analysis. All of this will be performed on a provided dataset and the expected result will be a program that solves the problem and a report with your findings.

The assignment will have these components:

1. Identify the closest related protein-coding sequences in the E. coli genome for each of the input protein sequences given in the fasta file (protein_fragments.fa). This should be done by downloading the E. coli protein-coding sequences (from this file – genes_e_coli_new.fa), creating a BLAST database from them and performing a local BLAST search against this database. Importantly, we have protein fragments that need to be matched against DNA sequences, so we have to do something to make this compatible. (Here’s a tutorial on running blast locally.) The result of this step should be a comma-separated list containing: input sequence id, best matching E. coli gene id and the associated e-value. (6 pts)

2. For each of the identified E. coli genes, you should find the associated promoter DNA sequence in the given file (~~proms_e_coli_new.fa,~~ proms_e_coli_fixed.fa). Please note that the input file contained sequences from groups A and B, we speculate (based on the empirical evidence), that these groups of genes should have different regulatory mechanisms. We need to identify 10 sequence motifs present in the promoters associated with group A and 10 motifs associated with group B, independently of each other. All motifs should be of length 15 . The result of this step should be two sets of 10 different motif position-specific matrices in a .pfm format. You can do this with a consensus-like method you implement yourself, or you can use the MEME-suite and process its output. (6 pts)

3. Given the two sets of motifs, we would like to select only the motifs that are specific to group A or group B. For this, we need to compute the occurrences of the motifs from .pfm in promoters from both groups. For this purpose, we should scan the promoter sequences of genes from groups A and B, and obtain a number of positions in these promoter sequences that have a log-odds score higher than 0. Now, given the combined lengths of the promoters in each group and the total number of “hits” in each group, we can use the binomial test to identify the motifs that are significantly enriched in one of the groups against in the other. The result of this step should be a list of all 20 motifs with associated, number of hits in the promoters from group A, and promoters from group B and the associated p-values for enrichment in group A and group B respectively. (4 pts)

All of these results and your implementation, should be described in a short report. (4 pts)

Altogether, you should submit your solution (as a commented .py file, results of tasks 1,2,3 and a report in PDF) in an e-mail to me with [WBO] in the subject line before June 9th.

WBO 11 – sequence motifs

Polish version below

Today in the lecture we were discussing sequence motids. Slides are here: wyk11-motifs-eng

Today’s practicals:

1. Download the promoter sequences from the E.coli genome: ecoli_proms.fa

2. Read about the Bio.motifs module and its implementation of position specific matrices in its documentation

3 . (Homework) Implement the simplified consensus method for motif finding (based on slide 20) and use it on the E.coli promoter sequences for motifs of length 3,4,5,6. The algorithm works as follows

based on any position in the first promoter sequence initialize a PWM
for every following step:
choose the position in one of the remaining sequences that is best matching the current PWM
modify your PWM, by adding the word at the chosen position to your PWM and remove this sequence from further considerations
Perform this procedure for all possible starting points in the first sequence and return a few (chosen by the user, by default 5) resulting PWMs that have the highest information content
Possible improvements can be considered: permuting the input sequence set and/or skipping some sequences that have a particularly bad match

4. Compare your results with the known E. coli promoter motifs

5. Use the MEME motif finding suite on your sequences and compare to your previous results and the real motifs

Homework: This week again you can earn additional 2 points for a homework. This time it is the assignment 3 (consensus method). As usual, please send it before the next week lecture as an attached python (.py) file by e-mail with “[WBO]” in the subject.

Polska wersja: Dziś rozmawialiśmy o reprezentacji i wyszukiwaniu motywów sekwencyjnych: Slajdy są tu: wyk11-motifs-eng

Zadania na dziś:

1. Pobierz listę promotorów e_coli w pliku fasta: ecoli_proms.fa

2. Zapoznaj się z modelem macierzy motywów z modułu Bio.motifs i jego dokumentacją

3. Zaimplementuj prostą, zachłanną metodę consensus dla wyszukiwania motywów, zastosuj ją do promotorów e.coli i obejrzyj wyniki dla motywów długości 3,4,5,6. Schemat algorytmu jest następujący:

zaczynając od dowolnego miejsca w pierwszej sekwencji skonstruuj na jego podstawie model PWM
dla każdej kolejnej sekwencji:
wybieraj sekwencję najlepiej pasującą do obecnego modelu PWM
Po wybraniu sekwencji skonstuuj nowy model PWM, wzbogacony o tę sekwencję
Po wykonaniu tej procedury dla wszystkich punktów startowych z sekwencji 1. zwróć kilka (np. 5) najlepszych (w sensie zawartości informacyjnej) macierzy

Powtórz tę procedurę dla kilku różnych permutacji pliku wejściowego.

4. Porównaj wyniki ze znanymi motywami promotorowymi w E. coli

5. Wykonaj analizę tych samych promotorów przy pomocy MEME i porównaj wyniki z otrzymanymi z consensusa

Praca domowa: W tym tygodniu można dostać 2 punkty za przesłanie do mnie swojej implementacji metody consensus. Jak zwykle mają Państwo tydzień czasu, proszę o wysłanie pliku .py mailem z kodem [WBO] w temacie.