OrthoFinder2算法解读

1）基因树的构建，

• 寻找序列之间的同源关系（相似度的计算）：diamond，blast，mmseq2
• 基因树的推断：DendroBLAST，其他的参考config.json

2）物种树的构建

• STAG：从无根gene tree推断无根物种树
• STRIDE：从无根物种树推断有根物种树

Note：

• STAG：Species Tree from All Genes (STAG)

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

一点点背景知识

什么是正交组（Orthogroup）？

一个祖先物种经历过一次物种事件之后，分别形成了物种A和物种B，

祖先物种gene分化得到geneA和geneB —— geneA和geneB，称为正交群（Orthogroup）

Note：paralogs，也可以和orthologs一起，存在于同一个Orthogroup中

参考资料，

[1] 来自lakeseafly的简书笔记
[2] https://en.wiktionary.org/wiki/orthogroup

为什么要研究正交群？

• 正交群中的所有基因都来自<font color='yellow'>单个祖先基因</font> —— 正交群中的所有基因都有类似的序列和功能

  由于基因重复和丢失在进化中经常发生，<u>1 vs 1 的直系同源物非常少见</u>，

  因此，分析正交群的所有情况（1 vs 1，1 vs more，more vs more）可以得到有关进化历史的信息。

• “真正识别直系同源物的唯一方法，是构建系统发育树” —— lakeseafly

软件安装

conda create -n orthofinder2 -c bioconda orthofinder -y

软件使用

以mouse、人类、斑马鱼、青蛙（frog？）、日本河豚、果蝇的蛋白质序列为例

1、数据下载

需要注意的是，

• OrthoFinder2要求输入的序列为氨基酸（蛋白质）序列（CDS对应的序列，即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

将上述命令保存为脚本并运行，

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

2、提取最长转录本

这一步，主要是为了减少OrthoFinder2的运行时间，可以自己写脚本试试。

3、运行OrthoFinder2

orthofinder -f primary_transcripts/

OrthoFinder2结果解读

1、多少基因在Orthogroup中

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

文件：Comparative_Genomics_Statistics/Statistics_Overall.tsv

2、每个物种中有多少基因在Orthogroup中

文件：Comparative_Genomics_Statistics/Statistics_PerSpecies.tsv

3、物种树（species tree）

文件：Species_Tree/SpeciesTree_rooted.txt（自行计算了bootstrap值）

可选工具：

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

需要注意的是，

• 当物种树推断错误时，并不会对正交组（Orthogroup）的推断产生影响，但是可能会对gene tree推断得到的orthologue组成产生影响（e.g. 基因的复制事件）

• 当已有推断正确的物种树，可以使用-s/-ft参数进行指定

• OrthoFinder2默认情况下，使用STAG算法（多序列比对的替代，一般具有计算得到更高的bootstrap values）进行分析，该算法使用每个位点的gene tree来计算bootstrap值（并非标准的bootstrap计算方法）；

  如果进行分析时，<u>想要使用标准bootstrap计算方法</u>，可指定-M msa（用于<font color='yellow'>串联</font>所有的序列进行比对），用于物种树的推断结果（结果中即可看到100% bootstrap值）

4、同源基因（Orthologues）

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

文件夹：Orthologues/，其下包含了对应物种的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，对应人类中的：ENSG00000205022、ENSG00000100836、ENSG00000258643

5、基因树（gene tree）

文件夹：

• Gene_Trees/

• Resolved_Gene_Trees/（与上述文件夹的结果有一定区别，更加严格）

  针对该文件夹中的基因树，OrthoFinder2执行了“Duplication-Loss-Coalescence analysis”

需要注意的问题：

• gene tree没有bootstrap values

关于基因复制事件

文件：

• Gene_Duplication_Events/SpeciesTree_Gene_Duplications_0.5_Support.txt：放入Dendroscope进行可视乎
• Gene_Duplication_Events/Duplications.tsv
• Comparative_Genomics_Statistics/
  • Duplications_per_Orthogroup.tsv
  • Duplications_per_Species_Tree_Node.tsv

解读Duplications.tsv的几个点：

• 当超过50%以上的物种都对应的基因以及重复产生的基因，认为该基因重复事件是可信的

关于正交群（Orthogroups）：orthologs，paralogs，etc.

文件夹：

• Orthogroups/Orthogroups.tsv，一行一个orthogroup，一列为一个物种（按），如下所示

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

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

• Orthogroup_Sequences/，包含对应正交群内的序列

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

OrthoFinder2分析的注意事项

来自原作者的总结。

分析上的建议

• 在进行几个物种的比较基因组分析时，不推荐加入外类群的氨基酸序列进行OrthoFinder2分析。原因是，可能会影响正交群的分化时间推断。
• 想要保证推断同源基因对的准确性，需要保证一定的物种数量 —— 需要保证所选的物种亲缘关系不是太远（you should break up long branches with intermediate species），<font color='yellow'>最低4个物种，6-10个最优</font>。
• 推荐使用每一个gene的最长转录本进行分析（选择这样做的一个原因，也是因为这样做有利于节省计算资源）。

如果基因组注释的不好，怎么办？

两句话，

• OrthoFinder2对出现missing genes的情况，还是能保证一个相对准确的结果，所以不用担心
• 当一个物种超过100000条转录本时，就很耗费计算资源了！

分析细节上的修改：config.json

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

MCL inflation factor的含义？

这个是OrthoFinder2从OrthoMCL继承过来的参数，-I。

其默认数值为1.5，修改它的意义在哪里？

• 当OrthoFinder2鉴定得到的蛋白质序列，感觉关系不是特别大，

  比如在蛋白质功能上天差地别，又或者说蛋白质变异过大，即需要将-I设置得大一些（e.g. 5）

数值大的-I，意味着什么？

-I增大，代表了增大了聚类分析中的颗粒性（granularity），<u>可以得到更加小范围的聚类结果</u>（也就是相当于将聚类结果设置得更加严格）

Note：增加了granularity，也就是相当于将聚类中心点增多，

比如原本能够将3个点聚类在一起，但是增加了granularity之后，就得到了2和1的聚类，3独自聚类。

以下为我的猜想：

猜想1（错误）：

反推到Orthogroup的鉴定上，假设当前的输入的蛋白质序列来自四个物种，在默认I值的基础上，只有当4个序列都满聚类条件时，才可以被认定为Orthogroup，

当I值增大以后，当出现3个满足条件的序列，其也可以归为一个Ortholog

猜想2：提升granularity，即允许更多聚类结果的呈现（增加了颗粒点）

[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相关概念】

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

额外

哪里可以下载到植物的最长转录本？

JGI网站下的.protein_primaryTranscriptOnly.fa文件