关键词 blastp，Hmmer，hmmsearch，hmmscan，Pfam，NCBI CDD
本文对 4 种兰花（Apostasia shenzhenica、Phalaenopsis equestris、Dendrobium catenatum、Gastrodia elata）的基因组使用 Blastp 和 Hmmer 两个基于不同算法的软件，筛选出包含 NB-ARC 结构域的蛋白。
简要步骤 Apostasia shenzhenica、Phalaenopsis equestris、Dendrobium catenatum 的蛋白序列数据可以从 NCBI 数据库获取，accession = PRJNA310678、PRJNA389183、PRJNA262478。Gastrodia elata 的蛋白序列数据可以从 Genome WareHouse Database 获得（http://bigd.big.ac.cn/gwh/Assembly/129/show）。从 Pfam 网站获取 NB-ARC 结构域（Pfam accession = PF00931）的 seed 序列信息（PF00931_seed.txt）和 hmm 模型信息（NB-ARC.hmm）。使用 Hmmer 软件的 hmmsearch 功能，利用已知的 NB-ARC 结构域的 hmm 模型检索 4 种兰花蛋白数据库中符合模型的蛋白质，总共筛选出 254 个蛋白质。使用 Blastp 软件，以 PF00931_seed.txt 文件中序列信息为模板，比对 4 种兰花蛋白数据库中的所有蛋白质，筛选出包含 seed 序列的蛋白质，总共筛选出 247 个蛋白质。将 Blastp 和 Hmmer 的结果整合并去冗，得到 265 个蛋白质。然后根据蛋白 ID 提取蛋白序列信息，将包含 265 个蛋白的序列信息的 fasta 文件提交到 NCBI 的保守结构域数据库（CDD，Conserved Domains Database，https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi）中进行检索，查看各蛋白所包含的全部结构域信息，进一步验证筛选出蛋白的准确性。上述内容主要参考文献《Genome- Wide Analysis of the Nucleotide Binding Site Leucine-Rich Repeat Genes of Four Orchids Revealed Extremely Low Numbers of Disease Resistance Genes》，Front. Genet.，08 January 2020。
实验结果总结：
在无法使用本地检索的情况下或为了简便，可以只使用 Blastp 网页筛选包含特定结构域的蛋白，虽然结果会有损失，但损失不大。Hmmer 和 Blastp 联合结果可以丰富单一算法的结果。详细步骤 一、准备数据
1.1、通过 PF00931 可以在 Pfam 上查看 NB-ARC domain 的详细内容。Pfam 通过 9 个 NB-ARC domain 序列作为 seed 建模，利用 HMM 算法在蛋白数据库中搜索，筛选出包含与模型吻合的 domain 的蛋白，总共比对出了 51786 条序列。在 Pfam 侧边栏的 Curation & model 板块下可以下载 Pfam 比对使用的 HMM 模型文件 NB-ARC.hmm，在侧边栏的 Alignments 板块下可以下载训练使用的 seed 序列 PF00931_seed.txt。
1.2、隐马可夫模型在预测时需要已知状态转移矩阵，Pfam 使用 9 个 seed 作为训练集，计算各种氨基酸之间的转换概率，得到状态转移矩阵，存放在 NB-ARC.hmm 文件中。NB-ARC.hmm 文件的主要内容即为氨基酸变换矩阵，矩阵总共 252 行，每行代表 NB-ARC domain 中的一个氨基酸位点。
二、使用 Blastp 和 Hmmer 筛选蛋白2.1、Hmmer 软件使用的是 HMM 算法，它既支持输入蛋白质序列，将输入的蛋白质训练作为训练集得到 HMM 模型，然后利用 HMM 模型对数据库搜索比对出包含 NB-ARC 结构域的蛋白（phmmer）；还支持直接输入已训练好的 HMM 模型（hmmsearch）。这里我使用从 Pfam 网站上下载的 NB-ARC 结构域的 HMM 模型（NB-ARC.hmm）作为输入，阈值设定为 1 e − 4 1e^{-4} 1e−4。若使用 网页版 hmmer，以网站上提供的数据库检索（Restrict by Taxonomy 参数输入 4 种兰花的物种名），总共得到 135 条结果；若使用 本地 hmmer 软件，以从 NCBI 和 GWD 下载的 4 种兰花的蛋白组作为数据库进行本地检索，得到 254 条结果，相比使用网页版，结果得到了极大的丰富。
Usage: hmmsearch [options] <hmmfile> <seqdb>$ hmmsearch -o genome/Apostasia_shenzhenica/Apostasia_shenzhenica_prot.align.hmm --tblout align_result/Apostasia_prot.align.hmm -E 0.0001 NB-ARC.hmm genome/Apostasia_shenzhenica/GCA_002786265.1_ASM278626v1_protein.fasta-o 完整结果输出路径，–tblout 简要结果输出路径
2.2、使用 Blastp 进行蛋白-蛋白的序列比对，输入文件为从 Pfam 下载的 9 个 seed 序列（PF00931_seed.txt），Blastp 使用的是默认参数（阈值为 0.05，在 BLAST 的新版本中将将默认阈值下调到了 0.05）。若使用 网页版 Blastp，蛋白数据库使用的是 Non-redundant protein sequences (nr)，物种的搜索范围为 4 种兰花，总共得到 235 个结果。（先下载 Blastp 的比对结果，格式为 Hit Table (csv)，去冗余后得到 235 个结果。若使用本地 blastp 软件（blastp -h 为简要说明，blastp -help 为详细说明），以从 NCBI 和 GWD 下载的 4 种兰花的蛋白组作为数据库进行本地检索，总共得到 247 非冗余结果。
由于 Gastrodia elata 的蛋白 fasta 文件中 protein_id 无法被 Blastp 识别，需先使用脚本更改 Gastrodia elata 的蛋白 fasta 文件中 protein_id；
将 Phalaenopsis equestris 的蛋白组数据（GCF_001263595.1_ASM126359v1_protein.fasta）构建为 blastp 可识别的蛋白数据库（-dbtype prot）（以 Phalaenopsis equestris 为例）
$ makeblastdb -in GCF_001263595.1_ASM126359v1_protein.fasta -dbtype prot -out Phalaenopsis_equestris_prot -parse_seqids利用 blastp 将文件（PF00931_seed.txt）中的序列逐个比对到 Phalaenopsis equestris 的蛋白数据库上
$ blastp -db Phalaenopsis_equestris_prot -query PF00931_seed.txt -out Phalaenopsis_equestris_prot.align.blastp -outfmt 6 -evalue 0.05PF00931_seed.txt 文件中序列里的 “.” 表示占位符，不代表任何氨基酸，真实的蛋白序列是不包含 “.” 的序列，blastp 会自动过滤掉序列中的 “.” 。
可以发现，使用下载的兰花蛋白组构建本地数据库的优点在于：
保证数据库是最新最全的，使结果更为准确、丰富；保证匹配到的蛋白 ID 一致，在使用网站匹配的结果中 Blastp 的蛋白 ID 与 Hmmer 不一致，且我没有找到转换方法。注意！
参考论文将 HMM 算法和 BLASTP 算法得到的结果合并去冗余后，使用 Pfam 主页（pfam.xfam.org）顶栏的 SEARCH 板块中 Batch search 进行线上批量分析，E-value=1e-4，查看去冗结果中各蛋白所命中的结构域，然后筛选出命中NB-ARC结构域的蛋白。上述步骤与 2.1 是冗余的！因为 Pfam 挖掘序列中结构域使用的是 Hmmer 软件的 hmmscan 功能（www.ebi.ac.uk/Tools/hmmer/search/hmmscan），将输入的蛋白质序列应用 Pfam 库中所有结构域的 HMM 模型，匹配出蛋白质序列中符合 E-value 的结构域，然后再筛选出命中 NB-ARC 的蛋白质。这与直接用 NB-ARC 的 hmm 模型遍历所有蛋白，筛选出命中的蛋白质本质是一样的，所以结果也是一样的。另一方面，由于存在于blastp结果中但没有出现在 hmmer结果中的蛋白无论如何也不会在 hmmscan 中命中 NB-ARC，所以论文作者实际上通过网上比对又将 2.2 中 blastp 的结果删去了，逻辑十分不合理。我认为可以解释的理由为：作者并不理解 “使用 Pfam 在线分析结构域” 的原理。综上，我没有执行论文中 “使用Pfam在线分析结构域” 的操作。
PS：一开始我假设作者是考虑到可能存在蛋白虽然命中了 NB-ARC，但此区段同时还命中了与 NB-ARC 相似的结构域，且NB-ARC 不是最佳命中。即有些蛋白质因为包含与 NB-ARC 相似的结构域而被 NB-ARC 的 HMM 模型筛选出来。上述情况的排除只有通过将 2.1 筛选出来的蛋白质，利用 hmmscan 与 Pfam 库中全部 HMM 模型比对后，才能实现。如果符合我的假设，那么 “使用Pfam在线分析结构域” 步骤应该在 “使用BLASTP算法挖掘包含NB-ARC结构域” 步骤之前。所以根据论文步骤，作者应该没有考虑到上述情况。
2.3、将 Hmmer 和 Blastp 结果合并去冗余后得到 265 个结果。可以发现，其与 Blastp 网页结果（235）相比增加了 12.8%，与 Blastp 本地结果（247）相比增加了 7.3%，与 Hmmer 本地结果（254）相比增加了 4.3%。在无法使用本地检索的情况下或为了简便，可以只使用 Blastp 网页筛选包含特定结构域的蛋白，虽然结果会有损失，但损失不大。
3、根据筛选出的蛋白 ID 提取蛋白的序列信息。将蛋白序列信息上传至 NCBI 的保守域数据库（Conserved Domains Database，CDD）（https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi），分析每个蛋白所带有的全部结构域，进一步验证筛选出蛋白的准确性。总共 265 条序列中只有 1 条序列没有找到任何结构域，剩下 264 条序列都命中了 NB-ARC 结构域，说明 Hmmer 和 Blastp 联合结果可以丰富单一算法的结果。
PS：CDD 是综合了 NCBI-curated、Pfam、SMART、COG、PRK 和 TIGRFAMs 六个数据库得到的，总共包含了 55570 个结构域的模型信息（PSSMs），远远大于 Pfam-A（V32.0）中包含的 16212 个结构域模型。