【SubPhaser-多倍体亚基因组分型流程解读】

写本篇文章，主要目的是从tmp文件和软件运行信息解读亚基因组分型分析。

进入tmp文件夹之后，其实就可以看到对应的文件：

如何查看什么步骤产生了怎么样的结果文件？

上述在进行SubPhaser试运行的时候，使用了nohup命令，该软件调用了什么软件、产生了什么结果文件等信息，都是记录在最终的nohup.out。

1、参数配置

从截图中，我们可以得到很多的信息

分析所使用的k是多少（默认情况下，k=15）
min_fold是多少
min_freq是多少
lower_count是多少
LTRfinder所使用的参数
LTRharvest所使用的参数
所使用cpu数
所使用的内存大小

软件先将基因组按染色体划分，结果保存于：/opt/biosoft/SubPhaser/example_data/Arabidopsis_suecica_tmp/Arabidopsis_suecica_chromosomes

2、Kmer计数

在此步骤成功生成.histo之后，会生成一个.ok文件

# e.g.
# 10.fasta                  # 10号染色体序列文件
# 10.fasta_15.fa            # 使用k=15切割，并对15mer进行频数统计的结果
# 10.fasta_15.fa.ok         # ok文件
10.fasta_15.histo           # Jellyfish结果文件
10.fasta_15.jf.ok           # ok文件

根据上述结果，构建了一个以Kmer类型为行，染色体ID为列的矩阵，每一个单元代表该类型的Kmer在该条染色体中的占比：

鉴定亚基因组特异性Kmer的几个重要阈值：

该Kmer在全基因组范围内的频数，需要超过200
该Kmer在A homoeolog中的频数要求至少是B homoeolog中的2倍，即判断为A中的特异性Kmer

3、聚类

作者对于聚类的描述如下，

使用k-means聚类方法，将染色体组聚集到某个类中（bootstrap默认情况下为1000，并且每次bootstrap只抽取原数据的50%进行聚类分析）
使用PCA对phasing结果是否成功进行评价（我没有仔细考究，但是我猜的使用方差来衡量）

the k-Means algorithm is used to cluster chromosomes into N groups (subgenomes) and perform 1,000 bootstrap resampling of 50% of these k-mers to proceed with analyses of the sampling distribution and statistical inference.

过程Output信息如下：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-brdMc0dh-1652088221884)(https://upload-images.jianshu.io/upload_images/24361169-7f8cfdadf4a40b6f.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)]

Output信息中，还包含了绘制PCA图所使用的脚本：Arabidopsis_suecica_k15_q200_f2.kmer.mat.R

结果文件：Arabidopsis_suecica_k15_q200_f2.kmer_pca.pdf

图示：

最终得到的亚基因组划分结果为：Arabidopsis_suecica_k15_q200_f2.chrom-subgenome.tsv

>#chrom  subgenome  bootstrap
1  SG1  100
2  SG1  100
3  SG1  100
4  SG1  100
5  SG1  100
6  SG2  100
7  SG2  100
9  SG2  100
8  SG2  100
10  SG2  100
13  SG2  100
11  SG2  100
12  SG2  100

4、统计检验

该部分的统计检验有2个目的：

检验该Kmer是否在对应亚基因组中呈现特异性

使用方法：student’s t-test
检验该Kmer是否在某一区域内呈现富集状态

使用方法：Fisher’s exact test（原理与GO富集分析相同，给忘记的同学们提个醒）

Output关键信息提取：Consistent with subgenome assignment: 248 (91.85%); potential exchange: 10 (3.70%)。

即，经统计检验之后，对应亚基因组某些部分或许与phasing结果不同，可能存在染色体重排等。

5、亚基因组特征分析

Output信息如下：

SubPhaser一些其他功能，即鉴定基因组中的哪些特征是在某一亚基因组中呈富集状态。
e.g. 转座子，gene，内含子，外显子，转录本
同时，SubPhaser还可以使用LTR-RTs（实际上是调包）来估计异源多倍体形成时间。

By default, the software identifies and analyzes subgenome-specific long terminal repeats retrotransposons (LTR-RTs) to calculate their insertion time and hence estimate the time boundaries from ancestor differentiation to allohybridization.

该分析所使用的软件：LTRharvest v1.6.1 & LTRfinder v1.07，同时使用TEsorter v1.3.0降低LTR-RTs检测的假阳性。

这边有一个非常重要的背景知识，即多倍体化之后，会激活LTR-RTs的活动：

Subgenome-specific LTR-RTs are considered to be actively inserted only when diploid progenitors have evolved as independent species (without exchanging LTR-RTs),

最终再将特异性Kmer回帖到鉴定出的LTR-RTs序列，鉴定出特异性的LTR-RTs序列。

关于异源多倍体化时间鉴定

（1）分化时间估计
使用来自不同亚基因组的LTR-RTs序列，使用Jukes-Cantor 69模型，对序列分化时间进行估计。

插入时间计算公式：T=K2rT = \frac{K}{2r}T=2rK，

r=1.3×10−8substitutionsperyearr = 1.3×10^{-8} substitutions per yearr=1.3×10−8substitutionsperyear
K，LTR之间的分化时间

（2）构建系统发育树：使用MAFFT进行多序列比对 —— IQ-tree构建系统发育树 —— ggtree可视化