「基因組學(xué)」使用CAFE進(jìn)行基因家族擴(kuò)張收縮分析

環(huán)境準(zhǔn)備

安裝軟件: 參考「基因組學(xué)」使用OrthoFinder進(jìn)行直系同源基因分析 安裝OrthoFinder，然后再安裝CAFE

wget https://github.com/hahnlab/CAFE/releases/download/v4.2/CAFE-4.2.tar.gz
cd CAFE
./configure 
make 
make install
# 如果沒有root權(quán)限
mkdir -p ~/opt/biosoft/CAFE-4.2/bin
make install prefix=~/opt/biosoft/CAFE-4.2

數(shù)據(jù)準(zhǔn)備: 一共會分析mouse, rat, cow, horse, cat, marmoset,macaque, gibbon, baboon, orangutan, chimpanzee, 和 human 12個物種的蛋白數(shù)據(jù)。既可以在Ensemble上下載，也可以從 https://share.weiyun.com/5ZIjBg8 (密碼：3jzdpm)下載。

解壓縮

tar xf twelve_spp_proteins.tar.gz
for i in `ls -1 twelve_spp_proteins/*.tar.gz`; do tar xf $i -C twelve_spp_proteins; done
rm twelve_spp_proteins/*.tar.gz

識別基因家族

識別物種內(nèi)和物種間的基因家族分為如下四步：

  1. 僅保留每個基因中有代表性的轉(zhuǎn)錄本，去除可變剪切和冗余基因
  2. 建立BLAST數(shù)據(jù)庫，使用blastp進(jìn)行 all-by-all 的比對
  3. 使用MCL基于blastp結(jié)果進(jìn)行聚類，基因序列相似的通常是一個基因家族
  4. 解析MCL的輸出結(jié)果，用作CAFE的輸入

將所有最長的轉(zhuǎn)錄本合并成單個文件

提取每個基因中最長的轉(zhuǎn)錄本，然后合并成單個文件。

python python_scripts/cafetutorial_longest_iso.py -d twelve_spp_proteins/
cat twelve_spp_proteins/longest_*.fa | seqkit rmdup - > makeblastdb_input.fa

All-by-all BLAST

先用makeblastdb建立BLAST數(shù)據(jù)庫

makeblastdb -in makeblastdb_input.fa -dbtype prot -out blastdb

之后用blastp進(jìn)行序列搜索，得到每個序列的相似序列

blastp -num_threads 20 -db blastdb -query makeblastdb_input.fa -outfmt 7 -seg yes > blast_output.txt &

-seg參數(shù)過濾低復(fù)雜度的序列(即氨基酸編碼為X), -num_threads線程數(shù)，此處設(shè)置為20.

這一步運(yùn)行的時間和你服務(wù)器性能有關(guān)，你可以讀完后續(xù)的內(nèi)容，再繼續(xù)分析。

使用MCL進(jìn)行序列聚類

現(xiàn)在，我們需要根據(jù)BLAST輸出中序列相似性信息尋找聚類。這些聚類將是后續(xù)用于CAFE分析的基因家族。聚類這一步將通過mcl處理。

使用shell命令將BLAST轉(zhuǎn)成MCL能夠識別的ABC格式

grep -v "#"  blast_output.txt | cut -f 1,2,11 > blast_output.abc

下一步，創(chuàng)建網(wǎng)絡(luò)文件(.mci)和字典文件(.tab)，然后進(jìn)行聚類

mcxload -abc blast_output.abc --stream-mirror --stream-neg-log10 -stream-tf 'ceil(200)' -o blast_output.mci --write-tab blast_output.tab &
# --stream-mirror: 為輸入創(chuàng)建鏡像，即每一個X-Y都有一個Y-X
# --stream-neg-log10: 對輸入的數(shù)值做-log10 轉(zhuǎn)換
# -stream-tf: 對輸入的數(shù)值進(jìn)行一元函數(shù)轉(zhuǎn)換，這里用到的是ceil(200)

根據(jù)mci文件進(jìn)行聚類, 其中主要調(diào)整的參數(shù)是-I, 它決定了聚類的粒度，值越小那么聚類密度越大，這個值沒有想象中的那么至關(guān)重要。一般設(shè)置為3，你也可以嘗試用其他值，然后比較結(jié)果。最終的目的是正確分析物種間的直系同源基因。

mcl blast_output.mci -I 3 
mcxdump -icl out.blast_output.mci.I30 -tabr blast_output.tab -o dump.blast_output.mci.I30

整理MCL的輸出結(jié)果

上一步MCL的輸出還不能直接用于CAFE，還需要對其進(jìn)行解析并過濾。

第一步是將原來的mcl格式轉(zhuǎn)成CAFE能用的格式

python python_scripts/cafetutorial_mcl2rawcafe.py -i dump.blast_output.mci.I30 -o unfiltered_cafe_input.txt -sp "ENSG00 ENSPTR ENSPPY ENSPAN ENSNLE ENSMMU ENSCJA ENSRNO ENSMUS ENSFCA ENSECA ENSBTA"

這里的"ENSG00" 是ENSEMBL編號中物種的標(biāo)識符。并且標(biāo)識符之前只能有一個空格，否則輸出結(jié)果就是下面這個情況

錯誤示范: 輸出

正確的輸出結(jié)果應(yīng)該是如下所示

正確示范

第二步，將那些基因拷貝數(shù)變異特別大的基因家族剔除掉，因?yàn)樗鼤斐蓞?shù)預(yù)測出錯。下面的腳本是過濾掉一個或多個物種有超過100個基因拷貝的基因家族，雖然不是特別的嚴(yán)格，但效果和根據(jù)拷貝數(shù)變異過濾類似

python python_scripts/cafetutorial_clade_and_size_filter.py -i unfiltered_cafe_input.txt -o filtered_cafe_input.txt -s

將原本的編號替換成有意義的物種名

sed  -i -e 's/ENSPAN/baboon/' -e 's/ENSFCA/cat/' -e 's/ENSBTA/cow/' -e 's/ENSNLE/gibbon/' -e 's/ENSECA/horse/' -e 's/ENSG00/human/' -e 's/ENSMMU/macaque/' -e 's/ENSCJA/marmoset/' -e 's/ENSMUS/mouse/' -e 's/ENSPPY/orang/' -e 's/ENSRNO/rat/' -e 's/ENSPTR/chimp/' filtered_cafe_input.txt
sed  -i -e 's/ENSPAN/baboon/' -e 's/ENSFCA/cat/' -e 's/ENSBTA/cow/' -e 's/ENSNLE/gibbon/' -e 's/ENSECA/horse/' -e 's/ENSG00/human/' -e 's/ENSMMU/macaque/' -e 's/ENSCJA/marmoset/' -e 's/ENSMUS/mouse/' -e 's/ENSPPY/orang/' -e 's/ENSRNO/rat/' -e 's/ENSPTR/chimp/' large_filtered_cafe_input.txt

物種樹推斷

構(gòu)建物種樹主要分為多序列聯(lián)配和系統(tǒng)發(fā)育樹推測兩步， 之后在已有進(jìn)化樹的基礎(chǔ)上構(gòu)建超度量樹用作CAFE輸入。

多序列聯(lián)配一般用的是單拷貝的直系同源基因，分別進(jìn)行多序列聯(lián)配之后然后合并成單個文件。接著用系統(tǒng)發(fā)育樹推測軟件進(jìn)行建樹，可選軟件有

  • 極大似然法: RAxML, PhyML, FastTree
  • 貝葉斯法: MrBayes

這里不展示具體過程，直接用已有的極大似然樹的結(jié)果(NEWICK格式)，保存為twelve_spp_raxml_cat_tree_midpoint_rooted.txt

(((cow :0.09289 ,( cat :0.07151 , horse :0.05727) :0.00398) :0.02355 ,(((( orang:0.01034 ,( chimp :0.00440 , human :0.00396) :0.00587) :0.00184 , gibbon:0.01331) :0.00573 ,( macaque :0.00443 , baboon :0.00422) :0.01431):0.01097 , marmoset :0.03886) :0.04239) :0.03383 ,( rat :0.04110 , mouse:0.03854) :0.10918);
已有進(jìn)化樹

推斷超度量樹

超度量樹(ultrametric tree)也叫時間樹，就是將系統(tǒng)發(fā)育樹的標(biāo)度改成時間，從根到所有物種的距離都相同。構(gòu)建方法有很多，比較常用的就是r8s.

這里用cafetutorial_prep_r8s.py構(gòu)建r8s的批量運(yùn)行腳本，然后提取超度量樹

python python_scripts/cafetutorial_prep_r8s.py -i twelve_spp_raxml_cat_tree_midpoint_rooted.txt -o r8s_ctl_file.txt -s 35157236 -p 'human,cat' -c '94'
r8s -b -f r8s_ctl_file.txt > r8s_tmp.txt
tail -n 1 r8s_tmp.txt | cut -c 16- > twelve_spp_r8s_ultrametric.txt

運(yùn)行CAFE

運(yùn)行CAFE有兩種模式，一種是CAFE的命令行模式，先執(zhí)行cafe進(jìn)行CAFE的shell, 然后在其中執(zhí)行命令。另一種是腳本模式，也就是你先把命令編輯完成，然后用cafe script_to_run.sh運(yùn)行。

估計出生-死亡參數(shù) \lambda

CAFE的主要功能就是根據(jù)給定的進(jìn)化樹和基因家族數(shù)估計一個或多個 birth-death(\lambda)參數(shù)。\lambda 參數(shù)描述的是基因出現(xiàn)或者消失的概率。

為整個樹估計單個\lambda

編輯cafetutorial_run1.sh。CAFE的命令不能有額外的空格出現(xiàn)在 tree后面的()中，以及l(fā)ambda 的 -t 后的()中，否則運(yùn)行時會無法正確解析文件導(dǎo)致報錯。

#!cafe
load -i filtered_cafe_input.txt -t 4 -l reports/log_run1.txt
tree ((((cat:68.710507,horse:68.710507):4.566782,cow:73.277289):20.722711,(((((chimp:4.444172,human:4.444172):6.682678,orang:11.126850):2.285855,gibbon:13.412706):7.211527,(macaque:4.567240,baboon:4.567240):16.056992):16.060702,marmoset:36.684935):57.315065):38.738021,(rat:36.302445,mouse:36.302445):96.435575)
lambda -s -t ((((1,1)1,1)1,(((((1,1)1,1)1,1)1,(1,1)1)1,1)1)1,(1,1)1)
report reports/report_run1

然后運(yùn)行如下命令

mkdir -p reports
cafe cafetutorial_run1.sh

結(jié)果統(tǒng)計

上一步運(yùn)行結(jié)束后的報告文件在reports/reportrun1.cafe,可以用已有的腳本分析哪些基因家族發(fā)生了擴(kuò)張或者搜索

python  python_scripts/cafetutorial_report_analysis.py -i reports/report_run1.cafe -o reports/summary_run1

在reports文件夾下會出現(xiàn)如下文件

  • summary_run1_node.txt: 統(tǒng)計每個節(jié)點(diǎn)中擴(kuò)張，收縮的基因家族數(shù)
  • summary_run1_fams.txt: 具體發(fā)生變化的基因家族

為高基因拷貝數(shù)的基因家族預(yù)測\lambda

之前過濾掉的高拷貝數(shù)變異的基因家族可以單獨(dú)進(jìn)行分析， 運(yùn)行命令如下

#!cafe
load -i large_filtered_cafe_input.txt -t 4 -l reports/log_run2.txt
tree ((((cat:68.710507,horse:68.710507):4.566782,cow:73.277289):20.722711,(((((chimp:4.444172,human:4.444172):6.682678,orang:11.126850):2.285855,gibbon:13.412706):7.211527,(macaque:4.567240,baboon:4.567240):16.056992):16.060702,marmoset:36.684935):57.315065):38.738021,(rat:36.302445,mouse:36.302445):96.435575)
lambda -l 0.00265952 -t ((((1,1)1,1)1,(((((1,1)1,1)1,1)1,(1,1)1)1,1)1)1,(1,1)1)
report reports/report_run2

為樹的不同部分預(yù)測多個\lambda

如果你認(rèn)為不同物種或者不同分支的基因家族進(jìn)化速率不同，那么可以讓CAFE預(yù)測多個\lambda值. 對lambda部分進(jìn)行調(diào)整， 相同數(shù)字表示\lambda相同，不同數(shù)字表示\lambda不同。

#!cafe
load -i filtered_cafe_input.txt -t 4 -l reports/log_run3.txt
tree ((((cat:68.710507,horse:68.710507):4.566782,cow:73.277289):20.722711,(((((chimp:4.444172,human:4.444172):6.682678,orang:11.126850):2.285855,gibbon:13.412706):7.211527,(macaque:4.567240,baboon:4.567240):16.056992):16.060702,marmoset:36.684935):57.315065):38.738021,(rat:36.302445,mouse:36.302445):96.435575)
lambda -s -t ((((3,3)3,3)3,(((((1,1)1,2)2,2)2,(2,2)2)2,3)3)3,(3,3)3)
report reports/report_run3

CAFE主要輸出文件格式

CAFE主要輸出內(nèi)容如下，下游分析所需信息需要通過對其解析獲得。

信息

Lambda是整個進(jìn)化樹的預(yù)測值

# IDs of nodes表示不同節(jié)點(diǎn)的編號，這里cat為0，horse為2，cat和horse所在的節(jié)點(diǎn)是1.

最后是每個基因家族的結(jié)果。以最開始的表示行為例，第一列對應(yīng)輸入基因家族的編號；第二列是Newick的進(jìn)化樹，cat_59中的59表示該基因家族在cat里有59個基因；第三列是Family-wide P-value，用于表明該基因家族是否是顯著性的擴(kuò)張或是收縮，這里是0.438，說明變化不明顯。在第三列的p值小于0.01時，第四列表明哪個分支的基因家族發(fā)生了變化，上圖中只有ID 11的基因家族有變化, 但是0,1,2,4分支并沒有變化。

參考資料

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容