blastn & blastp 寻找同源基因

今天接到一个任务，大致内容是在一个植物的全长转录组数据中找拟南芥的三个同源基因。简简单单的描述，我的想法也很简单，直接找基因的CDS序列做blastn比对就完事了，结果却没有那么顺利…记录一下踩的坑和解决办法。

1 blastn寻找同源基因

三个基因TAIR号是AT4G28590、AT2G43010和AT2G34640，从全长转录组测序报告中，我找到了非冗余的转录本序列文件CD-hit-est.fasta，首先第一步就是本地建核酸序列库。

1	makeblastdb -in CD-hit-est.fasta -dbtype nucl -input_type fasta -out Kc

因为给的是TAIR号，所以直接去TRIR官网查找相应基因的CDS序列做比对

手动创建query gene的fa序列文件

vim RCB_cds.fna
# 手动创建fa文件
>AT4G28590
ATGAGTTTCTTCGCTGTTGCTTGCTCCGCGCCAAGATCTTCTATGCTTCTCACCGGCTTGAATTCGAGCTTCTCTGATATGCATCGCAGCCCACTATTTGTTTTCCCGGTGACTATATCATCCCGGAGCGTGAAACGCTTCGCCGCTGTTTCGTCTGATTCCGTACTAGACCCTGAATCCAAAAATCAAACTCGGTCCCGTCGCAAAAATAAGGAAGCAGTTACGCCAATTGCTGAAACCGAGAACAATGAAAAGTTTCCGACAAAGGTCCCGCGTAAATCGAAGCGTGGGCGGCGGAGTGAAGCAGACGCTGTGGAAGATTACGTGAGAAGCTCCCTCGAGCGTACTTTCTCCACCATAAAGGAGCAGAATCCGGAGGTTTTTGAGAACAAGGAGAAGGCGAATTTCATCAAAGACAGAGGCGTTGATGAAGAAGAGGAAGAAGAAGAAGAGATGGTGGTGGAAGAGGAAGATCCAGATTGGCCAGTAGATACAGACGTTGGATGGGGAATCAAAGCTTCGGAGTATTTCGATACACATCCAATCAAAAACGTGGTTGGAGATGATGGGAGTGAGATTGATTGGGAAGGTGAGATTGATGATAGTTGGGTCAAGGAGATCAATTGTTTGGAATGGGAAAGCTTTGCTTTTCATCCTAGTCCACTCGTTGTCCTTGTATTCGAGCGATACAAAAGAGCTAGTGATAACTGGAAGACATTGAAGGAGCTTGAGAAAGCTATCAAAGTTTATTGGGATGCGAAAGATCGATTACCTCCACGGGCGGTTAAGATTGACCTGAACATCGAGACAGATTTGGCATATGCTCTTAAAGCTAAGGAATGCCCACAGATTCTCTTCTTACGCGGAAACCGGATTCTGTACAGGGAGAAAGACTTTCGCACGGCGGATGAATTGGTTCATATGATTGCGCATTTCTACTATAAAGCGAAGAGGCCTTCGTGTGTCGACAAGGCTAATGTAACCCCGTACTGTTAG

blastn比对

1	blastn -query RCB_cds.fna -out RCB_blastn_Kc.out -db Kc -outfmt 6 -evalue 1e-5 -num_threads 4

RCB_blastn_Kc.out是blast的m8格式输出文件，找到匹配长度最长的（也就是第一条）subject gene id，回到非冗余转录本，找到subject gene在哪行，最后找出转录本序列。

# 找到subject gene所在行（subject gene id中有所在行数，这里验证下）
cat CD-hit-est.fasta | grep -n "Kc-zong_1-10k_transcript/10791"
# 提取序列
awk 'NR>=10719 && NR<=10720' CD-hit-est.fasta

紧接着出现一个问题：AT2G43010和AT2G34640这两个基因无法通过blastn比对找到同源序列，evalue值不管放到多宽都比对不上。

因为这个植物在NCBI上没有参考基因组，我们课题组也只测了全长转录组而没有测基因组，所以当一开始没有比对出结果的时候，我一度怀疑是这种植物压根儿就没有这俩基因，或者这个样品叶片（测序的部位）在检测的时间点就没有转录相应的基因。

本地blast找不到同源基因，我又从近缘菊科植物开始折腾，思路是如果菊科有同源基因则寻找保守结构域，设计引物将CDS区域克隆出来。至今已发表的植物基因组可以从网站Plabipd（本站网址导航栏有收录）找到，这个网站很贴心地把物种种属关系也列了出来，可以很方便地找物种学名和近缘关系。

理想很丰满现实很骨感，我从菊科一级一级往上找，直到Eudicotyledoneae（真双子叶植物分支）才用blastn比对上同源基因，而且无一例外比对上的全是十字花科（拟南芥所在科）植物，根本不算近缘物种….无奈之下试了blast的其他功能，用氨基酸序列跑了一遍blastp，然后发现菊科也有序列可以比对上了！这才打开新世界的大门

2 blastp寻找同源基因

基于翻译阅读框对去冗余的全长转录本进行CDS预测（TransDecoder软件），结果以fasta格式保存，后续我会对这个文件验证一遍，先建蛋白库做blastp比对。

1	makeblastdb -in transdecoder.pep.fa -dbtype prot -input_type fasta -out nr_Kc

可以看到只有25128个编码蛋白基因，对于基因组大小在1G左右的菊科物种来说，这个基因数量过少。因此后续还需要对全长转录组数据再跑一遍验证一下，这个是后话。

通过TAIR号在TAIR官网查找蛋白序列，创建fa文件后进行本地blastp比对

1
2

blastp -query PIF4_pep.fna -out PIF4_blastp_nr_Kc.out -db nr_Kc -outfmt 6 -evalue 1e-5 -num_threads 4
blastp -query HMR_pep.fna -out HMR_blastp_nr_Kc.out -db nr_Kc -outfmt 6 -evalue 1e-5 -num_threads 4

注意下结果文件名写清楚什么基因，用的什么方法比对，比对的什么库。这个时候再查看各自的结果文件，发现有比对结果，再回到非冗余转录本文件找对应的cds序列。操作过程都一样，这里不再赘述了。

3 总结

找三条同源基因花了一整天的时间，主要原因还是对同源序列了解不够深刻。

同源就是有共同的进化祖先，序列相似性搜索可以通过检测过高的相似性来识别同源蛋白质或基因：当两个序列的相似性超过偶然的预期时，我们推断这两个序列存在同源性。当观察到过高的相似性时，这两个序列不是独立出现的，它们起源于一个共同的祖先。

通过算法进行序列对库比对的工具，比如blast等，是通过过高相似性来减少假阳性的结果。所以通过算法在统计学上找不到库里显著的匹配项，不代表这个物种中一定没有同源基因。

从这次blastn和blastp比对结果来看，核酸序列比对可能更不容易找到同源序列。其实也好理解，生物在进化的几亿年时间里，很难保证不同物种有高相似性的核酸序列。同个氨基酸有不同密码子（简并性），也能证明蛋白质一级结构才是对生物影响最大的，蛋白质序列相同，就会有相似结构和功能。因此，蛋白质序列也就是氨基酸序列，对相似性的搜索比核酸序列要敏感的多。