基因家族生信分析

1、页眉内容基因家族生信分析一、什么是基因家族概念：是来源于同一个祖先，有一个基因通过基因重复而产生两个或更多的拷贝而构成的一组基因， 他们在结构和功能上具有明显的相似性， 编码相似的蛋白质产物。划分：按功能划分 ：把一些功能类似的基因聚类，形成一个家族。按照序列相似程度划分: 一般将同源的基因放在一起认为是一个家族。1. 常见基因家族：WRKY基因家族 ：是植物前十大蛋白质基因家族之一，大量研究表明， WRKY 基因家族的许多成员参与调控植物的生长发育，形态建成与抗病虫。NBS-LRR抗病基因家族 ：是植物中最大类抗病基因家族之一。MADS-BOX基因家族： 是植物体内的重要转录因子，它们广泛地

2、调控着植物的生长、发育和生殖等过程。在植物中参与花器官的发育，开花时间的调节，在果实，根，茎，叶的发育中都起着重要的作用。热激蛋白 70 家族（HSP70）是一类在植物中高度保守的分子伴侣蛋白， 在细胞中协助蛋白质正确折叠。二、基因家族分析流程 ：利用蛋白保守域结构提取号在 Pfam 数据库提取其隐马尔科夫模型矩阵文件（ *.hmm）在数据库（ Ensemble 、JGI、NVBI）下载你所需要的物种的基因组数据（ *.fa,*.gff ）在虚拟机中 Bio-Linux 中的 hummsearch程序，用隐马尔科夫模型矩阵文件在蛋白序列文件中搜索含有该保守结构域的蛋白将蛋白序列导入 MEGA软

3、件构建进化树（可以阐明成员之间系统进化关系，从进化关系上揭示其多样性）利用 MEME搜索蛋白质的保守结构域利用 MEME搜索基因家族成员的 motif 可以揭示基因家族在物种内的多样化及其功能， 如果他们都含有相同的 motif 表明其功能具有相似性，如果部分家族成员含有其他不同的 motif, 很可能这些成员有其他特异功能，或者可以归分为一个亚族绘制基因染色体位置图从*.gff文件中抽取我们搜索到的基因位置信息，_v2.0/ 在线绘制基因染色体位置图通过染色体位置分布，可以了解基因主要分布字哪条染色体上，及是页眉内容否能形成基因簇（被认为是通过重组与错配促进基因交流）基因结构分析从 gff

4、文件中抽取基因的结构信息，绘制转录本结构图。计算串联重复基因的Ka,Ks1. 首先 将筛 选到 的基 因 的 cds 序列 进行 多序 列对比， 筛选 identity > 75% ，tength 大于对比的两条序列中较长的那条的长度的 75%，将筛选到的基因分别用 clustalw 进行比对，比对结果导入 KsKs_Calculster 计算 Ka，Ks、Ka/ks 比，计算核苷酸的非同义替代（ ka）与核苷酸的同义替代（ks）的平均速率。2. Ka/ks 比值 <1 表明：通过纯化选择降低了氨基酸变化的速率；比值=1 表示中性选择；比值 >1, 表明这些基因可能已经收到积

5、极选择，有利于适应性遗传，这些受正向选择的基因将作为以后的研究重点。软件的安装从图片中获得进入NCBI-blast 官网复制 blast-linux 版本的链接在 Linux 终端1. blast的安装#wget blast 链接#tar xvfz 文件名 解压缩文件# mv 解压缩文件 /root/local/app# mv 解压缩文件 blast页眉内容# vi .bashrc#在最后一行添加 export $PATH=/root/local/app/blast/bin:$PATH 并保存退出#source .bashrc 运行#blastp -version 查看是否安装成功。2humm

6、er的安装#yum install -y wget/安装 wget#wgethmmer 源码链接# vi .bashrc#(在最末端添加的语句 ) PATH=$PATH:/biosoft/ hmmer-#yum install -y gcc#./configure#make#make check#make install#which hmmsearch 查看是否安装成功。3.perl的安装#wget源代码链接# tar xvfz perl- 解压缩#./configure#make#make install 安装完成。3.bioperl的安装#wget -O -| bash#perlbrew

