Elegancki wiersz poleceń

Magda Mielczarek
Magda Mielczarek 27 lut, 7 minut czytania

 

System operacyjny Linux od lat jest jednym z najczęściej wykorzystywanych środowisk obliczeniowych m.in. przez naukowców. Może być on używany za pośrednictwem graficznego interfejsu użytkownika (GUI, z ang. Graphical User Interface), jednak używanie go z tak zwanego interfejsu wiersza poleceń (CLI z ang. Command-Line Interface) jest znacznie bardziej wydajne, co powoduje, że to właśnie CLI często stanowi domyślne środowisko pracy bioinformatyka. CLI umożliwia użytkownikowi integrację różnorodnych narzędzi (np. do analizy bioinformatycznej), automatyzację pracy, a także wydajną integrację z systemami kolejkowymi w środowiskach o wysokiej mocy obliczeniowej. Niekwestionowaną zaletą Linuxa jest także to, że zestaw narzędzi dostępnych z poziomu CLI jest rozwijany latami przez społeczność międzynarodową (zalety oprogramowania typu open-source) i dzięki temu jest bardzo rozbudowany. Dlatego właśnie większość czynności wykonywanych w wierszu poleceń można wykonać na wiele sposobów. CLI to potężne narzędzie ukryte za prostymi poleceniami pozwalającymi na szybką manipulację ogromnymi zbiorami danych np. biologii molekularnej takich jak np. całe genomy czy transkryptomy. Poniższy artykuł nie jest wprowadzeniem do CLI lecz ma na celu przybliżenie kilku popularnych narzędzi i trików przydatnych w przetwarzaniu danych sekwencji biologicznych.

 

Przykład 1. Pobieranie chromosomu myszy domowej.

 

Pierwszym przykładem może być wygodne pobieranie genomu referencyjnego. Genom referencyjny to pełna informacja genetyczna dla danego gatunku, która stanowi punkt odniesienia w badaniach naukowych. Genomy referencyjne są dostępne w wielu biologicznych bazach danych takich jak NCBI [1], Ensebml [2] czy USCS [3]. Genom referencyjny występuje w formacie FASTA , który charakteryzuje się tym, że zawiera opis sekwencji (czyli nagłówek, zaczynający się od znaku „>”) oraz z kolejnych linii zawierających sekwencję nukleotydową/białkową. Genom występuje zazwyczaj w wielu plikach tekstowych, w których przechowywane są m.in. chromosomy, genom mitochondrialny czy takie fragmenty genomu, które nie zostały (jeszcze) dopasowane do żadnego z chromosomów. Chociaż pliki można pobrać za pomocą przeglądarki internetowej, zwykle wygodniejsze i szybsze jest użycie polecenia wget, po którym następuje adres internetowy pliku, który chcemy pobrać (Przykład 1).

 

Przykład 2. Pobieranie pełnego garnituru chromosomów myszy domowej.

 

Polecenie wget pobiera plik w formacie FASTA dla chromosomu 1 myszy z bazy danych Ensembl. Wynik „Resolving…” pokazuje, że wget próbuje połączyć się z serwerem i ma zamiar rozpocząć pobieranie. Pasek postępu informuje na bieżąco o stanie procesu. Ostatni wiersz wyświetlonego komunikatu wskazuje, że plik został zapisany. Polecenie wget można zmodyfikować, aby zapisać więcej plików (Przykład 2). Poniższy przykład pobiera wszystkie chromosomy, z wyłączeniem innych elementów genomowych (np. mitochondrialnego DNA).

 

Zakres nazw chromosomów podany jest w nawiasach klamrowych ({}), podczas gdy [1-9XY] dopasowuje dowolny znak w nawiasach kwadratowych ([]). Sekwencja „1-9” odpowiada dowolnej cyfrze. Polecenie pobiera zatem autosomy od 1 do 9 oraz chromosomy płciowe X i Y. Znak „?” reprezentuje jeden dowolny znak. Dlatego pobierane są wszystkie autosomy, których nazwy rozpoczynają się od „1” (np. 13), po którym następuje dowolny pojedynczy znak. Warto zauważyć, że pliki z genomem referencyjnym są nazwane zgodnie z następującym wzorcem: (gatunek).
(wersja_genomu_odniesienia).(typ_cząsteczki).(chromosom).(nazwa_chromosomu).fa.gz, zatem plik Mus_musculus.GRCm38.dna.chromosome.1.fa.gz przedstawia sekwencję DNA całego pierwszego chromosomu myszy domowej, odpowiadającą wersji GRCm38 genomu referencyjnego.

Przykład 3. Wyświetlenie zawartości bieżącego katalogu.

 

Listę zawartości bieżącego katalogu (w tym przypadku tego, do którego pobraliśmy piki) wykonuje polecenie ls (Przykład 3).

