基因组重复序列注释全流程：RepeatModeler建库与RepeatMasker实战指南

拿到一个刚组装的基因组，第一步不是急着做基因注释，而是先把里面的重复序列全部“揪出来”。这个问题我吃过亏：早期做某个非模式物种基因组时，直接拿Repbase数据库跑了RepeatMasker，结果后来做基因结构预测时，假阳性高得离谱，一大半预测出来的“基因”其实都是转座子碎片。后来老老实实补了RepeatModeler从头建模这一步，再回头看，所有下游分析都稳了一个档次。

RepeatModeler和RepeatMasker这对组合，基本是基因组重复序列注释的事实标准。前者负责从头构建你当前物种的重复序列文库，后者负责用这个文库把基因组里所有重复位点标记出来、屏蔽掉。这套流程适用于刚组装好的基因组、需要做基因注释的物种、以及想系统了解重复序列组成的研究者。这篇文章我会把从环境安装、建模、注释到结果评估的完整过程拆开讲，包括那些文档里不会明说的坑。

1. 为什么注释重复序列必须放在基因组分析的第一站

很多做基因组项目的朋友容易低估重复序列的存在感。简单说，重复序列不是“基因组里的一点杂质”，而是基因组的主体构成。人类基因组中重复序列占比大约在45%-50%，植物里更夸张：玉米基因组超过80%，小麦甚至能到85%。你面对一个基因组时，里面一半以上可能都是转座子、反转座子、串联重复这些“非基因”的东西。

这些重复序列不是你做完基因注释后再“顺便处理”的边角料。它们会直接影响三件事：

第一，基因结构预测的准确性。无论是从头预测软件（Augustus、GlimmerHMM、SNAP），还是基于转录组证据的注释流程（BRAKER、MAKER），在识别基因模型时都极容易把转座子内部的序列误判为外显子。转座子编码区里有大量看似“规整”的开放阅读框，在不屏蔽重复序列的情况下，基因预测结果里会出现大量假基因，后续做功能注释和进化分析都是埋雷。

第二，比对结果的可靠性。RNA-seq reads比对回基因组时，如果重复序列没有屏蔽，reads会大量多重比对到不同位置的同源转座子拷贝上。主流比对工具对multimapping reads要么随机分配要么直接丢弃，最终导致基因表达量估计严重偏差。做变异检测也一样，重复区里的假阳性位点会让你后续验证实验白做一大堆。

第三，基因组结构分析的信号污染。做共线性分析、Hi-C互作图、重组率估算时，重复序列会制造大量虚假的同源信号。特别是对于物种内部或者近缘种之间的比较，转座子家族的扩增潮会在不同物种间留下“相同的印记”，不把它们识别出来，你在共线性区块里看到的可能只是重复序列的同源拷贝，而不是真正的直系同源区段。

1.1 模式物种可以用现成库，非模式物种必须从头建模

RepeatMasker本身支持直接指定物种，比如-species human或者-species maize，它会从自带的Dfam和Repbase库中提取对应物种的重复序列。这种方法在模式物种上很快，效果也不错。

但问题在于，Dfam和Repbase的物种覆盖是极度不均的。人和小鼠这些模式物种注释得很全，但你可能在做一个地方性的鱼类、一个经济昆虫或者一个野生植物，这些物种的重复序列在公共库里往往只有零星几条，甚至完全没有。用公共库去屏蔽一个远缘物种的基因组，会出现严重的漏报——大量物种特有的转座子新家族根本没有参考序列，自然不会被识别。

这就要靠RepeatModeler来做从头预测：它不依赖任何已有注释，直接从你的基因组序列里发现重复单元的家族结构和边界信息，构建一套属于你这个物种的“专属文库”。简而言之，RepeatModeler就是在你没有任何先验知识的情况下，让算法自己找出基因组里哪些序列是重复的、它们属于什么家族、边界在哪里。

1.2 从头建模的几类算法支柱

