Pobieranie sekwencji z bazy GenBank i praca w plikami FASTA przy pomocy linii komend

Dlaczego linia komend?

Pobieranie sekwencji przy użyciu narzędzi uruchamianych z linii komend jest w niektórych sytuacjach łatwiejsze, zwłaszcza gdy wiemy, które sekwencje (albo ich fragmenty) i/lub chcemy pobrać ich wiele.

Dostęp do GenBank-u przy pomocy narzędzi Entrez Direct

,,Entrez Direct'' to grupa programów opublikowanych przez NCBI umożliwiających dostęp do baz danych tej organizacji. Podręcznik do pakietu (wraz z instrukcją instalacji) pt. Entrez Direct: E-utilities on the UNIX Command Line można znaleźć pod adresem https://www.ncbi.nlm.nih.gov/books/NBK179288/. Nie będziemy oczywiście poznawać dogłębnie wszystkich narzędzi ale pokażę kilka przydatnych komend i sposobów ich wykorzystania.

Wyszukiwanie sekwencji - esearch

Do wyszukiwania informacji w bazach danych służy narzędzie esearch. Przy podstawowym użyciu, przyjmuje dwa argumenty:
    -db - określający bazę do której zostanie wysłane zapytanie
    -query - zapytanie
Wywołanie programu będzie wyglądało tak:
1
esearch -db identyfikator_bazy -query zapytanie
Copied!
Jeśli zapytanie składa się z wielu słów, należy je zamknąć w cudzysłowach.
Lista identyfikatorów wybranych baz danych znajduje się pod adresem https://www.ncbi.nlm.nih.gov/books/NBK25497/table/chapter2.T._entrez_unique_identifiers_ui/?report=objectonly. Jak widać jest ich wiele. Nas będzie interesowała teraz baza ,,Nucleotide'', która ma identyfikator nuccore lub nucleotide. To właśnie z niej wcześniej korzystaliśmy przez stronę internetową.
Warto też przejrzeć wbudowaną pomoc programu, gdzie znajdziemy nie tylko identyfikatory baz ale też opisy innych opcji programu. Uruchamiamy ją wpisując:
1
$: esearch -h
Copied!
Uruchom polecenie:
1
$: esearch -db nucleotide -query "atp6 Orobanche"
Copied!
Pokaże się wynik:
1
<ENTREZ_DIRECT>
2
<Db>nucleotide</Db>
3
<WebEnv>NCID_1_33234224_130.14.18.34_9001_1517069868_100176055_0MetA0_S_MegaStore_F_1</WebEnv>
4
<QueryKey>1</QueryKey>
5
<Count>17</Count>
6
<Step>1</Step>
7
</ENTREZ_DIRECT>
Copied!
Nie wydaje się zbyt przydatny, ale można zrozumieć, że znaleziono 17 wyników. Nie widać żadnych sekwencji. Zaraz pokażę jak uzyskać bardziej użyteczny rezultat.

Pobieranie sekwencji - efetch

