写这篇文章的主要目的是从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中间的两倍是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|imageView2/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、统计检验
统计检查有两个目的：
-
检验该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. 同时，SubPhaser还可以使用LTR-RTs(实际上是调包)估计异源多倍体的形成时间。
By default, the software identifies and analyzes subgenome-specific long terminal repeats retrotransposons (TR-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 = K 2 r T = \frac{K}{2r} T=2rK，
- r = 1.3 × 1 0 − 8 s u b s t i t u t i o n s p e r y e a r r = 1.3×10^{-8} substitutions per year r=1.3×10−8substitutionsperyear
- K，LTR之间的分化时间
（2）构建系统发育树：使用MAFFT进行多序列比对 —— IQ-tree构建系统发育树 —— ggtree可视化