RepeatModeler不是单一算法，它是一个流程，内部会调用多个算法模块协同工作。理解这些底层算法，有助于你判断运行结果是否可靠，而不是当一个只会敲命令的工具人。

• RECON：经典的重复序列家族识别算法，基于全基因组自我比对，利用多序列联配的信息来寻找重复家族。它的特点是敏感度高，但边界往往不够精确。
• RepeatScout：使用k-mer频率统计的方法，通过寻找高频k-mer来发现重复单元，速度快，适合发现高频重复。
• TRF（Tandem Repeats Finder）：专门负责串联重复序列（tandem repeats）的识别，比如着丝粒、端粒附近的大量短串联重复。
• LTR_retriever（通过-LTRStruct启用）：专门处理LTR反转座子。LTR反转座子在植物和部分动物基因组中占比极高，它的结构特征是两侧有长末端重复序列（LTR）、中间有编码区，这个模块会专门分析这类结构。

RepeatModeler 2.0以上的版本把这些工具的输出统一整合，输出一套非冗余的重复序列文库。要注意的是，它输出的文库是以家族为粒度的，也就是说，它给出的是一个家族里代表序列的“共识序列”，而不是基因组里每一个拷贝。同一家族不同拷贝之间存在序列差异，这些差异在后续RepeatMasker比对时会被容忍。

2. 环境准备：装RepeatModeler和RepeatMasker最容易踩的坑

先别急着敲命令，环境装不对，后面每一步都可能报错。这套工具链的依赖关系比较复杂，我建议用Conda或者Singularity容器管理，不要在系统Python环境里裸装，那样gcc版本、Perl模块和库文件冲突能让你折腾一整天。

2.1 核心依赖的来龙去脉

RepeatMasker本身是一个由Perl脚本驱动的流程，它的比对引擎默认是rmblast——这是NCBI BLAST的一个fork版本，专门针对重复序列注释优化过。除此之外，它还需要TRF来处理串联重复，以及famdb.py（Python脚本）来读取Dfam数据库。

RepeatModeler在RepeatMasker的基础上，又集成了RECON、RepeatScout等工具。4.0.0以上版本的RepeatModeler配合4.1.x版本的RepeatMasker，是目前最稳定的组合，建议尽量用这个版本组合。

这里有一个容易翻车的细节：RepeatMasker有自己的rmblast依赖版本，与系统里其他工具链可能冲突。

2.2 用Conda隔离环境实操

我的建议是在一个全新的conda环境中同时安装两者，让conda帮你处理依赖关系。

conda create -n repeat -c bioconda -c conda-forge repeatmodeler repeatmasker conda activate repeat

安装完成之后，先测试一下两个主要命令是否正常：

RepeatModeler -version RepeatMasker -version

如果能正常输出版本号，说明基础安装没问题。RepeatMasker第一次运行时需要设置Dfam数据库，或者在你手动指定-lib的情况下不依赖数据库。我们的场景是使用RepeatModeler生成的文库，所以不强制需要Dfam在线配置，但建议还是把Dfam装上，万一以后需要对比或者补充注释会用得上。

2.3 版本兼容和perl模块检查

RepeatMasker和RepeatModeler都是Perl脚本，对Perl内置模块有依赖。conda环境里一般会处理好，但如果你是从源码编译安装，务必检查JSON::PP、File::Which、Thread::Queue这几个常用模块。特别是执行时如果出现Can't locate JSON/PP.pm in @INC这类错误，说明Perl模块缺失。

另外重要的一点：不要把RepeatModeler/RepeatMasker跟其他版本敏感的比对软件放在同一个环境中使用。比如系统里如果装有另一个版本的rmblast或者ncbi-blast，PATH顺序会导致RepeatMasker调用到错误的比对引擎，最典型的报错是rmblast: error while loading shared libraries或者BLAST database format mismatch。此时检查一下当前环境里rmblast实际指向哪里：

which rmblast rmblast -version

