中文超碰在线,cliu久久久

====環(huán)境和測試數(shù)據(jù)準(zhǔn)備=====

需要先安裝OrthoFinder，這個參考前面一個帖子。然后今天我們主要來測試CAFE。

?安裝也比較簡單：

下載安裝包之后

https://github.com/hahnlab/CAFE/releases/download/v4.2.1/CAFE-4.2.1.tar.gz

cd CAFE

./configure

make

make install prefix=文件夾

注：在release文件夾中成功安裝了cafe，可以把它加入自己的環(huán)境變量。

測試數(shù)據(jù)準(zhǔn)備：一共下載了mouse, rat, cow, horse, cat, marmoset,macaque, gibbon, baboon, orangutan, chimpanzee, 和 human 12個物種的蛋白數(shù)據(jù)。（我是另外一個博主那里下載的：從https://share.weiyun.com/5ZIjBg8 (密碼：3jzdpm)下載。）

//解壓縮一下就可以了

tar xf twelve_spp_proteins.tar.gz

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

rm twelve_spp_proteins/*.tar.gz

===識別基因家族=====

識別物種內(nèi)和物種間的基因家族分為如下四步：（當(dāng)然也可以）

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 -in makeblastdb_input.fa-dbtype prot -out blastdb

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

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

第三步：使用MCL進(jìn)行序列聚類

根據(jù)BLAST輸出中序列相似性信息尋找聚類。這些聚類將是后續(xù)用于CAFE分析的基因家族。聚類這一步將通過mcl處理。使用shell命令將BLAST轉(zhuǎn)成MCL能夠識別的ABC格式（其實就是挑選三列，兩個ID和Evalue）。

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能用的格式。

pythonpython_scripts/cafetutorial_mcl2rawcafe.py? -i dump.blast_output.mci.I30 -o unfiltered_cafe_input.txt -sp "ENSG00 ENSPTR ENSPPY ENSPAN ENSNLE ENSMMUENSCJA ENSRNO ENSMUS ENSFCA ENSECA ENSBTA"

這里的"ENSG00" 是ENSEMBL編號中物種的標(biāo)識符。

unfiltered_cafe_input.txt文件如下所示：

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

pythonpython_scripts/cafetutorial_clade_and_size_filter.py -iunfiltered_cafe_input.txt -o filtered_cafe_input.txt –s

然后把ID換成物種名字：

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

第五步：物種樹推斷

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

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

極大似然法: RAxML, PhyML, FastTree

貝葉斯法: MrBayes

推斷超度量樹

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

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

pythonpython_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'

/gpfs03/home/jingjing/software/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

運行CAFE

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

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

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

mkdir -p reports

cafe cafetutorial_run1.sh

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

python? /gpfs03/home/jingjing/software/CAFE/script/Fulton_python_scripts/cafetutorial_report_analysis.py? -i reports/report_run1.cafe -o reports/summary_run1? ?（注意這些python程序要基于python2才行）

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

summary_run1_node.txt: 統(tǒng)計每個節(jié)點中擴(kuò)張，收縮的基因家族數(shù)

summary_run1_fams.txt: 具體發(fā)生變化的基因家族

看下CAFE的輸出結(jié)果：

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

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

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

cafe結(jié)果可視化

在網(wǎng)上搜cafe可視化發(fā)現(xiàn)并沒有什么資料，都是說自己進(jìn)行結(jié)果提取和畫圖。

這里有兩種方法，一是使用cafe自帶腳本計算出擴(kuò)張和收縮數(shù)目后自行繪圖，比如ggtree或者itol手動繪制等。

方法如下：

python cafetutorial_draw_tree.py -isummary_run1_node.txt -t '(((chimp:6,human:6):81,(mouse:17,rat:17):70):6,dog:93)' -d '(((chimp<3>,human<3>)<3>,(mouse<3>,rat<3>)<3>)<3>,dog<3>)' -o summary_run1_tree_rapid.png

其中cafetutorial_report_analysis.py生成結(jié)果文件中summary_run1_node.txt中包含每個節(jié)點擴(kuò)張和收縮的數(shù)目。cafetutorial_draw_tree.py簡單的對結(jié)果進(jìn)行繪圖，以上參數(shù)在resultfile.cafe中找。默認(rèn)繪制擴(kuò)張基因數(shù)目，可以添加-y Contractions來繪制收縮基因數(shù)目，然后根據(jù)這些使用其他工具繪制更美觀的圖片。

第二種方法是唯一一個工具直接對cafe結(jié)果進(jìn)行可視化的工具，CAFE_fig（安裝的時候也碰到很多問題）。

安裝時其主頁現(xiàn)實安裝ete3 3.0.0b35版本，不用管它，安裝最新版本。他提示安裝低版本是英文有部分人會因為其ete3無法使用PyQt5的內(nèi)容所以必須要降級。但是如果選擇使用3.0.0b35版本代表必須將PyQt5降級到PyQt4，也比較麻煩。而且我自己測試時全部使用最新版并沒有出現(xiàn)問題，所以安裝時不需要降級。