Python - wprowadzenie

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 wynikowe BLAST itd.
  • łą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=True powoduje zatrzymanie ,,translacji’’ po natrafieniu na kodon STOP

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, lub tblastx. Jeśli szukamy sekwencji DNA dla zapytania DNA to wybieramy blastn.
  • database - baza, z którą się kontaktujemy. Domyślnie jest to nt, obejmująca m. in. GenBank
  • sequence - 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
____________________________________________________________
Last updated on 21 Apr 2021
Published on 21 Apr 2021