如果指向的版本不是RepeatMasker需要的版本，最简单的做法是回到独立conda环境（如上面的repeat环境），确保PATH里面优先是conda的bin目录。

提示：很多人在同一个环境里同时装了EDTA、LTR_retriever和RepeatMasker，这时候版本牵扯最麻烦。建议EDTA单独一个环境，RepeatModeler/RepeatMasker单独一个环境，两者通过文件交换数据，互不干扰。

3. RepeatModeler从头建模：跑通核心流程

当你的环境准备妥当后，就可以开始正式的建模了。但建库之前，有一步容易被忽略的预处理，会对结果质量产生重大影响。

3.1 输入基因组的预处理检查

RepeatModeler对输入fasta的要求比较严格。我现在每次跑之前都会用seqkit stats先看一眼：

seqkit stats genome.fa

重点检查三件事：

• 序列ID不能有空格、|、:等特殊字符。RepeatModeler在内部分发任务时，这些字符容易造成脚本解析异常，报错信息也不直观。
• 过短的序列建议过滤掉。我一般会把小于500bp的contig删掉，因为太短的序列很难被聚类为有意义的重复家族，还会增加比对噪声。可以用seqkit seq -m 500做长度过滤。
• 序列不能全是大写或全是小写之外的奇怪字符，比如N占比过多、有非ACGTN字符，这些会在建库时干扰统计。

一个稳妥做法是：

seqkit seq -m 500 -g -N genome.fa > genome.clean.fa

-g表示去除gap，-N表示去除含过多N的序列，具体阈值可以用--max-n参数调整。处理完之后，再用seqkit stats确认一下数据规模。

注意：这个“清理”并不是说把N全部删掉，而是把以N为主的低质量序列替换成一条干净的序列，避免repeat建模被垃圾序列带偏。

3.2 构建数据库并启动RepeatModeler

RepeatModeler的第一步是把fasta文件转换成它内部的BLAST数据库格式。命令如下：

BuildDatabase -name my_species genome.clean.fa

这一步会生成以my_species为前缀的BLAST数据库文件。之后启动建模主流程：

RepeatModeler -database my_species -pa 16 -LTRStruct

参数说明：

• -database：指定刚才构建的数据库名字，不是fasta文件。
• -pa：并行运行的线程数，既影响CPU核数，也决定内存占用。一般建议用16到32个线程，不要盲目开到上百，因为内存会成倍增长。
• -LTRStruct：启用LTR反转座子的结构分析模块。这一项我建议默认打开，特别是对于植物基因组，LTR占了重复序列的很大比重，不打开会遗漏大量LTR。

运行过程中，屏幕会实时打印当前阶段。RepeatModeler内部主要分3个阶段：全基因组自我比对和家族聚类（RECON、RepeatScout）、LTR结构识别、最终的整合分类。整个流程在哺乳动物规模基因组（约3Gb）上，以16线程运行，通常需要2到5天的时间；如果遇到重复序列异常丰富的基因组（如小麦、玉米），时间可能更长。没有捷径可走，耐心等就行。

3.3 建模后的初始文库与处理逻辑

运行结束之后，会在输出目录（默认是RM_xxx）里生成以下关键文件：

• consensi.fa.classified：已经被分类到具体重复家族的共有序列。
• consensi.fa.unclassified：无法明确分类到已知家族的序列。
• round-*.fam：中间过程文件，一般不用管。

其中consensi.fa.unclassified非常重要，它里面往往藏着一些这个物种特有的、公共库完全未知的新家族序列。正确地处理这些未知序列，往往能显著提升后续RepeatMasker的发现率。

实际操作中，我会把两类合在一起去冗余后再用：

cat consensi.fa.classified consensi.fa.unclassified > my_species.raw.lib cd-hit-est -i my_species.raw.lib -o my_species.lib -c 0.8 -n 8 -M 0 -T 16