7、install-cpanm#/root/perl5/perlbrew/bin/cpanm Bio:Perl三、具体操作 ：1.保守域结构分析页眉内容下载蛋白保守结构域文件、cds、cDNA、gff注释文件和隐马尔科夫矩阵模型。以拟南芥为例：下载完成后，需要将文件传到Linux 系统上进行分析：打开虚拟机输入 ip a 将虚拟机 IP 连接到 Xshell 上，在 Xshell 上进行操作，将文件通过 xftp （同样需要连接 IP ）传到 Linux 系统上，然后进行解压。( 一个命令解压多个文件 :gunzip*.gz)接下来用 hummsearch寻找含有该蛋白保守结构域的蛋白及核酸序列安

8、装 hummsearchyum install -y wget#wgethmmer 源码链接/安装wget# vi .bashrc#(在最末端添加的语句 ) PATH=$PATH:/biosoft/ hmmer-#yum install -y gcc#./configure#make#make check#make install#which hmmsearch 查看是否安装成功。解压文件移动到 APP目录下面在 app 目录下面新建文件夹mkdir hmmer将删除安装包打开文字编辑器vi /.bashrc在文字编辑器里最后一行添加以上内容安装好 wget如果 make check 出现错误

9、百度用以下方法解决出现 complete 安装完成#source /.bashrc#which hmmsearch至此 hmmer安装完成。页眉内容虚拟机操作：1.导入下载好的文件；2.hmmsearch -cut_tc domtblout NB-ARC.txt NB-ARC.hmm Arabidopsis_thaliana.TAIR10可以用 editplus 打开 .txt 文件3. perl domain_xulie.pl 结果文件 .txt 蛋白序列文件domain.fa 1e-204.clustalw 进行多序列比对，得到aln 文件和 dnd 文件。5.hmmbuild 拟南芥特异

10、的hmm 模型文件domain.aln6.hmmsearchcut_tcdomtbloutresult.txtnewhmm 文件 蛋白质序列文件7.在 Excel 中，根据特定的evalue 进行筛选，并对第一列进行去重复，得到第一列去重复的 id，保存为 id.txt8. 用 perl 脚本根据 id 提取序列Per get_fa_by_id.pl id.txt 蛋白序列 wenjain >结果输出文件可以根据的得到的序列文件进行后续的构建进化树、motif 分析等。2.搜索基因家族成员的MOTIF2.1 需要准备的文件1. 拟南芥 NBS 基因蛋白质序列2. 蛋白保守结构域的隐马尔科

11、夫模型矩阵文件2.2 MOTIF 的搜索使用 meme 软件命令：meme nbs_pep.fa -protein -oc nbs_motif -nostatus -maxsize 600000 -mor anr-nmotifs 10 -minw 6 -maxw 50搜索结果存放在nbs_motif 文件夹中。文件夹中的eps 文件可以用AI打开编辑，可以另存为png或 jpg 格式，也可打开网页版，也可用tbtools软件打开，下载motif在基因上的位置信息。3. 绘制基因在染色体上的位置图3.1 需要准备的文件1. 拟南芥 NBS 基因 id2. 拟南芥的注释文件（ gff3 文件）页眉

12、内容3. 拟南芥基因组长度4.1 在线绘图工具： MapGene2Chrom4.2 samtools faidx 拟南芥 . 可得到拟南芥 . 该文件包括各个染色体，染色体长度。4.3 对基因的 id 文件在 Excel 中进行分列，去重复处理。4.4 使用处理过的 id 文件，对拟南芥的注释文件进行筛选使用 perl 脚本得到基因在染色体上的位置。命令： perlget_gene_gff.pl -in1 基因的 id 文件-in2 拟南芥 gff3 文件 -out 新文件名称4.5 新文件存放的是基因在染色体上的位置4.6 在在线文件 MapGene2chrom 中，将基因在染色体上的位置信

