1

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

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

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

我认为我越接近工作情况，情况如下：

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

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

我的快速阅读示例如下：

>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

当我查看此页面时，似乎我的问题是关于修复阈值，但显然我到目前为止还没有成功。

感谢您的任何帮助。

4

2 回答 2

0

这段代码对我有用（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% 阈值，不可能被解释为蛋白质。

于 2014-07-02T20:01:12.487 回答
0

一旦我在爆破肽方面遇到问题，似乎这是一个正确选择参数的问题。我花了很长时间才弄清楚它们实际上应该是什么（各种网站上的不一致和稀缺的数据，包括在这方面的 NCBI 文档中相当复杂）。我知道您对爆破核苷酸序列感兴趣，但据说您会在查看下面的代码时找到您的解决方案。特别注意参数为filter,composition_based_statistics和. 就我而言，它们似乎至关重要。word_sizematrix_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)
于 2014-02-11T21:28:13.650 回答