这里用了cd-hit-est来做冗余去除，-c 0.8表示一致性80%以上的序列视为冗余将其合并。这里之所以用80%而不是更高阈值，是因为同一重复家族的不同拷贝本身就有序列分歧，如果阈值设得过高（比如0.95），去冗余效果会很差，最终文库体积过大，重复家族会被拆得七零八落，影响下游注释的整洁度。

另外别忘了看分类情况。可以用grep统计一下库里各类重复序列的数量：

grep -c "LTR" my_species.lib grep -c "DNA transposon" my_species.lib grep -c "Unknown" my_species.lib

Unknown比例如果超过30%，通常说明部分序列在RepeatClassifier中没能匹配到已知蛋白结构域——这些序列有可能是新家族，也有可能是组装错误造成的异常序列，建议保留但后续关注。

4. RepeatMasker注释：用自家库把基因组完整扫一遍

文库构建完成之后，终于可以进入真正的全基因组注释环节。这一步做的事情很简单：用你建好的文库，对基因组中每一个位置进行比对，判断它属于哪个重复家族，并标记出精确的坐标边界。

4.1 屏蔽策略：软屏蔽还是硬屏蔽

RepeatMasker在标记重复序列之后，会生成一个“屏蔽后”的基因组序列文件。屏蔽有两种方式：

• 硬屏蔽：把重复序列的碱基全部替换成N。
• 软屏蔽：把重复序列的碱基改成小写（原来是大写），保留原始序列信息。

我强烈建议使用软屏蔽策略（RepeatMasker的-xsmall参数）。原因很简单：虽然屏蔽是为了防止下游工具误用重复序列，但硬屏蔽之后序列信息就彻底丢失了——万一某个基因注释流程需要在重复区域中重新比对特异reads，或者你需要手动检查某个被屏蔽位点是否真的是转座子，小写序列仍然保留了完整的碱基信息。绝大多数现代注释流程（MAKER、BRAKER）都能正确处理软屏蔽的fasta文件。

4.2 完整命令示例与参数说明

RepeatMasker -lib my_species.lib -pa 16 -xsmall -gff -dir repeat_mask_out genome.clean.fa

参数逐项说明：

• -lib：指定你自己的文库，这是与“依赖公共库”模式最大的区别。
• -pa：并行线程数，这个参数跟RepeatModeler里的含义一致。
• -xsmall：软屏蔽，保留小写碱基。
• -gff：输出GFF3格式的注释文件，这个后面做可视化或者作为证据输入给基因注释工具时非常有用。
• -dir：指定输出目录。RepeatMasker默认会把结果输出到输入文件所在目录，建议用这个参数保持输出文件集中整洁。

跑完之后，输出目录里会有以下几类文件，这是解读结果的核心：

4.3 如何高效解读.out和.tbl文件

.tbl文件是全局扫描结果的一页纸总结：

================================================== file name : genome.clean.fa sequences : 45 total length : 2500000000 bp (2.50 Gb) GC level : 40.12 % bases masked : 1200000000 bp ( 48.00 %) ==================================================

这里的主要指标是bases masked的百分比。对于哺乳动物基因组，40%-50%左右是一个正常水平；对于植物，60%-80%之间都很常见；如果你的物种本身已知重复含量极高，却没有达到相应的标准，那就要回头检查你的文库质量。

.out文件每一行代表一个“repeat match”，格式类似BLAST输出。列依次是：SW score、序列名、起始位置、终止位置、重复类别、重复家族名、分歧度等。比如典型的行：

2151 0.1 0.0 0.0 chr1 1000 2500 (p) LTR/Gypsy RLC_fam1 1000 2500 3e-20

解读时有一个关键指标：分歧度（divergence），就是格式里的0.1位置。它表示当前拷贝与文库共识序列的差异百分比。如果大量匹配的分歧度都很高（比如>15%），说明这些拷贝是进化上比较古老的插入事件，或者是同一超家族里已经分化出多个亚家族。如果你发现集中高分歧序列占比异常大，可能需要回到重复家族分类阶段做更深层的亚家族划分。

