使用Agora/ANGEs构建祖先核型

作者: 深山夕照深秋雨OvO | 来源:发表于2023-12-14 19:11 被阅读0次

Agora iOS SDK-开始聊天
Agora iOS SDK-多人聊天
Agora SDK使用小结
单子叶祖先 7条染色体
声网(Agora)视频通话-小程序
分子进化理论（4）
Agora
agora
音视频开发
声网 Agora Unity SDK 全功能版正式上线

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的输入文件
输出文件长这样: