CAFEを使って遺伝子重複イベントを評価する

#CAFEとは
複数種のゲノムレベルデータを基に、系統樹上で、遺伝子数が有意に変化した遺伝子ファミリーを検出します。
系統樹の枝長（分岐年代）に対して、遺伝子数の変動が有意に大きいファミリーをとってこれます。
特定の遺伝子ファミリーの重複イベントが対応位置で獲得された形質と関連があるかもという議論を行ったり、
これらの生物種は他の近縁分類群と比べて、非常に遺伝子の増減が大きい/少ないグループであるなどと議論します

#公式WebPage

CAFE v4.0
https://hahnlab.github.io/CAFE/
チュートリアルの中のpythonスクリプトが便利なので使っていく（チュート専用みたいな名前だが、十分汎用的）
チュートリアル用のファイル群は
https://iu.app.box.com/v/cafetutorial-files　
からダウンロードする

r8s
https://sourceforge.net/projects/r8s/

#系統樹作成
公式とは違ってorthofinderとIQ-treeを使っている。
公式のチュートはMCLとRAxML/PhyMLだが問題ないはず

##iq_treeによるMLツリーの作成

いかを参考におこなう

つくられたファイルはiqtree.contree とする

##r8sによる分岐年代ツリーへの変換

ダウンロードしたチュートリアル用のファイルから
python_scripts/cafetutorial_prep_r8s.py
を使って、r8sの実行用設定ファイルを系統樹情報から作成する(python はv2系)

r8s.sh

python cafetutorial_prep_r8s.py \
　-i iqtree.contree \            #"予め求めておいた種の系統樹（枝長付き）Newick形式"
  -o r8s_ctl_file.txt \          #(精製されるr8sの実行ファイルの名前を入れる）
  -s 10000000 \                   #"系統樹作成に使った形質数（塩基座位数/AA数）”
　-p ’sp1 , sp2’ \                #分岐年代推定に使う化石記録などがある分岐ノードを指定する　OTUを２つ指定すれば共通祖先を系統樹上から見つけてそこを指定していると判断してくれる
　-c ’100’ \　#化石記録などからえた分岐年数を指定する。複数わかっているノードがあるときは後述

形質数は"*.igtree"ファイルから
Input data: -- taxa with -- partitions and --- total sites (0% missing data)
のところから持ってきた

作られたr8s_ctl_file.txtファイルの中身を確認する

r8s_ctl_file.txt

