一、OrthoFinder的安装
$conda install -y orthofinder
使用之前记得查看帮助文档:
$orthofinder
OrthoFinder version 2.3.3 Copyright (C) 2014 David Emms
SIMPLE USAGE:
Run full OrthoFinder analysis on FASTA format proteomes in <dir>
orthofinder [options] -f <dir>
Add new species in <dir1> to previous run in <dir2> and run new analysis
orthofinder [options] -f <dir1> -b <dir2>
OPTIONS:
-t <int> Number of parallel sequence search threads [Default = 64]
-a <int> Number of parallel analysis threads [Default = 1]
-M <txt> Method for gene tree inference. Options 'dendroblast' & 'msa'
[Default = dendroblast]
-S <txt> Sequence search program [Default = diamond]
Options: blast, mmseqs, blast_gz, diamond
-A <txt> MSA program, requires '-M msa' [Default = mafft]
Options: muscle, mafft
-T <txt> Tree inference method, requires '-M msa' [Default = fasttree]
Options: iqtree, raxml-ng, fasttree, raxml
-s <file> User-specified rooted species tree
-I <int> MCL inflation parameter [Default = 1.5]
-x <file> Info for outputting results in OrthoXML format
-p <dir> Write the temporary pickle files to <dir>
-1 Only perform one-way sequence search
-X Don't add species names to sequence IDs
-n <txt> Name to append to the results directory
-o <txt> Non-default results directory
-h Print this help text
WORKFLOW STOPPING OPTIONS:
-op Stop after preparing input files for BLAST
-og Stop after inferring orthogroups
-os Stop after writing sequence files for orthogroups
(requires '-M msa')
-oa Stop after inferring alignments for orthogroups
(requires '-M msa')
-ot Stop after inferring gene trees for orthogroups
WORKFLOW RESTART COMMANDS:
-b <dir> Start OrthoFinder from pre-computed BLAST results in <dir>
-fg <dir> Start OrthoFinder from pre-computed orthogroups in <dir>
-ft <dir> Start OrthoFinder from pre-computed gene trees in <dir>
LICENSE:
Distributed under the GNU General Public License (GPLv3). See License.md
CITATION:
When publishing work that uses OrthoFinder please cite:
Emms D.M. & Kelly S. (2015), Genome Biology 16:157
If you use the species tree in your work then please also cite:
Emms D.M. & Kelly S. (2017), MBE 34(12): 3267-3278
Emms D.M. & Kelly S. (2018), bioRxiv https://doi.org/10.1101/267914
二、OrthoFinder的使用
OrthoFinder需要蛋白编码的Fasta（fa）文件作为输入，所以首先我们准备好要比对的物种的蛋白序列（至少要2个），存入一个文件夹。
我这里有9个需要比对的物种的fa文件：
1.典型用法：
$orthofinder -f orthofinder
其中：
-f 指定蛋白序列所在的文件夹（这里我的文件夹名为orthofinder）
此命令仅需指定输入文件夹，其余为默认参数。
2.进阶用法：
在OrthoFinder中我们可以通过不同的参数指定比对过程中使用的方法和其他参数设置，下面列出常用的几个参数：
常用参数：-t 序列搜索使用的线程数 （默认值为16）
-a 序列分析使用的线程数 （默认值为1）
-M 基因树推断方法（默认为dendroblast）可选：dendroblast ，msa
-S 序列比对使用的程序 （默认为Blast）可选：blast, mmseqs, blast_gz, diamond（推荐使用diamond，比对速度快）
-A 多序列联配方式，该选项仅当 -M msa 选项时才有效（默认为mafft）可选：muscle, mafft
-T 建树方式，该选项仅当 -M msa 选项时才有效 （默认为fasttree）可选：iqtree, raxml-ng, fasttree, raxml
-s 输入特定的根物种树
-I 设定MCL的膨胀系数（默认为1.5）
我做的研究，一般数据较多较大，比对我选择Diamond，建树选择fasttree,多序列比对msa，线程数由于设备限制，设置为16
$orthofinder -f orthofinder -M msa -S diamond -t 16 -a 16
3.高级用法：
通过上面的参数，我们可以做一个从蛋白序列到同源基因建树的流程，但是仍然有一些参数无法设置，比如建树过程中的bootstrap。这时候，我们就要修改程序的config.json文件。
如果你是用anaconda安装的orthofinder，那么config.json文件在你安装这个软件的conda环境里面的bin文件里：
首先我们打开config.json文件
{
"__comment": "Variable names that can be used:",
"__comment": "INPUT : The full path of the input filename (fasta file of sequences for and msa method, multiple sequence alignment fasta file for tree method)",
"__comment": "BASENAME : Just the filename without the directory path (a number of methods use this to name the output file automatically, see MergeAlign command for an example)",
"__comment": "PATH : Path to the directory containing the input file",
"__comment": "OUTPUT : The full path of the user specified output filename",
"__comment": "BASEOUTNAME : Just the filename without the directory path (of the output filename)",
"__comment": "IDENTIFIER : A name generated by OrthoFinder to uniquely identify the orthogroup (a number of methods use this to name the output file automatically, see RAxML command for an example). Not applicable for program_type search.",
"__comment": "DATABASE : For the search program_type, for use in the search_cmd. The full path of the database to search against",
"muscle":{
"program_type": "msa",
"cmd_line": "muscle -in INPUT -out OUTPUT"
},
"raxml":{
"program_type": "tree",
"cmd_line": "raxmlHPC-AVX -m PROTGAMMALG -p 12345 -s INPUT -n IDENTIFIER -w PATH > /dev/null",
"ouput_filename": "PATH/RAxML_bestTree.IDENTIFIER"
},
"raxml-ng":{
"program_type": "tree",
"cmd_line": "raxml-ng --msa INPUT --model LG+G4 --seed 12345 --threads 1",
"ouput_filename": "INPUT.raxml.bestTree"
},
"iqtree":{
"program_type": "tree",
"cmd_line": "iqtree -s INPUT -pre PATH/IDENTIFIER > /dev/null",
"ouput_filename": "PATH/IDENTIFIER.treefile"
},
"diamond":{
"program_type": "search",
"db_cmd": "diamond makedb --in INPUT -d OUTPUT",
"search_cmd": "diamond blastp -d DATABASE -q INPUT -o OUTPUT --more-sensitive -p 1 --quiet -e 0.001 --compress 1"
},
"blast_gz":{
"program_type": "search",
"db_cmd": "makeblastdb -dbtype prot -in INPUT -out OUTPUT",
"search_cmd": "blastp -outfmt 6 -evalue 0.001 -query INPUT -db DATABASE | gzip > OUTPUT.gz"
},
"mmseqs":{
"program_type": "search",
"db_cmd": "mmseqs createdb INPUT OUTPUT.fa ; mmseqs createindex OUTPUT.fa /tmp",
"search_cmd": "mmseqs search PATH/mmseqsDBBASENAME DATABASE.fa OUTPUT.db /tmp/tmpBASEOUTNAME --threads 1 ; mmseqs convertalis PATH/mmseqsDBBASENAME DATABASE.fa OUTPUT.db OUTPUT"
}
}
文件开头介绍了修改方法，你可以添加任何你想在流程中使用的软件，只要按照其格式进行命令修改就行。这里主要说一下对其中缘由参数的修改。
如要修改比对过程中的E-value值，那么在相应比对命令里修改数值，blast_gz 里的 -evalue ，diamond 里的 -e 。
如要在iqtree建树过程中增加bootstrap， 则在iqtree的"cmd_line":中添加
-bb 1000 （iqtree的超快bootstrap）或 -b 1000（传统bootstrap）
raxml-ng类似，你可以添加任何你想修改的参数，大家可以在熟悉各软件的的参数后，根据自己的需求更改。
三、OrthoFinder的结果解读
程序运行结束后，会在你指定的蛋白序列文件夹中生成一个Results_*** 文件夹，其中***是运行日期
(1) Results Files: Orthogroups
Orthogroups.tsv 用制表符分隔的文件，每一行是直系同源基因组对应的基因
Orthogroups.txt 类似于Orthogroups.csv，只不过是OrhtoMCL的输出格式
Orthogroups_UnassignedGenes.tsv 格式同Orthogroups.tsv，只不过是物种特异性的基因
Orthogroups.GeneCount.tsv 格式同Orthogroups.tsv, 只不过不再是基因名信息，而是以基因数
对于Orthogroups.GeneCount.tsv的结果，我们可以提取出来
我们选出各个物种中基因数大于0的基因家族，首先看物种1
我们不要第一行，然后看物种1，也就是$2，选出大于0的，然后我们需要的是基因家族编号，也即是第一列
$sed '1d' Orthogroups.GeneCount.csv |awk '$2 >0 {print $1}' >1.txt
(2) Results Files: Gene Trees and Species Tree
Gene_Trees: 每个直系同源基因基因组里的基因树
Recon_Gene_Trees：使用OrthoFinder duplication-loss coalescent 模型进行发育树推断
Potential_Rooted_Species_Trees: 可能的有根物种树
SpeciesTree_rooted.txt: 从所有包含STAG支持的直系同源组推断的STAG物种树
SpeciesTree_rooted_node_labels.txt: 同上，只不过多了一个标签信息，用于解释基因重复数据。
（3）Results Files: Orthogroup Statistics
Orthogroups_SpeciesOverlaps.csv： 不同物种间的同源基因的交集
SingleCopyOrthogroups.txt： 单基因拷贝组的编号
Statistics_Overall.csv：总体统计信息
Statistics_PerSpecies.csv：分物种统计信息
四： 高级用法
（1）添加新物种到之前的分析（previous_orthofinder_directory指的是包含“SpeciesIDs.txt”的目录）
$orthofinder -b previous_orthofinder_directory -f new_fasta_directory
（2）从之前的分析中移除物种
从输出目录下找到工作目录“WorkingDirectory”中的“SpeciesIDs.txt”文件，在要移除的物种那一行最前面加上一个“#”并保存，然后运行（previous_orthofinder_directory指的是包含“SpeciesIDs.txt”的目录）：
$orthofinder -b previous_orthofinder_directory
（3）同时添加和删除物种
编辑好“SpeciesIDs.txt”后，运行：
$orthofinder -b previous_orthofinder_directory -f new_fasta_directory
（4）更多高级功能请阅读官方文档，主要包括“Inferring MSA Gene Trees”、并行计算、单独运行BLAST、使用预先计算的BLAST结果以及回归检测。
