基因家族擴增與收縮分析

基因家族擴增與收縮分析

基因家族的擴張和收縮分析一般會使用orthoMCL進行同源基因識別，然後選擇直系同源基因進行物種樹構建，最後使用CAFE對聚類結果進行基因家族的擴張和收縮分析

直系同源基因鑑定

直系同源基因鑑定網上一般給出了兩個軟體，一個是orthoMCL，一個是orthofinder。orthoMCL雖然很長時間沒有進行過維護更新了，但大家進行基因家族擴張和收縮分析是依然經常性的使用，而orthofinder是16年出現的新軟體，本身使用和安裝起來更加方便，我也是比較推薦這個

提取最長轉錄本

進行mcl聚類之前，首先需要挑選每個基因最長的轉錄本形成一個fa文件，fa文件中是胺基酸序列，不是鹼基序列，這裡強調一下。然後cafe官方的例子中提供了一個腳本cafetutorial_longest_iso.py，但是這個腳本比較針對他提供的範例fa文件，對於正常分析時的fa文件是不適用的。我在網上找過很多的提取最長轉錄本的腳本，但是都不具有廣泛的適用性，最主要的原因還是幾個資料庫間對於基因組的注釋文件的規則有些不同，而且對於自己組裝形成的基因組也有不同，一般自己組裝的基因組gff文件是由EVM生成的。而有的生成出來一個基因具有多個轉錄本，而有的基因組注釋並不完全，本身有scaffold構成，一個mrna就是一個轉錄本，沒有重複。所以在提取最長轉錄本時需要根據實際情況進行提取。這裡我的基因組來源有三個，一個是ncbi，一個是ensembl，一個是evm。所以腳本內寫了三個模式，不同的gff使用不同的模式。主要思路是對gff進行解析，判斷其基因對應的所有mrna，以及所有mrna對應的cds。然後對所有cds進行長度計算，循環遍歷每個基因的所有轉錄本並進行長度比較，保留最長的那條寫入pep文件。
需要觀察注釋的就是解析gff文件這一步，比如對於ensembl注釋而言

image.png

判斷到mrna之後可以獲取到其parent所屬的gene是什麼，然就從下面的cds可以獲得這個cds對應的parent的transcript是什麼。而且cds的id就是蛋白id，所以可以使用蛋白序列進行提取，更加方便。
而對於ncbi注釋而言，道理是相同的

image.png

不過由於本身id前有部分前綴需要去掉，所以處理時記得注意一下
腳本參考本人github的extract_L_protein.py文件。
同時有時候部分注釋會出現點問題，到時候還是要具體情況具體分析，一切根據注釋來。然後輸入的fa文件記得對蛋白id進行處理，去掉多餘的description部分，一切與注釋中的結果對應即可

orthofinder運行

orthofinder的運行比較簡單。安裝好後運行

orthofinder -M msa -T fasttree -f dataset

orthofinder參數並不複雜
-M 指定基因樹推斷方法，默認『dendroblast』，可換成『msa』
-T 物種樹推斷方法，默認『fasttree』，可選fasttree, raxml, raxml-ng, iqtree。
-f 表示fa文件夾路徑
另外還有一些參數大家可以自己去看--help，比較簡單。
fa文件夾就是之前提取的最長轉錄本的fa文件的文件夾，所有物种放在一起。同時記得將基因id前加上物種名，如human|ENSP00000473549。這樣，幫助聚類之後分別基因來源。如何添加更改物種名這裡不做贅述，較為簡單，用python進行修改即可。
另外我測試時因為未知原因無法使用raxml參數進行建樹，使用默認fasttree才可行。

最後結果會在dataset文件夾中生成一個命名類似Results_Jan15的文件夾。裡面結果比較多。具體內容參看orthofinder官方文檔，或者參考另一篇orthofinder文章。

構建時間尺度物種樹

在orthofinder構建成功物種樹之後，還不能直接進入cafe中使用，需要使用r8s估計時間尺度。
r8s使用前需要編輯一個input文件，內容和方法參考陳連福老師博客

這裡使用cafe在github的腳本cafetutorial_prep_r8s.py。可以直接生成一個input文件
寫法如下

python2 ~/tools/CAFE/python_scripts/cafetutorial_prep_r8s.py -i ~/dataset/OrthoFinder/Results_Jan15/Species_Tree/SpeciesTree_rooted.txt -o r8s_ctl_file.txt -s 2335369 -p 'goat,cow' -c 25

