蛋白编码基因的注释
输入文件的准备
注释要注意基因组数据完整性和连续性，评估组装质量，这决定着注释结果的质量和准确性
需要以蛋白质序列、表达序列标签 (EST)、全长转录本序列和 RNA-Seq 数据的形式等外部证据。这些与屏蔽重复序列的基因组比对，以帮助确定正确的基因结构
方法学
计算阶段
重复序列的识别
重复序列的识别可以通过RepeatModeler完成。最新版本，即 RepeatModeler2，由多个不同的重复序列发现程序组成，即 RepeatScout、RECON、LTRharvest 和 LTR_retriever 。 RepeatScout 和 RECON 分别用于调查高度保守的 TE 和不太保守的 TE 的基因组。 LTRHarvest 用于搜索长末端重复 (LTR) 转座子，LTR_retriever 用于过滤来自 LTRHarvest 的假阳性结果。
第一步：使用组装的基因组构建数据库，如下所示：
~/RepeatModeler-2.0/BuildDatabase -name genome genome.fa
BuildDatabase命令采用以下参数：我们要创建的数据库的名称-name，后跟输入文件genome.fa
Repeat Modeler2并没有用所有的scaffold/contig建数据库，它会进行采样，因此每次运行会有不同的结果。在这种情况下如果过滤掉短contig（如长度小于1000bp的contig），就可以避免这些短contig被采样。因此，建议在建库之前过滤掉短contig。
创建数据库后可以按如下命令运行 RepeatModeler2
~/RepeatModeler-2.0/RepeatModeler -database genome -LTRStruct
选项-database表示要分析序列的数据库名，即创建数据库BuildDatabase提供的名称。-LTRStruct让RepeatModeler2运行LTRHarvest 和 LTR_retreiver，并将结果与 RepeatScout 和 RECON 输出相结合。
RepeatModeler2成功运行后会产生两个输出文件：genome-families.fa和genome-families.stk。这两个文件包含相同的信息，即在基因组中从头识别的共有重复序列。但是，genome-families.fa是 FASTA 格式，而genome-families.stk是可以提交到 Dfam 数据库的Stockholm-formatted文件
重复序列的屏蔽
识别出重复序列后，就可以在基因组数据上屏蔽他们。RepeatModeler2识别出的重复序列可以用RepBase数据库补充。RepBase是最常用的DNA重复元件数据库。以下代码可以完成
perl ~/RepeatMasker/util/buildRMLibFromEMBL.pl RMRBSeqs.embl > RMRBSeqs.embl.fasta
cat genome-families.fa RMRBSeqs.embl.fasta > genome-families. RepBase.fa
然后就可以用RepeatMasker来屏蔽重复序列。以下是运行RepeatMasker的命令
~/RepeatMasker -lib genome-families.RepBase.fa -pa 16 -xsmall genome.fa
RepeatMasker命令输入以下文件：屏蔽重复序列所需的通过 RepeatModeler2 和 RepBase输出的重复序列组合库 (lib) 以及需要被执行屏蔽重复序列的组装基因组(genome.fa)。此外，-xsmall选项指示工具执行softmask，-pa指定要使用的内核数。完成后，RepeatMasker 输出三个文件：一个屏蔽重复序列的基因组组装文件 (genome.masked.fa)， 一个包含屏蔽序列列表 的文件(genome.out)，以及一个总结序列重复内容的文件 (genome.tbl)。
证据排列
Exonerate可以组装蛋白序列和表达序列标签（EST)。如以下命令所示
~/exonerate --softmasktarget --model protein2genome --showvulgar no --showalignment no --showquerygff no --showtargetgff yes --percent 80 --ryo "AveragePercentIdentity: %pi\n" --query protein.fasta --target genome.masked.fa > genome.protein.exonerate
~/exonerate --softmasktarget --model est2genome --showvulgar no --showalignment no --showquerygff no --showtargetgff yes -percent 80 --ryo "AveragePercentIdentity: %pi\n" --query est. fasta --target genome.masked.fa > genome.est.exonerate
--softmasktarget选项允许对目标序列 (genome.masked.fa) 进行softmask。 --model选项的选择取决于所使用的查询序列的类型：protein2genome发出信号 Exonerate 以查找与目标匹配率至少 50% 的蛋白质，而est2genome用于查找前10个匹配项到目标的EST。根据所需的输出格式设置选项--showvulgar、--showalignment、--showquerygff、--showtargetgff和--ryo。 --percent选项指的是百分比自我评分阈值。
至于RNA-Seq reads，可以通过三种方式处理：（1）它们可以与感兴趣的基因组比对，比对后的reads作为证据提供； (2) aligned reads可以组装成转录本作为证据； (3) raw reads可以从头组装，从头组装作为证据。第一种方法可以使用 Hisat2 执行，这是一种流行的拼接感知对齐器。下面显示了如何使用 Hisat2 v. 2.1.0 的示例：
hisat2-build -p 16 genome.masked.fa genome.masked.fa
hisat2 -x genome.masked.fa -p 16 -1 sampleA.1.fq -2 sampleA.2. fq -S sampleA.sam
hisat2-build命令用于为屏蔽重复序列的基因组建立索引；它采用以下参数：输入文件 genome.masked.fa、要创建的索引文件的名称（在本例中为 genome.masked.fa）和要使用的线程数（本例中为-p 16 ）。然后可以使用hisat2命令执行读取比对，该命令采用以下参数：索引文件的基本名称 -x、读取对的第一个配对 -1和读取对的第二个配对 -2。 -S选项指示 Hisat2 以 SAM 格式输出读取映射，然后可以使用 SAMtools将其转换为排序的索引BAM文件。下面显示了运行 SAMtools v. 1.9 的示例命令
samtools view -Sb sampleA.sam | samtools sort - > sampleA. sorted.bam
samtools index sampleA.sorted.bam
已经表明，第一种方法的基因预测在统计学上比涉及组装 转录本组装的方法（例如第二种和第三种方法）更稳健。此外，第二种方法可能会导致在转录本组装步骤中从相邻的、间隔很近的基因中读取 RNA-Seq 进行合并，但denovo组装可以避免这个问题。
组装转录本与屏蔽重复序列的基因组的比对仍然可以用作注释阶段的证据来源；这可以通过程序Program to Assemble Spliced Alignments (PASA)来实现。
首先，通过“seqclean”命令对转录本序列（transcripts.fasta）进行质量修整（去除 Ns、低质量碱基和 polyA 尾），如下所示：
~/PASApipeline/bin/seqclean transcripts.fasta
上面产生了两个文件，其中任何一个都可以用于下一步；这些文件包括修剪后的转录序列 transcripts.fasta.clean)和一个制表符分隔的文件，其中包含执行的修剪类型和序列中发生修剪的位置 transcripts.fasta.cln。然后可以按如下方式运行PASA v.2.4.1 比对组装流程
~/PASApipeline/Launch_PASA_pipeline.pl \
-c alignAssembly.config \
-C -R -g genome.masked.fasta \
-t all_transcripts.fasta.clean \
-T -u transcripts.fasta \
--ALIGNERS gmap
这里，gmap 用于转录比对。 -c参数用于指定包含存储输出的 SQLite 数据库或 MySQL 数据库路径的配置文件。 “-C”和“-R”选项分别用于创建数据库和运行对齐/组装流程。屏蔽重复序列程序集、干净的转录本序列和原始转录本序列分别通过“-g”、“-t”和“-u”参数作为输入文件提供。 “-T”参数向 PASA 发出信号，表明转录本序列已使用“seqclean”命令进行了修剪。完成时会生成几个文件；感兴趣的是 FASTA 格式 (PASA.assemblies.fasta) 和一般特征格式 3 (GFF3) 格式 (PASA.assemblies.gff3) 的 PASA 程序集，它们是注释阶段所需的。 GFF3 文件是一个制表符分隔的文件，如表 2中所述。
基因预测
基因预测可以是denovo，也可以基于证据。证据驱动的基因预测是基因预测的首选和更准确的方法。这种方法需要将外部证据来源纳入基因预测 。可用于此目的的流行工具是 BRAKER2。 BRAKER2通过RNA-Seq比对利用两个单独工具（即 GeneMark-EX 和 AUGUSTUS）的denovo基因预测功能进行训练和基因结构推断。 GeneMark-EX 是一种强大的自训练算法，可以学习物种特异性参数并识别真核基因组中的蛋白质编码基因；如果可用，它会使用外部证据来完善其模型。 GeneMark-EX 识别的初始预测基因集随后被 AUGUSTUS 用作训练集，它通过整合外部证据产生最终的基因预测。当具有良好转录组覆盖率的 RNA-Seq 文库可用于感兴趣的物种时，可以使用以下 BRAKER2 命令：
perl ~/BRAKER/scripts/braker.pl \ --genome=genome.masked.fa \ --AUGUSTUS_BIN_PATH=~/Augustus/bin \ --AUGUSTUS_CONFIG_PATH=~/Augustus/config \ --GENEMARK_PATH=~/gm_et_linux_64 \ --BAMTOOLS_PATH=~/bamtools/bin \
--species=species1 \ --workingdir=Braker_output \ --bam=all.samples.bam \ --softmasking
“braker.pl”脚本采用以下参数：softmask程序集--genome、BAM 格式的单独或合并的 Hisat2 映射--bam、AUGUSTUS 路径-AUGUSTUS_BIN_PATH、GeneMark-EX- -GENEMARK_PATH)和其他附属程序，如 Bamtools --BAMTOOLS_PATH，运行 AUGUSTUS 所需的配置文件 --AUGUSTUS_CONFIG_PATH)和物种名称 -species；这是 AUGUSTUS 所必需的。此外，--softmasking选项向 BRAKER2 发出信号，表明输入程序集已被softmask。 BRAKER2
完成后输出三个文件：一个包含从 RNA-Seq 映射中提取的内含子信息并分别提供给 GeneMark-EX 和 AUGUSTUS 用于基因训练和预测的文件 (hintsfile.gff)，一个包含 GeneMark-EX 预测的基因的文件（genemark.gtf ）和一个包含 AUGUSTUS 预测的基因的文件 (augustus.gtf)。
注释阶段
一旦基因预测可用，下一步就是将它们与对齐的外部证据结合起来（见注释 9）。可以使用 EvidenceModeler (EVM)执行此任务。 EVM 将基因预测与外部证据的比对结合到加权共识基因结构中。 EVM v.1.1.1 需要以下输入文件：FASTA 格式的重复掩蔽程序集、GFF3 格式的证据（比对的蛋白质和转录序列和基因预测）以及证据权重文件。至于权重文件，它包括三列，即所使用证据的类别、类型和权重； “class”表示所用证据的类型，“type”包含所用证据的来源，“weight”是应用于每种证据类型的数值。图 2 显示了基于我们工作流计算阶段产生的证据的权重文件示例。 “ABINITIO_PREDICTION”指来自 BRAKER2 的 AUGUSTUS 和 GeneMark-EX 输出文件，“PROTEIN”指来自 Exonerate 的合并蛋白质和 EST 比对，“TRANSCRIPT”指 PASA 组件。在这里，我们使用 EVM 手册中推荐的权重值。典型 EVM 运行的第一步包括将输入文件划分为更小的数据集，如下所示：
ABINITIO_PREDICTION AUGUSTUS 1
ABINITIO_PREDICTION GeneMark_hmm 1
PROTEIN exonerate 1
TRANSCRIPT PASA 10
~/EvidenceModeler/EvmUtils/partition_EVM_inputs.pl \
--genome genome.masked.fasta \
--gene_predictions augustus.genemark.gff3 \
--protein_alignments exonerate.gff3 \
--transcript_alignments pasa.gff3 \
--segmentSize 100000 \
--overlapSize 10000 \
--partition_listing partitions_list.out
上述命令实质上是根据相应的 contigs 拆分 FASTA 和 GFF3 文件，较大的 contigs 进一步拆分为较小的重叠片段。 --genome选项指定 FASTA 格式的重复掩蔽装配。 --gene_predictions、--protein_alignments和--transcript_alignments文件分别指代 BRAKER2、Exonerate 和 PASA 输出文件。 --segmentSize设置为小于 1 Mb 以减少内存需求，--overlapSize设置为大于预期基因长度的长度以防止丢失完整的基因结构。 --partition_listing参数允许 EVM 总结在此步骤中创建的分区列表。第二步包括创建一组 EVM 命令，这些命令可以在第一步中的各个分区上运行，如下所示：
~/EvidenceModeler/EvmUtils/write_EVM_commands.pl \
--genome genome.masked.fasta \ --weights weights.txt \
--gene_predictions augustus.genemark.gff3 \
--protein_alignments exonerate.gff3 \
--transcript_alignments pasa.gff3 \
--output_file_name evm.out \
--partition partitions_list.out > commands.list
--output_file_name参数指定输出文件名，commands.list文件包含可以运行的命令，如下所示：
~/EvidenceModeler/EvmUtils/execute_EVM_commands.pl commands. list | tee run.log
每个命令的退出值保存在“run.log”文件中；退出值为零表示 EVM 命令已成功运行。最后，合并分区，解决冗余或不一致的预测。这可以使用以下命令完成：
~/EvidenceModeler/EvmUtils/recombine_EVM_partial_outputs.pl
--partitions partitions_list.out --output_file_name evm.out
原始evm.out文件以制表符分隔，详细说明了完全支持每个外显子结构的证据集。也可以使用以下命令将此文件转换为 GFF3 格式：
~/EvidenceModeler/EvmUtils/convert_EVM_outputs_to_GFF3.pl \
--partitions partitions_list.out \
--output evm.out \
--genome genome.masked.fasta
最后，可以用以下命令将每个分区的GFF3文件合并成一个文件。
find . -regex ".*evm.out.gff3" -exec cat {} \; > EVM.all.gff3
图3显示了所产生的GFF3文件的摘录。可以提取FASTA格式的CDS和注释基因的蛋白质，如下所示。
~/EVidenceModeler-1.1.1/EvmUtils/gff3_file_to_proteins.pl EVM.all.gff3 genome.fasta CDS > EVM.CDS.fasta ~/EVidenceModeler-1.1.1/EvmUtils/gff3_file_to_proteins.pl EVM.all.gff3 genome.fasta prot > EVM.prot.fasta
质控
一旦基因注释完成，检查注释的质量是很重要的。这个步骤不应该被忽视，因为它可能对后续的分析产生不利的影响。例如，习惯上使用一个物种的基因注释来协助密切相关物种的注释；如果存在不正确的注释，它们可以从一个项目延续到另一个项目。幸运的是，有几种方法可以检查基因注释的质量。例如，可以用BUSCO v.4.0在基因组水平和基因水平上进行两种不同的分析，具体如下。
~/busco-master/src/busco/run_BUSCO.py -i genome.fa -m genome -o genome -l eudicotyledons_odb10 -c 16
~/busco-master/src/busco/run_BUSCO.py -i EVM.prot.fasta -m protein -o EVM.prot -l eudicotyledons_odb10 -c 16
参数-i定义了要分析的输入文件（基因组或蛋白质FASTA），参数-o定义了存储输出文件的文件夹，参数-m指定了运行BUSCO的模式（基因组、蛋白质和转录组）。通过-l参数指定包含种群特定信息的世系文件（eudicotyledons_odb10）。使用的CPU数量是通过-c参数指定的。每次运行的结果包括完整的、重复的、缺失的和片段的基因的数量和百分比。可以对两次运行的结果进行比较，以确定所报告的BUSCOs的数量和百分比的任何差异。理想情况下，这种差异是不存在的或最小的。与基因组水平相比，在基因水平上有更多的缺失和破碎的BUSCO基因，将表明基因注释不完整，需要进一步调查
另一种有用的方法是通过序列同源性对基因进行功能注释[1]。例如，鉴于同源基因在物种之间是保守的，我们期望它们具有高度相似的序列，并在不同的物种中发挥相同的作用。序列同源性可以通过对注释的基因进行BLAST搜索来评估，这些数据库包括NCBI非冗余数据库、SwissProt/UniProtKB和RefSeq。使用BLAST v.2.9.0的例子命令显示在下面。
blastp -db database -query EVM.prot.fasta -outfmt 6 -evalue 1e-3 -out EVM.prot.blastp
blastp命令需要以下参数：一个BLAST格式的数据库-db，一个输入文件-query，如何显示结果-outfmt,在这种情况下选择表格格式），期望值e-value阈值-evalue和输出文件的名称-out。一个典型的BLAST输出文件显示在图4中，表3中对各列进行了解释。
值得注意的是，有些基因在这些数据库中可能没有任何匹配的基因；这种物种特有的基因可能是真正的生物变异的结果，也可能是被错误地注释了。因此，进行进一步的检查是很重要的。
一种方法是挖掘基因的蛋白域，因为这些蛋白域意味着基因的生化或调节功能；缺乏这种信息的基因将被排除在下游分析之外[25]。蛋白质结构域信息可以用InterProScan等工具推断出来[26]。下面用InterProScan v.5.41-78.0举例说明如何做到这一点。
~/interproscan-5.41-78.0/interproscan.sh -appl TIGRFAM, SUPER FAMILY,CDD,PRINTS,Pfam,PANTHER,ProDom,PrositeProfiles,PrositePatterns,Coils -i EVM.prot.fasta -o EVM.prot.interpro.out -f gff3
命令 "interproscan.sh "接受以下参数：以逗号分隔的要运行的应用程序列表-appl，输入文件-i，输出文件-o，以及输出格式-f；本例为GFF3。一个典型的输出文件显示在图5中。另一种方法是检查CDS序列中是否存在起始（"ATG"）和终止（"TAG"、"TAA "和 "TGA"）密码子；如果CDS序列中缺少这些密码子，相应的基因可以被标记为不完整，并排除在下游分析之外。
本文翻译自《Plant Bioinformatics：Methods and Protocols》 David Edwards 图片和表格可前往查看
