基因家族分析
1 介绍
· 基因家族(Gene family),是来源于同一个祖先,由一个基因通过基因复制或者加 倍而产生两个或更多的拷贝而构成的一组基因,它们在结构和功能上具有明显的相似性,编 码相似的蛋白质产物。
基因家族的遗传进化
· 协同进化( Concerted Evolution ) 两个相互作用的物种在进化过程中发展的相互适应的共同进化。一个物种由于另一物种影响而发生遗传进化的进化类 型。例如植物由于病原菌所施加的压力而与抗性基因表现出协同进化关系。
· 无功能化( Degeneration ) 由于有害突变的(非同义突变,可变剪切突变等等)不断积累,导致基因功能丧失,例如一些假基因。
· 新功能化( Subfunctionalization ) 基因在复制的过程中通过突变,遗传漂变等等使得一些基因有了新的功能。
2 分析流程
首先获得想要预测基因家族的基因以及序列,通过pfam
获得已知蛋白保守结构域的隐马尔科夫模型,通过hmmsearch
构建已知蛋白序列的结构域模型,初步筛选相关的家族基因,之后通过bedtools
的getfasta
获得初步筛选的蛋白序列fasta信息,合并到一个fasta文件中并进入CDD
或者pfam
按照关联结构域对家族蛋白进行进一步筛选。
在获得了我们想要进行分析的家族蛋白之后我们就要开始进行后续的构树和可视化等分析,首先我们需要使用mega
对这一组蛋白序列进行clustalw
多序列比对 ,随后构建系统发育树,构树方法有邻接法(Neighbor-Joining, NJ),最大似然法(Maximum likelihood,ML),最大简约法(Maximum parsimony,MP)和贝叶斯法(Bayesian inference, BI)不同方法各有优劣。在获得进化树之后可以使用在线工具EvolView
对进化树进行美化,添加一些参数使其展示的内容更加丰富。
得到了进化树我们就需要对其中的家族蛋白进行基因结构分析,这里我主要使用TBtools
,首先下载对应物种的基因组fasta文件以及基因组注释gtf/gff3文件导入TBtools
,
之后使用Gene Structure Shower
模块根据自己喜好构建并修改基因结构图。(很厉害的软件,基本上图中所有东西都可以按照自己兴趣调节)
得到基因结构图以后我们就要开始对基因家族的motif进行分析,
• motif是蛋白质分子具有特定功能的或者作为一个独立结构域一部分相近的二级结构聚合体;
• 基因保守域结构主要通过MEME在线网站分析(http://meme-suite.org/);也可以通过SMART、MOTIF Search;
• MEME在线网站进行两种分析: 1、MEME分析(重在发掘新的motif,比较敏感,是由序列查找motif的过程); 2、Mast分析(重在确定motif的存在,较全面,是由motif查找domain的过程);
• MEME输出的图形无法直接导出,只能通过截图软件进行,得到图的分辨率不佳 • TBtools能识别MEME网站中输出的XML文件,并能导出矢量图,提高图形质量。
得到进化树以及MEME图之后我们可以对其进行合并展示。如果之后我们想要再进一步进行分析的话,我们可以从基因家族的染色体定位入手,绘制基因家族在染色体中的定位图并且可以分析其上下游的调控元件或者判断是否有选择压力作用于这个蛋白质的编码基因,从而了解基因进化速度。
染色体定位可以根据测序文件,确定各条染色体大小,基因存 在于那条染色体上,并确定该基因在染色体的哪个区域, 使用map gene to chromosome( http://mg2c.iask.in/mg2c_v2.0/)等软件进行绘制。
调控元件的分析可以使用TBtools
通过gff3/gtf注释文件使用Gtf/Gff3 Sequence Extractor
模块获取上下游指定长度bp的序列并认定这其中可能包括有该基因的调控元件并进行进一步分析。
3 结果展示
3.1 构建模型
首先获取目的基因所在物种的参考蛋白序列以及目的基因蛋白的隐马尔科夫模型,可在ncbi
和Pfam
获得,不再过多赘述,之后构建家族模型。
#一般寻找基因家族,都可以通过保守结构域来预测,从而找到物种的某一基因家族,从#而进行之后的分析。 这里就需要用到HMMER,来鉴定物种某一基因家族。
#在鉴定基因家族时,常用到的工具是hmmsearch,里面常用的算法有三种。一般我们使#用--cut_tc算法对隐 马可夫模型进行搜索,tc算法是使用pfam提供的hmm文件中trusted #cutoof的值进行筛选,相对比较可靠。
hmmsearch Hsp20.hmm protein.fa > out
##Nramp.hmm 是上一步下载到的模型文件 protein.fa是全基因组蛋白序列文件,out是重##定向的输出的文件,把质量比较高的基因家族候选基因筛选出来E-value < 1 × 10e20
out输出文件
根据初步结果获得蛋白我们获取其fasta序列,根据Sequence
列获得蛋白序列id我们使用bedtools
等工具获取其fasta序列,合并并保存在一个fasta序列中。
随后可以使用ncbi的CD-search工具鉴定蛋白质或者核酸序列内的保守结构域或功能单位 或者使用Pfam
的Batch search
。
主要看
from
和to
两列,代表了关联结构域的长度,按照一定阈值筛选,这里排除了500个氨基酸以下的蛋白。
3.2 构建系统发育树
导入前面得到的fasta文件之后,使用Align by clustalW
对fasta中的蛋白序列进行多序列比对(ctrl+a选取所有网格使用默认参数即可),之后导出并存为.meg格式。
之后将.meg格式重新导入回mega中,使用最大似然法构建系统发育树(使用默认参数即可)。
构树选项针对导出结果还可以使用Evolview
美化发育树结果。