Drugim narzędziem pakietu Entrez Direct z którego będziemy korzystać jest efetch służący, jak wskazuje nazwa, do pobierania danych.
Znając numer dostępowy sekwencji można ją pobrać w ten sposób:
1
efetch -db identyfikator_bazy -format format_wyniku -id numer_dostępowy
Copied!
Na przykład:
1
$: efetch -db nucleotide -id KU180461 -format fasta
Copied!
Wynik jest zgodny z oczekiwaniami:
1
>KU180461.1 Orobanche coerulescens clone 12 ATPase subunit 6 (atp6) gene, partial cds; mitochondrial
2
CTGCTAACTCTCAGTTTGGTCCTACTTCTGATTCATTTTGTTACTAAAAAGGGAGGAGGAAACTCAGTAC
3
CAAATGCTTGGCAATCCGTGGTAGAGTTTATTTATGATTTCGTGCTGAACCTGGTAAACGAACAAATAGG
4
GGGTCTTTCCGGAAATGTTAAACAAAAGTTTTTCCCTTGCATCTTGGTCACTTTTACTTTTTTGTTATTT
5
TGTAATCTTCAGGGTATGATACCTTATAGCTTCACAGTTACAAGTCATTTTCTCATTACTTTAGGTCTCT
6
CATTTTCTCTTTTTATTGGCATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTT
7
ATTACCCGCAGGAGTCCCACTGCCATTAGCACCTTTTTTAGTACTCCTTGAGCTAATTTCTTATTGTTTT
8
CGCGCATTAAGCTTAGGAATACGTTTATTTGCTAATATGATGGCCGGTCATAGTTTAGTAAAGATTTTAA
9
GTGGGTTCGCTTGGACTATGCTATGTATGAATGATCTTTTGTATTTTATAGGGGATCTTGGTCCTTTATT
10
TATAGTTCTTGCATTAACCGGTCTTGAATTAGGTGTAGCTATATCACAAGCTCATGTTTCTACGATCTCA
11
ATCTGTATTTAC
Copied!
Wynik oczywiście możemy zapisać w pliku używając znaku >:
1
efetch -db nucleotide -id KU180461 -format fasta > atp6.fasta
Copied!
Przejrzyj pomoc programu, którą wywołasz komendą:
1
efetch -h
Copied!
Poeksperymentuj z różnymi formatami wyniku.
Teraz sprobujemy użyć polecenia efetch do pobrania serii sekwencji, tych samych, które należało pobrać przy pomocy przeglądarki w poprzedniej części. Zapisz ich numery, po jeden w każdej linii w pliku atp6_numery.txt
1
KU180471
2
KU180476
3
FJ595983
4
KC879635
5
HQ593782
6
KU180469
7
KU180474
8
KU180475
9
KU180468
10
KU180466
11
AF095276
Copied!
Teraz napiszemy skrypt o nazwie pobierz_sekwencje.sh, który w pętli będzie kolejno odczytywał numery sekwencji i przekazywał je do polecenia efetch, które będzie je pobierało i zapisywało w pliku. Nazwę pliku z danymi i nazwę pliku wynikowego będziemy podawać jako argumenty przy wywoływaniu skryptu.
1
#!/bin/bash
2
3
# Nazwa pliku z danymi podana jako pierwszy argument
4
dane=$1
5
# Nazwa pliku wynikowego podana jako drugi argument
6
fileout=$2
7
8
# Tworzony jest plik wynikowy, jeśli istnieje zostaje nadpisany
9
echo > $fileout
10
11
# Pętla odczytująca numery sekwencji z pliku $dane,
12
# są one przechowywane w zmiennej nrGB
13
while read nrGB
14
do
15
# Informacje o obieraniu i zapisywaniu sekwencji
16
echo "Pobieram sekwencję $nrGB i zapisuję do pliku $fileout"
17
# Pobieranie i zapisywanie sekwencji
18
efetch -db nucleotide -id $nrGB -format fasta >> $fileout
19
done < $dane
Copied!
Następnie uruchamiamy skrypt:
1
$: ./pobierz_sekwencje.sh atp6_numery.txt atp6_sekwencje.fasta
2
3
Pobieram sekwencję KU180471 i zapisuję do pliku atp6_sekwencje.fasta
4
Pobieram sekwencję KU180476 i zapisuję do pliku atp6_sekwencje.fasta
5
Pobieram sekwencję FJ595983 i zapisuję do pliku atp6_sekwencje.fasta
6
Pobieram sekwencję KC879635 i zapisuję do pliku atp6_sekwencje.fasta
7
Pobieram sekwencję HQ593782 i zapisuję do pliku atp6_sekwencje.fasta
8
Pobieram sekwencję KU180469 i zapisuję do pliku atp6_sekwencje.fasta
9
Pobieram sekwencję KU180474 i zapisuję do pliku atp6_sekwencje.fasta
10
Pobieram sekwencję KU180475 i zapisuję do pliku atp6_sekwencje.fasta
11
Pobieram sekwencję KU180468 i zapisuję do pliku atp6_sekwencje.fasta
12
Pobieram sekwencję KU180466 i zapisuję do pliku atp6_sekwencje.fasta
13
Pobieram sekwencję AF095276 i zapisuję do pliku atp6_sekwencje.fasta
Copied!
Sprawdź zawartość pliku wynikowego.

efetch - pobieranie fragmentów sekwencji

