Biopython: descripción general de BLAST

BLAST significa Basic Local Alignment Search Tool. Encuentra regiones de similitud entre secuencias biológicas. Biopython proporciona el módulo Bio.Blast para hacer frente a la operación NCBI BLAST. Puede ejecutar BLAST en una conexión local o en una conexión a Internet.

Entendamos estas dos conexiones brevemente en la siguiente sección:

Corriendo por Internet

Biopython proporciona el módulo Bio.Blast.NCBIWWW para llamar a la versión en línea de BLAST. Para hacer esto, necesitamos importar el siguiente módulo:

>>> from Bio.Blast import NCBIWWW

El módulo NCBIWW proporciona la función qblast para consultar la versión en línea de BLAST, https://blast.ncbi.nlm.nih.gov/Blast.cgi. qblast admite todos los parámetros admitidos por la versión en línea.

Para obtener ayuda sobre este módulo, use el siguiente comando y comprenda las características:

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

Por lo general, los argumentos de la función qblast son básicamente análogos a los diferentes parámetros que puede establecer en la página web BLAST. Esto hace que la función qblast sea fácil de entender y reduce la curva de aprendizaje para usarla.

Conexión y búsqueda

Para comprender el proceso de conexión y búsqueda de la versión en línea de BLAST, hagamos una búsqueda de secuencia simple (disponible en nuestro archivo de secuencia local) contra el servidor BLAST en línea a través de Biopython.

Step 1 - Crea un archivo llamado blast_example.fasta en el directorio de Biopython y proporcione la siguiente información de secuencia como entrada

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

>sequence B ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattca
tattctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc

Step 2 - Importar el módulo NCBIWWW.

>>> from Bio.Blast import NCBIWWW

Step 3 - Abra el archivo de secuencia, blast_example.fasta usando el módulo 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- Ahora, llame a la función qblast pasando datos de secuencia como parámetro principal. El otro parámetro representa la base de datos (nt) y el programa interno (blastn).

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

blast_resultscontiene el resultado de nuestra búsqueda. Se puede guardar en un archivo para su uso posterior y también, analizar para obtener los detalles. Aprenderemos cómo hacerlo en la próxima sección.

Step 5 - La misma funcionalidad se puede hacer usando el objeto Seq también en lugar de usar todo el archivo fasta como se muestra a continuación -

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

Ahora, llame a la función qblast pasando el objeto Seq, record.seq como parámetro principal.

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

BLAST asignará un identificador para su secuencia automáticamente.

Step 6 - El objeto result_handle tendrá el resultado completo y se puede guardar en un archivo para su uso posterior.

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

Veremos cómo analizar el archivo de resultados en la sección posterior.

Ejecución independiente de BLAST

Esta sección explica cómo ejecutar BLAST en el sistema local. Si ejecuta BLAST en el sistema local, puede ser más rápido y también le permite crear su propia base de datos para buscar secuencias.

Conectando BLAST

En general, no se recomienda ejecutar BLAST localmente debido a su gran tamaño, el esfuerzo adicional necesario para ejecutar el software y el costo involucrado. Online BLAST es suficiente para propósitos básicos y avanzados. Por supuesto, en algún momento es posible que deba instalarlo localmente.

Considere que está realizando búsquedas frecuentes en línea, lo que puede requerir mucho tiempo y un gran volumen de red, y si tiene datos de secuencia patentados o problemas relacionados con la IP, se recomienda instalarlo localmente.

Para hacer esto, debemos seguir los pasos a continuación:

Step 1- Descargue e instale el último binario de blast usando el enlace proporcionado: ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

Step 2- Descargue y descomprima la base de datos más reciente y necesaria utilizando el siguiente enlace: ftp://ftp.ncbi.nlm.nih.gov/blast/db/

El software BLAST proporciona muchas bases de datos en su sitio. Vamos a descargamos alu.n.gz archivo desde el sitio base de datos de explosión y descomprimirlo en la carpeta de aluminio. Este archivo está en formato FASTA. Para utilizar este archivo en nuestra aplicación de explosión, primero debemos convertir el archivo del formato FASTA al formato de base de datos de explosión. BLAST proporciona la aplicación makeblastdb para realizar esta conversión.

Utilice el siguiente fragmento de código:

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

Al ejecutar el código anterior, se analizará el archivo de entrada, alu.ny creará la base de datos BLAST como varios archivos alun.nsq, alun.nsi, etc. Ahora, podemos consultar esta base de datos para encontrar la secuencia.

Hemos instalado BLAST en nuestro servidor local y también tenemos una base de datos BLAST de muestra, alun para cuestionarlo.

Step 3- Creemos un archivo de secuencia de muestra para consultar la base de datos. Cree un archivo search.fsa y coloque los siguientes datos en él.

>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

Los datos de la secuencia se recopilan del archivo alu.n; por lo tanto, coincide con nuestra base de datos.

Step 4 - El software BLAST proporciona muchas aplicaciones para buscar en la base de datos y usamos blastn. blastn application requires minimum of three arguments, db, query and out. db se refiere a la base de datos contra la búsqueda; query es la secuencia para hacer coincidir y outes el archivo para almacenar los resultados. Ahora, ejecute el siguiente comando para realizar esta consulta simple:

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

Ejecutar el comando anterior buscará y dará salida en el results.xml archivo como se indica a continuación (parcialmente datos) -

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

El comando anterior se puede ejecutar dentro de Python usando el siguiente código:

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

Aquí, el primero es un control de la salida de explosión y el segundo es la posible salida de error generada por el comando de explosión.

Dado que hemos proporcionado el archivo de salida como argumento de línea de comando (out = “results.xml”) y establecemos el formato de salida como XML (outfmt = 5), el archivo de salida se guardará en el directorio de trabajo actual.

Analizando el resultado de BLAST

Generalmente, la salida BLAST se analiza como formato XML utilizando el módulo NCBIXML. Para hacer esto, necesitamos importar el siguiente módulo:

>>> from Bio.Blast import NCBIXML

Ahora, open the file directly using python open method y use NCBIXML parse method como se indica a continuación -

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

Esto producirá una salida de la siguiente manera:

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)