-i 表示輸入的物種樹，我們直接使用生成的物種樹，路徑如上所示。然後orthofinder的樹節點間有1代表分叉節點如((cow:xx.000000,goat:xx.000000)1:x.xxxx，需要刪掉1，最後應該是((cow:xx.000000,goat:xx.000000):x.xxxx
-o 表示生成的input文件
-s 表示構建物種樹時所用的鹼基數量，這個由於是orthofinder構建好了，所以我們需要自己去看其進行比對的序列文件路徑為OrthoFinder/Results_Jan15/MultipleSequenceAlignments/SpeciesTreeAlignment.fa，裡面包括了所用的物種比對序列，計算其中一個物種的序列長度即可。其實這裡部分有點猜想的意思在裡面，因為文檔就說是進行序列比對的鹼基數，我不清楚是所有物種全加上還是其中一個物種的長度。但是觀察過範例的鹼基數後覺得應該是一個物種，如果是所有的話數字會比較大。
-p 制定某一個分叉節點的時間，如這裡制定了羊和牛的分叉時間為25，該數字可以到Timetree進行查詢。
然後生成r8s_ctl_file.txt之後如果不滿意自動生成的input，還能進行手動修改，比如添加更多分叉節點使結果更加精確，以及分叉節點進行重新命名等，具體參照上面推薦的博客進行修改
準備好input文件後運行r8s，如何安裝r8s這裡不進行贅述，不難。

~/tools/r8s1.81/src/r8s -b -f r8s_ctl_file.txt > r8s_tmp.txt

生成文件r8s_tmp.txt，最後為生成的物種樹，設定分離節點會變成字符，這裡需要手動將字符刪掉。
比如之前設置的分離節點((cow:xx.000000,goat:xx.000000)oatcow:x.xxxx，會多出oatcow，需要將oatcow刪掉變成((cow:xx.000000,goat:xx.000000):x.xxxx。

CAFE分析基因家族擴張和收縮

接下來使用CAFE分析基因家族的擴張和收縮

準備cafe.data文件

orthofinder的結果文件OrthoFinder/Results_Jan15/Orthogroups/Orthogroups.GeneCount.tsv就是cafe.data，但是需要稍微修改一下，講最後一列的total去掉，同時第一列添加desc，內容可以全為null
，cafe的格式要求就是一列desc，一列id，然後物種在該基因家族中的拷貝數。單拷貝基因家族就是該基因家族的所有物種拷貝數為1，之前挑選最長轉錄本就是因為如果不適用最長轉錄本會導致大部分基因的拷貝數不為1，使單拷貝基因家族數目不準。修改方法如下

##去掉total
awk 'OFS="\t" {$NF="" ;print $0}' Orthogroups.GeneCount.tsv > cafe.data.1
##增加desc列
sed 's/^/(null)\t&/g' "/home/mwshi/project/xunlu/mcl/compliantFasta/OrthoFinder/Results_Jan15/Orthogroups/cafe.data.1">/home/mwshi/project/xunlu/mcl/compliantFasta/OrthoFinder/Results_Jan15/Orthogroups/cafe.data
##手動把第一列的null改成desc，雖然不改也沒關係。不管用什麼方法，只要使格式符合cafe輸入要求即可

準備cafetutorial_run.sh文件

接下來編輯cafe的運行文件cafetutorial_run.sh，內容如下

#!cafe
load -i cafe.data -t 4 -l log_run1.txt -p 0.05
tree (((chimp:6,human:6):81,(mouse:17,rat:17):70):6,dog:93)
lambda -s -t (((1,1)1,(2,2)2)2,2)
report resultfile

load 以指定要分析的基因家族文件
tree 以指定系統發育樹的結構，使用之前生成的時間樹即可，記得去掉中間分支節點的字符，否則報錯。
lambda 設定不同分支的基因家族進化速率，相同數字表示演化速率相同，不同數字表示演化速率不同，這裡感覺數字並不代表速率本身，只是一個編號，我自己運行時由於不知道互相演化速率差異所以全部填1進行運行，軟體自身也會對演化速率進行估計。

report 表示返回的結果文件
編輯好後運行cafe

cafe cafetutorial_run1.sh

結果儲存在resultfile.cafe中。

cafe結果可視化

這一步對我非常困難，因為在網上搜cafe可視化發現並沒有什麼資料，都是說自己進行結果提取和畫圖。
這裡有兩種方法，一是使用cafe自帶腳本計算出擴張和收縮數目後自行繪圖，比如ggtree或者itol手動繪製等。
方法如下

python2 ~/tools/CAFE/python_scripts/cafetutorial_report_analysis.py -i resultfile.cafe -o summary_run1
python ~/tools/CAFE/python_scripts/cafetutorial_draw_tree.py -i summary_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生成結果文件中summary_run1_node.txt中包含每個節點擴張和收縮的數目。cafetutorial_draw_tree.py簡單的對結果進行繪圖，以上參數在resultfile.cafe中找。默認繪製擴張基因數目，可以添加-y Contractions來繪製收縮基因數目，然後根據這些使用其他工具繪製更美觀的圖片
第二種方法是我找到的唯一一個工具直接對cafe結果進行可視化的工具，CAFE_fig。
安裝時其主頁現實安裝ete3 3.0.0b35版本，不用管它，安裝最新版本。他提示安裝低版本是英文有部分人會因為其ete3無法使用PyQt5的內容所以必須要降級。但是如果選擇使用3.0.0b35版本代表必須將PyQt5降級到PyQt4，也比較麻煩。而且我自己測試時全部使用最新版並沒有出現問題，所以安裝時不需要降級。

python3 ~/tools/CAFE_fig/CAFE_fig.py resultfile.cafe -pb 0.05 -pf 0.05 --dump test/ -g pdf --count_all_expansions

直接運行即可，另外，由於本人物種樹時間長度不足，導致畫出的圖非常小。這裡提供兩種方法。第一種是對圖片的比例尺進行估計，比如想讓樹邊長5倍，則在cafe時對樹節點值*5，擴大長度，然後畫圖完成後使用ai對比例尺的值/5修正即可。第二種方法是我閱讀其代碼，發現他使用的ETE3進行繪圖，參照ETE3說明文檔之後只需要在CAFE_fig.py中對pixels_per_mya進行修正。源碼中是1.0，如果想讓長度變成5倍則更正為5.0即可

image.png

然後由於未知原因，我畫出的圖片沒有顏色，暫時還不清楚為什麼。
以上就是我弄出的基因家族擴張和收縮分析流程，因為本身在這一塊還屬於新手上路，可能整個流程還有很多缺陷需要彌補，希望各位大佬不吝賜教。然後我也不清楚是否還有更簡單，更高效的方法流程，如果有大佬願意指點指點，感激不盡。