Dzięki programowi efetch można też uzyskać fragment sekwencji podając miejsce pierwszego i ostatniego nukleotydu żądanego fragmentu. Służą do tego opcje: -seq_start i -seq_stop.
Na przykład wytnijmy sekwencję genu nad2 z genomu mitochondrialnego Arabidopsis thaliana o identyfikatorze NC_001284. Miejsca początku i końca genu odczytamy z wyniku wyszukiwania:
Zakres genu
Teraz przekażmy je do polecenia:
1
efetch -db nuccore -id NC_001284 -format fasta -seq_start 134071 -seq_stop 138153
Copied!
Otrzymamy żądany fragment genomu.
Teraz powróćmy do problemu pobieraniu sekwencji we fragmentach z którym zetknęliśmy się w lekcji poświęconej pobieraniu sekwencji z GenBank-u przy pomocy strony internetowej. Sekwencja genu nad2 w genomie Arabidopsis thaliana (nr. NC_001284) występowała tam w dwu, oddzielonych od siebie fragmentach, przy czym część kodująca (CDS) miała aż pięć części.
nad2 u _Arabidopsis thaliana_
W dodatku były to sekwencje komplementarne (i odwrócone). Teraz pokażę jak używając prostego skryptu pobrać żądane fragmenty i połączyć je w jedną sekwencję.
Można oczywiście po kolei wykonywać polecenia pobierania fragmentów z sekwencji, za pomocą narzędzia efetch jak pokazałem powyżej i dołączać wyniki do pliku FASTA.
1
efetch -db nuccore -id NC_001284 -format fasta -seq_start 327890 -seq_stop 328078 > nad2_A_thaliana.fasta
2
efetch -db nuccore -id NC_001284 -format fasta -seq_start 329735 -seq_stop 330306 >> nad2_A_thaliana.fasta
3
efetch -db nuccore -id NC_001284 -format fasta -seq_start 332945 -seq_stop 333105 >> nad2_A_thaliana.fasta
4
efetch -db nuccore -id NC_001284 -format fasta -seq_start 79740 -seq_stop 80132 >> nad2_A_thaliana.fasta
5
efetch -db nuccore -id NC_001284 -format fasta -seq_start 81113 -seq_stop 81297 >> nad2_A_thaliana.fasta
Copied!
Komendy możesz wpisywać i uruchamiać w terminalu, albo umieścić je w skrypcie (np. pobierz_odcinki.sh) i go uruchomić.
Otrzymany plik będzie miał taką zawartość:
1
>NC_001284.2:327890-328078 Arabidopsis thaliana mitochondrion, complete genome
2
TTAAAGATATGAACTGAGTGCCATTTGATGAGTAACTGAGAACAAAGGAGAGGGGTATAGCAAGGATGAA
3
GTAATGAAAAAGGAAGTCATTGCTAGTAGTAACGACTTATTACGATCCATTGGTTCATATAGAATCCATG
4
TCCTAGGTGTATCAAAAAACATTCTTTTCACTAAGCGTATATAATAAAA
5
>NC_001284.2:329735-330306 Arabidopsis thaliana mitochondrion, complete genome
6
ACGACCTATAACGCTAGTCACTACTCCCACTGGGGCTAGAAAGTAAGCCCCACAACCCAAAGCGGCGAAG
7
AACAAATAGAATTTGCTACAAAAGCCGGCTAACGGGGGTATTCCTGCGTATGAGAACATAGTAATGGAGA
8
AGGTAATAGCCGAAATAGGATTCGTTTTGGCTAGAGCGCCCAAATCCGCTATATATTTGACACGGGTTTG
9
CCGTAATGCTGAAACTATGGCGAATGCATCCATCGTCATTAATGCATAAATAAAGATACCAATTAGTAGT
10
GATTGAATTCCTTCTATGGTTCCACATGAGAAACCAGTACGAATATAACCTACATGTCCAATTGAACTAT
11
GAGCTAGAGGTCTTTTGACTTTCGTTTGGGCCATGGCGGCCAGTGCTCCTAAGATCATAGAAGCAATGCT
12
GCAGAAAAAGAAGATTTGTTGCAATGTAGCTCCATAGGAACCATAAATAGAAACACGTAAAATATTAGCA
13
GAAATAGAGATTTTAGGCGCAATAGAAAGGAATGCTGTAACCGGGGTGGGTGAACCCTCATAGATATCTG
14
GTGCCCACATAT
15
>NC_001284.2:332945-333105 Arabidopsis thaliana mitochondrion, complete genome
16
GAAAAGGAACTGCAGTGATCTTGAATAGGAATCCTACAGCGATAGACAGAATCCCCATAAAAATACCACT
17
AGATCGAGCACCAGTGATTTCGTATCCGGTCAAAATCTTGGCTAATTGATCGAAGTGGGTAGCTCCAGTA
18
GACCCATAGATCATGGAACAA
19
>NC_001284.2:79740-80132 Arabidopsis thaliana mitochondrion, complete genome
20
CCAAACAATAATATTCCAGAGGAAAATGCACCTAAGATCAAATATTTCGAGCCGGCTTCCGTGGAAAATT
21
CAGACTTTCTTTTTGATGCTGCGATTACATAAAAACATAAACTTTGAGGCTCAATAGCTAAATACATGGC
22
AATTAAATCATGAGCCGAGATCATAAAGAGCATACCGCGAGTAGGAAGTGGAATTAATACAATGAATTCA
23
AAAGCATCAAACCTCTCTTGGTCGGAAGAATCGAAACACATCGAAATGGTACCAGCCGTACTTAATAATA
24
GAAAGATTTGGCAGAAATATGTAAAATTGTCCCTCCTAAAAAGATTATTCCAGAATAAATGGGCAATAGT
25
TAGGAGAGGTGCGCCAGCGGCGAGCAGAAGCAAGGTTATTAGA
26
>NC_001284.2:81113-81297 Arabidopsis thaliana mitochondrion, complete genome
27
ACACTAAGTAATCCAAGCCAACCCACATTACTGGCTAACGGCGGATAATCATATTTCTTAGAGGTACTAA
28
ATACAACTCCATGAATGAGCAAAATGGAGGTTGCATTAATGATAAAGATCTCTGGGGAAACCGCTAAAAA
29
AAGATTGAACATGTGTGGGAGGATCCGAACGAATTCTGCTTTCAT
Copied!
Zamiast uruchamiać wielokrotnie tą samą komendę, można utworzyć plik tekstowy z miejscami początku i końca fragmentów sekwencji a następnie stworzyć skrypt, który w pętli po kolei je pobierze i umieści w pliku.
Najpierw utworzymy plik, w którym zapiszemy kolejne miejsca początku i końca fragmentów sekwencji oddzielone znakiem tabulatora (odcinki.tsv).
1
327890 328078
2
329735 330306
3
332945 333105
4
79740 80132
5
81113 81297
Copied!
Teraz skrypt:
1
#!/bin/bash
2
# Plik z danymi podany jako argument 1 przy uruchamianiu skryptu
3
data=$1
4
# Plik wyjściowy podany jako argument 2 przy uruchamianiu skryptu
5
fileout=$2
6
# Numer sekwencji w GenBank-u podany jako argument 3
7
nrGB=$3
8
9
# Tworzymy plik wejściowy, jeśli już istnieje to jego zawartość zostaje usunięta
10
echo > $fileout
11
12
# Pętla. Wartości start i stop odczytywane są z pliku $data
13
while read start stop
14
do
15
# Skrypt informuje, który fragment pobiera
16
echo "Pobieram fragment: $start - $stop z sekwencji $nrGB i zapisuję do pliku $fileout"
17
# Pobranie fragmentu i zapisanie do pliku
18
efetch -db nuccore -id $nrGB -format fasta -seq_start $start -seq_stop $stop >> $fileout
19
done < $data
Copied!
1
$: ./pobierz_odcinki_petla.sh odcinki.tsv nad2_A_thaliana.fasta NC_001284
2
3
Pobieram fragment: 327890 - 328078 z sekwencji NC_001284 i zapisuję do pliku nad2_A_thaliana.fasta
4
Pobieram fragment: 329735 - 330306 z sekwencji NC_001284 i zapisuję do pliku nad2_A_thaliana.fasta
5
Pobieram fragment: 332945 - 333105 z sekwencji NC_001284 i zapisuję do pliku nad2_A_thaliana.fasta
6
Pobieram fragment: 79740 - 80132 z sekwencji NC_001284 i zapisuję do pliku nad2_A_thaliana.fasta
7
Pobieram fragment: 81113 - 81297 z sekwencji NC_001284 i zapisuję do pliku nad2_A_thaliana.fasta
Copied!
Sprawdź plik wynikowy.

