OrthoFinder2：直系同源基因的尋找以及Orthogroup構(gòu)建

OrthoFinder2算法解讀

1）基因樹的構(gòu)建，

  • 尋找序列之間的同源關(guān)系（相似度的計(jì)算）：diamond，blast，mmseq2
  • 基因樹的推斷：DendroBLAST，其他的參考config.json

2）物種樹的構(gòu)建

  • STAG：從無根gene tree推斷無根物種樹
  • STRIDE：從無根物種樹推斷有根物種樹

Note：

  • STAG：Species Tree from All Genes (STAG)

3）使用“duplication-loss-coalescent model”，推斷orthogroup（ortholog + paralog）

一點(diǎn)點(diǎn)背景知識

什么是正交組（Orthogroup）？

一個(gè)祖先物種經(jīng)歷過一次物種事件之后，分別形成了物種A和物種B，

祖先物種gene分化得到geneA和geneB —— geneA和geneB，稱為正交群（Orthogroup）

Note：paralogs，也可以和orthologs一起，存在于同一個(gè)Orthogroup中

參考資料，

[1] 來自lakeseafly的簡書筆記
[2] https://en.wiktionary.org/wiki/orthogroup

為什么要研究正交群？

  • 正交群中的所有基因都來自<font color='yellow'>單個(gè)祖先基因</font> —— 正交群中的所有基因都有類似的序列和功能

    由于基因重復(fù)和丟失在進(jìn)化中經(jīng)常發(fā)生，<u>1 vs 1 的直系同源物非常少見</u>，

    因此，分析正交群的所有情況（1 vs 1，1 vs more，more vs more）可以得到有關(guān)進(jìn)化歷史的信息。

  • “真正識別直系同源物的唯一方法，是構(gòu)建系統(tǒng)發(fā)育樹” —— lakeseafly

軟件安裝

conda create -n orthofinder2 -c bioconda orthofinder -y

軟件使用

以mouse、人類、斑馬魚、青蛙（frog？）、日本河豚、果蠅的蛋白質(zhì)序列為例

1、數(shù)據(jù)下載

需要注意的是，

  • OrthoFinder2要求輸入的序列為氨基酸（蛋白質(zhì)）序列（CDS對應(yīng)的序列，即protein coding genes序列）
# human
wget -c http://ftp.ensembl.org/pub/release-106/fasta/homo_sapiens/pep/Homo_sapiens.GRCh38.pep.all.fa.gz
# mouse
wget -c http://ftp.ensembl.org/pub/release-106/fasta/mus_musculus/pep/Mus_musculus.GRCm39.pep.all.fa.gz
# zebrafish
wget -c http://ftp.ensembl.org/pub/release-106/fasta/danio_rerio/pep/Danio_rerio.GRCz11.pep.all.fa.gz
# frog
wget -c http://ftp.ensembl.org/pub/release-106/fasta/xenopus_tropicalis/pep/Xenopus_tropicalis.Xenopus_tropicalis_v9.1.pep.all.fa.gz
# fugu
wget -c http://ftp.ensembl.org/pub/release-106/fasta/takifugu_rubripes/pep/Takifugu_rubripes.fTakRub1.2.pep.all.fa.gz
# fruit fly
wget -c http://ftp.ensembl.org/pub/release-106/fasta/drosophila_melanogaster/pep/Drosophila_melanogaster.BDGP6.32.pep.all.fa.gz

將上述命令保存為腳本并運(yùn)行，

sh datadownload.sh 2>220704.datadownload.log &

2、提取最長轉(zhuǎn)錄本

這一步，主要是為了減少OrthoFinder2的運(yùn)行時(shí)間，可以自己寫腳本試試。

3、運(yùn)行OrthoFinder2

orthofinder -f primary_transcripts/

OrthoFinder2結(jié)果解讀

1、多少基因在Orthogroup中

一般要求80%的gene可以分配到orthogroup中，

文件：Comparative_Genomics_Statistics/Statistics_Overall.tsv

2、每個(gè)物種中有多少基因在Orthogroup中

