ONA – zadanie 3 (20 pkt) – przetwarzanie macierzy kontaktów

Naszym zadaniem jest stworzenie pakietu umożliwiającego obróbkę macierzy kontaktów chromosomowych. Będzie on obejmował 3 zasadnicze zadania:

1. Wczytywanie listy kontaktów dla wybranego chromosomu do macierzy o zadanej rozdzielczości (5pkt). Mamy dany plik wejściowy, który zawiera 4 kolumny:

chromosom1 pozycja1 chromosom2 pozycja2

Każdy taki wpis (a będą ich miliony) opisuje kontakt pojedyńczej nici DNA pochodzącej z chromosomu1 na pozycji1 z nicią pochodzącą z chromosomu2 na pozycji2. Mając do dyspozycji taki plik oraz podany numer chromosomu C i długość segmentu S powinniśmy skonstruować kwadratową macierz symetryczną M kontaktów (rozmiaru N*N, gdzie N= L/S, gdzie L oznacza długość chromosomu C) pomiędzy segmentami zadanego chromosomu. Np. jeśli chromosom ma długość L=1000000 par zasad, a S=1000, otrzymamy macierz M rozmiaru 1000*1000, a wartość M[i,j] będzie oznaczać ile wierszy opisywało kontakt pomiędzy fragmentem [1000*i:1000*(i+1)] a fragmentem [1000*j:1000*(j+1)].

2. Normalizacja macierzy (10pkt). Mając daną macierz kwadratową M o zadanym rozmiarze N*N, często potrzebujemy dokonać jej normalizacji. W tym przypadku interesuje nas rozwiązanie równania M[i,j]= B[i]*[B[j]*T[i,j]+epsilon, gdzie wektor B (rozmiaru N) opisuje multiplikatywne obciążenia poszczególnych segmentów, zaś od Macierzy T (rozmiaru N*N) oczekujemy, że będzie symetryczna i znormalizowana (sum(T[:,i])=1.0 dla każdego i). Daje nam to układ równać, który wymaga rozwiązania iteracyjnego, minimalizującego błąd (epsilon). Wynikiem tej procedury powinna być macierz T, wektor B i wartość epsilon, a dozwolonym parametrem powinna być maksymalna liczba iteracji i maksymalny Epsilon.

3. Wyświetlanie  macierzy kontaktów (5 pkt). Kiedy już mamy macierz kontaktów (M lub T) to chcielibyśmy móc obejrzeć ją na ekranie lub zapisać do pliku. Potrzebny nam będzie program, który na podstawie macierzy i dodatkowych parametrów wyświetla zadaną macierz metodą imshow(). Program powinien pozwalać na następujące opcje:

• określanie podmacierzy do wyświetlenia (tylko prostokąt o współrzędnych x_1:x_2,y_1:y_2)
• Wygładzanie mapy przy pomocy filtra Gaussowskiego o określonym promieniu
• logarytmiczna lub liniowa mapa kolorów
• wybór schematu kolorów
• Zmniejszanie rozdzielczości dużych map do wizualizacji o stały czynnik (jeśli mapa ma np. 10000×10000 można zadać zmniejszenie jej przed wizualizacją – np. zmniejszenie o czynnik 10 powinno uśrednić dane z kwadratów rozmiaru 10×10 i wyświetlić mapę rozmiaru 1000×1000)

Wynik wysyłamy w postaci pojedyńczego skryptu (jako niespakowany załącznik) na adres prowadzącego zajęcia z dopiskiem ONA-3-2016 w tytule do dnia 20. VI 2016.

UWAGA: Przykłądowy plik z danymi (dla 1 chromosomu o długości 1mln par zasad) można znaleźć tutaj – uwaga, plik jest duży…)

UWAGA 2: Kontakty chromosomowe są dość często wzbogacone w kontakty blisko “przekątnej”, czyli kontakty niedużego zasięgu. W związku z tym warto rozważyć wersję programu, która ignoruje wartości w macierzy na samej przekątnej.

Wskazówki do rozwiązania po pierwszych pytaniach:

1) Czy dla dużych macierzy M (np. 1000×1000) możemy myśleć o rozwiązaniu dokłądnym tego układu równań?

2) Czy normalizacja kolumnowo/wierszowa macierzy T nasuwa nam jakieś oszacowanie na przybliżone wartości B[i]?

3) Czy jeśli ktoś dałby nam wartość przybliżoną B, to umielibyśmy znaleźć krok iteracji poprawiający to oszacowanie przy pomocy macierzy M?

4) Epsilon, o którym wspominam w zadaniu jest pewnego rodzaju skrótem myślowym. Kiedy występuje w równaniu, mamy na myśli różnicę między prawą, a lewą stroną tego konkretnego równania (powiedzmy, mały epsilon). Kiedy pojawia się jako kryterium zakończenia iteracji (duży Epsilon), jest to już różnica między wszystkimi lewymi i prawymi stronami równań. W zagadnieniu najmniejszych kwadratów, typowo ten duży Epsilon to suma kwadratów małych epsilonów. W naszym przypadku, gdzie liczba równań, a więc i składników sumy zależy istotnie od parametrów zadanych przez użytkownika, warto rozważyć duży Epsilon będący maksimum z wartości bezwzględnych małych epsilonów.