Nazwy plików mają rozszerzenie .gz, co oznacza, że pliki są skompresowane. To czy pliki należy rozpakować zależy od tego do jakich celów będą używane. W rozpatrywanym przypadku celem jest połączenie wszystkich chromosomów w jeden plik (jak to często ma miejsce w przypadku bioinformatycznej analizy sekwencji całego genomu) i przechowywanie ich w rozpakowanej formie (Przykład 4).

 

Przykład 4. Połączenie plików za pomocą narzędzia zcat oraz przekierowanie wyniku do nowego pliku.

 

Przykład 5. Alternatywne połączenie plików za pomocą narzędzia zcat oraz przekierowanie wyniku do nowego pliku. Ten sposób zapewnia połączenie plików w odpowiedniej kolejności.

Pliki można połączyć za pomocą polecenia „zcat”, gdzie gwiazdka („*”) oznacza dowolne znaki. Domyślnie scalone pliki zostaną wypisane na ekranie, co biorąc pod uwagę rozmiar genomu myszy, nie jest pożądane. Zamiast tego znak „>” przekieruje dane wyjściowe do pliku Mus_musculus.GRCm38.dna.ref.fa, łącząc wszystkie pliki o nazwach rozpoczynających się od „Mus_musculus.GRCm38.dna.chromosome.” i kończąc na rozszerzeniu „.fa”. Domyślnie polecenie zcat (cat dla plików nieskompresowanych) łączy pliki w porządku alfabetycznym, więc poszczególne chromosomy zostaną uporządkowane w następujący sposób: 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1, 2, 3, 4, 5, 6, 7, 8, 9, X, Y. Aby uporządkować pliki od chromosomu 1 do 19 i umieszczając na końcu chromosomy płci, można użyć np. pętli for (Przykład 5).

W pierwszej iteracji pętli następuje wyświetlenie zawartości pliku Mus_musculus.GRCm38.dna.chromosome.1.fa, w kolejnej iteracji, zawartości pliku Mus_musculus.GRCm38.dna. chromosome.2.fa i tak dalej aż do plików zawierających chromosmy 19, X i Y. Informacja odnośnie nazwy chromosomu jest przechowywana w zmiennej chr. Dane wyjściowe pętli zostaną zapisane do pliku o nazwie M.musculus.GRCm38.ref.num.order.fa.

Warto sprawdzić czy pliki zostały połączone w żądany sposób, do tego zadania przydatne będą wyrażenie regularne będące narzędziem do dopasowywania wzorców w sekwencji znaków. I tak na przykład, polecenie grep wyszukuje w plikach wiersze zawierające podany wzorzec. Warto znaleźć zatem te wiersze M.musculus.GRCm38.ref.num.order.fa, które zawierają znak większości („>”) stanowiący charakterystyczny znak nagłówka dla każdej sekwencji. Należy pamiętać, że znaki specjalne, takie jak znak większości („>”), należy podawać w cudzysłowie (Przykład 6). W przeciwnym razie wskazuje przekierowanie danych wyjściowych, jak w poprzednim przykładzie (Przykład 5). W wyniku użytego polecenia wyświetlane są jednie te wiersze, które stanowią nagłówek i są one wyświetlane w takiej kolejności w jakiej występują w pliku.

Przykład 6. Wyświetlenie tylko tych linii pliku FASTA, gdzie występuje tylko nagłówek sekwencji.

 

Przykład 7. Po pobraniu pliku w formacie FASTQ i rozpakowaniu, następuje zliczenie sekwencji (7488403) i obliczenie ich średniej długości (50 pz).

 

Awk to język programowania, który służy do wyszukiwania i przetwarzanie wzorców w plikach lub strumieniach danych, i który jest domyślnie dostępny w linii poleceń bash. Poniższe polecenie (Przykład 7) podaje całkowitą liczbę odczytów, a także ich średnią długość obliczoną na podstawie pliku w formacie FASTQ. Format ten opisuje sekwencje biologiczną za pomocą 4 wierszy: (i) identyfikatora sekwencji, (ii) samej sekwencji nukleotydowej, (iii) separatora oddzielającego wiersze drugi i czwarty oraz (iv) jakości wyznaczonej dla każdego z nukleotydów. Ponieważ sekwencja jest przechowywana w drugim z czterech wierszy, można ją wyodrębnić za pomocą operatora modulo (% czyli reszta z dzielenia liczb całkowitych). NR to numer aktualnie czytanego wiersza. Za pomocą warunku NR%4==2 wyodrębniamy tylko wiersze, których pozostała część po podzieleniu przez 4 (liczba wierszy opisujących pojedynczą sekwencję) jest równa 2 (tj. wiersz zawierający sekwencję). Następnie liczba sekwencji jest zwiększana o 1 (inkrementacja) i ta operacja jest wyrażona poprzez reads++. Sumowane są też długości sekwencji poprzez readlen+=length ($0). Ostatnią rzeczą, jaką robi awk, jest wykonanie poleceń następujących po warunku END, po którym awk wypisuje liczbę sekwencji (print reads) oraz ich średnią (print readlen/reads) w danym pliku FASTQ. Należy pamiętać, że w przypadku plików skompresowanych należy je najpierw zdekompresować lub odczytać za pomocą zcat lub bzcat – w zależności od typu kompresji.

 

