比較基因組學分析目錄
1：單拷貝基因構建物種樹以及計算分化時間
2：基因家族收縮與擴張分析
3：特異節(jié)點富集分析
前言
上篇推文中介紹到比較基因組學分析常用套路的第一步仇祭，利用單拷貝基因構建具有分化時間的物種樹规惰，補充一點氛堕，對于跨度較大的物種浓体，可以選擇單拷貝基因的方法卡辰，比如此次分析使用的物種。對于目級或者科級水平來講茂腥，推薦使用共線性基因建樹圈浇。以十字花科為例，如果用單拷貝基因您旁，可能只有1000多組烙常，但是使用共線性可能有接近5000組。共線性基因建樹可以使用WGDI的方法鹤盒，這部分內(nèi)容我以后會探索一下蚕脏。
Sun et al 2022 WGDI
本篇推文主要講基因家族的收縮與擴張分析，使用的軟件是cafe5侦锯，2020年發(fā)表驼鞭，相較于cafe4來講操作更加方便并且新增了模型(Gamma)
1. 安裝
git clone https://github.com/hahnlab/CAFE5.git
cd CAFE5
./configure
make
安裝過程其實有點復雜，可能不同的服務器會出現(xiàn)不同的錯誤尺碰，這個請自行解決
2. CAFE5使用
輸入文件至少要兩個挣棕，一個是基因家族數(shù)目統(tǒng)計文件Genefamilies_Count.tsv，一個是樹文件tree.txt（帶有分化時間）亲桥，還可以增加一個lambda文件
2.1 主要參數(shù)
--fixed_alpha, -a
Alpha value of the discrete gamma distribution to use in category calculations. If not specified, the alpha parameter will be estimated by maximum likelihood.
--lambda_per_family, -b
Estimate lambda by family (for testing purposes only).
--cores, -c
Number of processing cores to use, requires an integer argument. Default=All available cores.
--error_model, -e
Run with no file name to estimate the global error model file. This file can be provided in subsequent runs by providing the path to the Error model file with no spaces (e.g. -eBase_error_model.txt).
--Expansion, -E
Expansion parameter for Nelder-Mead optimizer, Default=2.
--rootdist, -f
Path to root distribution file for simulating datasets.
--help, -h
Help menu with a list of all commands.
--infile, -i
Path to tab delimited gene families file to be analyzed - Required for estimation.
--Iterations, -I
Maximum number of iterations that will be performed in lambda search. Default=300 (increase this number if likelihood is still improving when limit is hit).
--n_gamma_cats, -k
Number of gamma categories to use. If specified, the Gamma model will be used to run calculations; otherwise the Base model will be used.
--fixed_lambda, -l
Value (between 0 and 1) for a single user provided lambda value, otherwise lambda is estimated.
--log_config, -L
Turn on logging, provide name of the configuration file for logging (see example log.config file).
--fixed_multiple_lambdas, -m
Multiple lambda values, comma separated, must be used in conjunction with lambda tree (-y).
--output_prefix, -o
Output directory - Name of directory automatically created for output. Default=results.
--poisson, -p
Use a Poisson distribution for the root frequency distribution. If no -p flag is given, a uniform distribution will be used. A value can be specified (-p10, or --poisson=10); otherwise the distribution will be estimated from the gene families.
--pvalue, -P
P-value to use for determining significance of family size change, Default=0.05.
--chisquare_compare, -r
Chi square compare (not tested).
--Reflection, -R
Reflection parameter for Nelder-Mead optimizer, Default=1.
--simulate, -s
Simulate families. Either provide an argument of the number of families to simulate (-s100, or --simulate=100) or provide a rootdist file giving a set of root family sizes to match. Without such a file, the families will be generated with root sizes selected randomly between 0 and 100.
--tree, -t
Path to file containing newick formatted tree - Required for estimation.
--lambda_tree, -y
Path to lambda tree, for use with multiple lambdas.
--zero_root, -z
Include gene families that don't exist at the root, not recommended.
其實主要用的就是-i -p -k -y -t這些參數(shù)
2.2 輸入文件準備
2.2.1. Genefamilies_Count.tsv
制表符分隔的基因家族計數(shù)文件穴张，通常用OrthoMCL,
OrthoFinder等軟件獲取計數(shù)信息。
示例格式
Desc Family ID human chimp orang baboon gibbon macaque marmoset rat mouse cat horse cow
ATPase ORTHOMCL1 52 55 54 57 54 56 56 53 52 57 55 54
(null) ORTHOMCL2 76 51 41 39 45 36 37 67 79 37 41 49
HMG box ORTHOMCL3 50 49 48 48 46 49 48 55 52 51 47 55
(null) ORTHOMCL4 43 43 47 53 44 47 46 59 58 51 50 55
Dynamin ORTHOMCL5 43 40 43 44 31 46 33 79 70 43 49 50
......
....
..
DnaJ ORTHOMCL10016 45 46 50 46 46 47 46 48 49 45 44 48
我們首先利用OrthoFinder的Orthogroups.GeneCount.tsv文件生成符合要求的輸入文件
cp Results_May02/Orthogroups/Orthogroups.GeneCount.tsv CAFE/
awk 'OFS="\t" {$NF=""; print}' Orthogroups.GeneCount.tsv > tmp && awk '{print "(null)""\t"$0}' tmp > cafe.input.tsv && sed -i '1s/(null)/Desc/g' cafe.input.tsv && rm tmp
查看文件格式
Desc Orthogroup Aof.pro Ath.pro Atr.pro Cba.pro Cri.pro Csa.pro Csu.pro Kle.pro Mpo.pro Nco.pro Osa.pro Ppa.pro Smo.pro Tpl.pro Vca.pro Vvi.pro Zma.pro
(null) OG0000000 145 112 95 5 372 129 3 1 2 217 126 16 206 419 4 177 117
(null) OG0000001 9 4 3 1691 9 96 2 56 2 4 21 0 2 5 3 2 0
(null) OG0000002 32 117 62 1 92 117 2 0 20 81 119 77 40 193 5 107 161
(null) OG0000003 37 104 54 3 89 76 4 5 10 74 144 22 47 134 8 79 154
(null) OG0000004 73 104 51 4 40 80 2 10 12 76 87 33 22 136 5 97 135
(null) OG0000005 28 46 36 11 37 47 0 3 50 81 81 32 48 120 0 54 73
(null) OG0000006 41 43 74 6 38 57 0 4 25 57 52 19 33 155 0 87 40
(null) OG0000007 58 52 60 0 18 42 0 0 12 50 56 17 57 99 1 82 52
(null) OG0000008 38 57 26 7 52 47 4 6 19 40 59 43 20 29 1 41 80
(null) OG0000009 46 57 26 1 25 46 1 2 11 52 65 29 13 50 1 48 87
生成之后還需要剔除不同物種間拷貝數(shù)差異過大的基因家族两曼，否則會報錯，有內(nèi)置腳本可以使用玻驻，我在運行的時候需要去掉第一行才能使用
python ~/soft/CAFE5/tutorial/clade_and_size_filter.py -i cafe.input.tsv -o gene_family_filter.txt -s
笨方法
awk 'NR==1 || $3<100 && $4<100 && $5<100 && $6<100 && $7<100 && $8<100 && $9<100 && $10<100 && $11<100 && $12<100 && $13<100 && $14<100 && $15<100 && $16<100 && $17<100 && $18<100 && $19<100 {print $0}' cafe.input.tsv >gene_family_filter.txt
最后的文件格式悼凑，保證第一行的物種名字與進化樹的一致即可
Desc Orthogroup Aof Ath Atr Cba Cri Csa Csu Kle Mpo Nco Osa Ppa Smo Tpl Vca Vvi Zma
(null) OG0000020 26 37 23 4 35 28 0 1 24 28 43 24 27 47 0 35 42
(null) OG0000021 49 41 31 7 30 31 8 2 7 26 49 15 11 31 0 36 45
(null) OG0000022 27 25 31 0 27 34 2 1 23 25 46 18 27 44 1 39 45
(null) OG0000024 37 40 27 0 22 30 1 11 9 33 38 18 25 43 0 37 39
(null) OG0000029 28 26 23 2 24 25 1 2 5 32 34 31 17 35 1 30 40
(null) OG0000030 23 30 16 1 23 27 1 1 27 26 26 16 15 49 1 28 35
(null) OG0000031 28 36 26 3 27 23 8 1 3 17 37 10 18 38 3 34 28
(null) OG0000032 18 16 25 0 24 19 0 5 4 25 36 6 38 49 1 38 35
(null) OG0000033 19 17 20 0 12 16 4 6 18 35 42 4 23 39 3 45 28
(null) OG0000035 17 37 17 8 28 24 2 2 5 30 41 8 6 26 2 32 37
(null) OG0000036 22 15 17 3 22 19 7 12 13 20 24 20 30 38 3 35 22
(null) OG0000039 14 27 36 0 34 24 0 2 2 12 41 4 47 13 0 37 26
(null) OG0000040 15 30 9 1 19 35 0 2 11 25 26 19 12 48 0 39 27
2.2.2. tree.txt
本步驟直接使用mcmctree生成的FigTree.tre文件修改一下即可使用
grep "UTREE 1 =" FigTree.tre | sed -E -e "s/\[[^]]*\]//g" -e "s/[ \t]//g" -e "/^$/d" -e "s/UTREE1=//" > tree.txt
2.3 運行CAFE5
cafe5 -i gene_family_filter.txt -t tree.txt -o out -c 1
## 如果使用Gamma模型與泊松分布
cafe5 -i gene_family_filter.txt -t tree.txt -o out -c 1 -k 2 -p ##注意-k可以調(diào)，一般為2-5
結果還是報錯
報錯
查了一下解決方法璧瞬，可以將lambda值調(diào)低(0.0001)再進行計算
解決方案
重新運行
cafe5 -i gene_family_filter.txt -t tree.txt -o out -c 16 -l 0.0001 -k 2 -p
當大家出現(xiàn)這種錯誤時户辫，可以嘗試這種解決方案，單單刪除這些家族是沒用的嗤锉。
3.4 CAFE5輸出結果
3.4.1 結果文件
Gamma_asr.tre ## 每個基因家族的樹文件
Gamma_branch_probabilities.tab ## 每個分支計算的概率
Gamma_category_likelihoods.txt
Gamma_change.tab ## 每一個基因家族在每個節(jié)點的收縮與擴增數(shù)目
Gamma_clade_results.txt ##每個節(jié)點基因家族的擴增/收縮數(shù)目
Gamma_count.txt ## 每一個基因家族在每個節(jié)點的數(shù)目
Gamma_family_likelihoods.txt
Gamma_family_results.txt ## 基因家族變化的p值和是否顯著的結果
Gamma_results.txt ## 模型渔欢，最終似然值，最終Lambda值等參數(shù)信息瘟忱。
我們主要用的文件有Gamma_asr.tre(主要對應后面表格中的節(jié)點)奥额、Gamma_change.tab(看哪些基因家族在哪個節(jié)點發(fā)生變化)苫幢、Gamma_clade_results.txt(體現(xiàn)在樹上，每個節(jié)點基因家族的收縮/擴增數(shù)目)垫挨、Gamma_family_results.txt(顯著擴增/收縮的基因家族)
3.4.2 每個節(jié)點基因家族收縮/擴增數(shù)目的體現(xiàn)
其實有繪圖腳本韩肝，但是很久沒有更新，可能不適用于CAFE5九榔，我們可以自己畫
將基因家族的擴增/收縮數(shù)目體現(xiàn)在樹上哀峻，需要兩個文件，Gamma_asr.tre哲泊，Gamma_clade_results.txt
cat Gamma_clade_results.txt
#Taxon_ID Increase Decrease
Mpo<21> 232 1298
Ppa<20> 2231 371
<31> 134 65
<25> 949 220
<23> 134 209
Atr<13> 516 922
Kle<29> 493 741
<12> 245 56
<28> 314 340
<4> 118 176
Cri<17> 1669 287
<22> 214 184
<19> 445 93
<16> 291 352
Osa<3> 579 572
Aof<6> 935 840
<8> 142 112
Csa<1> 326 834
<7> 640 138
Tpl<15> 1147 395
<14> 204 273
Nco<11> 631 559
Zma<2> 1776 232
Vvi<5> 958 433
<10> 413 112
<9> 345 66
Smo<18> 842 1315
Cba<24> 744 1664
<30> 23 17
Ath<0> 1090 291
Csu<27> 305 1560
Vca<26> 438 1040
less Gamma_asr.tre
BEGIN TREES;
TREE OG0000021 = ((Kle<29>*_2:820.007,(((Mpo<21>*_7:428.285,Ppa<20>*_15:428.285)<23>_12:70.3982,((Cri<17>_30:404.796,(Tpl<15>_31:308.175,(Atr<13>_31:208.47,(Nco<11>*_26:176.909,(((Osa<3>*_49:45.2652,Zma<2>_45:45.2652)<7>_45:60.4652,Aof<6>*_49:105.73)<9>*_44:25.8965,(Vvi<5>_36:101.378,(Csa<1>*_31:83.0358,Ath<0>*_41:83.0358)<4>_36:18.3421)<8>_36:30.2489)<10>*_36:45.2824)<12>_31:31.5605)<14>_31:99.7053)<16>*_30:96.6205)<19>*_27:48.1939,Smo<18>_11:452.99)<22>*_13:45.694)<25>*_12:190.787,Cba<24>_7:689.47)<28>*_7:130.537)<31>_5:177.864,(Csu<27>*_8:874.03,Vca<26>*_0:874.03)<30>_5:123.841)<32>_5
可以看到少了個32
我們可以利用Gamma_change.tab文件去找一下
經(jīng)查看確實沒有剩蟀，忽略這一部分
將兩個文件的nodeid對應即可繪圖
基因家族收縮和擴張
3.4.3 其他整理
相較于CAFE4，這些結果并沒有直接體現(xiàn)顯著擴張/收縮的基因家族切威，或者我們想找一下某個節(jié)點具體擴張的基因育特，可以結合目前拿到的輸出文件進行進一步整理
cat Gamma_family_results.txt |grep "y"|cut -f1 >p0.05.significant
#提取顯著擴張或收縮的orthogroupsID
grep -f p0.05.significant Gamma_change.tab > Gamma_p0.05change.tab
#顯著擴張/收縮的基因家族在每個節(jié)點的收縮與擴增數(shù)目
cat Gamma_p0.05change.tab | cut -f1,2 | grep "+[1-9]" | cut -f1 > node0significant.expand
#Ath顯著擴張的orthogroupsID
wc -l node0significant.expand
#Ath顯著擴張的基因家族數(shù)目
cp ../../Results_May02/Orthogroups/Orthogroups.tsv ./
grep -f node0significant.expand Orthogroups.tsv | cut -f3 | sed "s/ /\n/g" | sed "s/\t/\n/g" | sed "s/,//g" | sort | uniq > node0significant.expand.genes
#提取Ath顯著擴張的基因，方法一
cp ../../Results_May02/Orthogroups/Orthogroups.txt ./
grep -f node0significant.expand Orthogroups.txt |sed "s/ /\n/g" | grep "Ath" | sort | uniq > node0significant.expand.genes
#提取Ath顯著擴張的基因牢屋，方法二
4 結語
總體來說且预，個人感覺CAFE5要比CAFE4方便很多，但是輸出結果的可視化方面還需要加強烙无，下一篇推文將簡單實踐某個節(jié)點的基因功能富集分析(GO/KEGG),類似于分析下圖的紅色節(jié)點锋谐。
image.png
如果覺得有幫助歡迎在博客打賞