.gff文件则是做下游可视化和筛选的利器，比如你可以在IGV中直接加载，看到每类转座子在染色体上的分布密度。

5. 质量评估与结果修正：屏蔽率合理才算真完成

RepeatMasker跑完，不代表你的重复注释就算完成了。我见过不少同学拿到tbl文件看到“bases masked 45%”就心满意足地去做下一步了，但很多时候这个数字是假象——要么漏了很多，要么把一些非编码保守区域也屏蔽掉了。以下几个关键检查点供你参考。

5.1 先跟近缘物种比对一下屏蔽率

每个类群都有自己的“参考范围”。我通常的做法是去查一下该物种近缘种（最好同属或者同科）已发表的基因组重复序列占比。比如你注释一个禾本科植物，水稻的重复序列占比约40%-50%，玉米约80%——但如果你的物种跟玉米亲缘较远、组装又是Scaffold级别，那最终报出85%的屏蔽率就需要警惕了，有可能是组装把大量未定位序列堆叠在一起造成了虚假重复。

如果你的物种屏蔽率显著低于近缘种，最可能的原因是RepeatModeler建库时遗漏了某些高拷贝家族。这时可以做一个快速检验：用RepeatMasker自带的-species模式比对Dfam数据库看看是否有明显更高的覆盖率：

RepeatMasker -species maize -pa 16 -xsmall genome.clean.fa

如果公共库的屏蔽率明显高于你的自建库，说明你的文库确实不够全。但也不必气馁，因为公共库的屏蔽率高有时也是“假阳性”——它把某些低复杂度区域也标成了转座子。更合理的做法是比较两者结果的交集，看看你的库到底漏掉了哪些真实重复。

5.2 Unknown分类过多时的对策

RepeatMasker输出的注释里，如果Unknown类占比过高（比如超过总屏蔽区域的20%），说明有很多序列虽然被比对到了重复区域，但无法被归类到任何已知重复家族。这有两种解读：

• 你的物种确实存在大量新家族。这在一些基因组独特、转座子谱系特殊的物种中很常见。
• RepeatClassifier的分类效果不好。RepeatClassifier依赖RepeatMasker自带的分类器数据，如果这些数据本身覆盖不足，则原本是已知家族的序列也可能被归为Unknown。

针对第二种情况，一个有效操作是把Unknown序列提取出来，手动和已知蛋白结构域数据库比对一遍。比如用InterProScan或者hmmscan搜索转座子蛋白结构域（如反转录酶、整合酶、转座酶），如果能命中已知结构域，就可以手动把分类标签补上。实在没有工具的，至少在嵌套注释里保留Unknown，不要随手删掉，因为它们也是基因组真实重复的一部分。

5.3 用注释结果反过来优化建模

重复序列注释是一个可以迭代优化的过程。我一般会把RepeatMasker输出中那些“高可信但未分类”的序列再次加入到库中，重新跑一轮RepeatModeler：

• 从.out文件中提取长度较长（>500 bp）且匹配到Unknown的区间。
• 将这些区间序列与已经建立的文库合并。
• 对这个合并库再做一次cd-hit去冗余。
• 用去冗余后的库重新RepeatMasker。

这个迭代过程通常能多找回5-10%的注释区域。不过也要适可而止，两三轮之后提升就开始边际递减，再往后加的多数是ambigous的低复杂度噪声。

6. 实战中的性能瓶颈与调优

最后聊聊性能问题。RepeatModeler和RepeatMasker在大基因组版本上真的是“吞机器兽”，尤其是RepeatModeler的自我比对阶段，内存和CPU的消耗很凶。这一部分我整理了实战中比较常见的瓶颈和应对经验。

6.1 内存和线程的合理分配

RepeatModeler在默认情况下，每个-pa线程大约需要2-4GB内存。一个3Gb基因组，16线程运行，高峰期内存可能冲到64GB以上。如果服务器内存不足，优先减少线程数量，而不是降低任务并行度。因为RepeatModeler内部的并行策略往往是按序列块分配的，线程太多但内存不够时，反而会造成整机swap崩溃。

