4
$\begingroup$

I have found a blog post with a script that I would like to use for my current research project: link

The script is incredibly fast and produces a smooth conservation plot. In the blog post, the author mentions that it would be totally possible to change the remote blast function to a local database. Instead of making a local database, I have opted to try to change the function to still use a remote blast, but with different parameters.

I found the following resource for learning about the NCBIWWW module: link

Using this, I rewrote the line in the script:

blast_handle = NCBIWWW.qblast(blast_method, "nr", search_gi)

to:

blast_handle = NCBIWWW.qblast(program="blastp", 
                                      database="refseq_protein", 
                                      sequence=search_gi, 
                                      entrez_query=" txid2 [ORGN] OR txid4751 [ORGN] ", 
                                      expect=20000.0, 
                                      alignments=10000, 
                                      descriptions=10000)

Strangely, I get the exact same graph as an output even with the different parameters. Am I doing something about this search incorrectly?

Thanks.

Update:

I got my BLAST to match the webpage. my current parameters are:

blast_handle = NCBIWWW.qblast(program="blastp", database="refseq_protein", sequence="P02686.3", word_size=6, expect=20000.0, hitlist_size=10000, gapcosts="11 1", matrix_name="BLOSUM62", filter="F", genetic_code=1, threshold=21, composition_based_statistics=2, alignments=10000, descriptions=10000, entrez_query=" txid2 [ORGN] OR txid4751 [ORGN] ")

With this new blast_handle, the script no longer works. I get the following error codes:

line 191, in main(sys.argv[1])

line 27, in main blast_rec = ncbi_manager.remote_blast(protein_gi, "blastp")

line 184, in remote_blast return rec_it.next()

line 575, in parse raise ValueError("Your XML file was empty")

ValueError: Your XML file was empty

$\endgroup$
  • 2
    $\begingroup$ Have you checked the hits you're getting back from BLAST in your different methods? If your query matches the same few high-scoring sequences regardless of algorithm, you'll then be creating an MSA of exactly the same sequences each time (maybe +/- some random seeding noise) which would mean your graph was the same. $\endgroup$ – Joe Healey Jul 18 '17 at 21:30
  • $\begingroup$ I just checked the edited BLAST with my intended parameters and compared it to the results using the NCBI web BLAST with the same parameters. The xml file that I read the qblast into only showed 18 hits, while the web page found 1088 hits. Clearly something isn't working... I'm going to keep reading up on this. $\endgroup$ – bluescholar1212 Jul 19 '17 at 23:31
  • $\begingroup$ update: added hitlist_size parameter to my blast call. Still not getting same results as the blast from the web page... $\endgroup$ – bluescholar1212 Jul 20 '17 at 16:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.