构建某一基因家族的进化树

2020-04-07 本文已影响0人巩翔宇Ibrahimovic

写在前面
绘制某一基因家族的进化树，这是基因组学研究过程中不可缺少的一部分。划分成不同的基因家族的依据一般是其具有某一相同的结构域，下面我将以水稻中GUB_WAK_bind类家族绘制一下进化树。参考中国农大徐明良老师于2015年发表于NG上的Zm-WAK文章。

软件安装

HMMER3.3

MEGA-X(windows)

HMMER3.3 官方地址：http://hmmer.org/download.html

#HMMER3.3支持conda安装
conda install -c bioconda hmmer

1.准备文件

这里我刚开始选用的是RAP-DB上的文件，但是后来发现它上面的注释文件对基因的简称并不友好，搜索各大数据库网站后，发现ENSEMBL PLANT网站上的注释文件符合我的要求。


#下载基因组文件
wget -c https://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_genome.fasta.gz
#下载CDS文件
https://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_cds_2020-03-24.fasta.gz
#下载GFF文件
wget -c https://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_representative_transcript_exon_2020-03-24.gtf.gz
#下载全基因组蛋白质文件，ensemble plant里面的注释能看到CGSNL Gene Name(也就是我们平时说的简称)
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-46/fasta/oryza_sativa/pep/Oryza_sativa.IRGSP-1.0.pep.all.fa.gz
#解压缩及重命名
gunzip -f
mv oldname newname

2.准备hmm文件

1.登录网站：[http://pfam.xfam.org/](http://pfam.xfam.org/)

2.关键字搜索"WAK"，选择**GUB_WAK_bind**

3.点击 **curation & model**，拉到最下面，选择**download**。

3.进行HMM搜索

hmmsearch GUB_WAK_bind.hmm protein.fasta > GUB_WAK_bind_out

4.寻找基因ID对应的蛋白质序列

# 截取id号
vim GUB_WAK_bind_out
# 获取id号所在的行号，然后再用sed命令截取行，再用grep命令将id号匹配并重定向。
在vim命令模式下，输入“：set nu”
# sed命令截取，并用管道符直接输入给grep，匹配重定向到id文件
sed -n '17,149p'  GUB_WAK_bind_out | grep -o "Os.*" | cut -d ' ' -f 1 > id_all_gub_wak
#寻找id对应的简称,从而人为去除不必要的基因,并提取需要的列
grep -f id_all_gub_wak protein.fasta | grep '>.*' | cut -d ' ' -f 1,4,8,9,10,11,12 > res.txt
#根据res，人工选择只有wak注释的基因，共60个
得到id_wak_only.txt
# 利用samtools工具来进行序列提取
# 首先建立索引文件
samtools faidx protein.fasta
# 再将id好作为输入，之后在重定向
# 参考链接：https://www.biostars.org/p/49820/
xargs samtools faidx protein.fa < id_wak_only.txt> GUB_WAK_protein
# 得到的序列文件是含有回车符的，我利用一个perl单行命令将fasta格式的多行序列变成单行的fasta格式序列，链接：http://www.biotrainee.com/thread-291-1-1.html
perl -pe '/^>/ ? print "\n" : chomp' GUB_WAK_protein | tail -n +2 > GUB_WAK_protein_out

5.寻找基因组DNA序列(可选)

>#根据基因id来抓取注释文件里面的位置信息
grep -f response2_id.txt transcript.gtf | cut -f 1,4,5,9| sed 's/gene.*\s"//'| sed 's/";//' > gene.bed
#建立基因组文件的索引
samtools faidx rice_genome.fasta
#根据名称和位置定位到基因组序列上
bedtools getfasta -fi rice_genome.fasta -bed gene.bed -name > gene_seq

最终得到的这个GUB_WAK_protein_out文件，拿到MEGA里面去绘制进化树,记得要修改后缀名为.fasta。
首先是导入文件进行序列的比较，然后再选edit里面绘制进化树，最后是退出来，进行界面的极大似然进化树的绘制，记得要选择bootstrap值(我选择的是1000),在test里面。
下一篇我们来谈谈进化树的分类及美化问题。