Danaus plexippus genome annotation

Funannotate v1.5.3 docker image was used to train Augustus v3.2.3, predict gene models, and perform functional annotation. As input for optimizing the performance of Augustus v3.2.3, funannotate used 2,404 PASA v2.3.3 gene models. To obtain this training gene model set, transcripts were de novo assembled with Trinity v2018-2.8.3 under settings --SS_lib_type RF, using all poly(A) RNA-seq paired reads after adapter removal with Trimmomatic v0.32. These transcripts were aligned to the genome under PASA using BLAT v36, obtaining a first set of gene models. The 500 longest non-redundant ORFs associated with the PASA gene models were used to train TransDecoder v5.2.0. Then the gene models were selected according to their abundance as estimated by Kallisto v0.44.0 under settings --rf-stranded using the Trinity normalized reads. Ultimately, BRAKER v2.0.3b trained Augustus with the retained gene models.<br> For gene prediction, funannotate aligned mRNAs and proteins from the previous annotation (official gene set 2, OGS2) with minimap v2.14-r883 under settings -ax splice --cs -u b -G 3000, and Diamond blastx v0.8.22, respectively. Protein alignments were further refined by funannotate, including 3 kb upstream and downstream of the region of alignment, and subsequently executing Exonerate v2.4.0. Additionally, funannotate parsed the introns supported by alignments of poly(A) RNA-seq reads generated with HISAT v2.1 under settings --rna-strandness RF --max-intronlen 10,000. This combination of hints (protein alignments, transcript alignments, and intron locations) was used by Augustus to predict a second set of 16,756 gene models. Of them, 9,695 were dubbed as highly supported, i.e. had more than 90% of their model supported either by intron hints, transcript alignments, or protein alignments. GeneMark-ET v4.35, under settings --max_intron 3,000 --soft_mask 2,000, was also run independently to predict a third set of gene models but only relying on intron hints.<br> The PASA, Augustus highly supported, Augustus not highly supported, and GeneMark prediction sets were combined by EVidenceModeler, assigning them 10, 5, 1, and 1 relative weights, respectively. The predictions were further filtered by removing genes shorter than 50 aa in length, or that had high sequence similarity (diamond blastp --sensitive --evalue 1e-10) to the repeat database included in funannotate, or that had more than 90% of the model intersecting regions masked by RepeatMasker. The filtered set of gene models was updated in order to include UTR information by two executions of the PASA annotation comparison using the Trinity transcripts and filtering gene models according to transcripts per million as calculated by Kallisto. Alternative transcripts were only kept if they were at least 10% as highly expressed as the most highly expressed transcript per gene.<br> Non-coding genes were annotated with the following tools: tRNA genes, tRNAscan-SE v.2.0; rRNA genes, RNAmmer v.1.2; and for a variety of other RNA genes, Infernal v1.1.1. Specifically, for miRNA-encoding genes, we used BLASTn to locate the most recent annotation of these genes. Lastly, FEELnc classified lncRNAs from the transcripts assembled by StringTie v1.3.2d , and considering protein-coding predictions described above.

本研究使用Funannotate v1.5.3容器镜像完成了Augustus v3.2.3的模型训练、基因结构预测与功能注释工作。为优化Augustus v3.2.3的预测性能，Funannotate以2404个PASA v2.3.3基因模型作为训练输入集。为获取该训练基因模型集，研究人员首先使用Trimmomatic v0.32去除测序接头，随后利用所有经处理的poly(A) RNA-seq双端测序读段，结合参数--SS_lib_type RF，通过Trinity v2018-2.8.3进行从头转录组组装。将组装得到的转录本通过BLAT v36比对至参考基因组，并借助PASA流程得到首批基因模型集。选取与PASA基因模型相关联的500条最长非冗余开放阅读框（Open Reading Frame, ORF），用于训练TransDecoder v5.2.0。随后基于Kallisto v0.44.0（参数--rf-stranded）对Trinity归一化后的测序读段进行定量，根据基因模型的表达丰度进行筛选，最终保留的基因模型被用于BRAKER v2.0.3b训练Augustus。
在基因预测环节，Funannotate分别使用minimap v2.14-r883（参数-ax splice --cs -u b -G 3000）与Diamond blastx v0.8.22，将来自既往注释版本（官方基因集2，OGS2）的mRNA与蛋白质序列进行比对。Funannotate还对蛋白质比对结果进行了精细化处理：将比对区域上下游各3kb纳入考量范围，并执行Exonerate v2.4.0进行进一步校验。此外，研究人员利用HISAT v2.1（参数--rna-strandness RF --max-intronlen 10000）比对poly(A) RNA-seq读段，获取支持的内含子信息作为注释线索。上述三类注释线索（蛋白质比对、转录本比对与内含子位置信息）被一并输入Augustus，得到第二组共16756个基因模型。其中9695个被定义为高支持度模型，即其模型结构中超过90%的区域可被内含子线索、转录本比对或蛋白质比对所验证。此外，研究人员还独立运行了GeneMark-ET v4.35（参数--max_intron 3000 --soft_mask 2000），仅基于内含子线索得到第三组基因模型。
研究人员通过证据模型整合器（EVidenceModeler）整合PASA基因模型集、高支持度Augustus模型集、低支持度Augustus模型集与GeneMark预测模型集，并分别为其赋予10、5、1、1的相对权重。随后对整合后的预测结果进行过滤：移除长度小于50个氨基酸的基因，移除与Funannotate内置重复数据库存在高度序列相似性（通过diamond blastp --sensitive --evalue 1e-10评估）的基因，以及移除超过90%模型区域被RepeatMasker屏蔽的基因。为向过滤后的基因模型集添加非翻译区（Untranslated Region, UTR）信息，研究人员基于Trinity组装的转录本，通过两次PASA注释比对流程完成更新，并依据Kallisto计算的每百万转录本数（Transcripts Per Million, TPM）对基因模型进行筛选。仅保留那些表达丰度不低于该基因最高表达转录本10%的可变剪接转录本。
非编码基因的注释通过以下工具完成：转运RNA（tRNA）基因使用tRNAscan-SE v.2.0；核糖体RNA（ribosomal RNA, rRNA）基因使用RNAmmer v1.2；其余各类RNA基因则使用Infernal v1.1.1进行注释。其中，针对microRNA（miRNA）编码基因，研究人员通过BLASTn定位该类基因的最新注释结果。最终，FEELnc基于StringTie v1.3.2d组装得到的转录本，并结合前述蛋白质编码基因预测结果，对长链非编码RNA（long non-coding RNA, lncRNA）进行分类注释。