仅作个人笔记 – wgdi推断祖先核型

0.仅作个人笔记使用
意思就是只考虑我能不能看懂
具体每个参数的含义参照官方文档, 这里全部用默认参数
https://wgdi.readthedocs.io/en/latest/Introduction.html

我的研究物种是A, B, C, D
其系统发生关系是:

正好是两个姐妹物种

1.挑选参考基因组 - 要求是组装质量过关/对于研究物种来说是外群
(按照我问到的是, 不同科的物种; 但是我在我的研究类群中发现是不同属的更好; 也就是物种E, F, G; 比如我选了sppE)

2.输入文件的准备
A, B, C, D以及参考基因组/外群 sppE的基因注释文件: .gff, 并提取最长转录本: .cds/.pep
iTools_Code/iTools Fatools getCdsPep -Ref A.genome.fa -Gff A.gff -OutPut 取个名

用到的脚本是软件开发者提供的
https://github.com/SunPengChuan/wgdi-example/tree/main/code

#4个物种+1个外群都是相同的处理
python 01.getgff.py A.gff A.gff2
python 02.gff_lens.py A.gff2 A A.gff3 A.length

#A.gff3和A.length是软件要求的输入文件, 见下图, 这里我们修改第一列
#确保每个物种的染色体名是独一无二的,
awk 'BEGIN{ FS="/t";OFS="/t" }{print "A"$1,$2,$3,$4,$5,$6,$7}'  A.gff3 > A.gff4
awk 'BEGIN{ FS="/t";OFS="/t" }{print "A"$1,$2,$3}'   A.length > A.length2

#对最长转录本进行重命名
python 03.seq_newname.py A.gff4 A.cds A.rename.cds
python 03.seq_newname.py A.gff4 A.pep A.rename.pep

#
A.gff4,  A.length2,   A.rename.pep,    A.rename.cds就是后面用到的输出文件

A.gff2

A.gff4

A.length2

A.rename.pep/稍微注意一下标注的地方

3. A, B, C, D分别和sppE做共线性点图以及拿到共线性块
   实际证明不同程度亲缘关系远近的都做一下; 按照我的例子是sppE/sppF/sppG都要做
   总共12个共线性点图
   另外sppE/sppF/sppG, 三者之间也要做, 共计3个共线性点图, 外群可以酌情增加或者减少

以上15个共线性点图就是祖先核型推断所用到的全部

以物种A(fonta)和sppE为例:

#后续配置文件中具体参数说明参加https://wgdi.readthedocs.io/en/latest/dotplot.html

makeblastdb -in sppE.rename.pep -dbtype prot -out sppE.pep
blastp -query A.rename.pep -out A.sppE.blastp -db sppE.pep -evalue 1e-5 -num_threads 30 -outfmt 6
#输出的A.sppE.blastp即后续文件

3.1 共线性点图的获取(只做这一步就有图了, 没有特殊要求不用往下做了)
vim A.sppE.dot.conf
#输入以下内容
[dotplot]
blast = A.sppE.blastp  #Blast的输出结果
blast_reverse = false
gff1 = A.gff4    #前面准备好的输入文件
gff2 = sppE.gff4    #前面准备好的输入文件
lens1 = A.length2     #前面准备好的输入文件
lens2 = sppE.length2     #前面准备好的输入文件
genome1_name =  A    #共线性点图的名字(纵坐标对应的物种)
genome2_name =  sppE   #共线性点图的名字(横坐标对应的物种)
multiple  = 1
score = 100
evalue = 1e-5
repeat_number = 10
position = order
ancestor_left = none
ancestor_top = none
markersize = 0.5
figsize = 10,10
savefig = A_sppE.dot.pdf   #输出文件

wgdi -d A.sppE.dot.conf
#A_sppE.dot.pdf即最后的结果
#单从画共线性点图来看, 比jcvi繁琐太多了

3.2 extracting collinearity
vim A.sppE.collinearity.conf
#输入以下内容
[collinearity]
gff1 = A.gff4    #前面准备好的输入文件
gff2 = sppE.gff4    #前面准备好的输入文件
lens1 = A.length2     #前面准备好的输入文件
lens2 = sppE.length2     #前面准备好的输入文件
blast = A.sppE.blastp  #Blast的输出结果
blast_reverse = false
multiple  = 1
process = 8
evalue = 1e-5
score = 100
grading = 50,40,25
mg = 40,40
pvalue = 0.2
repeat_number = 10
positon = order
savefile = A.sppE.collinearity  #输出文件


wgdi -icl A.sppE.collinearity.conf


