以下參照 hoptop的「基因組學(xué)」使用OrthoFinder進行直系同源基因分析這篇帖子試跑了一遍

1.軟件安裝

#用conda安裝orthofinder
conda create -n orthofinder orthofinder=2.2.7

#安裝CAFE
安裝二進制軟件底燎，最好有管理員權(quán)限

2.數(shù)據(jù)準(zhǔn)備

下載「基因組學(xué)」使用OrthoFinder進行直系同源基因分析 中 hoptop準(zhǔn)備的數(shù)據(jù)，總共有十二個物種的蛋白數(shù)據(jù)以及接下來會用到的Python腳本卜朗。

#解壓縮
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

3.識別基因家族：

①僅保留每個基因中有代表性的轉(zhuǎn)錄本熬的，去除可變剪切和冗余基因；
②建立BLAST數(shù)據(jù)庫拌喉，使用blastp進行 all-by-all 的比對速那；
③使用MCL基于blastp結(jié)果進行聚類，基因序列相似的通常是一個基因家族司光；
④解析MCL的輸出結(jié)果琅坡，用作CAFE的輸入。

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

#提取每個基因中最長的轉(zhuǎn)錄本
python python_scripts/cafetutorial_longest_iso.py -d twelve_spp_proteins/ 
#合并每個物種的最長轉(zhuǎn)錄本成為1個文件
cat twelve_spp_proteins/longest_*.fa | /data1/spider/ytbiosoft/seqkit rmdup - > makeblastdb_input.fa

②All-by-all BLAST

##我的blast工具安裝在miniconda的一個子環(huán)境里面残家，首先需要激活這個環(huán)境
source/data1/spider/miniconda3/bin/activate
conda activate bioinforspace

##先用makeblastdb建立BLAST數(shù)據(jù)庫
makeblastdb -in makeblastdb_input.fa -dbtype prot -out blastdb

##之后用blastp進行序列搜索榆俺，得到每個序列的相似序列
blastp -num_threads 20 -db blastdb -query makeblastdb_input.fa -outfmt 7 -seg yes > blast_output.txt &
#參數(shù)：
-seg參數(shù)過濾低復(fù)雜度的序列(即氨基酸編碼為X)
-num_threads線程數(shù)，此處設(shè)置為20

③使用MCL進行序列聚類
根據(jù)BLAST輸出中序列相似性信息尋找聚類，這些聚類將用于后續(xù)CAFE分析的基因家族茴晋。

##使用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)陪捷，然后進行聚類source /data1/spider/miniconda3/bin/activateconda activate orthofinder     #mcxload安裝在orthofinder環(huán)境中mcxload -abc blast_output.abc --stream-mirror --stream-neg-log10 -stream-tf 'ceil(200)' -o blast_output.mci -write-tab blast_output.tab &#參數(shù)：--stream-mirror: 為輸入創(chuàng)建鏡像，即每一個X-Y都有一個Y-X--stream-neg-log10: 對輸入的數(shù)值做-log10 轉(zhuǎn)換-stream-tf: 對輸入的數(shù)值進行一元函數(shù)轉(zhuǎn)換诺擅，這里用到的是ceil(200)

#根據(jù)mci文件進行聚類, 其中主要調(diào)整的參數(shù)是-I市袖，它決定了聚類的粒度，值越小那么聚類密度越大烁涌。一般設(shè)置為3苍碟。
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，還需要對其進行解析并過濾撮执。

##1.將原來的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"
##2.將那些基因拷貝數(shù)變異特別大的基因家族剔除掉微峰，因為它會造成參數(shù)預(yù)測出錯。以下腳本是過濾掉一個或多個物種有超過100個基因拷貝的基因家族抒钱。
python python_scripts/cafetutorial_clade_and_size_filter.py -i unfiltered_cafe_input.txt -o filtered_cafe_input.txt -s
##3.將原本的編號替換成有意義的物種名
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

4.物種樹推斷

構(gòu)建物種樹主要分為多序列聯(lián)配和系統(tǒng)發(fā)育樹兩步蜓肆，隨后，在已有進化樹的基礎(chǔ)上構(gòu)建超度量樹用作CAFE的輸入谋币。
①多序列聯(lián)配一般用的是單拷貝的直系同源基因仗扬，分別進行多序列聯(lián)配，然后合并成單個文件蕾额；
②使用系統(tǒng)發(fā)育樹推測軟件進行建樹（軟件： RAxML, PhyML, FastTree早芭； MrBayes），生成的系統(tǒng)發(fā)育樹用NEWICK格式保存為twelve_spp_raxml_cat_tree_midpoint_rooted.txt凡简；

cat 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);

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

#使用cafetutorial_prep_r8s.py構(gòu)建r8s的批量運行腳本筐眷，然后提取超度量樹
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'
/data1/spider/liupiao/biosoft/r8s1.81/src/r8s -b -f r8s_ctl_file.txt > r8s_tmp.txt &
tail -n 1 r8s_tmp.txt | cut -c 16- > twelve_spp_r8s_ultrametric.txt

5.運行CAFE

估計出生-死亡參數(shù)λ
CAFE的主要功能是根據(jù)給定的進化樹和基因家族數(shù)估計一個或多個λ參數(shù)黎烈。λ參數(shù)描述的是出現(xiàn)或消失的概率。
①為整個樹估計單個λ

##編輯腳本cafetutorial_run1.sh匀谣。注意：以下腳本的括號中不能出現(xiàn)任何空格照棋。

#!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

mkdir -p reports
/data/biosoft/cafe cafetutorial_run1.sh   ##運行腳本

##結(jié)果統(tǒng)計
#上一步運行結(jié)束后的報告文件存放在reports/reportrun1.cafe,可以用已有的腳本分析哪些基因家族發(fā)生了擴張或進行搜索
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é)點中擴張，收縮的基因家族數(shù)
summary_run1_fams.txt: 具體發(fā)生變化的基因家族

②為高基因拷貝數(shù)的基因家族預(yù)測λ
之前過濾掉的高拷貝數(shù)變異的基因家族可以進行單獨的分析武翎，運行命令如下：

#!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

#運行腳本
/data/biosoft/cafe cafetutorial_run2.sh

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

#!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

#運行腳本
/data/biosoft/cafe cafetutorial_run3.sh

CAFE主要輸出文件格式

Lambda是整個進化樹的預(yù)測值
IDs of nodes表示不同節(jié)點的編號趴捅，這里cat為0，horse為2霹疫，cat和horse所在的節(jié)點是1
最后是每個基因家族的結(jié)果拱绑。以最開始的表示行為例，第一列對應(yīng)輸入基因家族的編號丽蝎；第二列是Newick的進化樹猎拨，cat_59中的59表示該基因家族在cat里有59個基因；第三列是Family-wide P-value屠阻，用于表明該基因家族是否是顯著性的擴張或是收縮红省，這里是0.438，說明變化不明顯栏笆。在第三列的p值小于0.01時类腮，第四列表明哪個分支的基因家族發(fā)生了變化，上圖中只有ID 11的基因家族有變化, 但是0,1,2,4分支并沒有變化蛉加。

亟待解決的問題：怎么將cafe輸出結(jié)果可視化？缸逃？针饥？

參考
01 「基因組學(xué)」使用OrthoFinder進行直系同源基因分析
02 「基因組學(xué)」使用CAFE進行基因家族擴張收縮分析