文件：Comparative_Genomics_Statistics/Statistics_PerSpecies.tsv

3、物種樹（species tree）

文件：Species_Tree/SpeciesTree_rooted.txt（自行計(jì)算了bootstrap值）

可選工具：

  • Dendroscope —— 勾選 “Interpret as edge labels”
  • ETE Toolkit tree viewer

需要注意的是，

  • 當(dāng)物種樹推斷錯(cuò)誤時(shí)，并不會對正交組（Orthogroup）的推斷產(chǎn)生影響，但是可能會對gene tree推斷得到的orthologue組成產(chǎn)生影響（e.g. 基因的復(fù)制事件）

  • 當(dāng)已有推斷正確的物種樹，可以使用-s/-ft參數(shù)進(jìn)行指定

  • OrthoFinder2默認(rèn)情況下，使用STAG算法（多序列比對的替代，一般具有計(jì)算得到更高的bootstrap values）進(jìn)行分析，該算法使用每個(gè)位點(diǎn)的gene tree來計(jì)算bootstrap值（并非標(biāo)準(zhǔn)的bootstrap計(jì)算方法）；

    如果進(jìn)行分析時(shí)，<u>想要使用標(biāo)準(zhǔn)bootstrap計(jì)算方法</u>，可指定-M msa（用于<font color='yellow'>串聯(lián)</font>所有的序列進(jìn)行比對），用于物種樹的推斷結(jié)果（結(jié)果中即可看到100% bootstrap值）

4、同源基因（Orthologues）

OrthoFinder2分析的主要功能之一 —— “find the orthologue of a gene you're interested in”

文件夾：Orthologues/，其下包含了對應(yīng)物種的gene的同源基因

e.g. Orthologues/Orthologues_Drosophila_melanogaster/Drosophila_melanogaster__v__Homo_sapiens.tsv由3列組成：1）Orthogroup，2）Drosophila_melanogaster，3）Homo_sapiens

e.g. 果蠅的FBgn0005648，對應(yīng)人類中的：ENSG00000205022、ENSG00000100836、ENSG00000258643

5、基因樹（gene tree）

文件夾：

  • Gene_Trees/

  • Resolved_Gene_Trees/（與上述文件夾的結(jié)果有一定區(qū)別，更加嚴(yán)格）

    針對該文件夾中的基因樹，OrthoFinder2執(zhí)行了“Duplication-Loss-Coalescence analysis”

需要注意的問題：

  • gene tree沒有bootstrap values

關(guān)于基因復(fù)制事件

文件：

  • Gene_Duplication_Events/SpeciesTree_Gene_Duplications_0.5_Support.txt：放入Dendroscope進(jìn)行可視乎
  • Gene_Duplication_Events/Duplications.tsv
  • Comparative_Genomics_Statistics/
    • Duplications_per_Orthogroup.tsv
    • Duplications_per_Species_Tree_Node.tsv

解讀Duplications.tsv的幾個(gè)點(diǎn)：

  • 當(dāng)超過50%以上的物種都對應(yīng)的基因以及重復(fù)產(chǎn)生的基因，認(rèn)為該基因重復(fù)事件是可信的

關(guān)于正交群（Orthogroups）：orthologs，paralogs，etc.

文件夾：

  • Orthogroups/Orthogroups.tsv，一行一個(gè)orthogroup，一列為一個(gè)物種（按），如下所示

    Orthogroup    Sp.A Sp.B Sp.C
OG0000000.fa  gene1   gene2   gene3

需要注意的是，<font color='yellow'>orthogroup可以包含ortholog以及paralog</font>。

  • Orthogroup_Sequences/，包含對應(yīng)正交群內(nèi)的序列

    # e.g.
OG0000000.fa
OG0000001.fa
OG0000002.fa
OG0000003.fa
OG0000004.fa
OG0000005.fa
OG0000006.fa
OG0000007.fa
OG0000008.fa
OG0000009.fa
# etc.

OrthoFinder2分析的注意事項(xiàng)

來自原作者的總結(jié)。