RepeatMasker的内存需求相对温和，但它的并行效率更高，-pa直接控制比对线程数。在大型基因组上，我一般设置16到32个线程，运行时间通常能控制在10-20小时。如果你机器核数很多，可以适当地把线程加到64，但再往上提升就有限了，因为I/O和数据库访问会成为瓶颈。

6.2 RepeatModeler卡在LTR检索环节的排查

LTR检索（-LTRStruct）阶段经常是运行时间最长的部分，也是最容易“看起来卡住”的部分。

我遇到过一次：某植物基因组跑了两天还停留在LTR Retriever阶段，而且屏幕上长期不刷新。排查了很久，最后发现是建库时的序列ID里包含了一个不常见的unicode字符，导致LTR_retriever在处理子序列集时反复匹配失败并重试。清理ID之后，这个问题就不再出现。

所以如果你也遇到长时间卡在某一步，先做两件事：

1. 用top查看是否有某个子进程还在消耗CPU，如果完全是0%且长时间不变，基本可以断定是“死锁”，而不是“在计算”。
2. 检查序列ID的可打印字符范围。规范的做法是只保留字母、数字、下划线和连字符，其他一律替换或删除。

6.3 超大基因组的分区运行策略

有些基因组动辄十几Gb甚至几十Gb（比如小麦、松树），单体跑RepeatMasker会非常耗时，而且一旦中途报错，全部重来代价极高。这时可以用RepeatMasker自带的-pa与-gff配合，将基因组按染色体分拆成多个部分分别运行。

实际操作中，我建议把序列分成几组，每组大约2-3Gb的fasta文件，分别在不同节点上运行RepeatMasker，最后再合并结果。合并时注意tbl文件的统计不可直接相加，需要将各部分的.out和.gff合并后重新统计一次，或者直接写个简单脚本按行拼接后重新汇总。.masked文件可以按染色体拼接回去，只要保证每条序列唯一即可。

6.4 结果文件的归档和可重复性

把重复序列注释跑完之后，有几样东西一定要保存下来，不然后面写论文方法部分会抓瞎：

• 重复序列文库my_species.lib：这是你最重要的资产，不仅本物种后续需要，近缘物种做注释时也可以直接用。
• 建库的准确命令和参数：RepeatModeler的版本号、运行参数、随机种子（如果有）都要记录，最好写进一个README或Makefile。
• 最终.out和.tbl文件：这些是论文表格和图形的数据来源。
• 关键的中间日志：哪个步骤运行了多久、有没有出现警告，都会成为后面排查问题和回答审稿人提问的依据。

我习惯给每个基因组的重复注释建一个固定目录结构：

RepeatAnnotation/ ├── 01_clean_fasta/ ├── 02_repeatmodeler/ │ └── my_species.lib ├── 03_repeatmasker/ │ ├── genome.clean.fa.masked │ ├── genome.clean.fa.out │ ├── genome.clean.fa.tbl │ └── genome.clean.fa.out.gff ├── 04_summary/ └── README.md

这样即使过了半年再回来找结果，也不会对着文件名一头雾水。

回到我在最开头提到的那次教训——拿到一个没有现成文库的基因组时，别再图省事直接拿公共库硬跑。老老实实RepeatModeler建库、RepeatMasker注释、检查屏蔽率、迭代优化一轮，整个流程走下来虽然要花几天时间，但换来的是下游所有分析都站在一个干净的地基上。特别是基因注释结果如果被审稿人质疑假基因过多时，你可以理直气壮地说：重复序列屏蔽是完整走完的，有生成文库和注释结果的完整记录。

这套流程还有一个额外的好处：你建出来的my_species.lib，后续在做近缘物种比较基因组学、转座子进化分析、甚至群体遗传学的reads过滤时都能反复利用，相当于一次投入、多处收益。所以第一次花时间把库做好、做全，绝对不亏。