#NEXUS
begin trees;
tree nl_tree = [&R] (sp4:0.2,(sp3:0.5,(sp1:0.4, sp2:0.8));
begin rates;
blformat nsites= 10000000 lengths=persite ultrametric=no;
collapse;
mrca node_name sp1 sp2;   
fixage taxon=node_name age=100;
divtime method=pl algorithm=tn cvStart=0 cvInc=0.5 cvNum=8 crossv=yes;
describe plot=chronogram;
describe plot=tree_description;
end;

分岐年代の補正ポイントが複数あるときは
mrc行とfixageをコピーして追加する。
node_nameがノードの名前、２つめを指定するときは別の名前をつけてfixage taxon=で参照する。
ML treeを使うときは　lengths=persite
MP　のときはlengths=total
ルートがうまく決まらない場合はrerootで指定する

以上をもとに例えばこんな感じに変更する

r8s_ctl_file.txt

#NEXUS
begin trees;
tree nl_tree = [&R] (sp4:0.2,(sp3:0.5,(sp1:0.4, sp2:0.8));
begin rates;
blformat nsites= 10000000 lengths=persite ultrametric=no;
collapse;
mrca node_name1 sp1 sp2;   
fixage taxon=node_name1 age=100;
mrca node_name2 sp4 sp3;   
fixage taxon=node_name2 age=10;
divtime method=pl algorithm=tn cvStart=0 cvInc=0.5 cvNum=8 crossv=yes;
describe plot=chronogram;
describe plot=tree_description;
reroot taxon=sp1;
end;

r8sの実行と作られた分岐年代補正が入った系統樹の取得

r8s_exe.sh

r8s -b -f r8s_ctl_file.txt > r8s.tmp
tail -n 1 r8s.tmp | cut -c 16- > twelve_spp_r8s_ultrametric.txt

#CAFE本体の実行

各ファミリーの遺伝子数を,種ごとにまとめた表（orthofinder であれば　Orthogroups.GeneCount.csvファイル)を使い、以下の様にtsv形式でまとめる
（ヘッダーはDescriptionとIDが固定(?)、その次はタクソン名, famiA,Bは任意のファミリー名）
全てのタクソンで遺伝子数が同じファミリーは、エラーがでるので除く。
RTaseなどのトランスポゾン関連の遺伝子も、過剰に増減するので、除いたほうが良いと思う（研究対象でなければ）。

ortho_num.tsv

Description 	 ID	sp1 	sp2	sp3 	 sp4
FamiA	FamiA	10 	20	10 	 0
FamiB	FamiB	30 	25	10 	 20
FamiC	FamiC	10 	20	30 	 20
.....

以下のように実行ファイルを準備する

run.sh

#! cafe
load -i ortho_num.tsv -t 2 -l ortho_num.tsv.log
tree (((sp1:20, sp2:20):20,sp3:40),10,sp4:50)
lambda -s -t (((1.1)1,1)1,1)
report ortho_num.tsv.rep

load行の　-i には作成したファミリーの遺伝子数ファイル　-t CPU数 -l 出力するログファイル名　を入れる
この他に -p (有意な増減とするときのp-value, Default 0.01) などが入れられる。
tree行にはr8sで作った分岐年代情報付きtreeを入れる
lambda行には 最初に解析するときは　-s を指定し、-t には上のようにNewick でtreeを表して、全ノードの場所に１を入れたものを入れる
-s を指定すると birth-death parameter λ　（与えられたデータセットでの平均的(？)なファミリーの増減のパラメータ）を新規に推定する。
-t　ではλが同じだと考えられるノードを数字で表す。（ツリー全体で同じ遺伝子の増減率だと考えられれば１で埋める。何らかの情報から違うと期待されるクレードがあればそこを２で埋める　例：(((1,1)1,(2,2)2)2,2)　）

変動の大きいファミリーの処理

RUNの前に変動の大きいファミリー（遺伝子数が多いファミリー）を別ファイルにする

-sでλを推定をする場合は、例外的なファミリーは入っていない方が良いので、
変動の大きいファミリーは取り除いて、ラムダなどを計算し
後で、求めたλを指定して大きなファミリーを計算したほうが良い。
目安は１００以上の遺伝子数が１つ以上のOTUで見られるような遺伝子ファミリーは除いてから解析する。

これをしてくれるのが、
python_scripts/cafetutorial_clade_and_size_filter.py
スクリプトである。

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

大きなファミリーが入っていない
filtered_cafe_input.txt
と
大きいファミリーしか入っていない
large_filtered_cafe_input.txt
ができる。

まず
filtered_cafe_input.txtでλを求める。

ちょっと実行ファイルを変更して

run.sh

#! cafe
load -i filtered_cafe_input.txt -t 2 -l ortho_num.tsv.log　<--変更
tree (((sp1:20, sp2:20):20,sp3:40),10,sp4:50)
lambda -s -t (((1.1)1,1)1,1)
report ortho_num.tsv.rep

実行は

cafe run.sh

lamdaはlog ファイルの下部や、つくられた*.rep.cafeファイルの２行目の最後に書いてある。

このlamda値をつかって別ファイルにされた大きいファミリー
large_filtered_cafe_input.txt
について計算する。時間がかかるのとコアダンプすることも多いので注意。-t数を工夫するなどする。

large_run.sh

#! cafe
load -i large_filtered_cafe_input.txt -t 2 -l large_ortho_num.tsv.log
tree (sp4:50,(sp3:20,(sp1:10, sp2:10):20)
lambda -l 0.001 -t ((1,1)1)1)
report large_ortho_num.tsv.rep

lambda　-l には求めておいたλの値を入れる。

結果の可視化

上記のレポートファイルは分かりづらいので表や系統樹にしてしまうほうが楽

表にするには
python_scripts/cafetutorial_report_analysis.py

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

cafecore.py
cafecore.pyc
が同じディレクトリにないとうごかないので注意.
また、internal nodeの名前が入っていると動かない
（r8sでの分岐年代の基準として指定したノードの名前が持ち越されているとエラーになる）
ので、その部分をレポートファイルから除くこと。
（今回であれば、node_name1　node_name2 をsedで除けば良い。単に全テキストから除けばOK)

図にするには
cafetutorial_draw_tree.py
を使う。
標準だどpng（解せぬ）なのでpdfなどにしたいときは67行目当たりを以下のように変える(pythonを讃えよ）

if output_file:
#    fig.savefig(output_file, format='png', bbox_inches='tight', dpi=300)
    fig.savefig(output_file, format='pdf', bbox_inches='tight')

実行方法

python /cafetutorial_draw_tree.py \
    -i summary_run1_node.txt \
    -t  "(sp4:50,(sp3:20,(sp1:10, sp2:10):20)"   \ #r8sで得たnewick形式の枝長つき系統樹”
　　 -d  "(sp4<1>,(sp3<2>,(sp1<3>:, sp2<4>)<5>)"  \   #”上記cafetutorial_report_analysis.pyなどで得たノードIDの入ったnewick(*fams.txtに入っている)" 
　　-o "出力ファイルパス” \
    -y Rapid \ #"出力する情報を指定する　Expansions/Contractions/Rapid"