基因组组装完成后，跑完三代组装、Hi-C挂载、纠错这些大工程，很多人会下意识松一口气——但如果你直接把基因组丢给BRAKER或者Augustus做基因结构预测，接下来大概率会收获一堆结构错乱的基因模型。这不是组装的问题，而是绕过了重复序列注释这一步。RepeatModeler和RepeatMasker就是处理这个环节的两件标准工具：前者从目标基因组里从头挖掘重复家族并建模，后者用模型库把基因组里的重复序列找出来、标记好。这篇文章会从原理讲起，把建库、注释、结果解读到报错排查的完整流程都过一遍，适合刚组装完基因组、准备做下游注释的初学者，也适合想回头补做重复注释的老手。
1. 重复序列注释：组装完拿到基因组后的第一件事
很多人对重复序列的理解停留在“垃圾DNA”阶段，觉得不注释也能硬着头皮往下走。实际上，真核生物的基因组里转座子、串联重复、rDNA这些序列的占比远超预期：人类的基因组大约45%是重复序列，玉米和小麦这类植物基因组动辄80%到85%以上。你组装出来的每一个contig或scaffold里，都有大量几乎一模一样的拷贝散落在不同位置。如果不先做注释和屏蔽，后续所有基于序列比对和特征统计的分析都会被它们搅浑。
我见过最典型的翻车现场，发生在基因结构预测环节。基因预测软件会扫描开放阅读框、剪接位点、编码潜能等特征，而转座子编码区在序列特征上非常像真实基因，特别是Gypsy和Copia这类LTR反转座子，自带完整的gag、pol、env结构域，和宿主基因的编码特征高度相似。结果就是预测软件把转座子区域当成了基因，输出一堆假基因模型；反过来，真实基因如果内部插入了重复序列，又会被当成基因间区切掉，造成外显子遗漏。最后做功能注释的时候，BLAST比对结果一片混乱，蛋白结构域注释全都对不上号。
重复序列的影响范围还不止基因注释。做全基因组比对时，重复区域会产生大量假的共线性匹配，直接影响系统发育分析和比较基因组学的结论；做群体变异检测时，转座子插入多态性会干扰SNP/InDel的calling，导致假阳性变异满天飞；就连最基础的GC含量统计和染色体可视化，重复序列不处理，结果图都画不干净。所以我的经验是：组装完成后的第一件事，优先做重复注释和屏蔽，它直接决定了后续所有分析的地基稳不稳。
2. 工具分工与运行环境：RepeatModeler和RepeatMasker各自解决什么问题
2.1 一个负责“造库”，一个负责“用库”
这两件工具的名字看起来很像，容易让人混淆。我的理解方式很简单：RepeatModeler是“造库”的，RepeatMasker是“用库”的。
RepeatModeler做的是从头（de novo）重复建模。它不需要任何物种预先存在的重复序列数据库，而是直接对输入的基因组序列进行自我比对和搜索，通过内部整合的RECON、RepeatScout、TRF等算法，把基因组里反复出现的相似序列片段聚类、组装成一个个consensus序列，也就是重复家族的“模板”。如果加了-LTRStruct参数，它还会调用LTR_retriever，专门针对LTR反转座子做更精细的结构预测。跑完之后你会得到一个包含若干条representative consensus的库文件，这个库就是整个重复注释的核心资产。
RepeatMasker则拿着这个库，用RMBlast比对引擎去扫描整个基因组，把每一处和库中序列同源的位置都找出来，标记它的坐标、方向、家族分类、变异程度，然后根据你的要求把这段序列硬屏蔽成N，或者软屏蔽成小写字母。输出文件默认会给四个：屏蔽后的基因组序列、逐条注释表、汇总统计表、GFF注释文件。整条pipeline的产物专业、规范，也是目前被同行普遍接受的标准格式。
2.2 安装依赖和数据库准备
安装方面，最省心的方式是直接用conda建一个干净的环境：
conda create -n repeat -c bioconda repeatmodeler repeatmasker rmblast -y conda activate repeat装完之后建议先检查版本，确认RMBlast已经被RepeatMasker正确识别：
RepeatMasker -v如果输出里显示RMBlast不可用，说明环境变量PATH没有指向正确的rmblast目录，需要手动配置。这个坑我踩过不止一次，conda装完经常因为版本冲突导致RMBlast版本不匹配，建议装完后用一个小测试序列立刻验证。
数据层面，RepeatMasker默认会调用自带的Dfam数据库，Dfam是开放的，覆盖了相当一部分真核生物的重复家族。Repbase是另一个重要的商业数据库，覆盖物种更广、注释更精细，但需要注册并签署学术使用协议，下载后放到RepeatMasker的Libraries目录下。如果你研究的物种不是模式生物，我的建议是：重点依赖RepeatModeler自己生成的从头库，Dfam和Repbase作为补充，不要只依赖现成数据库——非模式物种的重复序列往往与库中已知家族差异很大，只跑RepeatMasker会把大量真实重复漏掉。
3. 完整实操流程：从建库到注释
3.1 预处理：先别急着建库
拿到基因组文件后，先别急着跑BuildDatabase，花十分钟把输入文件清理干净，能省下后面几天的时间。
第一步，查看基本统计信息：
seqkit stats genome.fa第二步，过滤掉过短的contig。RepeatModeler对碎片化输入很敏感，一堆几百bp的短contig会严重拖慢RECON的聚类过程，而且组装出的consensus质量也不高。一般我会用长度过滤条件，比如只保留1kb以上的序列：
seqkit seq -g -m 1000 -o genome.clean.fa genome.fa第三步，检查序列名。RepeatModeler和RepeatMasker都要求序列名简洁唯一，不能有空格、竖线、括号等特殊字符。装配工具生成的fasta头部有时会带额外的描述信息，这些都要清掉：
# 检查重复的序列名 seqkit seq -n genome.clean.fa | sort | uniq -d # 去掉序列名中第一个空格之后的内容（在ID后加空格，后面写描述） seqkit replace -p ' .*' -r '' genome.clean.fa -o genome.clean.fixed.fa如果你组装时用了线粒体、叶绿体序列，最好先把细胞器基因组分离出去。RepeatModeler不会区分核基因组和细胞器基因组，细胞器里大量的重复结构会被当成核内重复家族建模，白白增加后续的干扰。
3.2 BuildDatabase：把fasta转成BLAST数据库
RepeatModeler的第一步是构建一个属于自己的序列数据库格式：
BuildDatabase -name my_species_db genome.clean.fixed.fa这里的-name参数是数据库前缀，注意别带点和空格。命令跑完后，目录里会出现my_species_db.nhr、my_species_db.nin等格式文件，这就是RMBlast能识别的数据库文件。这个步骤通常很快，但如果你的基因组非常大（比如几十Gb），磁盘IO会成为瓶颈，建议放到SSD或NVMe盘上跑。
有一个细节：BuildDatabase默认只建核甘酸数据库，不需要额外指定BLAST类型；但如果你的基因组文件里有非常长的scaffold，它会在内部自动切分索引，这个不用管。
3.3 RepeatModeler从头建模：这是整个流程里最耗时的一步
接下来是重头戏：
nohup RepeatModeler -database my_species_db -pa 16 -LTRStruct > repeatmodeler.log 2>&1 &-pa参数控制并行线程数，一般按CPU核数的一半到三分之二设置即可，不是越大越好。RECON阶段存在内存瓶颈，线程开太猛容易直接把内存顶爆。我之前在一台128核、512GB内存的服务器上跑一个2.5Gb的植物基因组，-pa开64，结果RECON阶段直接内存溢出；降到32才稳定下来。
-LTRStruct参数强烈建议加上。它会调用LTR_retriever专门处理LTR反转座子，对LTR富集的基因组效果提升非常明显。如果没有这个参数，LTR的完整结构会被拆得七零八落，库的质量会打折扣。
运行时间方面，小基因组（500Mb以下）可能几小时就能跑完，到了一两个Gb级别的基因组往往要跑好几天，超大基因组跑几周都很正常。中间过程中间，工作目录里会出现RM_xxx.xxxxxx这样的临时目录，千万别删，最后的结果就藏在里面。日志要经常看：
tail -f repeatmodeler.log如果发现日志停在某一轮刷屏很久，可以先看看是不是在跑RECON的某个batch；如果长时间没有任何输出，再检查内存、磁盘是否够用。
跑完之后，进入RM_开头的目录，重点看这几个文件：
ls RM_* consensi.fa consensi.fa.classifiedconsensi.fa是全部重复家族consensus序列；consensi.fa.classified是经过RepeatClassifier分类后的版本。分类后的版本里，每条序列名后会带上类似#LTR/Gypsy、#DNA/hAT、#LINE/L1这样的分类标签，直接用这个文件作为RepeatMasker的库就行。
3.4 处理unknown家族：分类结果不理想时怎么办
RepeatClassifier会把无法明确归类的序列标成Unknown，这些unknown在后续的RepeatMasker输出里也会原样保留。如果unknown占比过高（比如超过20%），说明库质量不佳，或者这个物种的重复序列太特化、太老化。
我的处理方案是：把RepeatModeler输出的consensi库和Dfam的已知库合并，再重新跑一遍RepeatMasker；合并后的库会让一些原本被标成unknown的短序列重新找到同源归属。另一种做法是拿unknown序列去NCBI的nr库里做BLASTx，看它们是否编码逆转录酶、转座酶这类结构域，以此反推家族类型。这个过程比较费人工，但注释出来的结果更精细。
3.5 RepeatMasker注释：最后一步，也最需要细心
拿到合格的库之后（假设已经合并处理好，命名为final_lib.fa），新建输出目录并运行：
mkdir -p mask_output RepeatMasker -lib final_lib.fa -pa 32 -gff -xsmall -dir mask_output genome.clean.fixed.fa解释一下关键参数：
-lib：指定重复库文件，和-species参数互斥。用自己从头建的库就选它。-gff：输出GFF3格式的注释文件，下游做基因结构注释、可视化都用得到。-xsmall：软屏蔽。把重复序列区域转成小写字母，而不是替换成N。软屏蔽的好处是保留了序列的原始信息，做基因预测时，部分软件能利用小写区域来提高注释准确度，比如BRAKER的默认流程就支持软屏蔽基因组。-pa：并行线程数。注意RepeatMasker是按输入序列拆分任务的，如果你的输入文件里只有一条超长scaffold，开再高的并行度也只有一个线程在干活。这时候可以先用seqkit把scaffold按固定窗口切分，或者在组装时尽量保证序列数量足够多。
跑完之后，mask_output目录下会出现以下文件：
ls mask_output genome.clean.fixed.fa.masked genome.clean.fixed.fa.out genome.clean.fixed.fa.tbl genome.clean.fixed.fa.out.gff到这一步，重复序列注释的核心流程就结束了。但拿到文件只是开始，结果解读才是重头戏。
4. 结果文件解读：别只盯着masked基因组
4.1 四个输出文件，各有各的用途
很多人跑完RepeatMasker之后只拿.masked文件去跑下游流程，把.out和.tbl文件晾在一边。这其实损失了大量信息。我把四个文件的用途整理成了表：
| 文件 | 内容 | 主要用途 |
|---|---|---|
| .masked | 屏蔽/软屏蔽后的基因组序列 | 基因预测、序列比对等下游分析的输入 |
| .out | 逐条重复注释明细（每条位置、方向、家族、变异率） | 精细分析，如转座子插入时间推断、家族分布统计 |
| .tbl | 汇总统计表（各类重复占比、条数） | 快速查看重复序列总体含量，出柱状图、饼图 |
| .gff | GFF3格式的结构化注释 | 基因组浏览器可视化、纳入标准pipeline |
其中.tbl文件我建议第一步就看，它能快速告诉你这个基因组的重复含量总体情况：
================================================== total length: 123456789 bp GC level: 38.22 % bases masked: 65432100 bp ( 53.06 % ) ==================================================bases masked这一栏的百分比就是一个基因组“重复程度”的直接指标。如果这个值低得离谱，比如人类基因组跑出来只有5%，那基本可以断定库没选对或参数有问题。
4.2 .out文件里最容易看错的几列
.out文件是制表符分隔的明细表，每一行代表一个重复序列片段。列的数量不少，但最核心的几个含义如下：
| 列 | 含义 |
|---|---|
| SW score | Smith-Waterman对比得分，越高表示匹配越强 |
| perc div | 该拷贝与consensus之间的突变率（divergence） |
| perc del / perc ins | 缺失率和插入率 |
| query sequence | 基因组上的序列名 |
| query begin / end / left | 该片段在基因组上的起止位置及左侧剩余长度 |
| matching repeat | 匹配到的重复家族名称 |
| repeat begin / end / left | 该片段在consensus序列上的对应位置 |
比较容易混淆的是perc div这列，很多人以为是“序列相似度”，实际上它是差异度，数字越大表示这个拷贝和consensus差异越大，通常也意味着插入时间更古老。结合LTR反转座子5'端和3'端LTR的差异率，可以估算转座子的爆发时间，这是转座子演化分析的标准做法。
还有一个细节：匹配到互补链的时候，方向会在matching repeat或query位置信息中标出C（complementary），分析时要注意区分正负链。负链的重复经常会被忽略，如果再往下游提取序列，没考虑方向会提取出反向互补序列，导致后续分析全部跑偏。
4.3 用自带工具快速出图
RepeatMasker安装目录的util子目录里有一堆现成的脚本，不用自己造轮子：
# 将.out转为GFF3（如果之前没加-gff参数） util/rmOutToGFF3.pl genome.clean.fixed.fa.out > repeat.gff3 # 按家族分类汇总统计 util/processRepeats.pl genome.clean.fixed.fa # 计算divergence分布，输出可读的文件 util/divergence.pl genome.clean.fixed.fa.out > divergence.txtprocessRepeats.pl跑完会在同目录下生成一堆以div、bp、cnt开头的统计文件，配合R语言就能画出一张重复序列divergence分布图。这张图能直观反映重复序列在基因组演化时间轴上的分布规律，是论文里常用的图。
5. 参数调优和并行加速：大基因组怎么跑才高效
5.1 RepeatModeler的断点续跑和资源控制
RepeatModeler最让人头疼的一个问题就是时间太长，而且一旦中途挂掉，前面跑的东西全废。好在它支持断点续跑：只要RM_开头的目录还在，直接用同样的命令再启动一次，它会自动检测到已有的运行记录并继续。
举一个实际的例子：
# 第一次运行挂了之后，不加-nohup这些花活，直接重新跑 RepeatModeler -database my_species_db -pa 16 -LTRStruct -recoverDir RM_12345.abcdefg-recoverDir参数后面跟的就是之前的临时目录名，这个参数救过我很多次。另外，建议在运行前把当前工作目录的磁盘剩余空间检查一遍，RepeatModeler中间产生的临时文件体量可能达到基因组的10倍以上，特别是RECON阶段的中间文件，硬盘不够会很尴尬。
5.2 RepeatMasker的并行度陷阱
前面提过，RepeatMasker的并行是按输入序列数量来分的。如果你输入的是一个染色体级别的fasta文件，里面只有几十条超长序列，那-pa 32实际只会有几十个任务在跑，绝大多数核都在摸鱼。更常见的问题是内存争抢：每个线程同时扫描一条超大scaffold时，内存占用飙升，服务器卡到SSH都连不上。
解决办法有两个方向。一是把超长序列按窗格切分，比如每10Mb一段，切完之后再跑。切的时候注意保留坐标信息，或者用专门的工具（如seqkit sliding）生成带坐标的窗口序列，方便后续把结果拼回去。二是把RepeatMasker分染色体提交，每条染色体一个作业，互不干扰。
还有一个参数细节：如果用了-gff，会额外产生GFF输出的计算开销；如果只需要统计信息，可以不加-gff。另外，-html参数会生成网页报告，但代价是运行时间变长，服务器上跑大批量数据时一般不建议加。
5.3 库的质量决定了注释的上限
我一直强调，RepeatMasker只是执行者，真正的灵魂在于库。如果你的库里面只有几十条consensus，而物种基因组实际有几千个重复家族，那RepeatMasker再努力也只能注释出冰山一角。所以在继续往下分析之前，可以做一个快速的质量体检：
# 统计库里的家族数量 grep -c ">" final_lib.fa # 统计分类后的家族类型分布 grep ">" final_lib.fa | sed 's/.*#//' | sort | uniq -c | sort -k1 -rn如果分类类型明显偏少，说明RepeatModeler的聚类不够充分，可以先检查输入的基因组是否包含足够的重复拷贝信息。组装质量太差、重复序列被压成极少几段时，RepeatModeler很难建立高质量模型。这时候回头优化组装比硬调参数更有用。
6. 高频报错排查与避坑指南
整个流程里，你会遇到各种莫名其妙的报错。下面这张表是我在实际项目中积累的，每一条都踩过或帮别人排查过。
| 现象 | 可能原因 | 解决办法 |
|---|---|---|
| BuildDatabase报duplicate sequence names | fasta里存在重复序列名 | 用seqkit seq -n检查，再用seqkit rename或replace改名 |
| RepeatModeler在RECON阶段内存溢出 | 线程开太多，基因组过大 | 降低-pa，限制单批次内存；给作业脚本加内存上限（如-Xmx） |
| RepeatModeler跑完了但consensi.fa是空的 | 输入基因组重复含量极低，或序列碎片化太严重 | 检查输入文件是否只有几条短序列；过滤短contig并重跑 |
| RepeatMasker .tbl显示0% masked | 库文件路径错误，或库与物种差异太大 | 确认-lib参数生效，先跑一个已知重复富集的染色体测试 |
| RMBlast报段错误（segmentation fault） | 输入fasta里有非法字符，或线程过多 | 清理fasta，只用ACGTN字符；降低-pa重新跑 |
| 输出文件里unknown比例超过30% | 库质量差，或物种重复序列特化严重 | 合并Dfam/Repbase库重跑，或人工分类unknown |
| 运行日志长时间不更新 | 内存不足导致进程卡死，或作业被系统OOM killer杀掉 | 检查dmesg、作业系统日志，降低-pa或增加内存重跑 |
| 软屏蔽后小写区域被下游软件当成小写N | 引用了不支持软屏蔽的注释工具 | 检查工具文档；部分流程需要显式开启soft-masking支持 |
这里我特别想展开说一个隐蔽的坑：RepeatMasker用的RMBlast对fasta里的非法字符非常敏感。看起来一模一样的fasta文件，可能因为里面混入了IUPAC模糊碱基以外的字符就报错。建议在预处理阶段强制把所有碱基统一成ACGTN，其余一律替换掉：
seqkit seq -w 0 genome.clean.fixed.fa | tr 'RYMKSW' 'NNNNNN' > genome.clean.strict.fa另一个值得注意的地方，是RepeatModeler和RepeatMasker对路径的依赖。它们内部会调用相对路径，如果你在别的目录里直接指定绝对路径运行，有时会找不到依赖文件。我的习惯是：专门建一个项目目录，所有输入、数据库、输出全部放在同一个目录层级下，避免路径混乱。
7. 实操中的个人体会
我自己的习惯是把重复注释做成一个可复用的流程，而不是每次手动敲命令。即使只是单人项目，也强烈建议用Snakemake或Nextflow把BuildDatabase、RepeatModeler、RepeatMasker、结果统计这些步骤串起来，这样后续换基因组、换参数、或者拿到更新版本的库时，一条命令就能重跑整个流程。
跑大型基因组的过程中，还有一个容易被忽略的点：别在登录节点上直接运行。RepeatModeler和RepeatMasker都是长时间、大内存的作业，务必通过作业调度系统（SLURM/PBS）提交。日志文件也要定期检查，包括磁盘空间和内存占用，很多失败其实都是可以提前预判的。
从开始跑RepeatModeler到最终得到一份满意的重复注释结果，往往需要多次迭代，每次都要认真看tbl、看divergence分布、抽检具体区域。注释结果的好坏不只影响一篇论文，更会沉淀成后续所有分析的基础资源。把这些步骤走扎实，后面做基因家族分析、比较基因组、群体遗传的时候，你会省下大量返工的时间。