13、息文件复制到，input1 框中，在 input2 中粘入 samtools得到的 fai 文件。4.绘制转录本的结构图4.1 需要准备的文件1. 拟南芥 NBS 基因转录本 id（通过家族成员鉴定得到的蛋白2.拟南芥基因的注释文件（gtf 文件）id 文件）3.在线绘图工具： Gene Structure Display Server 2.0http： /4.2 具体方法1. 准备 gtf 文件：输入命令： gffreadgff3 注释文件-T-o输出文件（gtf文件）2.editplus 打开 gtf 文件，去除 ” transcript:”3. 使用 perl 脚本提取拟南芥转录本结构信

14、息：命令：perlget_gtf.pl-in1拟南芥转录本 id 文件 -in2 gtf 文件 -out 输出文件（ nbs_gtf.txt）4.通过在线绘图工具，进行绘图。5.筛选出串联重复基因5.1 准备文件1.拟南芥 NBS 基因 CDS 序列串联重复基因筛选标准【 （a）length of alignable sequence covers>75% of longer页眉内容gene,and(b) similarity of aligned regions >75%】参考文献： Extent of gene duplication in the genomes of Dro

15、sophila, nematode, andyeast.2.由于筛选时产生的文件较多，因此创建新的目录：mkdir 新目录3.用 editplus 打开家族成员的id 文件，对转录本id 进行处理，使一个基因只拿一个转录本。4.把 id 复制到 Excel，首先排序处理， 然后进行分列， 然后以第一列删除重复值。最后将第一列和第二列进行合并。将处理好的id 导入 Linux 。5.使用 perl 脚本提取 cds 序列：命令： perlget_fa_by_id.plid文件 拟南芥cds序列文件> cds.fa6.使用 blast 软件筛选串联重复基因6.1建立目标序列的数据库： mak

16、eblastdb-in cds.fa -dbtype nucl -title cds.fa6.2进行多序列比对： blastn -query cds.fa-db cds.fa -evalue 1e-20 -outfmt9-out result.txt6.3用 editplus 打开6.4得到 cds 序列的长度，使用 samtools工具建立索引：命令 ：samtools faidxcds.fa6.5 用 perl 脚本对result.txt 进行筛选， perl KAKS_SHAIXUAN.pl-in1-in2result.txt-outcleanresult.txt6.6 用 editpl

17、us 打开，将内容复制到Excel ，在id 后插入一列用公式：if（ A1>B1,A1&B1,B1&A1 ）。然后全选，以第 C 列删除重复值。并保存到新的文件中，并导入到 Linux 中7. 计算串联重复基因的 KaKs 。7.1 准备文件1.串联重复基因的CDS 序列文件7.2 将成对的串联重复序列保存在一个文件中。、方法一：复制需要找到的序列的id，在 editplus 中按 ctrl+F 搜索，找到后复制粘贴到一个文件中。方法二：首先将成对的id 保存在同一个文件，导入到Linux 中，在 Linux 中，利页眉内容用 perl 脚本提取序列： perlget_

18、fa_by_id.pl新的 id 文件 cds.fa 文件 >id1.fa7.3 计算 KaKS1.计算之前需要使用CLUSTAW 对序列进行比较。可获得id1.aln。2.使用 KaKs_calculator 工具将 id1.aln 文件转换成 id1.axt 文件命令： axtvenvertorid1.alnid1.axt3.计算 KaKs，输入命令：KaKs_calculstor-iid.axt-oid1_kaks.txt4.如果报错，则把两条序列长度保持一致。依此将所有的串联重复基因对，进行计算。四基因家族成员的鉴定（未知隐马尔科夫模型）1. 鉴定测略?在 NCBI 数据库中尽量

19、多下载几个物种的需要鉴定的蛋白保守结构域序列，以及所要研究物种的所有蛋白序列?在虚拟机中本地建库，并进行blast建库命令： makeblastdb -in 研究物种的蛋白序列文件-dbtype prot -title 库名称?进行序列比对命令： blastp -query 下载的多个物种序列文件-db 库名称 -evalue 1e-10 -outfmt6 -out 结构域 .blast?使用 sed 命令去除表头和结尾得到新的new 结构域 . blast?Awk print$1结new构域 . Blast |less 可查看打印的结果? Awk print$1结构new域 . Blast

