GlimmerHMM安装与使用-生信工具24

GlimmerHMM

01 概述

GlimmerHMM是一种基于广义隐马尔科夫模型（GHMM）的新型基因预测工具。虽然该基因预测工具符合GHMM的总体数学框架，但它还结合了从GeneSplicer程序中改编的剪接位点模型。可变长度的特征状态（例如外显子、内含子、基因间区域）是使用Nth-order插值马尔科夫模型（IMM）实现的，如Delcher等人1999年所描述的，N=8。目前，GlimmerHMM的GHMM结构包括每个相位的内含子、基因间区域和四种类型的外显子（起始外显子、内部外显子、终止外显子和单独外显子），如下图所示。

GlimmerHMM在预测DNA序列中的基因时使用了一些假设。主要假设如下：

1. 每个基因的编码区以起始密码子ATG开始（但可以预测部分基因），
2. 除了最后一个密码子，基因中没有在读框内的终止密码子，
3. 每个外显子与前一个外显子保持一致的读码框。

这些约束显著提高了计算最优基因模型的效率，因为它们限制了GHMM算法的搜索空间。另一方面，系统无法检测到真正的移码。

02 下载与安装

wget -c https://ccb.jhu.edu/software/glimmerhmm/dl/GlimmerHMM-3.0.4.tar.gz
tar -zxvf GlimmerHMM-3.0.4.tar.gzhttps://ccb.jhu.edu/software/glimmerhmm/man.shtml   #官网及操作手册

03 训练特定生物的GlimmerHMM

3.1 需要特定生物的训练

我们早期实现的基因预测工具GlimmerM在为水稻训练时，明显比为拟南芥训练时产生的基因模型更优（Pertea和Salzberg 2002）。这对于GlimmerHMM也是有效的。虽然当前的GlimmerHMM训练应该对密切相关的生物体效果良好，但用户应该使用我们的训练程序重新训练系统以适应其他真核生物，以提高基因预测的准确性。对于GlimmerHMM已经训练的所有物种，我们期望随着这些物种DNA序列数量的增加，通过重新训练系统，性能会进一步提高。 

3.2 数据收集

首先，需要注意的是，彻底收集一个良好的训练数据集是训练任何基因预测工具之前的一个步骤。这一步非常重要，因为用于训练的数据质量与最终基因预测工具的准确性成正比。作为任何特定物种的基因预测工具，GlimmerHMM需要一个包含尽可能多的该物种基因组完整编码序列的训练数据集。通过调查公共数据库，可以找到所有有实验室证据支持的目标生物体的已知基因。如果这些基因数量足够多，它们将构成一个良好的训练数据集。不幸的是，这种情况很少见，所以需要使用其他方法来构建可靠的数据集。根据我们的经验，训练数据集通常是通过使用BLAST等程序将长ORF（超过500个碱基）与非冗余蛋白质序列数据库进行比对，以便将已知蛋白质映射到基因组中来构建的。 

3.3 训练GlimmerHMM

要训练GlimmerM，应运行以下命令： 

trainGlimmerHMM <mfasta_file> <exon_file> [optional_parameters]`<mfasta_file>`和`<exon_file>`分别是包含已知基因外显子坐标的多FASTA文件和文件。`<mfasta_file>`是一个包含训练序列的多FASTA文件，格式如下：>seq1
AGTCGTCGCTAGCTAGCTAGCATCGAGTCTTTTCGATCGAGGACTAGACTT
CTAGCTAGCTAGCATAGCATACGAGCATATCGGTCATGAGACTGATTGGGC
>seq2
TTTAGCTAGCTAGCATAGCATACGAGCATATCGGTAGACTGATTGGGTTTA
TGCGTTA

`<exon_file>`是一个包含外显子坐标的文件，坐标相对于`<mfasta_file>`中包含的序列；不同的基因用空行分隔；格式如下：

seq1 5 15
seq1 20 34seq1 50 48
seq1 45 36seq2 17 20

在此示例中，`seq1`有两个基因：一个在正链上，另一个在互补链上。这里可以找到fasta和外显子文件的真实示例。

如果没有足够的数据用于训练剪接位点，训练过程将失败并退出，并显示警告信息。在这种情况下，用户应该收集更多具有内含子的已知基因，然后再次尝试训练过程。如果没有足够的数据用于训练基因的翻译起始和终止位点，用户需要收集更多基因。默认情况下，假定至少有50个带有标准起始密码子（ATG）和终止密码子（TAA/TAG/TGA）的基因。

训练程序可以接受的可选参数包括：

-i i1,i2,...,in	isochores to be considered for training (e.g. if two isochores are desired with 0-40% GC content and 40-100% then the option should be: -i 0,40,100; default is set to -i 0,100 )
-f val	val = average value of upstream UTR region if known
-l val	val = average value of downstream UTR region if known
-n val	val = average value of intergenic region if known
-v val	val = value of flanking region to be considered around genes (default=200)
-b val	val = build 1st or 2nd order markov model for the splice sites (default=1)

- `-i i1,i2,...,in`：用于训练的等温区（例如，如果希望两个等温区的GC含量为0-40%和40-100%，则选项应为：`-i 0,40,100`; 默认设置为`-i 0,100`）- `-f val`：`val` = 如果已知，则为上游UTR区域的平均值- `-l val`：`val` = 如果已知，则为下游UTR区域的平均值- `-n val`：`val` = 如果已知，则为基因间区的平均值- `-v val`：`val` = 考虑在基因周围的侧翼区域值（默认=200）- `-b val`：`val` = 构建1阶或2阶马尔科夫模型用于剪接位点（默认=1）

3.4 修改训练期间获得的参数