分析上的建議

  • 在進(jìn)行幾個(gè)物種的比較基因組分析時(shí)，不推薦加入外類群的氨基酸序列進(jìn)行OrthoFinder2分析。原因是，可能會影響正交群的分化時(shí)間推斷。
  • 想要保證推斷同源基因?qū)Φ臏?zhǔn)確性，需要保證一定的物種數(shù)量 —— 需要保證所選的物種親緣關(guān)系不是太遠(yuǎn)（you should break up long branches with intermediate species），<font color='yellow'>最低4個(gè)物種，6-10個(gè)最優(yōu)</font>。
  • 推薦使用每一個(gè)gene的最長轉(zhuǎn)錄本進(jìn)行分析（選擇這樣做的一個(gè)原因，也是因?yàn)檫@樣做有利于節(jié)省計(jì)算資源）。

如果基因組注釋的不好，怎么辦？

兩句話，

  • OrthoFinder2對出現(xiàn)missing genes的情況，還是能保證一個(gè)相對準(zhǔn)確的結(jié)果，所以不用擔(dān)心
  • 當(dāng)一個(gè)物種超過100000條轉(zhuǎn)錄本時(shí)，就很耗費(fèi)計(jì)算資源了！

分析細(xì)節(jié)上的修改：config.json 

find ~ -name "config.json"
vim PATH/config.json

MCL inflation factor的含義？

這個(gè)是OrthoFinder2從OrthoMCL繼承過來的參數(shù)，-I。

其默認(rèn)數(shù)值為1.5，修改它的意義在哪里？

  • 當(dāng)OrthoFinder2鑒定得到的蛋白質(zhì)序列，感覺關(guān)系不是特別大，

    比如在蛋白質(zhì)功能上天差地別，又或者說蛋白質(zhì)變異過大，即需要將-I設(shè)置得大一些（e.g. 5）

數(shù)值大的-I，意味著什么？

-I增大，代表了增大了聚類分析中的顆粒性（granularity），<u>可以得到更加小范圍的聚類結(jié)果</u>（也就是相當(dāng)于將聚類結(jié)果設(shè)置得更加嚴(yán)格）

Note：增加了granularity，也就是相當(dāng)于將聚類中心點(diǎn)增多，

比如原本能夠?qū)?個(gè)點(diǎn)聚類在一起，但是增加了granularity之后，就得到了2和1的聚類，3獨(dú)自聚類。

以下為我的猜想：

猜想1（錯(cuò)誤）：

反推到Orthogroup的鑒定上，假設(shè)當(dāng)前的輸入的蛋白質(zhì)序列來自四個(gè)物種，在默認(rèn)I值的基礎(chǔ)上，只有當(dāng)4個(gè)序列都滿聚類條件時(shí)，才可以被認(rèn)定為Orthogroup，

當(dāng)I值增大以后，當(dāng)出現(xiàn)3個(gè)滿足條件的序列，其也可以歸為一個(gè)Ortholog

猜想2：提升granularity，即允許更多聚類結(jié)果的呈現(xiàn)（增加了顆粒點(diǎn)）

[1] https://www.biostars.org/p/169145/#:~:text=OrthoMCL%20uses%20an%20inflation%20of,numbers.
[2] http://micans.org/mcl/man/mcl.html#opt-I
[3] OrthoFinder2 Issues
[4] Hierarchical clustering

參考資料

[1] https://davidemms.github.io/orthofinder_tutorials/running-an-example-orthofinder-analysis.html

[2] https://davidemms.github.io/orthofinder_tutorials/exploring-orthofinders-results.html

[3] https://github.com/davidemms/OrthoFinder#what-orthofinder-provides

[4] https://github.com/davidemms/OrthoFinder#orthogroups-orthologs--paralogs

拓展閱讀資料

別人的博客寫的太好了

[1] https://blog.csdn.net/sinat_41621566/article/details/112320002s【OrthoFinder2相關(guān)概念】

[2] https://blog.csdn.net/sinat_41621566/article/details/112320002

額外

哪里可以下載到植物的最長轉(zhuǎn)錄本？

JGI網(wǎng)站下的.protein_primaryTranscriptOnly.fa文件

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容