基因家族分析（1）数据准备及软件安装

欢迎关注Bioinfor 生信云微信公众号！

基因组数据的准备

进行全基因水平的基因家族鉴定之前，需要准备好一套基因组数据，基因组数据可以从公共数据库下载，也可以根据基因组文献提供的地址到指定网站进行下载。一套完整的数据至少包括如下内容：

1. 基因组序列文件，fasta 格式
2. 基因结构注释文件，gff3 格式
3. 所有蛋白质序列文件，fasta 格式
4. 所有 cds 序列文件，fasta 格式
   由于不同基因组数据库存储的数据格式及命名有各自规范，很多数据下载后不能直接用来做基因家族分析，需要提前处理好，处理原则如下：
   1.去除所有可变剪切，一个编码基因保留一个转录本
   2.去除 gff3 文件中非编码基因及重复序列等信息下面介绍 JGI、Ensembl 和 NCBI 下载的数据的处理

JGI/phytozome 数据处理

phytozome 是一个收录植物基因组数据的网站，数据整理比较规范，已经提供了去除可变剪切的 cds 和 protein 序列文件。只有 gff3 文件需要过滤处理。

示例数据为拟南芥数据，下载于phytozome13

Athaliana_167_TAIR10.cds_primaryTranscriptOnly.fa #cds序列
Athaliana_167_TAIR10.gene_exons.gff3 #基因结构文件
Athaliana_167_TAIR10.protein_primaryTranscriptOnly.fa #蛋白文件
Athaliana_167_TAIR9.fa #基因组文件

## 提取最长转录本基因ID
awk '$1 ~ /^>/ {print $1}' Athaliana_167_TAIR10.cds_primaryTranscriptOnly.fa |   sed 's/^>//' > Ath_mRNA.id

## gff3文件相对于cds和蛋白序列文件，ID部分多了.TAIR10的字符串需要去除掉
sed 's/\.TAIR10//g' Athaliana_167_TAIR10.gene_exons.gff3 > Ath_1.gff3

## 基于mRNA id对gff3文件进行过滤
perl gff_filter_bymRNAID.pl   Ath_1.gff3    Ath_mRNA.id   geneID_mrnaID.table   Ath_final.gff3

## 重命名蛋白序列和cds序列文件及基因组名称，方便后续使用
mv Athaliana_167_TAIR10.cds_primaryTranscriptOnly.fa   Ath.cds.fa
mv Athaliana_167_TAIR10.protein_primaryTranscriptOnly.fa  Ath.pep.fa
mv Athaliana_167_TAIR9.fa Ath.genome.fa

Ensembl 数据处理

Ensembl 数据库可以下载动物、植物、真菌、细菌等物种基因组数据。数据格式规范，mRNA ID 和 cds ID 基本一致，但没有去除可变剪切的版本，需要自己进行手动处理

示例数据为拟南芥数据，下载自 ensembl.

Arabidopsis_thaliana.TAIR10.47.gff3 # 基因结构文件  
Arabidopsis_thaliana.TAIR10.dna.toplevel.fa # 基因组序列文件 
Arabidopsis_thaliana.TAIR10.cds.all.fa # cds序列文件
Arabidopsis_thaliana.TAIR10.pep.all.fa # 蛋白序列文件

# 去除gff3文件中ID部分多余字符
cp Arabidopsis_thaliana.TAIR10.47.gff3  Ath.gff3.tmp #复制一份
sed -i 's/=gene:/=/g' Ath.gff3.tmp
sed -i 's/=transcript:/=/g'  Ath.gff3.tmp
sed -i 's/=CDS:/=/g'  Ath.gff3.tmp

# 基于gff3提取最长cds序列ID，并过滤gff3文件

perl gff_filter_longest.pl    Ath.gff3.tmp   Ath_gene_mrna_cds.ids  Ath_final.gff3

# 提取最长cds ID列表
awk '{print $3}' Ath_gene_mrna_cds.ids > Ath_mRNA.id

##基于最长cds ID信息提取cds和蛋白质序列文件
seqtk subseq Arabidopsis_thaliana.TAIR10.cds.all.fa Ath_mRNA.id  >  Ath.cds.fasta

seqtk subseq Arabidopsis_thaliana.TAIR10.pep.all.fa Ath_mRNA.id  >  Ath.pep.fasta

#  基因组文件重命名
mv Arabidopsis_thaliana.TAIR10.dna.toplevel.fa Ath.genome.fasta

NCBI及其它数据库的处理有需要的请私信我

#没有cds和蛋白序列的情况
如果没有 cds 和蛋白序列，可以基于 gff 和基因组序列文件使用 gffread进行提取.

gffread Ath_final.gff3 -g Ath.genome.fasta -x Ath.cds.fasta #提取cds序列
gffread Ath_final.gff3 -g Ath.genome.fasta -y Ath.pep.fasta #提取蛋白序列

软件安装

conda安装

用conda安装比对、结构域预测、motif鉴定、进化树构建、多序列比对结果过滤、fasta序列处理工具等等
blast
hmmer
meme
fasttree
trimal
seqkit
gffread
McscanX
JCVI

R包的安装

Peptides #蛋白质等电点和分子量的统计
seqlogo #绘制seqlogo图
pheatmap #绘制热图
msa #多序列比对的R包

windows软件

染色体核型图mapchart
进化树构建 mega

在线软件

进化树美化 https://itol.embl.de/
motif 预测meme MEME - Submission form (meme-suite.org)
基因结构绘制 Gene Structure Display Server 2.0 (gao-lab.org)
顺式作用元件预测 PlantCARE, a database of plant promoters and their cis-acting regulatory elements (ugent.be)