Sekwencjonowanie DNA w technologiach nowej generacji prowadzi do odczytania milionów krótkich sekwencji DNA zwanych odczytami. Często potrzeba odnaleźć pozycje takich odczytów w genomie, czyli zmapować je do genomu. Proces mapowania odczytów (sekwencji DNA) do genomu referencyjnego polega (w dużym uproszczeniu) na przypisaniu każdemu odczytowi współrzędnych genomowych n, m wskazujących odpowiednio chromosom oraz numer pary zasad w sekwencji tego chromosomu z genomu referencyjnego. Te współrzędne oznaczają miejsce z którego dany odczyt najprawdopodobniej został wygenerowany (chromosom oraz początkową pozycję na tym chromosomie).
Jednym z parametrów wykorzystywanych w analizach genomowych jest tzw. pokrycie. Pokrycie jest zdefiniowane następująco: c = (k,n,x,y), gdzie k określa liczbę odczytów zmapowanych na dany fragment, n oznacza numer chromosomu, natomiast x i y to współrzędne genomowe określające odpowiednio początek i koniec (włącznie z tymi pozycjami) fragmentu chromosomu – segmentu.
Zmapowane odczyty są przechowywane w plikach. Każda linia pliku z odczytami zawiera 3 liczby oddzielone tabulacją:
- chromosom
- współrzędna chromosomowa – miejsce zmapowania odczytu
- długość odczytu
Dla uproszczenia przyjmujemy, że odczyty w plikach są posortowane najpierw względem chromosomu, potem współrzędnej, a następnie długości odczytu. Wiele odczytów może mapować się w jedno miejsce.
Czasem może nas interesować jakaś złożona funkcja operująca na pokryciu z więcej niż jednego eksperymentu. Jeżeli dysponujemy dwoma plikami (A i B) ze zmapowanymi odczytami możemy obliczyć np. iloczyn pokryć. Iloczyn pokryć A z B jest zdefiniowany następująco: cm = (km,n,x,y), gdzie km = kA * kB.
Polecenie:
Napisz w Python’ie następujące funkcje:
- (1 pkt)
read_reads(plik)wczytującą (i zwracającą) odczyty z pliku - (3 pkt)
compute_coverage(odczyty)obliczającą pokrycie na podstawie odczytów. - (1 pkt)
write_coverage(pokrycie, plik_wyjsciowy)zapisującą w pliku o nazwieplik_wyjsciowyobliczone pokrycie w 4 kolumnach oddzielonych tabulacją: chromosom, początek segmentu, koniec segmentu, pokrycie. - (5 pkt)
multiply_coverage(plikA, plikB, plik_wyjsciowy), zapisującą w pliku o nazwieplik_wyjsciowyiloczyn pokryć odczytami z plików o nazwachplikAiplikB.
Komentarz:
- Aby wyjściowy plik z pokryciem nie był zbyt duży, fragmenty genomu o takim samym pokryciu, które nachodzą lub bezpośrednio sąsiadują ze sobą (kolejne pozycje), powinny być połączone w jeden segment. Należy rozważyć różne przypadki rozmieszczenia odczytów w segmencie (patrz przykłady).
- Powyższe funkcje powinny znaleźć się w pojedynczym pliku
zad1_[nr indeksu].py.Funkcje mogą wywoływać siebie nawzajem. W pliku z rozwiązaniem mogą znajdować się też dodatkowe funkcje, ale powinny być one poniżej funkcji wymaganych w rozwiązaniu. - Kod powinien być wyposażony w komentarze.
- Zakładamy, że podczas wykonania, potrzebne pliki (
plikA.txtiplikB.txt) znajdują się w bieżącym katalogu.
Rozwiązania wysyłamy do 29. listopada (do godziny 23:55) e-mailem do prowadzących laboratoria z kopią do wykładowcy, używając w tytule tagu [WDI-1-2015]. Aby uzyskać punkty za zadanie trzeba swoje rozwiązanie przedstawić osobiście prowadzącemu lub prowadzącej laboratoria.
Wszelkie pytania dotyczące treści zadania prosimy umieszczać w komentarzach do tego postu.
W przykładach wszystkie odczyty są zmapowane na tym samym chromosomie. Czy nasz program ma uwzględniać też możliwość, że są na różnych chromosomach?
Tak, trzeba uwzględnić to, że odczyty mogą być zmapowane na różne chromosomy. Wtedy trzeba policzyć pokrycie dla każdego z uwzględnionych chromosomów.
Czy funkcja read_reads() może zwracać odczyty w postaci listy list ??
i czy ta funkcja jako parametr przyjmuje zmienną w postaci: plik= ‘nazwa_pliku’ czy może zmienną w postaci
plik=open(‘nazwa_pliku’,’r’) ??
1. To, co dokładnie zwraca funkcja read_reads jest Państwa decyzją. Może to być na przykład lista list. Ważne, żeby był to ten sam typ danych, którego używa funkcja compute_coverage.
(Analogicznie z pokryciem – wybór typu danych, który będzie reprezentował pokrycie pozostawiamy Państwu.)
2. Przyjmijmy, że we wszystkich funkcjach, które przyjmują jako argument ‘plik’ (‘plik’, ‘plik_wyjsciowy’ itp.), chodzi o napis = nazwę pliku, np. ‘odczytyA’.
W przykładach pokrycia mają takie same przedziały tzn. np
pokrycie A pokrycie B
1 1 2 3 1 1 2 1
1 2 4 4 1 4 6 1
1 4 6 3
w tym przykładzie (1) podobnie jak w przykładzie 2 “przedziały” się pokrywają tzn w pokryciu A występują
przedziały na pozycjach 1 2, 4 6, tak jak w pokryciu B, czy może być tak sytuacja ze te przedziały nie nachodzą na siebie tzn np:
pokrycie A:
1 1 2 1
1 2 5 1
pokrycie B:
1 1 3 1
1 3 7 1
Tak, może się zdażyć taka sytuacja. Trzeba rozważyć różne warianty “nakładania” się pokryć (np. sytuacje w których powstaje więcej segmentów o różnych pokryciach).
W opisie zadania pokrycie jest zdefiniowane następująco: c = (k,n,x,y), natomiast w przykładach została zamieniona kolejność: c=(n,x,y,k). Rozumiem, że rozwiązując zadanie mamy zachować kolejność taką jak w treści zadania(a nie w przykładach)?
Zgadza się, w przykładach omyłkowo została zamieniona kolejność na (n,x,y,k). Zostało to już poprawione. Należy zachować kolejność: (k,n,x,y).
Dot. przykładów: jak to jest, że jeśli wszystkie odczyty w pliku A zostały zmapowane na chromozomie 1, to pokrycia tych odczytów mogą zostać zmapowane na różnych chromozomach?
Przepraszam za podwójny post, jednak po obejrzeniu przykładów, zastanawiam się czy pokrycia nie są raczej zapisywane w formie c = (n,x,y,k)
Czy program ma przewidywać że pokrycie w pliku A będzie dotyczyć np chromosomów :1,2,4 a w pliku B:1,2,5
Pokrycie jest określone 4 parametrami: (k,n,x,y). Jednym z tych parametrów (n) jest chromosom, więc tak należy uwzględnić chromosom w plikach z pokryciami (proszę porównać z przykładami). Jedna uwaga: w przykładach pliki z odczytami A i B zawierają odczyty pochodzące z jednego chromosomu. Poprawny program musi również uwzględniać sytuacje gdy w plikach z odczytami znajdują się odczyty pochodzące z więcej niż jednego chromosomu.