Łączymy efetch, sed, tr i rev

Jak widać poszczególne fragmenty występują oddzielnie. Teraz powinniśmy połączyć je w jedną sekwencję. Można to zrobić w edytorze tekstu ale od czego mamy linię komend i narzędzia sed, tr i rev?
Zadanie wykonamy w pięciu etapach:
    1.
    Usuniemy linie z opisami sekwencji
    2.
    Połączymy fragmenty w jedną sekwencję
    3.
    Przekształcimy ją w sekwencję komplementarną
    4.
    Odwrócimy sekwencje
    5.
    Dodamy opis w pierwszej linii pliku
Zastanów się, najpierw ogólnie, jak te zadania zrealizować używając programów sed, tr i rev?
    1.
    Usunąć linie zaczynające się od znaku >: tr lub sed
    2.
    Usunąć wszystkie znaki końca linii: \n: tr
    3.
    Zmienić na sekwencję komplementarną: tr
    4.
    Odwrócić sekwencję: rev
    5.
    Wprowadzić do pierwszej linii pliku opis sekwencji
Teraz przekształćmy pomysł w konkretny kod. Utwórz skrypt polacz_odcinki.sh:
1
#!/bin/bash
2
3
# Nazwa pliku z sekwencjami
4
surowy=nad2_A_thaliana.fasta
5
# Plik wyjściowy
6
file=polaczony_$surowy
7
8
# Kopiujemy plik z sekwencjami do pliku wyjściowego,
9
# który będziemy modyfikować
10
cp $surowy $file
11
12
# Usuwanie linii z opisami, czyli linii zaczynających się od znaku >
13
sed -i '/^>/d' $file
14
15
# Łączymy wszystkie linie w jedną, usuwając oznaczenia końca linii
16
# Wynik umieszczamy w pliku tymczasowym
17
cat $file | tr -d "\n" > $file.tmp
18
19
# Zmieniamy na sekwencję komplementarną, wynik umieszczamy
20
# w drugim pliku wyjściowym
21
22
cat $file.tmp | tr "ACGT" "TGCA" > $file.tmp2
23
# Odwracamy sekwencję, wynik umieszczamy w pliku wyjściowym
24
# zmieniając jego zawartość
25
26
cat $file.tmp2 | rev > $file
27
28
# Wstawiamy pierwszą linię z opisem i znak nowej linii
29
sed -i '1s/^/>NC_001284_Arabidopsis_thaliana_nad2\n/' $file
30
31
# Usuwamy pliki tymczasowe
32
rm $file.tmp*
Copied!
Można też skonstruować ,,jednolinijkowca'':
1
cat nad2_A_thaliana.fasta | sed '/^>/d' | tr -d "\n" | tr "ACGT" "TGCA" | rev | sed '1s/^/>NC_001284_Arabidopsis_thaliana_nad2\n/' > polaczony_RC_nad2_A_thaliana.fasta
Copied!

