拿到一个刚组装的基因组第一步不是急着做基因注释而是先把里面的重复序列全部“揪出来”。这个问题我吃过亏早期做某个非模式物种基因组时直接拿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来发现重复单元速度快适合发现高频重复。TRFTandem 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.pyPython脚本来读取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这几个常用模块。特别是执行时如果出现Cant locate JSON/PP.pm in INC这类错误说明Perl模块缺失。另外重要的一点不要把RepeatModeler/RepeatMasker跟其他版本敏感的比对软件放在同一个环境中使用。比如系统里如果装有另一个版本的rmblast或者ncbi-blastPATH顺序会导致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 构建数据库并启动RepeatModelerRepeatModeler的第一步是把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.libUnknown比例如果超过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默认会把结果输出到输入文件所在目录建议用这个参数保持输出文件集中整洁。跑完之后输出目录里会有以下几类文件这是解读结果的核心文件名内容genome.clean.fa.masked屏蔽后的基因组序列小写部分就是重复序列genome.clean.fa.out表格形式的详细重复注释一个区间一行是主要的分析数据源genome.clean.fa.tbl汇总统计表直接给出总屏蔽率、各类重复占比等关键数字genome.clean.fa.out.gffGFF3格式版本的注释文件方便导入IGV等基因组浏览器genome.clean.fa.polyout低复杂度序列注释一般在解释结果时不用重点看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之后这个问题就不再出现。所以如果你也遇到长时间卡在某一步先做两件事用top查看是否有某个子进程还在消耗CPU如果完全是0%且长时间不变基本可以断定是“死锁”而不是“在计算”。检查序列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过滤时都能反复利用相当于一次投入、多处收益。所以第一次花时间把库做好、做全绝对不亏。