Month: May 2020

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 int8uint8. 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:

  1. read_image_greyscale (image_file) – zwracający wymiary obrazka oraz ciąg bitów wartości na podstawie pliku w formacie png (2 pkt)
  2. write_image_greyscale (x,y,bit_vector,output_image_file) – zapisujący obrazek do formatu png w skali szarości (2pkt)
  3. 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)
  4. 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)
  5. 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)
  6. 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)
  7. 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)
  8. 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:

Coronaviruses_004_lores

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)

  1. 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,).
  2. Zobacz jak wyniki tego typu zadań wyglądają w Wolfram Alpha
  3. Zarejestruj się w systemie SageMath i tam spróbuj rozwiązać te same problemy
  4. (*) 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:

seireqn

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.

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?