Łączymy esearch i efetch

Dotychczas używaliśmy dwu narzędzi: służący do wyszukiwania esearch i pobierający sekwencje efetch oddzielnie. Ale możemy je razem połączyć za pomocą potoku tak aby wyszukiwanie połączyć z pobieraniem.
1
esearch -db nucleotide -query "atp6[All Fields] AND Orobanche[Organism]" | efetch -format fasta > atp6_Orobanche.fasta
Copied!

samtools faidx- pobieranie sekwencji z pliku FASTA

Pobierz plik http://ggoralski.pl/files/filogenetyka-data/Orobanchaceae-trnL-trnF-aligned.fasta
Sprawdź jego zawartość. Używaliśmy go już we wcześniejszych ćwiczeniach, zawiera wyrównane sewkencje trnL-trnF.
Teraz sprawdź nazwy sekwencji używając komendy grep:
1
$: grep ">" Orobanchaceae-trnL-trnF-aligned.fasta
2
3
>KY484464_O._teucrii
4
>KY484493_O._flava
5
>KY484489_O._mayeri
6
>KY484471_O._kochii
7
>KY484474_O._elatior
8
>KU238865_O._coerulescens
9
>KY484502_P._ramosa
10
>KY484503_P._purpurea
11
>KX524675_Lindenbergia_siniaca
Copied!
Teraz będziemy chcieli pobrać z pliku sekwencję o nazwie KX524675_Lindenbergia_siniaca. Można to zrobić oczywiście używając edytora tekstu lub programu do pracy z plikami FASTA, ale czasem wygodniej to zrobić używając linii komend. Użyjemy do tego zadania programu faidx, który jest częścią pakietu samtools.
Pakiet samtools zawiera wiele przydatnych programów do pracy z sekwencjami, przede wszystkim w innych formatach niż FASTA, ale faidx jest dość przydatne w codziennej pracy także w tym formacie.
Strona domowa projektu znajduje się pod adresem: http://www.htslib.org/. Tam też można przeczytać manuale programów.
Na Debianie można zainstalować pakiet komendą:
1
sudo apt-get install samtools
Copied!
W innym przypadku, można pobrać pliki źródłowe ze strony i skompilować wg. znajdującej się tam instrukcji.
Uruchamianie narzędzi z pakietu samtools wygląda następująco:
1
samtools nazwa_narzedzia opcje
Copied!
W przypadku faidx będzie to wyglądało tak:
1
samtools faidx opcje
Copied!
Opcjami, które użyjemy będzie kolejno nazwa pliku, z którego pobierzemy sekwencję i nazwa sekwencji. Uruchom:
1
samtools faidx Orobanchaceae-trnL-trnF-aligned.fasta KX524675_Lindenbergia_siniaca
Copied!
Na ekranie wyświetli się żądana sekwencja, oczywiście można ją zapisać w pliku używając >.
Zauważ, że podajemy pełną nazwę sekwencji (bez znaku >), wpisanie jej fragmentu nie zadziała. Może to być pewnym mankamentem, na przykład gdy nazwy są długie albo gdy chcemy pobrać wiele sekwencji, których nazwy mają jakąś część wspólną (np. nazwę rodzaju.). W taki przypadku pomocne będzie polecenie grep.
Uruchom:
1
$: grep "O." Orobanchaceae-trnL-trnF-aligned.fasta
2
3
>KY484464_O._teucrii
4
>KY484493_O._flava
5
>KY484489_O._mayeri
6
>KY484471_O._kochii
7
>KY484474_O._elatior
8
>KU238865_O._coerulescens
Copied!
W ten sposób otrzymujemy pełne nazwy wszystkich sekwencji należących do Orobanche. Teraz powinniśmy je przekazać do polecenia faidx, ale najpierw musimy pozbyć się znaków >:
1
$: grep "O." Orobanchaceae-trnL-trnF-aligned.fasta | sed 's/>//'
2
3
KY484464_O._teucrii
4
KY484493_O._flava
5
KY484489_O._mayeri
6
KY484471_O._kochii
7
KY484474_O._elatior
8
KU238865_O._coerulescens
Copied!
Czas użyć ich przy wywołaniu faidx. Polecenie:
1
grep "O." Orobanchaceae-trnL-trnF-aligned.fasta | sed 's/>//' | samtools faidx Orobanchaceae-trnL-trnF-aligned.fasta
Copied!
nie zadziała. Trzeba wykorzystać xargs:
1
grep "O." Orobanchaceae-trnL-trnF-aligned.fasta | sed 's/>//' | xargs samtools faidx Orobanchaceae-trnL-trnF-aligned.fasta > Orobanche-trnL-trnF.fasta
Copied!
Sprawdź zawartość pliku wynikowego.
Przy okazji pobierania sekwencji, można je oczywiście zmodyfikować. Na przykład możemy zmienić skrót O. na pełną nazwę rodzajową Orobanche i usunąć wszystkie znaki - oznaczające brak nukleotydu. Użyjemy w tym celu poecenia sed. Przy usuwaniu - pamiętaj o dodaniu g, który pozwala na wykonaniu polecenia na wszystkich znakach w linii.:
1
grep "O." Orobanchaceae-trnL-trnF-aligned.fasta | sed 's/>//' | xargs samtools faidx Orobanchaceae-trnL-trnF-aligned.fasta | sed 's/O\./Orobanche/' | sed 's/-//g' | > Orobanche-trnL-trnF.fasta
Copied!
Last modified 1yr ago