如何在rBlast中运行blastp？蛋白比对无匹配结果问题求助

rBlast完全支持blastp蛋白序列比对，帮助文档仅标注blastn属于文档更新滞后的问题，已有大量用户成功使用该功能完成蛋白比对，你遇到的无匹配结果问题可按以下步骤排查解决：

• 路径匹配校验：你当前代码中构建数据库时使用的是变量fa（对应listF[i]遍历到的文件），但初始化blast对象时硬编码了固定的数据库文件路径，首先确认你所用的查询序列确实属于你初始化blast时指定的数据库，避免建库和比对使用的文件不一致。
• 建库有效性校验：单独运行makeblastdb语句，检查控制台是否输出建库成功的提示，确认数据库对应目录下生成了.pin、.psq、.phr三个后缀的库文件，且文件大小正常，排除输入fasta格式错误、非法字符导致的建库失败问题。
• 放宽比对阈值：predict.BLAST默认使用1e-10的e-value阈值，同时自带比对长度、一致性过滤规则，短序列或同源性较低的序列容易被过滤。可显式调低阈值测试，示例如下：

# 显式指定宽松阈值运行比对
cl <- predict(bl, seq, evalue = 1, args = "-outfmt 6")

如果调整阈值后能输出结果，再根据你的研究需求调整到合理阈值范围即可。

• 序列有效性校验：打印你读取的seq对象，确认读取到的是合法的氨基酸序列，无终止密码子、特殊符号等异常。也可将查询序列导出后用本地命令行blastp直接比对同一数据库，确认序列本身确实存在匹配，排除序列本身的问题。
• 版本更新：旧版本rBlast存在blastp参数传递bug，可升级到最新版本后重试：

# Bioconductor版本升级
BiocManager::update("rBlast")

参考可运行代码 #

# 确认数据库路径，避免遍历变量出错
db_path <- "Trich_prot_fasta/Tri5640_1_GeneModels_FilteredModels1_aa.fasta"
# 构建蛋白库时添加-parse_seqids参数保留序列ID信息
makeblastdb(db_path, dbtype = "prot", args = "-parse_seqids")
# 初始化blastp对象
bl <- blast(db = db_path, type = "blastp")
# 读取查询序列
seq <- readAAStringSet("NDRkinase/testSeq.txt")
# 放宽阈值运行比对
cl <- predict(bl, seq, evalue = 1e-5)
# 输出比对结果
print(cl)

内容的提问来源于stack exchange，提问作者SKay