20、> id.txt?cat id.txt|sort|uniq > idd.txt去重复?Perlget_fa_by_id.plidd.txt 去重复 蛋白序列 wenjain >结果输出文件?在 Pfam 或者 NCBI 的 cdd 中搜索检查是否有相关蛋白结构域。之后再进行 motif分析一些命令及软件应用说明（参考一些视频资料）hmmsearch使用说明用途：利用蛋白保守结构域的隐马尔科夫模型搜索蛋白序列中具有该保守结构域的蛋白用法： hmmsearch -cut_tc -domtblout result.txt*.hmm页眉内容说明： result.txt是输出的结果文件

21、，*.hmm在 pfam数据库下载的模型，hmmbuild使用说明用途：利用clustalw比对生成的aln文件构建蛋白保守结构域的隐马尔科夫模型用法：hmmbuildnew.hmmdomain.aln说 明 ： new.hmm是结果文件也就是构建的蛋白保守结构域的隐马尔科夫模型，domain.aln是 clustalw比对生成的aln 文件domain_xulie.pl脚本使用说明用途：提取hmmsearch搜索结果中蛋白序列中保守结构域的序列，用于构建新的物种特异的蛋白保守结构域的隐马尔科夫模型用 法 ： perldomain_xulie.pl（ 脚 本 不 在 使 用 目 录 下 要 写

22、 全 路 径 ） hmmoutfiledomain.fastaE-value说明： hmmoutfile是 hmmsearch搜索结果文件domain.fasta是结果存放文件也就是蛋白序列中保守结构域的序列，E-value是提取序列时设定的 E 值get_fa_by_id.pl 使用说明用途：通过 ID号获取其相应的基因或蛋白序列用法： perl perl get_fa_by_id.pl id.txtcds.fastat>id_cds.fasta说明： id.txt 是包含你的 ID 的文件， cds.fasta 是你丛数据库中下载的包含所有cds 序列的文件， id_cds.fast

23、a是输出文件内容是 ID 对应的序列samtools faidx用途：提取fasta文件信息用法： samtools faidx *.fa说明：输入文件是 fasta 文件，自动生成输出目录结果的 fai 文件第一列是你输入的 fasta 文件的 ID 第二列是其序列长度Gffread使用说明用途：将基因组注释文件 gff3 转化成基因的注释文件 gtf 用法： gffread my.gff3 -T -o my.gtf说明： my.gff3是输入文件 基因组注释文件，my.gtf是输出文件是基因的注释文件Get_gene_gff.pl用途：想要绘制基因的染色体位置图必须要拿到基因在染色体上的具

24、体信息，该脚本就是从总的 gff 文件抽去你需要的基因的信息如：所在染色体，起始终止位置等信息用法：PerlGet_gene_gff.pl-in1 gene_id.txt-in2my.gff3 -out gene_location.txt说明：gene_id.txt是第一个输入文件基因的ID文件， my.gff3是第二个输入文件是物种基因组所有蛋白序列， gene_location.txt是结果输出文件页眉内容Get_gtf.pl脚本使用说明用途：从基因注释文件gtf 文件中提取转录本的结构信息用法： perl get_gtf.pl-in1 id.txt-in2 gene.gtf-outstr

25、ucture.txt说明： id.txt是第一个输入文件是id 文件 ， gene.gtf是第二个输入文件是基因注释文件gtf 文件structure.txt是结果输出文件存放着转录本的结构信息KaKs_shaixuan.pl使用说明用途：在多序列比对结果文件中筛选identity大于 75% ，比对上的序列长度大于对比的两条序列中最长序列的长度的75%用法： perl KaKs_shaixuan.pl -in1 cds.fai -in2 result.txt -out shaixuan.txt说明： cds.fai是 samtools faidx对 cds 的序列文件fasta作用，生成的文件， resu