基因组组装完成后跑完三代组装、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 21 -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汇总统计表各类重复占比、条数快速查看重复序列总体含量出柱状图、饼图.gffGFF3格式的结构化注释基因组浏览器可视化、纳入标准pipeline其中.tbl文件我建议第一步就看它能快速告诉你这个基因组的重复含量总体情况 total length: 123456789 bp GC level: 38.22 % bases masked: 65432100 bp ( 53.06 % ) bases masked这一栏的百分比就是一个基因组“重复程度”的直接指标。如果这个值低得离谱比如人类基因组跑出来只有5%那基本可以断定库没选对或参数有问题。4.2 .out文件里最容易看错的几列.out文件是制表符分隔的明细表每一行代表一个重复序列片段。列的数量不少但最核心的几个含义如下列含义SW scoreSmith-Waterman对比得分越高表示匹配越强perc div该拷贝与consensus之间的突变率divergenceperc 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位置信息中标出Ccomplementary分析时要注意区分正负链。负链的重复经常会被忽略如果再往下游提取序列没考虑方向会提取出反向互补序列导致后续分析全部跑偏。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 namesfasta里存在重复序列名用seqkit seq -n检查再用seqkit rename或replace改名RepeatModeler在RECON阶段内存溢出线程开太多基因组过大降低-pa限制单批次内存给作业脚本加内存上限如-XmxRepeatModeler跑完了但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分布、抽检具体区域。注释结果的好坏不只影响一篇论文更会沉淀成后续所有分析的基础资源。把这些步骤走扎实后面做基因家族分析、比较基因组、群体遗传的时候你会省下大量返工的时间。