比较基因组分析，首先需要选定要比较的物种，然后下载这些物种的基因组数据。

基因组数据准备

• 数据准备
  1.基因结构注释文件， gff/gtf格式
  2.蛋白序列文件，fasta格式
  3.cds序列文件，fasta格式
• 数据来源
  1.自己组装的基因组
  2.数据库下载
• 处理原则
  1.一个基因保留一个转录本
  2.ID要一致

数据库

• GI-Phytozome下载植物基因组
  https://phytozome-next.jgi.doe.gov/

提供去除可变剪切的cds和pep
Gff/gtf文件需要过滤

• Ensembl 基因组数据库
  植物： http://plants.ensembl.org/index.html， 79 个记录
  后生动物：http://metazoa.ensembl.org/index.htm ，112个记录
  原生生物：http://protists.ensembl.org/index.html ，237个记录
  真菌：http://fungi.ensembl.org/index.html ， 1,014 个记录
  细菌： http://bacteria.ensembl.org/index.html, 44,048 个记录

蛋白序列和cds序列均包含所有转录本，需要过滤
gff文件包含所有转录本，需要过滤

• NCBI下载基因组
  https://www.ncbi.nlm.nih.gov/genome/
  有些物种只有基因组序列，没有上传基因注释结果
  基因组的染色体名称为NCBI的accession编号
  蛋白序列和cds序列均包含所有转录本，需要过滤
  蛋白序列ID和基因ID无法直接对应，需要借助gff文件

• 物种特异数据库
  拟南芥： https://www.arabidopsis.org/
  苦荞： http://www.mbkbase.org/Pinku1/
  水稻： http://rice.plantbiology.msu.edu/
  ...

由于可变剪切的存在，一条基因仅保留一条最长的转录本。蛋白、CDS序列ID需要和注释文件一致。

提取最长转录本

Phytozome 提供最长转录本数据，可以直接下载。当使用Ensembl 下载的数据时，可以参考以下脚本提取最长cds。

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

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

## 基于gff3提取最长cds序列ID，并过滤gff3文件
perl ./gff_longest.pl Ath.gff3 Ath_id Ath_longest.gff3

## 提取最长的cds ID
awk '{print $2}' Ath_id > Ath_longest_id

## 基于最长的cds提取cds和蛋白质序列文件
seqtk subseq Arabidopsis_thaliana.TAIR10.cds.all.fa  Ath_longest_id > Ath_longest.cds.fasta

seqtk subseq   Arabidopsis_thaliana.TAIR10.pep.all.fa Ath_longest_id > Ath_longest.pep.fasta

## 如果没有cds和蛋白文件，也可以基于过滤后gff3文件从genome里提取cds并翻译成蛋白
gffread Ath_longest.gff3 \
-g Arabidopsis_thaliana.TAIR10.dna.toplevel.fa \ #基因组序列
-x Ath_longest.cds.fasta \ #输出cds序列
-y Ath_longest.pep.fasta #输出蛋白序列

欢迎关注Bioinfor生信云！