Biopython - przegląd BLAST

BLAST oznacza Basic Local Alignment Search Tool. Znajduje regiony podobieństwa między sekwencjami biologicznymi. Biopython dostarcza moduł Bio.Blast do obsługi operacji NCBI BLAST. Możesz uruchomić BLAST w połączeniu lokalnym lub przez połączenie internetowe.

Zrozummy w skrócie te dwa połączenia w następnej sekcji -

Bieganie przez Internet

Biopython dostarcza moduł Bio.Blast.NCBIWWW do wywoływania wersji online BLAST. Aby to zrobić, musimy zaimportować następujący moduł -

>>> from Bio.Blast import NCBIWWW

Moduł NCBIWW zapewnia funkcję qblast do sprawdzania wersji online BLAST, https://blast.ncbi.nlm.nih.gov/Blast.cgi. qblast obsługuje wszystkie parametry obsługiwane przez wersję online.

Aby uzyskać pomoc dotyczącą tego modułu, użyj poniższego polecenia i zapoznaj się z funkcjami -

>>> help(NCBIWWW.qblast) 
Help on function qblast in module Bio.Blast.NCBIWWW: 
qblast(
   program, database, sequence, 
   url_base = 'https://blast.ncbi.nlm.nih.gov/Blast.cgi', 
   auto_format = None, 
   composition_based_statistics = None, 
   db_genetic_code =  None, 
   endpoints = None, 
   entrez_query = '(none)', 
   expect = 10.0, 
   filter = None, 
   gapcosts = None, 
   genetic_code = None, 
   hitlist_size = 50, 
   i_thresh = None, 
   layout = None, 
   lcase_mask = None, 
   matrix_name = None, 
   nucl_penalty = None, 
   nucl_reward = None, 
   other_advanced = None, 
   perc_ident = None, 
   phi_pattern = None, 
   query_file = None, 
   query_believe_defline = None, 
   query_from = None, 
   query_to = None, 
   searchsp_eff = None, 
   service = None, 
   threshold = None, 
   ungapped_alignment = None, 
   word_size = None, 
   alignments = 500, 
   alignment_view = None, 
   descriptions = 500, 
   entrez_links_new_window = None, 
   expect_low = None, 
   expect_high = None, 
   format_entrez_query = None, 
   format_object = None, 
   format_type = 'XML', 
   ncbi_gi = None, 
   results_file = None, 
   show_overview = None, 
   megablast = None, 
   template_type = None, 
   template_length = None
) 
   
   BLAST search using NCBI's QBLAST server or a cloud service provider. 
   
   Supports all parameters of the qblast API for Put and Get. 
   
   Please note that BLAST on the cloud supports the NCBI-BLAST Common 
   URL API (http://ncbi.github.io/blast-cloud/dev/api.html). 
   To use this feature, please set url_base to 'http://host.my.cloud.service.provider.com/cgi-bin/blast.cgi' and 
   format_object = 'Alignment'. For more details, please see 8. Biopython – Overview of BLAST
   
https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE = BlastDocs&DOC_TYPE = CloudBlast 
   
   Some useful parameters: 
   
   - program blastn, blastp, blastx, tblastn, or tblastx (lower case) 
   - database Which database to search against (e.g. "nr"). 
   - sequence The sequence to search. 
   - ncbi_gi TRUE/FALSE whether to give 'gi' identifier. 
   - descriptions Number of descriptions to show. Def 500. 
   - alignments Number of alignments to show. Def 500. 
   - expect An expect value cutoff. Def 10.0. 
   - matrix_name Specify an alt. matrix (PAM30, PAM70, BLOSUM80, BLOSUM45). 
   - filter "none" turns off filtering. Default no filtering 
   - format_type "HTML", "Text", "ASN.1", or "XML". Def. "XML". 
   - entrez_query Entrez query to limit Blast search 
   - hitlist_size Number of hits to return. Default 50 
   - megablast TRUE/FALSE whether to use MEga BLAST algorithm (blastn only) 
   - service plain, psi, phi, rpsblast, megablast (lower case) 
   
   This function does no checking of the validity of the parameters 
   and passes the values to the server as is. More help is available at: 
   https://ncbi.github.io/blast-cloud/dev/api.html

Zwykle argumenty funkcji qblast są w zasadzie analogiczne do różnych parametrów, które można ustawić na stronie internetowej BLAST. To sprawia, że ​​funkcja qblast jest łatwa do zrozumienia, a także zmniejsza krzywą uczenia się, aby z niej korzystać.

Łączenie i wyszukiwanie

Aby zrozumieć proces łączenia i wyszukiwania wersji online BLAST, przeprowadźmy proste wyszukiwanie sekwencji (dostępne w naszym lokalnym pliku sekwencji) na serwerze online BLAST poprzez Biopython.

Step 1 - Utwórz plik o nazwie blast_example.fasta w katalogu Biopython i podaj poniższe informacje o sekwencji jako dane wejściowe

Example of a single sequence in FASTA/Pearson format: 
>sequence A ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc 

>sequence B ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattca
tattctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc

Step 2 - Importuj moduł NCBIWWW.

>>> from Bio.Blast import NCBIWWW

Step 3 - Otwórz plik sekwencji, blast_example.fasta przy użyciu modułu python IO.

>>> sequence_data = open("blast_example.fasta").read() 
>>> sequence_data 
'Example of a single sequence in FASTA/Pearson format:\n\n\n> sequence 
A\nggtaagtcctctagtacaaacacccccaatattgtgatataattaaaatt 
atattcatat\ntctgttgccagaaaaaacacttttaggctatattagagccatcttctttg aagcgttgtc\n\n'

Step 4- Teraz wywołaj funkcję qblast przekazującą dane sekwencji jako główny parametr. Drugi parametr reprezentuje bazę danych (nt) i program wewnętrzny (blastn).

>>> result_handle = NCBIWWW.qblast("blastn", "nt", sequence_data) 
>>> result_handle 
<_io.StringIO object at 0x000001EC9FAA4558>

blast_resultszawiera wynik naszego wyszukiwania. Można go zapisać w pliku do późniejszego wykorzystania, a także przeanalizować, aby uzyskać szczegółowe informacje. Dowiemy się, jak to zrobić w następnej sekcji.

Step 5 - Tę samą funkcjonalność można wykonać za pomocą obiektu Seq, zamiast używać całego pliku fasta, jak pokazano poniżej -

>>> from Bio import SeqIO 
>>> seq_record = next(SeqIO.parse(open('blast_example.fasta'),'fasta')) 
>>> seq_record.id 
'sequence' 
>>> seq_record.seq 
Seq('ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatat...gtc', 
SingleLetterAlphabet())

Teraz wywołaj funkcję qblast przekazując obiekt Seq, record.seq jako główny parametr.

>>> result_handle = NCBIWWW.qblast("blastn", "nt", seq_record.seq) 
>>> print(result_handle) 
<_io.StringIO object at 0x000001EC9FAA4558>

BLAST automatycznie przypisze identyfikator do twojej sekwencji.

Step 6 - obiekt result_handle będzie miał cały wynik i można go zapisać w pliku do późniejszego wykorzystania.

>>> with open('results.xml', 'w') as save_file: 
>>>   blast_results = result_handle.read() 
>>>   save_file.write(blast_results)

Zobaczymy, jak przeanalizować plik wynikowy w dalszej części.

Uruchamianie samodzielnego BLASTa

Ta sekcja wyjaśnia jak uruchomić BLAST w systemie lokalnym. Jeśli uruchomisz BLAST w systemie lokalnym, może być szybszy, a także pozwala na stworzenie własnej bazy danych do przeszukiwania sekwencji.

Łączenie BLAST

Ogólnie, lokalne uruchamianie BLAST-a nie jest zalecane ze względu na jego duży rozmiar, dodatkowy wysiłek potrzebny do uruchomienia oprogramowania i związane z tym koszty. Online BLAST jest wystarczające do podstawowych i zaawansowanych celów. Oczywiście czasami może być konieczne zainstalowanie go lokalnie.

Weź pod uwagę, że często wyszukujesz w Internecie, co może wymagać dużo czasu i dużego obciążenia sieci, a jeśli masz zastrzeżone dane sekwencyjne lub problemy związane z IP, zaleca się zainstalowanie go lokalnie.

Aby to zrobić, musimy wykonać poniższe kroki -

Step 1- Pobierz i zainstaluj najnowszy plik binarny blast za pomocą podanego linku - ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

Step 2- Pobierz i rozpakuj najnowszą i niezbędną bazę danych, korzystając z poniższego linku - ftp://ftp.ncbi.nlm.nih.gov/blast/db/

Oprogramowanie BLAST zapewnia wiele baz danych w swojej witrynie. Miejmy pobrać alu.n.gz plik ze strony bazy podmuch i rozpakować go do folderu alu. Ten plik jest w formacie FASTA. Aby użyć tego pliku w naszej aplikacji do obróbki strumieniowej, musimy najpierw przekonwertować plik z formatu FASTA na format bazy danych do obróbki strumieniowej. BLAST zapewnia aplikację makeblastdb do wykonania tej konwersji.

Użyj poniższego fragmentu kodu -

cd /path/to/alu 
makeblastdb -in alu.n -parse_seqids -dbtype nucl -out alun

Uruchomienie powyższego kodu sparsuje plik wejściowy alu.n i utworzy bazę danych BLAST jako wiele plików alun.nsq, alun.nsi itd. Teraz możemy przeszukać tę bazę danych, aby znaleźć sekwencję.

Zainstalowaliśmy BLAST na naszym lokalnym serwerze, a także mamy przykładową bazę danych BLAST, alun zapytać przeciwko temu.

Step 3- Utwórzmy przykładowy plik sekwencji, aby wysłać zapytanie do bazy danych. Utwórz plik search.fsa i umieść w nim poniższe dane.

>gnl|alu|Z15030_HSAL001056 (Alu-J) 
AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCT 
TGAGCCTAGGAGTTTGCGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAA 
AGAAAAAAAAAATAGCTCTGCTGGTGGTGCATGCCTATAGTCTCAGCTATCAGGAGGCTG 
GGACAGGAGGATCACTTGGGCCCGGGAGTTGAGGCTGTGGTGAGCCACGATCACACCACT 
GCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAACAAATAA 
>gnl|alu|D00596_HSAL003180 (Alu-Sx) 
AGCCAGGTGTGGTGGCTCACGCCTGTAATCCCACCGCTTTGGGAGGCTGAGTCAGATCAC 
CTGAGGTTAGGAATTTGGGACCAGCCTGGCCAACATGGCGACACCCCAGTCTCTACTAAT 
AACACAAAAAATTAGCCAGGTGTGCTGGTGCATGTCTGTAATCCCAGCTACTCAGGAGGC 
TGAGGCATGAGAATTGCTCACGAGGCGGAGGTTGTAGTGAGCTGAGATCGTGGCACTGTA
CTCCAGCCTGGCGACAGAGGGAGAACCCATGTCAAAAACAAAAAAAGACACCACCAAAGG 
TCAAAGCATA 
>gnl|alu|X55502_HSAL000745 (Alu-J) 
TGCCTTCCCCATCTGTAATTCTGGCACTTGGGGAGTCCAAGGCAGGATGATCACTTATGC 
CCAAGGAATTTGAGTACCAAGCCTGGGCAATATAACAAGGCCCTGTTTCTACAAAAACTT 
TAAACAATTAGCCAGGTGTGGTGGTGCGTGCCTGTGTCCAGCTACTCAGGAAGCTGAGGC 
AAGAGCTTGAGGCTACAGTGAGCTGTGTTCCACCATGGTGCTCCAGCCTGGGTGACAGGG 
CAAGACCCTGTCAAAAGAAAGGAAGAAAGAACGGAAGGAAAGAAGGAAAGAAACAAGGAG 
AG

Dane o sekwencji są zbierane z pliku alu.n; dlatego pasuje do naszej bazy danych.

Step 4 - Oprogramowanie BLAST udostępnia wiele aplikacji do przeszukiwania bazy danych i używamy blastn. blastn application requires minimum of three arguments, db, query and out. db odnosi się do bazy danych przeciwko wyszukiwaniu; query to sekwencja do dopasowania i outto plik do przechowywania wyników. Teraz uruchom poniższe polecenie, aby wykonać to proste zapytanie -

blastn -db alun -query search.fsa -out results.xml -outfmt 5

Uruchomienie powyższego polecenia spowoduje wyszukanie i wyświetlenie danych wyjściowych w results.xml plik jak podano poniżej (częściowo dane) -

<?xml version = "1.0"?> 
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" 
   "http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd">
<BlastOutput> 
   <BlastOutput_program>blastn</BlastOutput_program> 
   <BlastOutput_version>BLASTN 2.7.1+</BlastOutput_version> 
   <BlastOutput_reference>Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb 
      Miller (2000), "A greedy algorithm for aligning DNA sequences", J 
      Comput Biol 2000; 7(1-2):203-14.
   </BlastOutput_reference> 
   
   <BlastOutput_db>alun</BlastOutput_db> 
   <BlastOutput_query-ID>Query_1</BlastOutput_query-ID> 
   <BlastOutput_query-def>gnl|alu|Z15030_HSAL001056 (Alu-J)</BlastOutput_query-def>
   <BlastOutput_query-len>292</BlastOutput_query-len> 
   <BlastOutput_param> 
      <Parameters> 
         <Parameters_expect>10</Parameters_expect> 
         <Parameters_sc-match>1</Parameters_sc-match> 
         <Parameters_sc-mismatch>-2</Parameters_sc-mismatch> 
         <Parameters_gap-open>0</Parameters_gap-open> 
         <Parameters_gap-extend>0</Parameters_gap-extend> 
         <Parameters_filter>L;m;</Parameters_filter> 
      </Parameters> 
   </BlastOutput_param> 
   <BlastOutput_iterations> 
      <Iteration> 
         <Iteration_iter-num>1</Iteration_iter-num><Iteration_query-ID>Query_1</Iteration_query-ID> 
         <Iteration_query-def>gnl|alu|Z15030_HSAL001056 (Alu-J)</Iteration_query-def> 
         <Iteration_query-len>292</Iteration_query-len> 
         <Iteration_hits> 
         <Hit>
            <Hit_num>1</Hit_num> 
            <Hit_id>gnl|alu|Z15030_HSAL001056</Hit_id> 
            <Hit_def>(Alu-J)</Hit_def> 
            <Hit_accession>Z15030_HSAL001056</Hit_accession> 
            <Hit_len>292</Hit_len>
            <Hit_hsps> 
               <Hsp>
                 <Hsp_num>1</Hsp_num> 
                  <Hsp_bit-score>540.342</Hsp_bit-score> 
                  <Hsp_score>292</Hsp_score>
                  <Hsp_evalue>4.55414e-156</Hsp_evalue> 
                  <Hsp_query-from>1</Hsp_query-from>
                  <Hsp_query-to>292</Hsp_query-to> 
                  <Hsp_hit-from>1</Hsp_hit-from> 
                  <Hsp_hit-to>292</Hsp_hit-to> 
                  <Hsp_query-frame>1</Hsp_query-frame>
                  <Hsp_hit-frame>1</Hsp_hit-frame>
                  <Hsp_identity>292</Hsp_identity>
                  <Hsp_positive>292</Hsp_positive> 
                  <Hsp_gaps>0</Hsp_gaps> 
                  <Hsp_align-len>292</Hsp_align-len>
                  
                  <Hsp_qseq>
                     AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCTTGAGCCTAGGAGTTTG
                     CGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAAAGAAAAAAAAAATAGCTCTGCTGGTGGTGCATG
                     CCTATAGTCTCAGCTATCAGGAGGCTGGGACAGGAGGATCACTTGGGCCCGGGAGTTGAGGCTGTGGTGAGCC
                     ACGATCACACCACTGCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAACAAATAA
                  </Hsp_qseq> 

                  <Hsp_hseq>
                     AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCTTGAGCCTAGGA
                     GTTTGCGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAAAGAAAAAAAAAATAGCTCTGCT
                     GGTGGTGCATGCCTATAGTCTCAGCTATCAGGAGGCTGGGACAGGAGGATCACTTGGGCCCGGGAGTTGAGG
                     CTGTGGTGAGCCACGATCACACCACTGCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAAC
                     AAATAA
                  </Hsp_hseq>

                  <Hsp_midline>
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||
                  </Hsp_midline>
               </Hsp> 
            </Hit_hsps>
         </Hit>
         ......................... 
         ......................... 
         ......................... 
         </Iteration_hits> 
         <Iteration_stat> 
            <Statistics> 
               <Statistics_db-num>327</Statistics_db-num> 
               <Statistics_db-len>80506</Statistics_db-len> 
               <Statistics_hsp-lenv16</Statistics_hsp-len> 
               <Statistics_eff-space>21528364</Statistics_eff-space> 
               <Statistics_kappa>0.46</Statistics_kappa> 
               <Statistics_lambda>1.28</Statistics_lambda> 
               <Statistics_entropy>0.85</Statistics_entropy>
            </Statistics>
         </Iteration_stat>
      </Iteration> 
   </BlastOutput_iterations>
</BlastOutput>

Powyższe polecenie można uruchomić w Pythonie za pomocą poniższego kodu -

>>> from Bio.Blast.Applications import NcbiblastnCommandline 
>>> blastn_cline = NcbiblastnCommandline(query = "search.fasta", db = "alun", 
outfmt = 5, out = "results.xml") 
>>> stdout, stderr = blastn_cline()

Tutaj pierwszy z nich jest uchwytem do wyjścia wybuchu, a drugi jest możliwym wyjściem błędu wygenerowanym przez polecenie wybuchu.

Ponieważ dostarczyliśmy plik wyjściowy jako argument wiersza poleceń (out = „results.xml”) i ustawiliśmy format wyjściowy na XML (outfmt = 5), plik wyjściowy zostanie zapisany w bieżącym katalogu roboczym.

Analiza wyniku BLAST

Generalnie wynik BLAST jest analizowany jako format XML przy użyciu modułu NCBIXML. Aby to zrobić, musimy zaimportować następujący moduł -

>>> from Bio.Blast import NCBIXML

Teraz, open the file directly using python open method i use NCBIXML parse method jak podano poniżej -

>>> E_VALUE_THRESH = 1e-20 
>>> for record in NCBIXML.parse(open("results.xml")): 
>>>     if record.alignments: 
>>>        print("\n") 
>>>        print("query: %s" % record.query[:100]) 
>>>        for align in record.alignments: 
>>>           for hsp in align.hsps: 
>>>              if hsp.expect < E_VALUE_THRESH: 
>>>                 print("match: %s " % align.title[:100])

To da następujący wynik -

query: gnl|alu|Z15030_HSAL001056 (Alu-J) 
match: gnl|alu|Z15030_HSAL001056 (Alu-J) 
match: gnl|alu|L12964_HSAL003860 (Alu-J) 
match: gnl|alu|L13042_HSAL003863 (Alu-FLA?) 
match: gnl|alu|M86249_HSAL001462 (Alu-FLA?) 
match: gnl|alu|M29484_HSAL002265 (Alu-J) 

query: gnl|alu|D00596_HSAL003180 (Alu-Sx) 
match: gnl|alu|D00596_HSAL003180 (Alu-Sx) 
match: gnl|alu|J03071_HSAL001860 (Alu-J) 
match: gnl|alu|X72409_HSAL005025 (Alu-Sx) 

query: gnl|alu|X55502_HSAL000745 (Alu-J) 
match: gnl|alu|X55502_HSAL000745 (Alu-J)