Orthofinder2鉴定单拷贝基因并构建进化树

实验目的:比较基因组学,从基因组获取单拷贝基因,并构建有根的进化树,并估计分化时间,最终画出有分化时间的进化树
实验流程:
orthofinder2获取单拷贝基因——muscle多序列比对——Gblocks寻找保守位点——iqtree2构建进化树——MCMCTREE计算分歧时间——
(1)单拷贝基因的查找 Orthofinder;
(2)多序列比对 Muscle / MAFFT / ClustalW / T-coffee, 本实验使用Muscle,准确度最高。 mafft的准确性真的低
(3)调取保守区域,并收尾连接,形成supergene ,使用软件Gblocks  
(4)进化树构建 RaxML,PhyML或Mrbayes,本实验使用iqtree2
构建树的教程:https://www.yuque.com/wusheng/gw7a9p/mcc73y
(5)分化时间分析 divergence time,本实验使用mcmctree 是PAML 软件包的一个功能。分化时间进化树可视化使用MCMCTreeR
(6)基因扩张收缩分析 CAFE
(7)基因正选择codeML PAML中一个程序


  1. 提取所有基因组每个基因的最长的转录本
    使用getLongProteins.py脚本获取,最后一个参数是指定相同的基因后的分隔符
    getLongProteins.py input.pep.fa out.pep.fa .
    checkFa.py *.pep.fa检测fa文件序列中是否有其他字符,如果有的话会输出其他的字符,需要提前删除这些其他字符。
    2.使用Orthofinder进行分析
    把上面处理后的蛋白文件放入OrthoFinder/primary文件夹中,或者是其他名称也行,注意修改下面的-f后文件夹名称。一般尽量把蛋白序列命名为abbr.pep.fa这种格式。

第一次运行使用

orthofinder -t 16 -a 16 -f OrthoFinder/primary
-t 16 -a 16是指定使用的cpu进程数量为16
-f是指定输入的蛋白文件的位置

如果你分析完成了,想再添加新的物种,只需要新建一个文件夹run3,把新加的物种蛋白序列放入run3。再运行下面的命令即可

orthofinder -t 16 -a 16 -b OrthoFinder/Results_Jul08 -f run3
-b参数是指定之前分析的输出的文件夹
-f是指定新增的蛋白序列的文件夹

如果是想删除一个物种,同时新增一些物种,命令还是下面这样

orthofinder -t 16 -a 16 -b OrthoFinder/Results_Jul08 -f run3
但是注意,需要进入 OrthoFinder/Results_Jul08/WorkingDirectory这个文件夹中,找到SpeciesIDs.txt这个文件,在需要删除的物种前加上#符号即可。新增的物种蛋白序列还是放入run3

orthofinder输出的结果目录文件结构说明

Species_Tree 目录 里面是物种的进化树
SpeciesTree_rooted.txt 基于所有蛋白序列构建的进化树(可能不是真实的进化关系,因为涉及到染色体倒位等结构变异)
SpeciesTree_rooted_node_labels.txt 和上面一样的树,有label
Single_Copy_Orthologue_Sequences 里面是单拷贝直系同源基因的fa文件,每个fa文件里是每一个物种的一条序列。
Orthogroups目录
Orthogroups_SingleCopyOrthologues.txt 单拷贝基因的OG开头的编号的列表文件。总行数就是每个物种里单拷贝基因的数量。
Comparative_Genomics_Statistics目录
OrthologuesStats_one-to-one.tsv 这个文件里可以看到每个基因组和其他基因组比较,1对1的直系同源基因的数量。

如果singlecopy的基因数量比较少,这时候需要使用参数过滤,例如:在70%的物种里是单拷贝的基因的列表,这样能够多一些。使用getSingleCopy.py Results_XXXXX/Orthogroups/Orthogroups.GeneCount.tsv all 0.7即可获取对应的基因家族的名称。
从Results_XXXXX/Orthogroup_Sequence里找到对应家族的序列。

3.使用singlecopy2tree.bash脚本完成比对和保守位点提取
bash singlecopy2tree.bash OrthoFinder/Results_Jul08/Single_Copy_Orthologue_Sequences orthofinder_inputpath

  1. 使用iqtree2构建单拷贝基因的进化树
##使用iqtree2构建进化树(此处需要16个cpu,其实少一点也可以,速度很快)
iqtree2 -s all.fa -T 16 -B 1000 -m MFP --bnni -o ${outgroup}

5.进化树进化时间估计
http://www.timetree.org/查询两个物种之间的分化时间。输入拉丁学名,即可查询到。

image.png

一般情况下,使用的是CI后面的两个值,对应图中的3.9和53.1。这个例子不太好,差的有点太大了。最好选择两个物种之间估计的分化时间差距不大的。注意不是两个物种分化时间不大,是CI后的括号内的值的差距小。
PAML使用的是CI后面的值,R8S使用的是Median Time值。
PAML 的mcmctree用于计算分化时间(速度比较慢,但是比较准)
R8S也可以用于计算分化时间(速度快,准确率没有mcmctree高)

#R8S安装方法
wget -q https://sourceforge.net/projects/r8s/files/r8s1.81.tar.gz \
  && tar -zxvf r8s1.81.tar.gz \
  && cd r8s1.81/src \
  && sed -i 's|/usr/include/sys/errno.h||' Makefile.linux \
  && sed -i 's/continuousML.o //' Makefile.linux \
  && sed -i 's/continuousML.o:/#continuousML.o:/' Makefile.linux \
  && make -f Makefile.linux

OrthoFinder/Results_Jul08/Orthogroups/Orthogroups.GeneCount.tsv:每一行对应一个基因家族在该基因组对应的基因数目;

获取所有的直系同源基因的ID列表

cat name.list|tr ">" "\t"  >SingleGeneID.tsv
awk '
{
    for (i=1; i<=NF; i++) {
        if (NR == 1) {
            transposed[i] = $i
        } else {
            transposed[i] = transposed[i] "\t" $i
        }
    }
}
END {
    for (i=1; i<=NF; i++) {
        print transposed[i]
    }
}' SingleGeneID.tsv > Specs.singlecopy.id
tail -1 Specs.singlecopy.id >header
sed -i '$ d' Specs.singlecopy.id
cat header Specs.singlecopy.id > Species.singlecopy.tsv
rm -rf SingleGeneID.tsv Specs.singlecopy.id header 

最终Species.singlecopy.tsv是所有的单拷贝基因的列表,第一行是每个物种的名称,每一列是对应的基因的ID。每一行之间是单拷贝直系同源基因。

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

推荐阅读更多精彩内容