这一步是可选的，但许多情况下，手动调整系统参数可以提高预测的准确性。这部分是因为某些生物体的系统训练所需参数可能未知（例如基因间区的平均大小），所以可能需要尝试不同的参数值。我们注意到，类似于Majoros和Salzberg 2004中描述的梯度上升程序，在参数空间中搜索，通常可以优化训练的准确性。

选择适当的阈值来确定特定位点检测的假阴性-假阳性比率可能是一个挑战性问题，因为真实位点数量较少，而假阳性数量却非常多。这需要一系列决策，以最大限度地提高识别任务的准确性，而不显著损失灵敏度。在GlimmerHMM的训练过程中，通过检查假阳性率的权衡来选择调用序列为真实剪接位点的阈值。系统创建一个排序的阈值列表，调整评分函数，以便错过1, 2, 3个等真实位点。默认阈值选择为假阳性率下降不到1%的得分。为了在设置信号阈值时提供更大的灵活性，我们的训练程序允许用户查阅特定于训练数据的假阳性和假阴性率，并设置自己的阈值。

运行trainGlimmerHMM后，可以查看日志文件以找到GlimmerHMM的一些参数的默认值。训练目录中的config.file文件为每个等温区指定使用的配置文件。对于默认情况，当没有考虑等温区时，只有一行信息：

C+G <= 100 train_0_100.cfg

所有参数的更改都可以在train_0_100.cfg文件中完成。每个参数在.cfg文件的一行中描述。以下表格中描述了可以更改的参数：

line in the .cfg file	Description
acceptor_threshold value	acceptor site threshold value; the false negative/false positive rates for different acceptor thresholds can be consulted from the false.acc file.
donor_threshold value	donor site threshold value; the false negative/false positive rates for different donor thresholds can be consulted from the false.don file.
ATG_threshold value	start codon threshold value; the false negative/false positive rates for different start codon thresholds can be consulted from the false.atg file.
Stop_threshold value	stop codon threshold value; the false negative/false positive rates for different stop codon thresholds can be consulted from the false.stop file.
split_penalty value	a factor penalizing the gene finder's tendency to split genes.
intergenic_val value	the minimum intergenic distance espected between genes.
intergenic_penalty value	a penalty factor for genes situated closer than the intergenic_val value.
MeanIntergen value	the average size of intergenic regions.
BoostExon value	a factor to increase the sensitivity of exon prediction.
BoostSplice value	a factor to increase the score of good scoring splice sites.
BoostSgl value	a factor to increase the predicted number of single exon genes.
onlytga value	set this value to 1 if TGA is the only stop codon in the genome.
onlytaa value	set this value to 1 if TAA is the only stop codons in the genome.
onlytag value	set this value to 1 if TAG is the only stop codons in the genome.

- `acceptor_threshold value`：- 受体位点阈值；不同受体阈值的假阴性/假阳性率可以在`false.acc`文件中查阅。- `donor_threshold value`：- 供体位点阈值；不同供体阈值的假阴性/假阳性率可以在`false.don`文件中查阅。- `ATG_threshold value`：- 起始密码子阈值；不同起始密码子阈值的假阴性/假阳性率可以在`false.atg`文件中查阅。- `Stop_threshold value`：- 终止密码子阈值；不同终止密码子阈值的假阴性/假阳性率可以在`false.stop`文件中查阅。- `split_penalty value`：- 惩罚基因预测工具倾向于拆分基因的一个因子。- `intergenic_val value`：- 基因之间预期的最小基因间距。- `intergenic_penalty value`：- 基因间距小于`intergenic_val`值时的惩罚因子。- `MeanIntergen value`：- 基因间区的平均大小。- `BoostExon value`：- 增加外显子预测灵敏度的因子。- `BoostSplice value`：- 增加高评分剪接位点分数的因子。- `BoostSgl value`：- 增加预测单外显子基因数量的因子。- `onlytga value`：- 如果TGA是基因组中唯一的终止密码子，则将此值设置为1。- `onlytaa value`：- 如果TAA是基因组中唯一的终止密码子，则将此值设置为1。- `onlytag value`：- 如果TAG是基因组中唯一的终止密码子，则将此值设置为1。

04 运行 GlimmerHMM

程序GlimmerHMM需要两个输入：一个FASTA格式的DNA序列文件和一个包含程序训练文件的目录。如果未指定训练目录，默认情况下使用当前工作目录。目前，该程序没有实现任何选项。

05 参考文献

Delcher, A.L., Harmon, D., Kasif, S., White,O. and Salzberg, S.L. Improved microbial gene identification with GLIMMER Nucleic Acids Research, 27:23 (1999), 4636-4641.

Majoros, W.H., Pertea, M., and Salzberg, S.L. TigrScan and GlimmerHMM: two open-source ab initio eukaryotic gene-finders Bioinformatics 20 2878-2879.

Majoros,W.M. and Salzberg, S.L. (2004) An empirical analysis of training protocols for probabilistic gene finders. BMC Bioinformatics 5:206.

Pertea, M., X. Lin, et al. (2001). "GeneSplicer: a new computational method for splice site prediction." Nucleic Acids Res 29(5): 1185-90.

Pertea, M. and S. L. Salzberg (2002). "Computational gene finding in plants." Plant Molecular Biology 48(1-2): 39-48.

Pertea, M., S. L. Salzberg, et al. (2000). "Finding genes in Plasmodium falciparum." Nature 404(6773): 34; discussion 34-5.

Salzberg, S. L., M. Pertea, et al. (1999). "Interpolated Markov models for eukaryotic gene finding." Genomics 59(1): 24-31.

Yuan, Q., J. Quackenbush, et al. (2001). "Rice bioinformatics. analysis of rice sequence data and leveraging the data to other plant species." Plant Physiol 125(3): 1166-74.