基因家族成员鉴定
就基因家族工作做一简单介绍
基本思路
数据准备
确定好研究的基因家族后(比如: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