Analiza bioinformatyczna jest nierozerwalnym elementem całego procesu badania sekwencji DNA. Praca bioinformatyków jest więc kluczowa w realizacji usługi sekwencjonowania. Aby efektywnie przeprowadzić analizę bioinformatyczną potrzebne są odpowiednie narzędzia i rozwiązania. Jednym z takich rozwiązań jest potokowe przetwarzanie danych.  

Czym jest przetwarzanie potokowe?

Potoki przetwarzania danych są sposobem na sekwencyjne przetwarzanie informacji. Zadania takie jak sortowanie, czy wyszukiwanie duplikatów , są podzielone na osobne bloki, które łączą się ze sobą szeregowo. Właśnie stąd wiemy jaka jest prawidłowa kolejność wykonania wszystkich zadań.

Do czego służy przetwarzanie potokowe?

Potokowe przetwarzanie danych w codziennej pracy bioinformatyka wykorzystywane jest przede wszystkim po to, aby zaoszczędzić czas i zwiększyć efektywność pracy. Przykładowo, jeśli chcemy wyświetlić tylko te linie z pliku, które rozpoczynają się od “>”, a następnie posortować je, nie musimy stosować komendy wyszukania odpowiednich linii i zapisywać je do osobnego pliku, a dopiero później sortować ten plik. Zamiast tego możemy jednorazowo wyszukać określone linie i skierować je potokiem bezpośrednio do komendy sortującej.

Jak używać potoków w konsoli?

Potok danych w systemach opartych o UNIX przekierowujemy w konsoli za pomocą znaku “|”. Jest on łącznikiem między dwoma blokami komend. Zamiast tworzyć plik pośredni:

grep "^>" file_in.txt > file_tmp.txt
sort file_tmp.txt > file_out.txt

możemy przekazać dane potokiem:

grep "^>" file_in.txt | sort > file_out.txt

Tym sposobem nasze polecenie jest jednolinijkową komendą, z której jasno wynika przepływ danych, bez tworzenia pośrednich plików, o których należy pamiętać, aby je usunąć, gdy nie będą już potrzebne.

Komenda łączona za pomocą potoków może być dowolnie długa:

command1 | command2 | command3 | command3 | … | commandX

Zalety potokowego przetwarzania danych

• Przy przetwarzaniu potokowym dane są przechowywane w szybkiej pamięci RAM, dzięki temu unikamy tworzenia plików pośrednich – ich zapis na dysku i ponowny odczyt przez inny program powoduje spowolnienie przetwarzania ze względu na konieczność użycia dysków twardych, szczególnie przy przetwarzaniu dużych wolumenów danych. 
• Zdecydowanie bardziej wygodne i szybsze jest sporządzanie całych skryptów i poleceń. Zajmują one często tylko jedną linię kodu.
• Takie rozwiązanie daje również większą kontrolę nad tworzonymi plikami i przepływem danych, a także sprawia wrażenie czystego i przejrzystego kod.

Oczywiście zdarzają się również sytuacje, w których potok danych nie jest możliwy. Takim przykładem może być konieczność obróbki danych przez program, który kilkukrotnie wczytuje dane wejściowe. W trakcie przetwarzania potokowego dane przekazywane są do programu tylko jeden raz, w czasie rzeczywistym.

Przykład wykorzystania przetwarzania potokowego

Jako przykład przedstawiamy analizę danych o taksonomii wpisów w bazie danych sekwencji genomowych NCBI. Z pliku https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxcat_readme.txt możemy wyczytać, iż w bazie przechowywane są genomy organizmów przynależnych do następujących kategorii/królestw:

  A = Archaea
  B = Bacteria
  E = Eukaryota
  V = Viruses and Viroids
  U = Unclassified
  O = Other

Pierwsze spojrzenie na dane: ściągamy plik z taksonomicznymi ID, rozpakowujemy go i wybieramy 5 losowych wierszy spośród pierwszych 5000.

URL="https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxcat.zip"
curl -s "$URL" | zcat | head -n 5000 | shuf | head -n 5
E       3095    3095
E       8737    8737
B       1066    1066
B       447     66973
E       8420    8420

Dane mają postać Kategoria SpeciesTaxonomicID TaxonomicID. Jeżeli SpeciesTaxonomicID jest równy TaxonomicID to wpis dotyczy gatunku, w innym przypadku wpis dotyczy organizmu będącego podgatunkiem. Przykładowo, wiersz B 447 66973 dotyczy bakterii Legionella bozemanii serogroup 1 podgatunku bakterii gatunku Fluoribacter bozemanae.

Zliczymy teraz ile organizmów odnotowanych jest taksonomicznie w każdej z kategorii/w każdym królestwie. Podobnie jak poprzednio ściągamy plik, rozpakowujemy, przepisujemy tylko pierwszą kolumnę, sortujemy, zliczamy ile jest kolejnych takich samych wierszy i ostatecznie sortujemy w odwrotnej kolejności by wypisać począwszy od kategorii o największej liczbie organizmów. Wynik przekierowujemy jednocześnie na ekran i do pliku komendą tee.

curl -s "$URL" | zcat \
    | awk '{print $1}' | sort -T /dev/shm | uniq -c | sort -rn \
    | tee all_stats.txt
1338805 E
 495917 B
 203340 V
  15530 O
  12701 A
    958 U

Opcja -T /dev/shm na systemach Linux powoduje, iż przy sortowaniu dużych plików powstające pliki tymczasowe umieszczane są w pamięci RAM – daje to mniejsze obciążenie dysku i szybsze sortowanie kosztem zużycia RAM.

Znane są już dane ilościowe gatunków wraz z podgatunkami, w następnym kroku chcielibyśmy zliczyć ile gatunków występuje w każdej kategorii. W tym celu możemy zastosować filtr grep -e '\s([0-9]+)\s+\1', pozostawi on jedynie te wiersze, dla których druga i trzecia kolumna są identyczne.

Dodatkowo, w poniższej komendzie używamy tee (duplikacja strumienia) oraz przekierowania bash strumienia do kolejnego polecenia >(polecenie_potokowe|...) tak, aby wykonać dwa obliczenia: poprzednie oraz kolejne z filtrem jedynie dla gatunków. W ten sposób przy dużych plikach możemy wykonywać kilka obliczeń jednocześnie odczytując duży plik jedynie raz.

curl -s "$URL" | zcat \
    | tee >( awk '{print $1}' | sort -T /dev/shm | uniq -c | sort -rn \
           | tee all_stats.txt >&2 ) \
    | grep -e '\s\([0-9]\+\)\s\+\1' \
    | awk '{print $1}' | sort -T /dev/shm | uniq -c | sort -rn \
    | tee species_stats.txt
1338805 E
 495917 B
 203340 V
  15530 O
  12701 A
    958 U
1302432 E
 451723 B
  34725 V
  15530 O
  12341 A
    958 U

Wyniki zostały także zapisane w plikach all_stats.txt species_stats.txt. Możemy zaobserwować, iż największa różnica w liczbie odnotowanych podgatunków względem liczby odnotowanych gatunków jest w domenie wirusów. Przydatna jest tu komenda tail, która wypisuje zawartość plików wraz z ich nazwami.

tail -n +1 *stats.txt

==> all_stats.txt <==
1338805 E
495917 B
203340 V
 15530 O
 12701 A
   958 U

==> species_stats.txt <==
1302432 E
451723 B
 34725 V
 15530 O
 12341 A
   958 U