09 - Biopython
Biopython jak wskazuje nazwa jest zestawem narzędzi do obliczeń biologicznych w języku Python.
Strona domowa, zawierająca także dokumentację, znajduje się pod adresem https://biopython.org/.
Dostarczane moduły pozwalają m.in. na:
- obsługę plików z danymi zapisanymi w wielu formatach stosowanych w bioinformatyce, takich jak
FASTA,GenBank, pliki wynikoweBLASTitd. - łączenie się i pobieranie informacji z baz danych, m. in. NCBI (GenBank, BLAST)
- podstawową pracę z sekwencjami
- korzystanie z narzędzi
EMBOSS - itd…
W środowisku Anaconda można zainstalować Biopythona poleceniem:
conda install biopython
Inną możliwością jest instalacja przy pomocy pip:
pip3 install biopython
Teraz sprawdźmy, czy instalacja przebiegła pomyślnie:
import Bio
Bio.__version__
'1.78'
Sekwencje - podstawy
Bioinformatyka opiera się w dużej mierze na pracy z sekwencjami DNA, RNA i białek. Poprzednio, kiedy używaliśmy ich podczas nauki programowania, zapisywaliśmy je po prostu jako ciągi znaków:
sekwencja = 'AGGCATCGGTAGACCCAG'
print(type(sekwencja))
<class 'str'>
Biopython pozwala na tworzenie obiektów typu Seq, które umożliwiają działania specyficzne dla sekwencji nukleotydowych czy aminokwasowych.
Najpierw utwórzmy taki obiekt, używając podanej sekwencji znaków:
from Bio.Seq import Seq
sekw_DNA = Seq('AGGCATCGGTAGACC---CAG')
print(f'Sekwencja: {sekw_DNA}')
print(f'Typ: {type(sekw_DNA)}')
Sekwencja: AGGCATCGGTAGACC---CAG
Typ: <class 'Bio.Seq.Seq'>
Teraz sprawdźmy niektóre możliwości specyficzne dla sekwencji biologicznych:
print(f'Sekwencja: {sekw_DNA}')
print(f'Sekw. komplementarna: {sekw_DNA.complement()}')
print(f'Sekw. odwrócona komplementarna: {sekw_DNA.reverse_complement()}')
print(f'Sekw. bez "-": {sekw_DNA.ungap()}')
# Zmiana na RNA
sekw_RNA = sekw_DNA.transcribe()
print(f'Sekw. RNA: {sekw_RNA}')
print(f'Sekw. komplementarna RNA: {sekw_RNA.complement_rna()}')
print(f'Sekw. aminokwasów z DNA: {sekw_DNA.translate()}')
print(f'Sekw. aminokwasów z RNA: {sekw_RNA.translate()}')
print(f'Sekw. aminokwasów z DNA (STOP): {sekw_DNA.translate(to_stop=True)}')
print(f'Sekw. aminokwasów z RNA (STOP): {sekw_RNA.translate(to_stop=True)}')
Sekwencja: AGGCATCGGTAGACC---CAG
Sekw. komplementarna: TCCGTAGCCATCTGG---GTC
Sekw. odwrócona komplementarna: CTG---GGTCTACCGATGCCT
Sekw. bez "-": AGGCATCGGTAGACCCAG
Sekw. RNA: AGGCAUCGGUAGACC---CAG
Sekw. komplementarna RNA: UCCGUAGCCAUCUGG---GUC
Sekw. aminokwasów z DNA: RHR*T-Q
Sekw. aminokwasów z RNA: RHR*T-Q
Sekw. aminokwasów z DNA (STOP): RHR
Sekw. aminokwasów z RNA (STOP): RHR
Zauważ że:
- po transkrypcji nie otrzymujemy nici komplementarnej, ponieważ założenie jest takie, że wyjściowo mamy nić DNA kodującą.
- metoda
translate()zwróciła tę samą sekwencję aminokwasów z sekwencji DNA i RNA - argument
to_stop=Truepowoduje zatrzymanie ,,translacji’’ po natrafieniu na kodonSTOP
Na obiektach Seq możemy także wykonywać pewne czynności charakterystyczne dla łańcuchów znaków:
print(f'Sekwencja: {sekw_DNA}')
print(f'Długość: {len(sekw_DNA)}')
print(f'Liczba "-": {sekw_DNA.count("-")}')
print(f'Pierwsza zasada: {sekw_DNA[0]}')
print(f'Ostatnia zasada: {sekw_DNA[-1]}')
print(f'Pierwsze 3 zasady: {sekw_DNA[0:3]}')
print(f'Małe litery: {sekw_DNA.lower()}')
# Zmiana na łańcuch znaków
sekw_str = str(sekw_DNA)
print(f'Typ sekw_str: {type(sekw_str)}')
Sekwencja: AGGCATCGGTAGACC---CAG
Długość: 21
Liczba "-": 3
Pierwsza zasada: A
Ostatnia zasada: G
Pierwsze 3 zasady: AGG
Małe litery: aggcatcggtagacc---cag
Typ sekw_str: <class 'str'>
Sekwencje możemy też porównywać ze sobą
sekw_1 = Seq('AGCAA')
sekw_2 = Seq('GGCAT')
sekw_3 = Seq('AGCAA')
print(f'{sekw_1} == {sekw_2}: {sekw_1 == sekw_2}')
print(f'{sekw_1} == {sekw_3}: {sekw_1 == sekw_3}')
AGCAA == GGCAT: False
AGCAA == AGCAA: True
Warto wspomnieć, że listę wszystkich metod i atrybutów dostępnych dla obiektów danego typu, można uzyskać funkcją dir() (oczywiście warto też czytać właściwą dokumentację):
dir(Seq)
['__add__',
'__class__',
...
'back_transcribe',
'complement',
'complement_rna',
'count',
'count_overlap',
'encode',
'endswith',
'find',
'index',
'join',
'lower',
'lstrip',
'reverse_complement',
'reverse_complement_rna',
'rfind',
'rindex',
'rsplit',
'rstrip',
'split',
'startswith',
'strip',
'tomutable',
'transcribe',
'translate',
'ungap',
'upper']
Pliki w formacie FASTA i GB
Pobierz plik sekwencje_atp6.fasta do katalogu, w którym uruchamiasz kod.
Wczytajmy plik, w tym celu użyjemy metody parse(), podając jako argumenty nazwę pliku oraz format danych.
from Bio import SeqIO
sekwencje = SeqIO.parse('sekwencje_atp6.fasta', 'fasta')
Sprawdźmy jakiego typu jest obiekt sekwencje:
print(type(sekwencje))
<class 'Bio.SeqIO.FastaIO.FastaIterator'>
Jest to iterator (z iteratorem spotkaliśmy się już w lekcji poświęconej m.in. wyrażeniom regularnym), zatem wykorzystamy go w pętli, sprawdzimy jakie obiekty zawiera:
for s in sekwencje:
print(f'Typ: {type(s)}')
print(s)
print('---------------------------\n')
Typ: <class 'Bio.SeqRecord.SeqRecord'>
ID: AY007817_Daucus_carota
Name: AY007817_Daucus_carota
Description: AY007817_Daucus_carota
Number of features: 0
Seq('ACTTAACTGCTCGAAAGGCAGAGGCTGAGGCAACAGCTGACCACGGTTCTTTTC...AAT')
---------------------------
...
Typ: <class 'Bio.SeqRecord.SeqRecord'>
ID: KX524674_Lindenbergia_siniaca
Name: KX524674_Lindenbergia_siniaca
Description: KX524674_Lindenbergia_siniaca
Number of features: 0
Seq('CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTGTTACTAAAAAGGGA...TAC')
---------------------------
Jak widać, mamy do czynienia z obiektami typu SeqRecord, zaraz przyjrzymy się im bliżej.
Ponieważ obiekt sekwencje jest iteratorem, pobieranie elementów wiąże się z ich usuwaniem, zatem jeśli chcemy elementy wykorzystać wielokrotnie, można przekazać je do listy, choć najpierw musimy je ponownie wczytać:
sekwencje = SeqIO.parse('sekwencje_atp6.fasta', 'fasta')
sekwencje = list(sekwencje)
print(f'Typ: {type(sekwencje)}')
Typ: <class 'list'>
Teraz przyjrzyjmy się obiektowi SeqRecord:
sekw_1 = sekwencje[0]
print(sekw_1)
ID: AY007817_Daucus_carota
Name: AY007817_Daucus_carota
Description: AY007817_Daucus_carota
Number of features: 0
Seq('ACTTAACTGCTCGAAAGGCAGAGGCTGAGGCAACAGCTGACCACGGTTCTTTTC...AAT')
Jak można przeczytać w dokumentacji SeqRecord jest obiektem używanym do przechowywania sekwencji (w obiekcie typu Seq) oraz jej identyfikatora (ID), nazwy (name), opisu (Description) oraz ewentualnie dodatkowych cech i adnotacji.
Uzyskajmy teraz poszczególne, istotniejsze elementy z obiektu:
print(f'ID: {sekw_1.id}')
print(f'Nazwa: {sekw_1.name}')
print(f'Opis: {sekw_1.description}')
print(f'Sekwencja: {sekw_1.seq}')
ID: AY007817_Daucus_carota
Nazwa: AY007817_Daucus_carota
Opis: AY007817_Daucus_carota
Sekwencja: ACTTAACTGCTCGAAAGGCAGAGGCTGAGGCAACAGCTGACCACGGTTCTTTTCTTTAATCCCTTAACCTTCAATTAATGCCTGGTAAACATGTAGAAAAGACAGCTGAGAAGGCCTTTCTTTACTTTAGGGTCGCTCTGTTCATCCAACTCGAAATATCGTAAGAGAATAAAAAGAATCTTACGCCCAACAAGTCCCGTGCCTTTCTTGGTCGGACCAACCCAACCGGCTCTTTCCGACAAGTCTTTATGCATTTTTAGAGCAAGAAGCGGAACTATAGAGATGAAATTAAAAGTTTTTCTTACGATTACGCCCAACAGCCCACTTGAGCAATTTGCCATTCTCCCGTTGGTTCCTATGAATATAGGACACTTGTATTTCTCATTCACAAATCCCTCTTTGTTTATGCTACTCACTCTCAGTTTGGTCCTACTTTTGGTTCATTTTGTTACTAAAAACGGAGGAGGAAACTCAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTTATGATTTCGTGCCGAACCCGGTAAACGAACAAATAGGTGGTCTTTCCGGAAATGTTAAACAAAAGTTTTCCCCTCGCATCTCGGTCACTTTTACTTTTTCGTTCTTTCGTAATCCACAGGGTATGATACCTTATAGCTTCACAGTTACAAGTCATTTTCTCATTACTTTGGGCCTCTCATTTTCCATTTTTATTGGCATTACTCTCGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTTTCATTACCCGCAGGAGTCCCACTGCCGTTAGCACCTTTTTTAGTACTCCTTGAGCTAATCCCTCATTGTTTTCGCGCATTAAGCTCAGGAATACGTTTATTTGCTAATATGATGGCCGGTCATAGTTCAGTAAAGATTTTAAGTGGGTTCGCTTGGACTATGCTATGTATGAATGATCTTTTCTATTTCATAGGAAATCTGGGTCCTTTATTTATAGTTCTTGCATTAACCGGTCCGGAATTAGGTGTAGCTATATCACAAGCTCATGTTTCTACAATCTCAATCTGTATTTACTTGAATGATGCTACAAATCTCCATCAAAGTGAAGTGGTTATTTATAATTGAACAAAAGCGAGGAGCGAAGGCCGCAAT
Jak widać, ID, nazwa i opis zawierają tę samą wartość, jest to po prostu ciąg znaków po symbolu > w linii opisującej sekwencję.
Wczytajmy teraz plik w formacie genbank, który zawiera dużo więcej informacji na temat sekwencji.
Pobierz plik sekwencje_atp6.gb i zapisz w folderze, w którym uruchamiasz kod. Następnie uruchom kod:
from Bio import SeqIO
sekwencje_GB = SeqIO.parse('sekwencje_atp6.gb', 'gb')
sekwencje_GB = list(sekwencje_GB)
sekw_GB = sekwencje_GB[1]
print(sekw_GB)
ID: KU180476.1
Name: KU180476
Description: Centaurea scabiosa clone H7 ATPase subunit 6 (atp6) gene, partial cds; mitochondrial
Number of features: 3
/molecule_type=DNA
/topology=linear
/data_file_division=PLN
/date=10-MAY-2016
/accessions=['KU180476']
/sequence_version=1
/keywords=['']
/source=mitochondrion Centaurea scabiosa
/organism=Centaurea scabiosa
/taxonomy=['Eukaryota', 'Viridiplantae', 'Streptophyta', 'Embryophyta', 'Tracheophyta', 'Spermatophyta', 'Magnoliopsida', 'eudicotyledons', 'Gunneridae', 'Pentapetalae', 'asterids', 'campanulids', 'Asterales', 'Asteraceae', 'Carduoideae', 'Cardueae', 'Centaureinae', 'Centaurea']
/references=[Reference(title='Horizontal gene transfer found in Orobanche coerulescens', ...), Reference(title='Direct Submission', ...)]
/structured_comment=OrderedDict([('Assembly-Data', OrderedDict([('Sequencing Technology', 'Sanger dideoxy sequencing')]))])
Seq('CTGCTAACTCTCAGTTTTGTCCTACTTCTGATTCATTTTGTTACTAAAAAAGGA...TAC')
Jak widać, liczba informacji się znacznie zwiększyła.
Warto porównać zawartość poszczególnych atrybutów z informacjami w formacie genbank, w pliku, który został wczytany.
Użyjmy znów funkcji dir(), aby sprawdzić listę dostępnych metod i atrybutów (oczywiście pomocna może też być dokumentacja:
dir(sekw_GB)
['__add__',
'__bool__',
...
'annotations',
'dbxrefs',
'description',
'features',
'format',
'id',
'letter_annotations',
'lower',
'name',
'reverse_complement',
'seq',
'translate',
'upper']
Spróbujmy zatem wydobyć z obiektu niektóre dodatkowe dane, zapisane jako annotations:
print(sekw_GB.annotations)
{'molecule_type': 'DNA', 'topology': 'linear', 'data_file_division': 'PLN', 'date': '10-MAY-2016', 'accessions': ['KU180476'], 'sequence_version': 1, 'keywords': [''], 'source': 'mitochondrion Centaurea scabiosa', 'organism': 'Centaurea scabiosa', 'taxonomy': ['Eukaryota', 'Viridiplantae', 'Streptophyta', 'Embryophyta', 'Tracheophyta', 'Spermatophyta', 'Magnoliopsida', 'eudicotyledons', 'Gunneridae', 'Pentapetalae', 'asterids', 'campanulids', 'Asterales', 'Asteraceae', 'Carduoideae', 'Cardueae', 'Centaureinae', 'Centaurea'], 'references': [Reference(title='Horizontal gene transfer found in Orobanche coerulescens', ...), Reference(title='Direct Submission', ...)], 'structured_comment': OrderedDict([('Assembly-Data', OrderedDict([('Sequencing Technology', 'Sanger dideoxy sequencing')]))])}
Jak widać otrzymaliśmy słownik, zawierający informacje, które wyświetlały się po znakach /, w pliku źródłowym znajdują się w sekwencji FEATURES.
Spróbujmy uzyskać niektóre z nich:
print(f'molecule_type: {sekw_GB.annotations["molecule_type"]}')
print(f'date: {sekw_GB.annotations["date"]}')
print(f'organism: {sekw_GB.annotations["organism"]}')
print(f'taxonomy: {sekw_GB.annotations["taxonomy"]}')
molecule_type: DNA
date: 10-MAY-2016
organism: Centaurea scabiosa
taxonomy: ['Eukaryota', 'Viridiplantae', 'Streptophyta', 'Embryophyta', 'Tracheophyta', 'Spermatophyta', 'Magnoliopsida', 'eudicotyledons', 'Gunneridae', 'Pentapetalae', 'asterids', 'campanulids', 'Asterales', 'Asteraceae', 'Carduoideae', 'Cardueae', 'Centaureinae', 'Centaurea']
Skoro po odwołaniu się do atrybuty seq dostajemy obiekt typu Seq, możemy wykorzystać jego możliwości (dla przejrzystości wyświetlimy pierwsze 40 znaków):
print(f'sekwencja: {sekw_GB.seq[:40]}')
print(f'komplementarna: {sekw_GB.seq.complement()[:40]}')
print(f'odwr. komplem.: {sekw_GB.seq.reverse_complement()[:40]}')
sekwencja: CTGCTAACTCTCAGTTTTGTCCTACTTCTGATTCATTTTG
komplementarna: GACGATTGAGAGTCAAAACAGGATGAAGACTAAGTAAAAC
odwr. komplem.: GTAAATACAGATTAAGATCGTAAAAACATAAGCTTGTAAT
Skoro mamy listę sekwencji, możemy ją iterować uzyskując potrzebne informacje dla każdej, co można wykorzystać np. do utworzenia pliku w formacie csv/tsv:
plik = open('sekwencje.tsv', 'wt')
plik.write('Nr_GB:\tOrganizm\tDługość\tSekwencja\n')
for s in sekwencje_GB:
plik.write(f'{s.name}\t{s.annotations["organism"]}\t{len(s.seq)}\t{s.seq}\n')
plik.close()
Sprawdź zawartość pliku sekwencje.tsv.
Możemy też wygenerować plik w formacie FASTA formatując opisy sekwencji wg. uznania:
plik = open('sekwencje_eksport.fasta', 'wt')
for s in sekwencje_GB:
plik.write(f'>{s.name}_{s.annotations["organism"].replace(" ","_")}\n')
plik.write(f'{s.seq}\n')
plik.close()
Sprawdź zawartość pliku sekwencje_eksport.fasta.
Jeśli nie chcemy modyfikować opisu sekwencji, do uzyskania formatu FASTA, możemy użyć metody format(), zostanie wtedy użyte pole description:
for s in sekwencje_GB:
print(s.format("fasta"))
>KU180471.1 Artemisia campestris clone H5 ATPase subunit 6 (atp6) gene, partial cds; mitochondrial
CTGCTAACTCTCAGTTTGGTCCTACTTCTGATTCATTTTGTTACTAAAAAAGGAGGAGGA
AACTTAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTTATGATTTCGTGCTGAAC
CTGGTAAACGAACAAATAGGCGGTCTTTCCGGAAATGTTAAACAAAAGTTTTTCCCTTGC
ATCTTGGTCACTTTTACTTTTTTGTTATTTTGTAATCTTCAGGGTATGATACCTTATAGC
TTCACAGTTACAAGTCATTTTCTCATTACTTTAGGTCTCTCATTTTCTCTTTTTATTGGC
ATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTTATTACCCGCA
GGAGTCCCACTGCCGTTAGCACCCTTTTTAGTACTCCTTGAGCTAATTTCTTATTGTTTT
CGCGCATTAAGCTTAGGAATACGTTTATTTGCTAATATGATGGCCGGTCATAGTTTAGTA
AAGATTTTAAGTGGGTTCGCTTGGACTATGCTATGTATGAATGATATTTTGTTTTTTATA
GGGGATCTTGGTCCTTTATTTATAGTTCTTGCATTAACTGGTCTGGAATTAGGTGTAGCT
ATATTACAAGCTTATGTTTTTACGATCTTAATCTGTATTTAC
>KU180476.1 Centaurea scabiosa clone H7 ATPase subunit 6 (atp6) gene, partial cds; mitochondrial
CTGCTAACTCTCAGTTTTGTCCTACTTCTGATTCATTTTGTTACTAAAAAAGGAGGAGGA
AACTTAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTTATGATTTCGTGCTGAAC
CTGGTAAACGAACAAATAGGGGGTCTTTCCGGAAATGTTAAACAAAAGTTTTTCCCTTGC
ATCTTGGTCACTTTTACTTTTTTGTTATTTTGTAATCTTCAGGGTATGATACCTTATAGC
TTCACAGTTACAAGTCATTTTCTCATTACTTTAGGTCTCTCATTTTCTATTTTTATTGGC
ATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTTATTACCCGCA
GGAGTCCCACTGCCGTTAGCACCTTTTTTAGTACTCCTTGAGCTAATTTCTTATTGTTTT
CGCGCATTAAGCTTAGGAATACGTTTATTTGCTAATATGATGGCCGGTCATAGTTTAGTA
AAGATTTTAAGTGGGTTCGCTTGGACTATGCTATGTATGAATGATATTTTTTATTTTATA
GGGGATCTTGGTCCTTTATTTATAGTTCTTGCATTAACCGGTCTGGAATTAGGTGTAGCT
ATATTACAAGCTTATGTTTTTACGATCTTAATCTGTATTTAC
...
Biopython i bazy NCBI Entrez
Biopython umożliwia także kontaktowanie się i pobieranie danych z baz NCBI, dzięki czemu możemy np. wyszukiwać i pobierać sekwencje z bazy GenBank, czemu poświęcimy nieco więcej uwagi.
Entrez jest systemem, który pozwala na łączenie się z bazami danych NCBI, obejmującymi m.in. bazy sekwencji nukleotydów i genomów (GenBank), białek, publikacji naukowych itd.
Jeśli chcemy skontaktować się z bazą udostępnianą przez Entrez, powinniśmy znać jej nazwę.
Połączmy się zatem z Entrez, pobierzmy nazwy baz i wyświetlmy je:
from Bio import Entrez
# Powinniśmy określić PRAWDZIWY adres e-mail
Entrez.email = "twoj.prawdziwy@adres.e-mail"
# Pobranie informacji z Entrez
handle = Entrez.einfo()
# Parsowanie uzyskanej informacji do słownika
bazy = Entrez.read(handle)
# Zamykanie obiektu handle
handle.close()
print(bazy)
{'DbList': ['pubmed', 'protein', 'nuccore', 'ipg', 'nucleotide', 'structure', 'genome', 'annotinfo', 'assembly', 'bioproject', 'biosample', 'blastdbinfo', 'books', 'cdd', 'clinvar', 'gap', 'gapplus', 'grasp', 'dbvar', 'gene', 'gds', 'geoprofiles', 'homologene', 'medgen', 'mesh', 'ncbisearch', 'nlmcatalog', 'omim', 'orgtrack', 'pmc', 'popset', 'proteinclusters', 'pcassay', 'protfam', 'biosystems', 'pccompound', 'pcsubstance', 'seqannot', 'snp', 'sra', 'taxonomy', 'biocollections', 'gtr']}
Jak widać, otrzymaliśmy słownik z jednym kluczem (DbList) i listą wartości, które są nazwami baz danych.
Przy okazji zwróćmy uwagę na obiekt, który nazwaliśmy (i tak jest zwykle nazywany) handle.
print(type(handle))
<class 'http.client.HTTPResponse'>
Jest to obiekt, który ,,opakowuje’’ dane tekstowe, sformatowane w różny sposób, tak aby m.in. można ich było używać w sposób ustandaryzowany.
Więcej na ten temat można poczytać np. tu: http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec%3Aappendix-handles.
Zaleca się, żeby po pobraniu z niego danych, ,,zamykać" go: handle.close()
Możemy teraz pobrać opis, dla każdej bazy:
print(f'Baza{" ":<11s} | Opis')
print("-"*70)
for baza in bazy['DbList']:
# Pobieramy dane o bazie
handle = Entrez.einfo(db=baza)
dane_bazy = Entrez.read(handle)
handle.close()
print(f'{dane_bazy["DbInfo"]["DbName"]:<16s}| {dane_bazy["DbInfo"]["Description"]}')
Baza | Opis
----------------------------------------------------------------------
pubmed | PubMed bibliographic record
protein | Protein sequence record
nuccore | Core Nucleotide db
ipg | Identical Protein Groups DB
nuccore | Core Nucleotide db
structure | Three-dimensional molecular model
genome | Genomic sequences, contigs, and maps
annotinfo | Annotinfo Database
assembly | Genome Assembly Database
bioproject | BioProject Database
biosample | BioSample Database
blastdbinfo | BlastdbInfo Database
books | Books Database
cdd | Conserved Domain Database
clinvar | ClinVar Database
gap | dbGaP Data
gapplus | Internal Genotypes and Phenotypes database
grasp | grasp Data
dbvar | dbVar records
gene | Gene database
gds | GEO DataSets
geoprofiles | Genes Expression Omnibus
homologene | HomoloGene Database
medgen | Medgen Database
mesh | MeSH Database
ncbisearch | NcbiSearch Database
nlmcatalog | NLM Catalog Database
omim | OMIM records
orgtrack | Orgtrack Database
pmc | PubMed Central
popset | PopSet sequence record
proteinclusters | Protein Cluster record
pcassay | PubChem BioAssay Database
protfam | protfam DB
biosystems | BioSystems Database
pccompound | PubChem Compound Database
pcsubstance | PubChem Substance Database
seqannot | SeqAnnot Database
snp | Single Nucleotide Polymorphisms
sra | SRA Database
taxonomy | Taxonomy db
biocollections | Biocollections db
gtr | GTR Database
Spróbujmy teraz wyszukać sekwencje w bazie sekwencji nukleotydów. W tym celu utworzymy zapytanie, z użyciem nazw pól, które zostaną przeszukane. Nazwy pól używanych w bazach sekwencji, można znaleźć np. tu: https://www.ncbi.nlm.nih.gov/books/NBK49540/ Ograniczymy liczbę pobieranych rekordów do 5, żeby niepotrzebnie nie obciążać serwera.
from Bio import Entrez
# Powinniśmy określić PRAWDZIWY adres e-mail
Entrez.email = "twoj.prawdziwy@adres.e-mail"
# Zapytanie:
# organizm: Orobanche oraz gen: atp6
zapytanie = "Orobanche[ORGN] AND atp6[GENE]"
# Pobranie informacji z Entrez
handle = Entrez.esearch(db='nucleotide',
term=zapytanie,
# ustawienie formatu zwracanych identyfikatorów
# na "accession number"
idtype='acc',
# maksymalna liczba zwracanych rekordów
# (domyślnie 20)
retmax = 5
)
wyniki = Entrez.read(handle)
# Zamknij obiekt handle po użyciu.
handle.close()
print(wyniki["Count"])
print(wyniki["IdList"])
10
['KU180474.1', 'KU180473.1', 'KU180472.1', 'KU180470.1', 'KU180469.1']
Przyjrzyjmy się jeszcze całej zawartości obiektu (słownika) wyniki:
print(wyniki)
{'Count': '10', 'RetMax': '5', 'RetStart': '0', 'IdList': ['KU180474.1', 'KU180473.1', 'KU180472.1', 'KU180470.1', 'KU180469.1'], 'TranslationSet': [{'From': 'Orobanche[ORGN]', 'To': '"Orobanche"[Organism]'}], 'TranslationStack': [{'Term': '"Orobanche"[Organism]', 'Field': 'Organism', 'Count': '164147', 'Explode': 'Y'}, {'Term': 'atp6[GENE]', 'Field': 'GENE', 'Count': '129226', 'Explode': 'N'}, 'AND'], 'QueryTranslation': '"Orobanche"[Organism] AND atp6[GENE]'}
Pobierzmy teraz pierwszą z sekwencji w formacie FASTA:
handle = Entrez.efetch(db="nucleotide",
# Numer dost. pierwszej ze znalezionych sekwencji
id=wyniki["IdList"][0],
# Format pobranej sekwencji
rettype="fasta")
print(handle.read())
handle.close()
>KU180474.1 Orobanche grenieri clone G8 ATPase subunit 6 (atp6) gene, partial cds; mitochondrial
CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTGTTACTAAAAAGGGAGGAGGAAACTCAGTAC
CAAATGCTTTTCAATCCGTGTTAGAGCTTATTTATGATTTTGTGCCGAACCTGGTAAACGAACAAATAGG
TGGTCTTTCCGGAAATGTTAAACAAAAGTTTTTCCCTTGCATCTCGGTTACTTTTACTTTTTCGTTATTT
CGTAATCTGCAGGGTATGATACCTTATAGCTTCACAGTAACAAGTCATTTTATCGTTACTTTGGGTCTCT
CATTTTCTCTTTTTATTGGCATTACTATAGTGGGATTTCAAAAAAATGGGCTTCATTTTTTAAGCTTCTC
ATTACCCGCAGGAGTCCCACTGCCGTTAGCACCTTTTTTAGTACTCCTTGAGCTAATCCCTCATTGTTTT
CGCGCATTAAGCTTAGGAATACGTTTATTTGCTAATATGATGGCCGGTCATAGTTTAGTAAAGATTTTAA
GTGGGTTCGCTTGGACTATGCTATGTATGAATGATCTTTTCTATTTCATAGGGGATCCTGGTCCTTTATT
TATAGTTCTTGCATTAACCGGTCTTGAATTAGGTGTAGCTATATCACAAGCTCATGTTTCTACGATCTCA
ATCTGTATTTAC
Oczywiście możemy pobrać dane w formacie gb:
handle = Entrez.efetch(db="nucleotide",
id=wyniki["IdList"][0],
# Format pobranej sekwencji
rettype="gb")
print(handle.read())
handle.close()
LOCUS KU180474 642 bp DNA linear PLN 10-MAY-2016
DEFINITION Orobanche grenieri clone G8 ATPase subunit 6 (atp6) gene, partial
cds; mitochondrial.
ACCESSION KU180474
VERSION KU180474.1
KEYWORDS .
SOURCE mitochondrion Orobanche grenieri
ORGANISM Orobanche grenieri
Eukaryota; Viridiplantae; Streptophyta; Embryophyta; Tracheophyta;
Spermatophyta; Magnoliopsida; eudicotyledons; Gunneridae;
Pentapetalae; asterids; lamiids; Lamiales; Orobanchaceae;
Orobancheae; Orobanche.
REFERENCE 1 (bases 1 to 642)
AUTHORS Kwolek,D., Denysenko,M., Goralski,G., Cygan,M., Mizia,P.,
Piwowarczyk,R. and Joachimiak,A.J.
TITLE Horizontal gene transfer found in Orobanche coerulescens
JOURNAL Unpublished
REFERENCE 2 (bases 1 to 642)
AUTHORS Kwolek,D., Denysenko,M., Goralski,G., Cygan,M., Mizia,P.,
Piwowarczyk,R. and Joachimiak,A.J.
TITLE Direct Submission
JOURNAL Submitted (23-NOV-2015) Dept. of Plant Cytology and Embryology,
Institute of Botany, Jagiellonian University, ul. Gronostajowa 9,
pok. 2.33, Krakow 30-387, Poland
COMMENT ##Assembly-Data-START##
Sequencing Technology :: Sanger dideoxy sequencing
##Assembly-Data-END##
FEATURES Location/Qualifiers
source 1..642
/organism="Orobanche grenieri"
/organelle="mitochondrion"
/mol_type="genomic DNA"
/db_xref="taxon:320802"
/clone="G8"
/country="Georgia: Kortaneti"
/collected_by="Piwowarczyk, R."
/identified_by="Piwowarczyk, R."
gene <1..>642
/gene="atp6"
CDS <1..>642
/gene="atp6"
/codon_start=1
/product="ATPase subunit 6"
/protein_id="ANC68076.1"
/translation="LLTLSLVLLFVHFVTKKGGGNSVPNAFQSVLELIYDFVPNLVNE
QIGGLSGNVKQKFFPCISVTFTFSLFRNLQGMIPYSFTVTSHFIVTLGLSFSLFIGIT
IVGFQKNGLHFLSFSLPAGVPLPLAPFLVLLELIPHCFRALSLGIRLFANMMAGHSLV
KILSGFAWTMLCMNDLFYFIGDPGPLFIVLALTGLELGVAISQAHVSTISICIY"
ORIGIN
1 ctactcactc tcagtttggt cctacttttt gttcattttg ttactaaaaa gggaggagga
61 aactcagtac caaatgcttt tcaatccgtg ttagagctta tttatgattt tgtgccgaac
121 ctggtaaacg aacaaatagg tggtctttcc ggaaatgtta aacaaaagtt tttcccttgc
181 atctcggtta cttttacttt ttcgttattt cgtaatctgc agggtatgat accttatagc
241 ttcacagtaa caagtcattt tatcgttact ttgggtctct cattttctct ttttattggc
301 attactatag tgggatttca aaaaaatggg cttcattttt taagcttctc attacccgca
361 ggagtcccac tgccgttagc acctttttta gtactccttg agctaatccc tcattgtttt
421 cgcgcattaa gcttaggaat acgtttattt gctaatatga tggccggtca tagtttagta
481 aagattttaa gtgggttcgc ttggactatg ctatgtatga atgatctttt ctatttcata
541 ggggatcctg gtcctttatt tatagttctt gcattaaccg gtcttgaatt aggtgtagct
601 atatcacaag ctcatgtttc tacgatctca atctgtattt ac
//
Jeśli chcemy uzyskać poszczególne informacje, wynik należy nieco inaczej odczytać:
from Bio import Entrez
from Bio import SeqIO
Entrez.email = "twoj.prawdziwy@adres.e-mail"
handle = Entrez.efetch(db="nucleotide",
id=wyniki["IdList"][0],
rettype="gb")
# Taki sposób odczytu, pozwala odwoływać się do poszczególnych danych
pobrane = SeqIO.read(handle, "genbank" )
handle.close()
print(f'Numer: {pobrane.id}')
print(f'Organizm: {pobrane.annotations["organism"]}')
print(f'Sekwencja: {pobrane.seq[:40]}...')
print(f'Opis: {pobrane.description}')
Numer: KU180474.1
Organizm: Orobanche grenieri
Sekwencja: CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTG...
Opis: Orobanche grenieri clone G8 ATPase subunit 6 (atp6) gene, partial cds; mitochondrial
Wszystkie znalezione sekwencje możemy pobrać przekazując parametrowi id listę numerów dostępowych.
Warto w takim przypadku zapisać wyniki w pliku, żeby niepotrzebnie nie obciążać serwera kolejnymi, takimi samymi zapytaniami.
Później możemy je odczytać w celu dalszej analizy.
Połączmy wszystko w całość.
Przy pisaniu skryptów, nie zaszkodzi umieszczać linii kodu, które będą nas informowały, co obecnie robi skrypt, zwłaszcza, że pobieranie danych może zająć dużo czasu i czasem możemy przypuszczać, że np. skrypt się zawiesił.
from Bio import SeqIO
from Bio import Entrez
Entrez.email = "twoj.prawdziwy@adres.e-mail"
zapytanie = "Orobanche[ORGN] AND atp6[GENE]"
# Wyszukiwanie sekwencji
handle = Entrez.esearch(db='nucleotide',
term=zapytanie,
idtype='acc',
retmax = 5
)
wyniki = Entrez.read(handle)
print(f'Wyszukano {wyniki["Count"]} sekwencji.')
# Pobieranie sekwencji
print('Pobieram sekwencje.')
handle = Entrez.efetch(db="nucleotide",
id=wyniki["IdList"],
rettype="gb")
pobrane = handle.read()
handle.close()
plik = 'Pobrane-atp6.gb'
print(f'Zapisuję wyniki w pliku {plik}.')
plik_zapis = open(plik, 'w')
plik_zapis.write(pobrane)
# Odczyt danych
print(f'Odczytuję dane z pliku: {plik}.')
dane = SeqIO.parse(plik, 'gb')
sekwencje = list(dane)
print(f'Odczytano {len(sekwencje)} sekwencji.')
# Tu można wpisać dowolny kod analizujący sekwencje
# np:
for sekw in sekwencje:
print(f'{sekw.name} - {sekw.annotations["organism"]:<30s}: {sekw.seq[:10]}')
print('\nKoniec')
Wyszukano 10 sekwencji.
Pobieram sekwencje.
Zapisuję wyniki w pliku Pobrane-atp6.gb.
Odczytuję dane z pliku: Pobrane-atp6.gb.
Odczytano 5 sekwencji.
KU180474 - Orobanche grenieri : CTACTCACTC
KU180473 - Orobanche coerulescens : CTGCTAACTC
KU180472 - Orobanche cernua : CTACTCACTC
KU180470 - Orobanche elatior : CTACTCACTC
KU180469 - Orobanche alba subsp. alba : CTACTCACTC
Koniec
Biopython i BLAST
Kolejnym ważnym zastosowaniem Biopythona jest możliwość korzystania z narzędzia BLAST, który umożliwia wyszukanie sekwencji podobnych do sekwencji wysyłanej jako zapytanie.
Zacznijmy od prostego przykładu:
from Bio.Blast import NCBIWWW
# Sekwencja, która będzie zapytaniem
zapytanie = "\
CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTGTTACTAAAAAGGGAGGAGGAAACTCAGTAC\
CAAATGCTTTTCAATCCGTGTTAGAGCTTATTTATGATTTTGTGCCGAACCTGGTAAACGAACAAATAGG"
print("Wyszukuję dane...")
handle = NCBIWWW.qblast(program = "blastn",
database = "nt",
sequence = zapytanie,
hitlist_size = 10)
print("Dane pobrane.")
# Zapis do pliku
plik_zapis = open("Wynik-BLAST.xml", "w")
plik_zapis.write(handle.read())
plik_zapis.close()
handle.close()
Wyszukuję dane...
Dane pobrane.
Wyszukiwanie i pobieranie wyników zajmuje czasem dużo czasu, należy więc możliwie ograniczyć kontakt z serwerem i raczej zapisać wyniki w pliku, z którego następnie go odczytamy.
Wywołując funkcję qblast() ustawiliśmy cztery parametry:
program- określa program używany dla wyszukiwania, możliwe wartości:blastn,blastp,blastx,tblastn, lubtblastx. Jeśli szukamy sekwencji DNA dla zapytania DNA to wybieramyblastn.database- baza, z którą się kontaktujemy. Domyślnie jest tont, obejmująca m. in.GenBanksequence- sekwencja, dla której szukamy sekwencji podobnych, można też podac numer dostępowy GenBank.hitlist_size- maksymalna liczba wyników, domyślnie jest to 50.
Parametrów, które możemy ustawić jest znacznie więcej. Warto zajrzeć do dokumentacji: https://biopython.org/docs/1.75/api/Bio.Blast.NCBIWWW.html a także do przewodnika (choć dotyczy on główie interfejsu www), który znajduje się pod adresem: ftp://ftp.ncbi.nlm.nih.gov/pub/factsheets/HowTo_BLASTGuide.pdf.
Jak wspomniałem, czasem uzyskanie wyników zajmuje długi czas, ma to związek z przeciążeniami serwerów. Zwykle praca przebiega zdecydowanie sprawniej, jeśli się łączymy poza czasem pracy naukowców w USA i/lub w weekendy.
Jeśli udało się pobrać wyniki, przejrzyj zawartość pliku Wynik-BLAST.xml.
Teraz odczytajmy dane z pliku:
from Bio import SearchIO
wyniki = SearchIO.read('Wynik-BLAST.xml', 'blast-xml')
Po przekazaniu obiektu wyniki do funkcji print(), otrzymamy kilka informacji oraz tabelkę z krótkim opisem wyników.
print(wyniki)
Program: blastn (2.11.0+)
Query: No (140)
definition line
Target: nt
Hits: ---- ----- ----------------------------------------------------------
# # HSP ID + description
---- ----- ----------------------------------------------------------
0 1 gi|1025818475|gb|KU180474.1| Orobanche grenieri clone ...
1 1 gi|1025818471|gb|KU180472.1| Orobanche cernua clone G1...
2 1 gi|1025818461|gb|KU180467.1| Orobanche gracilis clone ...
3 1 gi|1025818465|gb|KU180469.1| Orobanche alba subsp. alb...
4 1 gi|1025818453|gb|KU180463.1| Orobanche picridis clone ...
5 1 gi|1025818467|gb|KU180470.1| Orobanche elatior clone 1...
6 1 gi|1025818457|gb|KU180465.1| Orobanche caryophyllacea ...
7 1 gi|1126299743|gb|KX524674.1| Lindenbergia siniaca isol...
8 1 gi|1025818463|gb|KU180468.1| Phelipanche purpurea subs...
9 1 gi|1025818459|gb|KU180466.1| Phelipanche ramosa clone ...
Sprawdźmy typ obiektu:
print(type(wyniki))
<class 'Bio.SearchIO._model.query.QueryResult'>
Jak widać, jest to obiekt typu QueryResult.
Przechowuje on m.in. wyniki dopasowań, pobierzmy pierwszy z nich, sprawdźmy, jakiego jest typu i przekażmy do funkcji print():
print(f'Typ: {type(wyniki[0])}')
print(wyniki[0])
Typ: <class 'Bio.SearchIO._model.hit.Hit'>
Query: No
definition line
Hit: gi|1025818475|gb|KU180474.1| (642)
Orobanche grenieri clone G8 ATPase subunit 6 (atp6) gene, partial cds;...
HSPs: ---- -------- --------- ------ --------------- ---------------------
# E-value Bit score Span Query range Hit range
---- -------- --------- ------ --------------- ---------------------
0 1.8e-63 253.76 140 [0:140] [0:140]
Jak widać, mamy do czynienia z obiektem typu Hit (zob. http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec132.
Zawiera wszystkie wyniki pochodzące z jednego elementu bazy, czyli jednej sekwencji o określonym numerze dostępowym.
W powyższym przykładzie nie mamy takich przypadków, ale w jednej docelowej sekwencji (zwłaszcza jeśli jest ona sekwencją genomu) może znaleźć się wiele dopasowań do wyszukiwanej sekwencji.
Takie dopasowanie nazywamy HSP (ang. high-scoring pair), które może się składać z kilku fragmentów (HSPFragment).
Sprawdźmy HSP z siódmego wyniku (Hit):
hsp = wyniki[6][0]
print(hsp)
Query: No definition line
Hit: gi|1025818457|gb|KU180465.1| Orobanche caryophyllacea clone 1 AT...
Query range: [0:140] (1)
Hit range: [0:140] (1)
Quick stats: evalue 9.3e-61; bitscore 244.74
Fragments: 1 (140 columns)
Query - CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTGTTACTAAAAAGGGAGGAGG~~~ATAGG
||||||||||||||||||||||||||||||||||||||||||||| |||||||||||||~~~|||||
Hit - CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTGTTACTCAAAAGGGAGGAGG~~~ATAGG
Query to dopasowana sekwencja (lub jej dopasowany fragment), który wysłaliśmy jako zapytanie, Hit to znaleziony, pasujący fragment nici DNA.
W powyższym podsumowaniu wyświetla się fragment obu nici.
Pomiędzy nimi mamy oznaczenia pokazujące wizualnie pasujące, lub nie, pary zasad.
Znak | oznacza, że w danym miejscu w obu sekwencjach znajduje się taka sama zasada.
Obie sekwencje, a także inne dane możemy pobrać.
hsp = wyniki[6][0]
print(type(hsp))
print(f'Nazwa: {hsp.hit_description}')
print(f'Query: {hsp.query.seq[:70]}')
print(f'Hit: {hsp.hit.seq[:70]}')
print(f'Hit ID: {hsp.hit.id}')
# Identyczne miejsca
print(f'Ident.: {hsp.ident_num}')
# Miejsc ze znakiem -
print(f'Gaps.: {hsp.gap_num}')
# Otrzymamy numer dostępowy GenBank
nr = hsp.hit.id.split('|')[3]
print(f'Nr. GB: {nr}')
<class 'Bio.SearchIO._model.hsp.HSP'>
Nazwa: Orobanche caryophyllacea clone 1 ATPase subunit 6 (atp6) gene, partial cds; mitochondrial
Query: CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTGTTACTAAAAAGGGAGGAGGAAACTCAGTAC
Hit: CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTGTTACTCAAAAGGGAGGAGGAAAGTCAGTAC
Hit ID: gi|1025818457|gb|KU180465.1|
Ident.: 138
Gaps.: 0
Nr. GB: KU180465.1
Sprawdźmy pierwszy (i jedyny) HSPFragment:
hsp_fr = wyniki[6][0][0]
print(type(hsp_fr))
print(f'Nazwa: {hsp_fr.hit_description}')
print(f'Query: {hsp_fr.query.seq[:70]}')
print(f'Hit: {hsp_fr.hit.seq[:70]}')
<class 'Bio.SearchIO._model.hsp.HSPFragment'>
Nazwa: Orobanche caryophyllacea clone 1 ATPase subunit 6 (atp6) gene, partial cds; mitochondrial
Query: CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTGTTACTAAAAAGGGAGGAGGAAACTCAGTAC
Hit: CTACTCACTCTCAGTTTGGTCCTACTTTTTGTTCATTTTGTTACTCAAAAGGGAGGAGGAAAGTCAGTAC
Mając numer sekwencji w której znaleziono dopasowanie, możemy pobrać ją w całości, a także dowiedzieć się więcej na temat organizmu.
Najpierw pobierzmy dane z bazy nucleotide:
from Bio import Entrez
from Bio import SeqIO
Entrez.email = "twoj.prawdziwy@adres.e-mail"
handle = Entrez.efetch(db="nucleotide",
id=nr,
retmode="text",
rettype="gb")
pobrane = SeqIO.read(handle, "genbank" )
handle.close()
Teraz pobierzmy poszczególne informacje:
print(f'Numer: {pobrane.id}')
print(f'Organizm: {pobrane.annotations["organism"]}')
print(f'Taksonomia: {pobrane.annotations["taxonomy"]}')
Numer: KU180465.1
Organizm: Orobanche caryophyllacea
Taksonomia: ['Eukaryota', 'Viridiplantae', 'Streptophyta', 'Embryophyta', 'Tracheophyta', 'Spermatophyta', 'Magnoliopsida', 'eudicotyledons', 'Gunneridae', 'Pentapetalae', 'asterids', 'lamiids', 'Lamiales', 'Orobanchaceae', 'Orobancheae', 'Orobanche']
Pobierzmy teraz wyniki dla nowej sekwencji i wypiszmy wybrane dane na temat wyników dopasowania.
Przy okazji pokażę jak można dane taksonomiczne pobrać z bazy taxonomy.
Najpierw pobierzmy dane z BLAST:
from Bio import Entrez
from Bio import SearchIO
from Bio.Blast import NCBIWWW
Entrez.email = "twoj.prawdziwy@adres.e-mail"
# Plik z wynikami wyszukwania:
plik = "Wynik_2-BLAST.xml"
# Sekwencja, która będzie zapytaniem
sekw = "GCGATGGCTTCAGTCAATGCAAAGCCTAGAATAGCATAGCCAAATAATTGTTTAGCTAAGGGCGT"
print("Wyszukuję dane...")
handle = NCBIWWW.qblast(program = "blastn",
database = "nt",
sequence = sekw,
hitlist_size = 10)
print("Dane pobrane.")
# Zapis do pliku
plik_zapis = open(plik, "w")
plik_zapis.write(handle.read())
plik_zapis.close()
handle.close()
Wyszukuję dane...
Dane pobrane.
Teraz odczytujemy wyniki i pobieramy dane organizmów u których znaleziono sekwencje:
def pobierz_nucleotide(nr):
"""Pobiera dane z bazy nucleotide"""
Entrez.email = "twoj.prawdziwy@adres.e-mail"
handle = Entrez.efetch(db="nucleotide",
id=nr,
retmode="text",
rettype="gb")
wynik = SeqIO.read(handle, "genbank" )
handle.close()
return wynik
def pobierz_taksonomie(nazwa):
"""Pobiera dane z bazy taxonomy"""
# Wyszukiwanie nazwy organizmu w bazie
handle = Entrez.esearch(db="taxonomy",
term=nazwa)
organizm = Entrez.read(handle)
handle.close()
# Pobieramy dane dla organizmu przekazując
# identyfikator
handle = Entrez.efetch(db="taxonomy" ,
id=organizm['IdList'] ,
retmode="xml" )
taksonomia = Entrez.read(handle)
handle.close()
return taksonomia
def drukuj_dane(hsp):
"""Drukuje dane"""
# Numer dostępowy GenBank
nr = hsp.hit.id.split('|')[3]
# Pobieranie sekwencji z bazy `nucleotide`
nucl = pobierz_nucleotide(nr)
# Pobieranie danych taksonomicznych
tax = pobierz_taksonomie(nucl.annotations['organism'])
print(f'Nazwa: {tax[0]["ScientificName"]}')
print(f'Systematyka: {tax[0]["Lineage"]}')
print(f'Numer: {nr}')
print(f'Opis sekw.: {hsp.hit_description}')
# Identyczne miejsca
print(f'Ident.: {hsp.ident_num}')
# Miejsc ze znakiem -
print(f'Gaps.: {hsp.gap_num}')
print('_'*60)
# Odczyt z pliku
znalezione = SearchIO.read(plik, 'blast-xml')
for hit in znalezione:
for hsp in hit:
drukuj_dane(hsp)
Nazwa: Ajuga ciliata
Systematyka: cellular organisms; Eukaryota; Viridiplantae; Streptophyta; Streptophytina; Embryophyta; Tracheophyta; Euphyllophyta; Spermatophyta; Magnoliopsida; Mesangiospermae; eudicotyledons; Gunneridae; Pentapetalae; asterids; lamiids; Lamiales; Lamiaceae; Teucrioideae; Ajuga
Numer: MT075725.1
Opis sekw.: Ajuga ciliata chromosome 1 mitochondrion, complete sequence
Ident.: 59
Gaps.: 0
...
____________________________________________________________
Nazwa: Caldesia oligococca
Systematyka: cellular organisms; Eukaryota; Viridiplantae; Streptophyta; Streptophytina; Embryophyta; Tracheophyta; Euphyllophyta; Spermatophyta; Magnoliopsida; Mesangiospermae; Liliopsida; Alismatales; Alismataceae; Caldesia
Numer: KU642324.1
Opis sekw.: Caldesia oligococca ATPase subunit 9 (atp9) gene, complete cds; mitochondrial
Ident.: 52
Gaps.: 0
____________________________________________________________