比較基因組分析粘昨，首先需要選定要比較的物種码耐，然后下載這些物種的基因組數(shù)據(jù)。
基因組數(shù)據(jù)準備
- 數(shù)據(jù)準備
1.基因結(jié)構(gòu)注釋文件脆荷， gff/gtf格式
2.蛋白序列文件，fasta格式
3.cds序列文件懊悯，fasta格式 - 數(shù)據(jù)來源
1.自己組裝的基因組
2.數(shù)據(jù)庫下載 - 處理原則
1.一個基因保留一個轉(zhuǎn)錄本
2.ID要一致
數(shù)據(jù)庫
- GI-Phytozome下載植物基因組
https://phytozome-next.jgi.doe.gov/
提供去除可變剪切的cds和pep
Gff/gtf文件需要過濾
- Ensembl 基因組數(shù)據(jù)庫
植物： 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序列均包含所有轉(zhuǎn)錄本，需要過濾
gff文件包含所有轉(zhuǎn)錄本捧毛，需要過濾
NCBI下載基因組
https://www.ncbi.nlm.nih.gov/genome/
有些物種只有基因組序列观堂，沒有上傳基因注釋結(jié)果
基因組的染色體名稱為NCBI的accession編號
蛋白序列和cds序列均包含所有轉(zhuǎn)錄本，需要過濾
蛋白序列ID和基因ID無法直接對應呀忧，需要借助gff文件物種特異數(shù)據(jù)庫
擬南芥： https://www.arabidopsis.org/
苦蕎： http://www.mbkbase.org/Pinku1/
水稻： http://rice.plantbiology.msu.edu/
...
由于可變剪切的存在师痕，一條基因僅保留一條最長的轉(zhuǎn)錄本。蛋白而账、CDS序列ID需要和注釋文件一致胰坟。
提取最長轉(zhuǎn)錄本
Phytozome 提供最長轉(zhuǎn)錄本數(shù)據(jù)，可以直接下載福扬。當使用Ensembl 下載的數(shù)據(jù)時腕铸，可以參考以下腳本提取最長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和蛋白質(zhì)序列文件
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 #輸出蛋白序列