Przykład 8. Polecenie head -n 8 wyświetla ośmiem pierwszych wierszy pliku w formacie FASTQC. Polecenie sed zmienia format oryginalnego pliku na plik w formacie FASTA, którego pierwsze osiem linii również zostaje wyświetlonych na ekranie.

 

Sed jest edytorem strumieniowym służącym do przekształcania tekstu. W tym przykładzie (Przykład 8) jest on używany do konwersji pliku FASTQ do FASTA poprzez zredukowanie liczby linii z czterech (identyfikator, sekwencja, separator oraz jakość sekwencji) do dwóch (identyfikator i sekwencja). Zatem, wzorzec 1~4 wybiera pierwszą linię z 4. W linii tej, początkowy znak „@” jest zastępowany przez „>”, ponieważ „s” przed wzorcem oznacza podstawienie. Po średniku, wzorzec 2~4p wskazuje, że należy wypisać co czwartą linię począwszy od drugiej (tj. linii zawierającej sekwencję). Efekt jest zapisywany w nowym pliku o nazwie SRR5078057_1.fasta ( znak „>” wskazuje na przekierowanie wyniku poleceń do pliku zamiast wyświetlenia go na ekranie monitora).

 

Przedstawione przykłady są jedynie namiastką tego, jak można usprawnić pracę z dużymi zbiorami danych biologii molekularnej. Ucząc się tej sztuki warto także sięgnąć nie tylko po podstawowe narzędzia powłoki bash, ale także narzędzia dedykowane biologom i bioinformatykom, takie jak np. bioawk [4]. Bioawk dodatkowo ułatwia przetwarzanie popularnych formatów biologicznych, takich jak BED [5], SAM [6], VCF [7], GFF [8], FASTA [9] oraz FASTQ [10]. Więcej przykładów obróbki danych wraz z przykładowymi plikami zostało przedstawionych w publikacji „Extraordinary Command Line: Basic Data Editing Tools for Biologists Dealing with Sequence Data” [11].

Podsumowując, obsługa linii poleceń bash jest jedną z najbardziej podstawowych umiejętności służącej edycji i manipulacji dużymi zbiorami danych biologicznych. Niekwestionowaną zaletą jest to, że dzięki wykorzystaniu wiersza poleceń możliwe jest szybkie przetwarzanie danych, ponieważ nie musi ono polegać na oprogramowaniu graficznym, które samo w sobie pochłania część czasu działania procesora i zazwyczaj wymaga ciągłej interakcji z użytkownikiem. Zaczynając od wykorzystania podstawowych poleceń, użytkownik może łatwo rozszerzyć umiejętności i zmodyfikować je do specyficznych potrzeb własnej analizy, nawet bez głębokiej wiedzy na temat programowania.

 

Źródła:

[1] www.ncbi.nlm.nih.gov

[2] www.ensembl.org

[3] www.genome.ucsc.edu

[4] www.github.com/lh3/bioawk

[5] Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics 2010; 26(6): 841-2; http://dx.doi.org/10.1093/bioinformatics/btq033.

[6] Li H, Handsaker B, Wysoker A, et al. The sequence alignment/map format and SAMtools. Bioinformatics 2009; 25(16): 2078-9; http://dx.doi.org/10.1093/bioinformatics/btp352.

[7] Danecek P, Auton A, Abecasis G, et al. The variant call format and VCFtools. Bioinformatics 2011; 27(15): 2156-8; http://dx.doi.org/10.1093/bioinformatics/btr330

[8] https://mblab.wustl.edu/GTF22.html

[9] Lipman D J, Pearson W R. Rapid and Sensitive Protein Similarity Searches. Science (80) 1985; 227(4693): 1435-41; http://dx.doi.org/10.1126/science.2983426.

[10] Cock PJ, Fields CJ, et al. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Res. 2010 Apr;38(6):1767-71; doi: 10.1093/nar/gkp1137

[11] Mielczarek M, Czech B, et al. Extraordinary Command Line: Basic Data Editing Tools for Biologists Dealing with Sequence Data. The Open Bioinformatics Journal 2020. 13 137-145; doi: 10.2174/1875036202013010137.

 

Podziel się: