短核苷酸序列的 Biopython 爆炸参数答案

【问题标题】：Biopython blast parameters for short nucleotidic sequences短核苷酸序列的 Biopython 爆炸参数
【发布时间】：2014-03-07 21:55:53
【问题描述】：

我正在尝试使用 NCBIWWW 通过 biopython 运行 blastn。
我在给定的示例文件上使用 qblast 函数。
我定义了一些方法，当我的 fasta 包含足够长的序列时，一切都像魅力一样。唯一失败的情况是当我需要爆炸来自 Illumina 测序的读数太短时。所以我想说这可能是因为提交作品时没有自动重新定义爆破参数。

我尽我所能接近blastn-short条件（参见here中的表C2）但没有任何成功。

看来我无法输入正确的参数。

我认为我越来越接近工作情况是：

result_handle = NCBIWWW.qblast("blastn", "nr",
                                fastaSequence,
                                word_size=7,
                                gapcosts='5 2',
                                nucl_reward=1,
                                nucl_penalty='-3',
                                expect=1000)

感谢您提供任何提示/建议以使其发挥作用。

我的 fasta 阅读示例如下：

>TEST 1-211670
AGACTGCGATCCGAACTGAGAAC

我得到的错误如下：

>ValueError: Error message from NCBI: Message ID#24 Error: Failed to read the Blast query: Protein FASTA provided for nucleotide sequence

当我查看this page 时，似乎我的问题在于修复阈值，但显然到目前为止我还没有成功。

感谢您的帮助。

【问题讨论】：

为什么不把这个问题放到 bioinformatics.stackexchange 上？

标签： python biopython blast ncbi

【解决方案1】：

一旦我在爆破肽方面遇到问题，并且似乎是正确选择参数的问题。我花了很长时间才弄清楚它们实际上应该是什么（各种网站上的不一致和稀缺的数据，包括在这方面的 NCBI 文档中相当复杂）。我知道您对爆破核苷酸序列感兴趣，但据说您会在查看下面的代码时找到您的解决方案。特别注意filter、composition_based_statistics、word_size 和matrix_name 等参数。就我而言，它们似乎至关重要。

blast_handle = NCBIWWW.qblast("blastp", "nr",
                              peptide_seq,
                              expect=200000.0,
                              filter=False,
                              word_size=2,
                              composition_based_statistics=False,
                              matrix_name="PAM30",
                              gapcosts="9 1",
                              hitlist_size=1000)

【讨论】：

感谢您的帮助。据我所知，矩阵类型仅在蛋白质对齐时适用；对于word_size，我应该有合适的（来自文档），而composition_based_statistics 和filter 我不知道。我都试过了，但似乎都没有帮助..

【解决方案2】：

此代码适用于我（Biopython 1.64）：

from Bio.Blast import NCBIWWW

fastaSequence = ">TEST 1-211670\nAGACTGCGATCCGAACTGAGAAC"

result_handle = NCBIWWW.qblast("blastn", "nr",
                               fastaSequence,
                               word_size=7,
                               gapcosts='5 2',
                               nucl_reward=1,
                               nucl_penalty='-3',
                               expect=1000)

print result_handle.getvalue()

也许您传递了错误的 fastaSequence。 Biopython 不会从 SeqRecords（或任何东西）到普通 FASTA 进行任何转换。您必须提供如上所示的查询。

Blast 确定序列是读取前几个字符的核苷酸还是蛋白质。如果它们在阈值以上的“ACGT”中，则为核苷酸，否则为蛋白质。因此，您的序列处于“ACGT”的 100% 阈值，不可能被解释为蛋白质。

【讨论】：