3.3 block的提取 #因为我是动物, 没有ks分析的需求所以跳过
vim A.sppE.blockinfo.conf
#输入以下内容
[blockinfo]
blast = A.sppE.blastp  #Blast的输出结果
gff1 = A.gff4    #前面准备好的输入文件
gff2 = sppE.gff4    #前面准备好的输入文件
lens1 = A.length2     #前面准备好的输入文件
lens2 = sppE.length2     #前面准备好的输入文件
collinearity = A.sppE.collinearity   #3.2的输出结果
score = 100
evalue = 1e-5
repeat_number = 10
position = order
ks = none   #没有跑ks, 所以是none
ks_col = ks_NG86
savefile = A.sppE.blockinfo.csv  #输出文件


wgdi -bi A.sppE.blockinfo.conf


3.4 Correspondence
vim A.sppE.correspondence.conf
#输入以下内容
[correspondence]
blockinfo = A.sppE.blockinfo.csv
lens1 = A.length2     #前面准备好的输入文件
lens2 = sppE.length2     #前面准备好的输入文件
tandem = true
tandem_length = 200
pvalue = 0.2
block_length = 5
multiple  = 1
homo = 0,1
savefile = A.sppE.blockinfo.corre.csv  #输出文件


wgdi -c A.sppE.correspondence.conf

4.祖先核型文件ancestor file的准备 - 最重要最耗时的一步
ref: https://github.com/SunPengChuan/wgdi-example/blob/main/Karyotype_Evolution.md

大体思路是分别用A, B, C, D 和 参考基因组/外群sppE作共线性点图, 找到4个物种共有的核型变化
比如ABCD的1号染色体和sppE的1号染色体共线性都特别好特别完整
那么可以断定这就是祖先核型的Anc1号染色体

几个要点:
基于共线性点图推导出核型变化
研究物种共有的变化, 即其共同祖先的变化
祖先核型数量已知;
如果未知, 并不能像ANGEs/Agora等等祖先核型软件一样自动帮你推一个核型数量; 这个时候需要自己根据共线性情况判断一个核型数量, 原则就是可以解释的通你从共线性点图中看到的所有核型变化

比如我发现共线性最好的是A物种的1号染色体(当然这里是fonta)
那么就在祖先核型文件中写下:
A1 1 1316 red 1
第一列是A物种的1号染色体, 第二列第三列是基因编号/顺序
第四列是颜色(祖先核型的Anc1画成红色), 第五列我没看懂
表示A物种的1号染色体即A1的第一个基因到第1316个基因属于祖先核型的1号染色体
祖先核型有几条染色体, ancestor file就有多少种颜色

软件开发者的例子

这里我做了很多共线性以及其他的分析确定了祖先核型的数量, 文件的写法和过程略去
(其实个人觉得wgdi只是用来做可视化的工具； 祖先核型的推断要求的输入文件是根据我们自己看共线性点图写出来的, 换句话我自己看共线性图的时候就已经推断出祖先核型了)

作者示例文件给的是用一个物种, 但实操之后这里用多个物种更为合理
我的文件第四列有30种颜色表示我的物种的祖先核型总共有30条染色体
每一个祖先染色体在现存的研究物种中共线性最好的染色体即第一列, 这里共线性最好是通过和多个外群做共线性分析确定的。
第一列是不同物种的不同染色体, 第二列第三列是这条染色体有多少基因, 第四列是祖先核型的颜色, 第五列不管

这里可以多提一点其他写法:
如果我推断出来一条祖先染色体在C物种断裂成了C1和C2两条染色体 (C物种有一次specific的断裂), 如果我想用C物种的信息构建祖先核型
写法是:
C1 1 1400 red 1
C2 1 1200 red 1

allspp.ancient.input.txt

把所有物种的 .gff4 和 .length2 和 .rename.pep    cat到一起
得到allspp.gff4、allspp.rename.pep、allspp.length2
再加上上文写出来的allspp.ancient.input.txt
就是所有的输入文件

#4.1 得到初步的祖先核型
vim 1.anc.kary.conf
#输入以下内容
[ancestral_karyotype]
gff = allspp.gff4
pep_file = allspp.rename.pep
ancestor = allspp.ancient.input.txt
mark = allspp.ANC
ancestor_gff =  allspp.ancestor.gff     #输出文件
ancestor_lens =  allspp.ancestor.lens   #输出文件
ancestor_pep =  allspp.ancestor.pep    #输出文件
ancestor_file =  allspp.ancestor.txt     #输出文件


wgdi -ak 1.4spalax.anc.kary.conf

#这个时候每个研究物种都要重复相同的流程
#我这里举一个物种的例子
#第一步是blastp比对
makeblastdb -in allspp.ancestor.pep -dbtype prot -out  allspp.ancestor
blastp  -query sppA.rename.pep -out sppA.ANC.blastp -db allspp.ancestor  -evalue 1e-5 -num_threads 60 -outfmt 6

