写本篇文章，主要目的是从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信息如下：
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模型，对序列分化时间进行估计。
插入时间计算公式：，
- K，LTR之间的分化时间
（2）构建系统发育树：使用MAFFT进行多序列比对 —— IQ-tree构建系统发育树 —— ggtree可视化
网友评论