0.忘了是在那篇文献看到的，用了Agora，再把Agora的输出结果喂给ANGEs
Agora和ANGEs都是构建祖先核型的软件
个人实测是Agora构建的结果太碎了, ANGEs运行速度又太慢(大于半个月)
所以用Agora将同源基因初步聚类, 以这些blocks喂给ANGEs能大幅提高运行速度
软件的安装参考官方文档
Agora: https://github.com/DyogenIBENS/Agora
ANGEs: https://cchauve.github.io/ANGES/index.html
这里仅作一个流程的记录, 所有疑问最好去官方的manual中寻找答案
1.Orthofinder (这里也可以用序列比对的结果而不是基因的结果, 但是亲测用基因的结果gap会比较少)
输入文件是各个物种的 .pep文件，以运行orthofinder，最好给各个物种的序列重新命名，以便区分不同物种，图方便物种树也用orthofinder的结果。当然自己重新构建物种树也可以
可以参考https://www.jianshu.com/p/336b65ca1b67
物种树的枝长不会影响Agora，ANGEs不清楚
这一步骤略
2.将OrthoFinder的结果处理为Agora要求的输入格式
Agora的输入文件有三类:
2.1 第一个是物种树
并对进化树的每个节点进行命名, 由于我的物种数量比较少, 所以命名是我手动加的
如图:
将三个节点分别命名为A0/A1/A2
同理对应的进化树的文本格式就是:
((sppA:0.1, sppB:0.2)A1:0.2, (sppC:0.3, sppD:0.4)A2:0.2)A0;
((sppA, sppB)A1, (sppC, sppD)A2)A0; 将枝长省略更为直观
2.2 第二个是每个末端节点(每个物种)对应的blocks的位置信息, 这里我们用的是每个同源基因的bed文件
前三列是位置信息, 第四列是正负链(1是正链/-1是负链), 第五列是基因名
gene.sppA.list
文件命名为 gene.sppA.list
这个文件很容易通过gff文件获取
比如对应我手里的这个gff文件:
sppA.gff
grep -w 'mRNA' sppA.gff > tmp1
awk 'BEGIN{ FS="\t";OFS="\t" }{print $1,$4,$5,$7,$9}' | awk -F ";" '{print $1}' | sed 's/ID=//g' > tmp2
sed 's/-/-1/g' tmp2 | sed 's/+/1/g' > gene.sppA.list
#对每个物种重复此操作
#为了方便管理, 把这类输入文件都专门放到一个GeneList的文件夹中
2.3 第三个是每个节点包括祖先节点对应的blocks, 这里我们用的是每个节点的同源基因
文件内容也很简单, 每一行都是一个基因家族(OG)所包含的基因，空格分隔符, 如下图:
orthologyGroups.A1.list
获取方法有两种
2.3.1 OrthoFinder的运行结果Phylogenetic_Hierarchical_Orthogroups文件夹中有按照节点划分的基因家族, 我在前文中命名为A1的节点, 被OrthoFinder命名为了N1节点. 也可以在前面命名节点的时候就按照OrthoFinder的来, 这样省去一些文本处理的工作
这个文件已经符合要求了, 只需要把无关的列删除, 并把分隔符转为空格分隔符即可
PS: 不知道 竖线 "|" 会不会产生影响, 这里 "物种名" + "|" + "基因名" 的命名格式是做其他分析用到的, 真实运行的时候是没有这个格式
PS: 只要确保不同物种的基因名独一无二即可
OrthoFinder的输出结果 - N1.tsv
sed 's/^\([^ \t]\+[ \t]\+\)\{3\}//' N1.tsv | sed 's/,/ /g' > tmp1 #去除文件的前三列并删除多余的逗号
sed 's/[ \t]\+/ /g' tmp1 | sed '1d' > orthologyGroups.A1.list #将多个tab/空格转换为一个空格, 并去除第一行
rm tmp1
sed -i '/^[[:space:]]*$/d' orthologyGroups.A1.list #删除空白行
sed -i 's/^[\t ]\+//' orthologyGroups.A1.list #删除行首多余的空格/tab
#对每个节点重复这个操作
#为了方便管理, 把这类输入文件都专门放到一个orthologyGroup的文件夹中
2.3.2 自己手动提一下(不建议用这个方法)
A1节点是sppA和sppB的祖先节点, 所以我们只需要把这两个物种的同源基因的信息放在一起即可
用到的是Orthofinder输出的Orthogroups/Orthogroups.tsv
这个文件内容其实差不多, 每一行都是一个基因家族(OG)以及这个基因家族对应的基因名
sppA和sppB在文件的第四列和第五列
Orthogroups.tsv
awk 'BEGIN{ FS="\t";OFS="\t" }{print $4,$5}' Orthogroups.tsv | sed 's/,/ /g' > tmp1 #提取sppA和sppB所在的列并删除多余的逗号
sed 's/[ \t]\+/ /g' tmp1 | sed '1d' > orthologyGroups.A1.list #将多个tab/空格转换为一个空格, 并去除第一行
rm tmp1
sed -i '/^[[:space:]]*$/d' orthologyGroups.A1.list #删除空白行
sed -i 's/^[\t ]\+//' orthologyGroups.A1.list #删除行首多余的空格/tab
#对每个节点重复这个操作
#同理如果是前文提到的A0节点也就是所有物种的共同祖先节点
#那就提取4个物种的信息
awk 'BEGIN{ FS="\t";OFS="\t" }{print $2,$3,$4,$5}' Orthogroups.tsv | sed 's/,/ /g' > tmp1 #提取四个物种的信息
#为了方便管理, 把这类输入文件都专门放到一个orthologyGroup的文件夹中
3.Agora
./software/Agora-master/src/agora-vertebrates.py \ #因为我的物种是脊椎动物, 所以用这个脚本
SpeciesTree_rerooted.txt \ #输入文件1
orthologyGroup/orthologyGroups.%s.list \ #输入文件3
GeneList/genes.%s.list \ #输入文件2
-workingDir=Agora.result #输出的文件夹
-target=A0 \ #目标节点也就是想看哪个节点的祖先核型
2> all.log #日志文件
虽然命令行指定了想看哪个节点的祖先核型, 但是Agora仍然会输出每个节点的祖先核型
可以先检查Agora的结果是否符合预期，如果符合预期就不用再往下做了
Agora也自带了结果的绘制与解读，详情请参阅官方文档
我使用Agora的效果并不好, 所以并没有过多的阅读相关部分
Agora.result/ancGenomes/vertebrates-workflow这个文件夹就是最后的结果
或者说可以用作ANGEs的输入文件
输出文件长这样:
Agora的输出文件 - ancGenome.A0.list.bz2
第一列就是构建的祖先核型，CAR_1表示这是祖先核型中最长的那一条染色体
第二列/第三列是祖先染色体上基因的顺序; 第四列是正负链
第六列即祖先染色体上的这个基因在现存物种中的同源基因
这个文件继续往下翻阅
会看到我的结果已经到了CAR_316, 祖先核型当然不可能有300条染色体
也没有说祖先染色体上才个位数的基因数量
所以这个软件在我的数据集上的结果并不理想
需要进一步使用ANGEs
网友评论