#4.2  collinearity
vim 2.gol.ANC.collinearity.conf
#输入以下内容
[collinearity]
gff1 = sppA.gff4      #最开始准备的输入文件
gff2 = allspp.ancestor.gff  #4.1的输出文件
lens1 = sppA.length2      #最开始准备的输入文件
lens2 = allspp.ancestor.lens    #4.1的输出文件
blast = sppA.ANC.blastp    #刚刚blast的比对结果
blast_reverse = false
multiple  = 1
process = 8
evalue = 1e-5
score = 100
grading = 50,40,25
mg = 40,40
pvalue = 0.2
repeat_number = 10
positon = order
savefile = sppA.ANC.collinearity   #输出文件


wgdi -icl 2.jud.ANC.collinearity.conf


#4.3 Blockinfo
vim 3.sppA.ANC.blockinfo.conf
#输入以下内容
[blockinfo]
gff1 = sppA.gff4      #最开始准备的输入文件
gff2 = allspp.ancestor.gff  #4.1的输出文件
lens1 = sppA.length2      #最开始准备的输入文件
lens2 = allspp.ancestor.lens    #4.1的输出文件
blast = sppA.ANC.blastp    #刚刚blast的比对结果
collinearity = sppA.ANC.collinearity   #4.2的输出结果
score = 100
evalue = 1e-5
repeat_number = 10
position = order
ks = none
ks_col = ks_NG86
savefile =  sppA.ANC.block.csv   #输出文件


wgdi -bi 3.sppA.ANC.blockinfo.conf


#4.4 Correspondence
vim 4.car.ANC.correspondence.conf
#输入以下内容
[correspondence]
blockinfo =  sppA.ANC.block.csv   #4.3的输出文件
lens1 = sppA.length2      #最开始准备的输入文件
lens2 = allspp.ancestor.lens    #4.1的输出文件
tandem = (true/false)
tandem_length = 200
pvalue = 0.2
block_length = 5
multiple  = 1
homo = 0,1
savefile = sppA.ANC.block.correspondence.csv    #输出文件


wgdi -c 4.jud.ANC.correspondence.conf


#4.5存疑, 询问得知这一步不是必需的, 如果只关注核型的演化, 这一步可以略过
#4.5 ancestral_karyotype_repertoire
vim 5.sppA.ANC.akr.conf
#输入以下内容
[ancestral_karyotype_repertoire]
blockinfo =  sppA.ANC.block.correspondence.csv  #4.4的输出文件
blockinfo_reverse = False
gff1 = sppA.gff4      #最开始准备的输入文件
gff2 = allspp.ancestor.gff  #4.1的输出文件
gap = 5
mark = ANC.on.car
ancestor = allspp.ancestor.txt     #4.1的输出文件
ancestor_new =  allspp.ancestor.on.sppA.txt   #sppA的核型结果
ancestor_pep =  allspp.ancestor.on.sppA.pep    #sppA的输出文件
ancestor_pep_new =  allspp.ancestor.on.sppA.pep_new    #sppA的输出文件
ancestor_gff =  allspp.ancestor.on.sppA.gff      #sppA的输出文件
ancestor_lens =  allspp.ancestor.on.sppA.lens     #sppA的输出文件


wgdi -akr 5.sppA.ANC.akr.conf


#4.6 因为是动物, 所以不需要运行Polyploid classification

#4.7 karyotype_mapping
vim 7.sppA.ANC.KaryMap.conf
#输入以下内容
[karyotype_mapping]
blast = sppA.ANC.blastp    #刚刚blast的比对结果
blast_reverse = false
gff1 = sppA.gff4      #最开始准备的输入文件
gff2 = allspp.ancestor.gff  #4.1的输出文件
score = 100
evalue = 1e-5
repeat_number = 5
ancestor_left = ancestor location file (Only one of ('left', 'top') can be reserved)
ancestor_top = allspp.ancestor.on.sppA.txt  #4.6的输出文件
the_other_lens = sppA.length2   #最开始准备的输入文件
blockinfo = sppA.ANC.block.correspondence.csv  #4.4的输出文件
blockinfo_reverse = false
limit_length = 5
the_other_ancestor_file =  ANC.on.sppA.km.txt   #输出文件


wgdi -km 7.sppA.ANC.KaryMap.conf


#4.8 生成图
vim 9.sppA.ANC.karyplot.conf
#输入以下内容
[karyotype]
ancestor = ANC.on.sppA.km.txt  #4.7的输出文件
width = 0.5
figsize = 10,6.18
savefig = sppA.Anc.Kary.pdf

如图所示