OrthoFinder寻找同源基因并建树

  • 提到寻找物种之间的同源基因,大家一般想到的是OrthoMCL (http://orthomcl.org/orthomcl/) ,OrthoMCL是现在用的最多的一款来找直系同源基因(Orthologs)以及旁系同源基因 (Paralog) 的软件。
  • 但是OrthoMCL 的最新版本是2013年7月公布的v2.0版本,已经很久没更新过了。根据官网的教程至少得十多步才能完成整个运行流程,包括Mysql数据库配置、修改OrthoMCL配置文件、转换序列格式、过滤、比对、解析结果和聚类等步骤,特别繁琐。

于是OrthoFinder(https://github.com/davidemms/OrthoFinder)这个软件应运而生,OrthoFinder是一个快速寻找同源基因,并可选择快速建树的软件,基于python3编写。该软件安装运行十分简单,且自2015年发布第一个版本,到目前仍在持续更新。

一、OrthoFinder的安装

1.通过anaconda安装(推荐)

anaconda应该是生信分析中必不可少的软件,具体在此不多赘述,anaconda安装教程移步(https://www.jianshu.com/p/62f155eb6ac5),在anaconda中安装OrthoFinder的命令为

conda install -y orthofinder

安装完成后输入orthofinder可查看其帮助文档即为成功


OrthoFinder version 2.2.7 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 = 16]
 -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 = blast]
                   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
 -n <txt>          Name to append to the 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

2.通过下载安装

1.从GitHub下载最新版本的下载OrthoFinder 下载(这里我以OrthoFinder-2.2.7.tar.gz为例)

2.在终端解压orthofinder压缩包:

tar xzf OrthoFinder-2.2.7.tar.gz

3.安装OrthoFinder的依赖程序:
这里需要安装的有:MCL, FastME and DIAMOND
具体程序下载和安装教程可在其官网查看:
MCL: http://micans.org/mcl/
FastME:http://www.atgc-montpellier.fr/fastme/binaries.php
DIAMOND:https://github.com/bbuchfink/diamond/releases

4.测试OrthoFinder:
用终端打开OrthoFinder所在文件夹,运行:

./orthofinder -f ExampleDataset -S diamond

5.加入orthofinder到环境变量(可选)

二、OrthoFinder的使用

OrthoFinder需要蛋白编码的Fasta文件作为输入,所以首先我们准备好要比对的物种的蛋白序列,存入一个文件夹。

1.典型用法:

orthofinder -f Dataset

其中:
-f 指定蛋白序列所在的文件夹(这里我的文件夹名为Dataset)
此命令仅需指定输入文件夹,其余为默认参数。

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)

有时候我们不需要运行整个从蛋白序列到建树的整个流程,我们能只需要运行其中的部分并得到其结果,那么就可以选择分步运行

分步运行:

停止选项

-op blast比对后停止运行

-og 推断同源组后停止运行

-os 同源组写入序列文件后停止运行(需要-M msa选项)

-oa 推断同源组序列比对后停止运行(需要-M msa选项)

-ot 推断同源组基因树后停止运行

开始选项

-b <blast结果文件夹 > 从blast结果开始运行

-fg <同源组结果文件夹> 从同源组结果开始运行

-ft <基因树文件夹> 从基因树文件开始运行

最后说一下我一般常用的参数,我做动物研究,一般数据较多较大,比对我选择Diamond,建树选择fasttree,多序列比对msa,线程数由于设备限制,设置为16

orthofinder -f Dataset -M msa -S diamond -t 16 -a 16 

可参照下图更直观理解:

参数图解

3.高级用法:

通过上面的参数,我们可以做一个从蛋白序列到同源基因建树的流程,但是仍然有一些参数无法设置,比如建树过程中的bootstrap。这时候,我们就要修改程序的config.json文件。

如果你是用anaconda安装的orthofinder,那么config.json文件在:

anaconda2/bin/config.json

如果你通过下载压缩包安装,config.json文件就在orthofinder目录里。

首先我们打开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_*** 文件夹,其中***是运行日期

直系同源组相关结果文件,将不同的直系同源基因进行分组:
  • Orthogroups.csv:用制表符分隔的文件,每一行是直系同源基因组对应的基因。
  • Orthogroups.txt: 类似于Orthogroups.csv,只不过是OrhtoMCL的输出格式
  • Orthogroups_UnassignedGenes.csv: 格式同Orthogroups.csv,只不过是物种特异性的基因
  • Orthogroups.GeneCount.csv:格式同Orthogroups.csv, 只不过不再是基因名信息,而是以基因数。
直系同源相关文件,分析每个直系同源基因组里的直系同源基因之间关系,结果会在Orthologues_***文件夹下,其中***是日期:
  • Gene_Trees: 每个直系同源基因基因组里的基因树
  • Recon_Gene_Trees:使用OrthoFinder duplication-loss coalescent 模型进行发育树推断
  • Potential_Rooted_Species_Trees: 可能的有根物种树
  • SpeciesTree_rooted.txt: 从所有包含STAG支持的直系同源组推断的STAG物种树
  • SpeciesTree_rooted_node_labels.txt: 同上,只不过多了一个标签信息,用于解释基因重复数据。
比较基因组学的相关结果文件:
  • Orthogroups_SpeciesOverlaps.csv: 不同物种间的同源基因的交集
  • SingleCopyOrthogroups.txt: 单基因拷贝组的编号
  • Statistics_Overall.csv:总体统计信息
  • Statistics_PerSpecies.csv:分物种统计信息

最后就可以利用以上结果,继续愉快地做分析啦~~~

参考文章:

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,761评论 5 460
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,953评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,998评论 0 320
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,248评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,130评论 4 356
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,145评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,550评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,236评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,510评论 1 291
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,601评论 2 310
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,376评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,247评论 3 313
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,613评论 3 299
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,911评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,191评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,532评论 2 342
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,739评论 2 335