就基因家族工作做一簡單介紹
基本思路
數據准備
確定好研究的基因家族后（比如：NBS，MADS-box etc.），下面就可以下載相關數據。
- 所研究物種的基因組序列； genome.fa
- 所研究物種蛋白序列；pep.fa
- 所研究物種gff文件
- 目標基因家族的隱馬科夫模型
- or RefSeq 對應基因家族的蛋白序列
對應基因組信息可根據發表文章中提供的路徑進行下載即可；
對於隱馬科夫模型可以從Pfam進行下載（比如：NBS）
點擊BROWSE后，進行搜索NBS即可；本次下載PF00931
下載RefSeq對應抗病序列如下所示；
或者下載RefSeq植物所有信息，根據注釋信息選取相關序列
見基因功能注釋
鑒定家族成員
-
基於保守結構模型進行鑒定
所用工具為hmmer
hmmbuild/hmmsearch/hmmscan/hmmalign 這幾個功能是主要用於蛋白質結構與分析和注釋的hmmer中小工具
軟件安裝
conda install -c bioconda hmmer=3.3
在鑒定基因家族時，常用到的工具是hmmsearch，里面常用的算法有三種。一般我們使用--cut_tc算法對隱馬可夫模型進行搜索，tc算法是使用pfam提供的hmm文件中trusted cutoof的值進行篩選，相對比較可靠。
以2015年發表文章Identification and distribution of the NBS-LRR gene family in the Cassava genome為例，解析其步驟。
分析流程如下圖所示
首先根據NB-ARC 模型，利用hmmersearch 篩選質量較高的基因（E-value < 1 × 10−20）。 而后使用clustalw2對高質量的序列進行多序列比對，多序列比對后，對這些置信的序列進行隱馬可夫模型的構建（使用hmmbuild），最后用新建的模型，重新對所需基因進行篩選。
第一步 篩選高質量NBS基因
## hmmersearh進行搜索
hmmsearch --cut_tc --domtblout NBS-ABC.out NBS-ARC.hmm \
Arabidopsis_thaliana.TAIR10.pep.all.fa
# --cut_tc: use profile's TC trusted cutoffs to set all thresholding
# --domtblout <f> : save parseable table of per-domain hits to file <f>
## 根據e-value 篩選高質量基因
grep -v "#" NBS-ABC.out|awk '($7 + 0) < 1E-20'|cut -f1 -d " "|sort -u > NBS-ARC_qua_id.txt
第二步：使用clustalw2進行多序列比較,構建新模型
比對具體步驟請看序列比對和構建進化樹（clustalw和phylip）
hmmbuild NBS-ARC.second.out NBS-ARC_qua.aln
hmmsearch --cut_tc --domtblout NBS-ARC.second.out \
NBS-ARC_qua.hmm ../Arabidopsis_thaliana.TAIR10.pep.all.fa
第三步：再次進行過濾
grep -v "#" NBS-ABC.second.out|awk '($7 + 0) < 1E-03' | cut -f1 -d " "|sort -u >final.NBS.list
-
基於同源注釋
簡而言之，利用已知的NBS基因作為ref，進行blastp比對
如果要取得可信度較高的結果，可講上述兩種方法得到的基因取交集即可
com m -12 gene1 gene1 >common.list
歡迎掃碼交流
