【如何拆分多倍体亚基因组？】SubPhaser软件

SubPhaser原文链接：https://nph.onlinelibrary.wiley.com/doi/10.1111/nph.18173

简单认识Subphaser

这款软件所针对的研究对象是谁？

neoallopolyploid（近多倍体）
homoploid hybrid（整倍体）—— 这个忘了的话，看看自己写的这篇文章：浅谈多倍体概念以及多倍体形成的几种途径

这个软件的工作流程是怎样的？

先寻找亚基因组特异性Kmer，再根据这些Kmer，将亚基因组分开。

（1）软件安装

git clone https://github.com/zhangrengang/SubPhaser.git
cd SubPhaser

# 创建subphaser虚拟环境 & 安装subphaser
conda env create -f SubPhaser.yaml
conda activate SubPhaser
python setup.py install

安装完，可以看看这个软件有多大，2.5G恐怕是我目前用到的依赖最多的软件了。

（2）测试

# 工作目录：/opt/biosoft/SubPhaser/example_data
# 使用拟南芥作为测试数据
nohup bash test_Arabidopsis.sh &

# 结果文件存储在：/opt/biosoft/SubPhaser/example_data/Arabidopsis_suecica_phase-results
# 转移至本地的存储位置：C:\Hongpu Chen\Data\【知识库】\【生物信息学】\【生物信息学-软件】\【工作目录】如何拆分亚基因组\Arabidopsis_suecica_phase-results

上述运行的脚本，内容如下：

prefix=Arabidopsis_suecica
# download genome
url=https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/019/202/805/GCA_019202805.1_ASM1920280v1/GCA_019202805.1_ASM1920280v1_genomic.fna.gz
[ ! -s ${prefix}_genome.fasta.gz ] && \
    wget $url -O ${prefix}_genomic.fna.gz -c && \
    mv ${prefix}_genomic.fna.gz ${prefix}_genome.fasta.gz

DT=`date +"%y%m%d%H%M"`
# run subphaser
options="-pre ${prefix}_" # to avoid conficts
subphaser -i ${prefix}_genome.fasta.gz -c ${prefix}_sg.config $options 2>&1 | tee ${prefix}.log.$DT

（3）测试数据结果解读

SubPhaser的结果，主要就是如下6个图

A：不同类型Kmer的频数统计图
B：不同类型Kmer的聚类结果（用heatmap & 聚类树来显示）
C：PCA
D：染色体圈图
E：亚基因组LTR-RTs的插入时间估计
F：对LTR/Gypsy成分进行1000次后验检验的系统发育树构建结果

其中，D图给出的信息很多，从外圈到内圈，分别代表：

1.使用k-means方法，得到的染色体分组信息
2.对应区段内，亚基因组特异性kmer的富集结果
3.对应区段内，经过标准化的亚基因组kmer组分估计
4-6.对应区段内，对应亚基因组特异性的kmer频数概率密度（划分为几个亚基因组，这边就对应有几圈）
7.亚基因组特异性LTR-RTs & 非特异性LTR-RTs的概率密度
8.不同区段之间的同源关系分析

D图解读

1.基于k-means聚类算法，将输入基因组，分成了<font color='yellow'>2</font>个亚基因组（黄色 & 蓝色部分）
2.对应区段的特异性kmer富集结果（我这边猜测的话，哪一种kmer显著富集就显示哪一种结果）
3.经过标准化的kmer占比
4-5.对应特异性kmer的概率密度结果
6.对应LTR-RTs的概率密度结果
7.同源关系

（4）从结果解读反推亚基因组分型分析

我所选用的测试物种是：Arabidopsis suecica，其配置文件如下：

# Arabidopsis_suecica_sg.config
1|CM032900.1    6|CM032905.1,7|CM032906.1
2|CM032901.1,3|CM032902.1   9|CM032908.1,8|CM032907.1,10|CM032909.1
4|CM032903.1,5|CM032904.1   13|CM032912.1,11|CM032910.1,12|CM032911.1

编写SubPhaser配置文件过程中，需要注意的几个点：

每一行为一同源染色体组（homoeologous chromosome set）
并不需要按染色体编号排序
使用空格分隔
需要重命名时，使用|
需要将多条染色体当作一条进行看待的时候，使用,进行连接

额外

如何推断scaffold/chromosome之间的同源关系？

之前没咋接触过基因组学方面的分析，因此开始学习SubPhaser这个软件之后，又转头去补充了一些基因组学、homology方面的知识。这几天一看Github的issues，张仁纲老师已经给出了上述问题对应的解决方案。

我这边浅浅做一个总结。

1、minimap2比对结果

基于minimap2对Fragaria vesca 和Fragaria ananassa的比对结果，可以设置config文件为：

1-1 1-2 1-3 1-4
2-1 2-2 2-3 2-4
3-1 3-2 3-3 3-4
4-1 4-2 4-3 4-4
5-1 5-2 5-3 5-4
6-1 6-2 6-3 6-4
7-1 7-2 7-3 7-4

图示：

2.HiC热图

根据HiC热图的结果，可以得到如下的config设置：

chr01a chr01b chr01c chr01d
chr02a chr02b chr02c chr02d
chr03a chr03b chr03c chr03d

参考资料

[1] https://github.com/zhangrengang/SubPhaser
[2] https://github.com/zhangrengang/SubPhaser/issues/4