Учебники

Биопайтон — обзор BLAST

BLAST расшифровывается как Basic Local Alignment Tool . Он находит области сходства между биологическими последовательностями. Biopython предоставляет модуль Bio.Blast для работы с NCBI BLAST. Вы можете запустить BLAST в локальном или интернет-соединении.

Давайте кратко разберемся в этих двух связях в следующем разделе —

Работает через интернет

Biopython предоставляет модуль Bio.Blast.NCBIWWW для вызова онлайн-версии BLAST. Для этого нам нужно импортировать следующий модуль —

>>> from Bio.Blast import NCBIWWW

Модуль NCBIWW предоставляет функцию qblast для запроса онлайн-версии BLAST, https://blast.ncbi.nlm.nih.gov/Blast.cgi . qblast поддерживает все параметры, поддерживаемые онлайн-версией.

Чтобы получить какую-либо справку об этом модуле, используйте приведенную ниже команду и поймите функции —

>>> 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

Обычно аргументы функции qblast в основном аналогичны различным параметрам, которые вы можете установить на веб-странице BLAST. Это облегчает понимание функции qblast, а также сокращает кривую обучения ее использованию.

Подключение и поиск

Чтобы понять процесс подключения и поиска в онлайн-версии BLAST, давайте выполним простой поиск последовательностей (доступен в нашем локальном файле последовательностей) на онлайн-сервере BLAST через Biopython.

Шаг 1 — Создайте файл с именем blast_example.fasta в каталоге Biopython и передайте приведенную ниже информацию о последовательности в качестве входных данных.

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

>sequence B ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattca
tattctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc

Шаг 2 — Импортируйте модуль NCBIWWW.

>>> from Bio.Blast import NCBIWWW

Шаг 3 — Откройте файл последовательности blast_example.fasta с помощью модуля ввода-вывода Python.

>>> 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'

Шаг 4 — Теперь вызовите функцию qblast, передавая данные последовательности в качестве основного параметра. Другой параметр представляет базу данных (nt) и внутреннюю программу (blastn).

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

blast_results содержит результат нашего поиска. Он может быть сохранен в файл для последующего использования, а также проанализирован для получения подробной информации. Мы узнаем, как это сделать, в следующем разделе.

Шаг 5 — Те же функциональные возможности могут быть выполнены с использованием объекта Seq, а не с использованием всего файла fasta, как показано ниже —

>>> 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())

Теперь вызовите функцию qblast, передав объект Seq, record.seq в качестве основного параметра.

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

BLAST автоматически назначит идентификатор вашей последовательности.

Шаг 6 — объект result_handle будет иметь весь результат и может быть сохранен в файл для дальнейшего использования.

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

Мы увидим, как разобрать файл результатов в следующем разделе.

Запуск автономного BLAST

В этом разделе объясняется, как запустить BLAST в локальной системе. Если вы запускаете BLAST в локальной системе, это может быть быстрее, а также позволяет создавать собственную базу данных для поиска по последовательностям.

Подключение BLAST

В общем, запуск BLAST локально не рекомендуется из-за его большого размера, дополнительных усилий, необходимых для запуска программного обеспечения, и связанных с этим затрат. Онлайн BLAST достаточно для базовых и продвинутых целей. Конечно, иногда вам может потребоваться установить его локально.

Предположим, что вы проводите частые поиски в Интернете, что может потребовать много времени и большого объема сети, и если у вас имеются запатентованные данные последовательности или проблемы, связанные с IP, рекомендуется установить его локально.

Для этого нам нужно выполнить следующие шаги —

Шаг 1 — Загрузите и установите последний двоичный файл Blast, используя данную ссылку — ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

Шаг 2 — Загрузите и распакуйте последнюю и необходимую базу данных, используя ссылку ниже — ftp://ftp.ncbi.nlm.nih.gov/blast/db/

Программное обеспечение BLAST предоставляет множество баз данных на своем сайте. Давайте загрузим файл alu.n.gz с сайта базы данных Blast и распакуем его в папку alu. Этот файл в формате FASTA. Чтобы использовать этот файл в нашем приложении Blast, нам нужно сначала преобразовать файл из формата FASTA в формат базы данных Blast. BLAST предоставляет приложение makeblastdb для этого преобразования.

Используйте приведенный ниже фрагмент кода —

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

Выполнение приведенного выше кода проанализирует входной файл alu.n и создаст базу данных BLAST как несколько файлов alun.nsq, alun.nsi и т. Д. Теперь мы можем запросить эту базу данных, чтобы найти последовательность.

Мы установили BLAST на наш локальный сервер, а также имеем образец базы данных BLAST, позволяющий выполнять запросы к нему.

Шаг 3 — Давайте создадим пример файла последовательности для запроса базы данных. Создайте файл search.fsa и поместите в него приведенные ниже данные.

>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

Данные о последовательности собраны из файла alu.n; следовательно, это соответствует нашей базе данных.

Шаг 4 — Программное обеспечение BLAST предоставляет множество приложений для поиска в базе данных, и мы используем blastn. Приложение blastn требует минимум трех аргументов: db, query и out. дБ относится к базе данных против поиска; запрос — это последовательность для сравнения, а файл для сохранения результатов. Теперь выполните следующую команду для выполнения этого простого запроса:

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

Выполнение вышеуказанной команды произведет поиск и выдаст вывод в файле results.xml, как показано ниже (частично данные) —

<?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>

Вышеприведенную команду можно запустить внутри питона, используя следующий код:

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

Здесь первый — дескриптор взрыва, а второй — возможный вывод ошибки, сгенерированный командой взрыва.

Поскольку мы предоставили выходной файл в качестве аргумента командной строки (out = «results.xml») и установили выходной формат как XML (outfmt = 5), выходной файл будет сохранен в текущем рабочем каталоге.

Разбор BLAST Результат

Обычно вывод BLAST анализируется в формате XML с использованием модуля NCBIXML. Для этого нам нужно импортировать следующий модуль —

>>> from Bio.Blast import NCBIXML

Теперь откройте файл напрямую с помощью метода открытия Python и используйте метод разбора NCBIXML, как показано ниже —

>>> 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])

Это произведет вывод следующим образом —