Python jest językiem ogólnego przeznaczenia, za pomocą którego możemy tworzyć programy, analizować dane, tworzyć i analizować bazy danych, tworzyć strony www i aplikacie webowe. Python ma też zastosowanie w tworzeniu gier wideo, machinne learning oraz big data. Ponieważ jest to jeden z najpopularniejszych języków programowania i jest to projekt typu open source, wielu programistów na całym świecie tworzy kod pyhona. Ten kod jest często udostępniany jako dodatkowe biblioteki. Biblioteki umożliwiają zwiększenie możliwości i potencjału języka Python, pozwalając na przyspieszenie i automatyzację wielu procesów. Dzięki temu coraz więcej osób korzysta z Pytona, a to przyciąga kolejnych developerów, którzy tworzą nowe biblioteki dla tego języka. To z kolei poszerza krąg jego użytkowników i tak następuje sprzężenie zwrotne, które powoduje, że Python staje się coraz popularniejszym językiem programowania.
Jedną z bibliotek do Pythona jest biblioteka Biopython, zawiera ona zbiór funkcji, które są bardzo przydatne do analizy danych sekwencyjnych. Za pomocą funkcji z Bioblioteki Biopython można między innymi:
Więcej szczegółowych informacji na temat biblioteki Biopython można znaleźć na stronie projektu. Na dzisiejszych ćwiczeniach dowiemy się, jak można pobrać sekwencje genomowe z bazy NCBI. Jak wykonać analizy blast w pythonie oraz w jaki sposób możemy analizować wyniki blast za pomocą bibliotecki python.
# import potrzebnych bibliotek i funkcji
import Bio
import pandas as pd
import numpy as np
from Bio.Blast.Applications import NcbiblastnCommandline
from Bio.Blast import NCBIXML
from Bio import SearchIO
from Bio import Entrez
W celu pobrania genomu wykorzystamy skrypt get_genome.sh, który wykorzystywaliśmy w poprzednim semestrze. Skrypt ten pozwala na pobranie genomu po numerze Biosample odpowiadający danemu genomowi. W jaki sposób możemy uzyskać informacje o biosample? Wykorzystamy do tego biobliotekę Biopython i funkcję Entrez.
# korzystając z entrez należy każdorazowo podać swoje dane
Entrez.email = "bartosz.kozak@upwr.edu.pl"
Najpierw wykorzystamy narzędzie esearch i przeszukamy bazę danych w celu odnalezienia interesujących nas danych. Bazą, którą będziemy przeszukiwać będzie baza "genomów", a zapytaniem będzie nazwa organizmu, dla którego chcemy pobrać informacje genomowe. W naszym przykładzie będzie to łubin wąskolistny.
handle = Entrez.esearch(db="genome", term="Lupinus angustifolius")
record = Entrez.read(handle)
handle.close()
record
lupin_id = record['IdList'][0]
lupin_id
Uzyskujemy informacje o numerze Id dla łubinu, następnie możemy wykorzystać tą informację do połączenia jej za pomocą polecenia elink z bazą sekwencji.
lupin_id
Entrez.email = "bartosz.kozak@upwr.edu.pl"
handle = Entrez.elink(dbfrom="genome", id=lupin_id,
linkname="genome_nuccore")
record = Entrez.read(handle)
sequence_id = record[0]['LinkSetDb'][0]['Link'][0]['Id']
sequence_id
Możemy uzyskać informację o sekwencji jednego z chromosomów łubinu wąskolistnego. Dzięki tej informacji będziemy mogli pobrać record odpowiadający id tej sekwencji, a z recordu wyciągnąć informację o id Biopróbki. Ta informacja jest nam potrzebna do użycia w skrypcie get_genome.sh
sequence_id
handle = Entrez.efetch(db="nucleotide", id=sequence_id,
rettype="gb", retmode="xml")
record = Entrez.read(handle)
biosample = record[0]['GBSeq_xrefs'][1]['GBXref_id']
biosample
Ćwiczenie 1:
Proszę zapisać numer biopróbki "biosample" do pliku tekstowego "bio_ids.txt"
Po utworzeniu pliku "bio_ids.txt" możemy pobrać ze strony internetowej skrypt get_genome.sh Następnie uruchamiamy skrypt podając jako parametr wejściowy nazwę pliku "bio_ids.txt". Skrypt automatycznie pobierze genom ze strony RefSeq.
Przed przystąpieniem do analiz Blast musimy utworzyć lokalną bazę danych dla programu Blast. W tym celu tworzymy katalog ~/Blast_DB. Następnie przenosimy do tego katalogu pobrany plik genomu. Jeżeli plik jest spakowany (rozszerzenie ".gz") musimy go rozpakować poleceniem gzip -d nazwa_pliku. Następnie wykorzystując program "makeblastdb" tworzymy bazę danych. Polecenie do utworzenia bazy wygląda następująco:
makeblastdb -in nazwa_pliku_genomu -pare_seqids -dbtype nuclĆwiczenie 2:
Proszę utworzyć w katalogu ~/Blast_DB bazę danych dla pobranego genomu łubinu wąskolistnego.
Po utworzeniu bazy danych możemy przejść do Pythona i wykonać blast. W przykładzie będziemy przeprowadzać blast sekwencji markerów Blast markerów DArT_SNP do genomu referencyjnego Łubinu wąskolistnego.
cl = NcbiblastnCommandline(query="/Dydaktyka/DArT.fasta",
db = "/home/bartek/Cw_genomy/genomy/GCA_002285895.2_ASM228589v2_genomic.fna",
out="test.tab",
outfmt=6)
# wyświetlamy komendę Blast
print(cl)
# wykonujemy komendę Blast
stdout, stderr = cl()
Wyniki Blast zostały zachowane jako plik formatu xml. Aby możliwa była łatwa analiza wyników przerobimy plik w formacie xml na plik tekstowy rozdzielony znakami tabulacji.
Plik w formacie tekstowym możemy wczytać do pythona jako obiekt typu pandas i wykonać jego dalszą analizę.
wyniki = pd.read_csv('test.tab', sep='\t',
header=None)
wyniki.columns = ['CloneId', 'NLL_PseChr', 'Ident', 'Length subiect', 'Diff', 'F', 'G',
'Length query', 'Start', 'Stop', 'Expect', 'Score']
wyniki.head()
len(wyniki)
Wyniki analizy blast pozwoliły na identyfikację prawie 7000 fragnentów genomu o zbliżonej sekwencji do markerów DArTseq. Jednakże w naszych analizach chemy wybrać tylko fragmenty najbardziej podobne sekwencyjne do markerów DArTseq. Markery te zostały stworzone dla łubinu wąskolistnego i mają długość 69 pz. Są to markery tupu SNP, dlatego dopuszczalna jest różnica w obrębie jednego nukleotydu między sekwencją markera, a sekwencją genomu referencyjnego. Wykorzystamy bibliotekę python do przefiltrowania wyników blast.
wyniki = wyniki[wyniki['Ident']>=98.5]
wyniki = wyniki[wyniki['F']==0]
wyniki =
wyniki[wyniki['Length subiect']==69].reset_index().drop('index',1)
len(wyniki)
Filtrowanie pozwoliło nam na ograniczenie wyników Blast do 2868 sekwencji. Oznacza to że okolo 60% sekwencji markerów zostało odnalezionych w genomie referencyjnym. Mając na uwadze, że genom referencyjny jest obecnie kompletny jedynie w około 50%, ten wynik nie pownienie dziwić. Wyniki blast możemy zapisać jako plik w formacie csv.
wyniki.to_csv('wyniki.csv', index=False)
Plikiki GFF zawierają informacje o lokalizacji poszczególnych genów w genomie. Jeżeli chcemy porównać sekwencje oraz lokalizację genomową jakiegoś genu (genów) pomiędzy dwoma organizmami, pliki GFF są bardzo użytecznym narzędziem. Pliki GFF możemy pobrać ze strony RefSeq, podobnie jak pobieraliśmy pliki z sekwencją genomu. Możemy w tym celu wykorzystać skrypt get_genome.sh. Musi on zostać jednak zmodyfikowany. Pliki GFF mają rozszerzenie 'gff', a nie fan.
Ćwiczenie 3:
Celem tego zadania jest pobranie pliku w formacie gff z bazy RefSeq dla Medicago truncatula. W tym celu należy pobrać informację o biosample dla medicago, a następnie zapisać ją do pliku tekstowego (np. medicago.txt). Następnie zmodyfikować skrypt get_genome.sh tak aby pobierał pliki z rozszerzeniem ".gff". Następnie należy uruchomić skypt i pobrać plik gff.
# wczytanie pliku GFF, selekcja mRNA
Medicago_GFF = pd.read_csv('/home/bartek/Cw_genomy/genomy/GCA_003473485.2_MtrunA17r5.0-ANR_genomic.gff', header=None, sep='\t|;', engine='python', skiprows=7)
Medicago_GFF.head()
Ćwiczenie 4:
Celem tego zadania jest pobranie pliku z sekencjami kodującymi CDS dla Medicago truncatula. Należy postępować analogicznie jak w przypadku ćwiczenia 3. Modyfikacja pliku get_genomu w tym przypadku polega na zmianie polecenia wget aby pobrane zostały pliki z "_cds_from_genomic.fna.gz"