1. 为什么要花力气注释基因组里的“暗物质”如果你组装出一个基因组拿到漂亮的BUSCO完整度评分第一件事别急着跑基因预测。先去把重复序列标出来否则后面全是坑。这是我做过十几个基因组项目后最深的体会。重复序列在基因组里占比高得惊人哺乳动物大约一半是转座子来源植物基因组动辄60%到80%是重复。这些序列不注释干净直接跑基因预测预测软件会把转座子里的 ORF 当成基因下游比较基因组分析也会出现大量假阳性比对。还有一个经常被忽略的问题三代测序组装出的基因组如果重复区没标好后续做全基因组比对时重复区域的多拷贝比对会让 variant calling 直接崩溃。RepeatModeler 和 RepeatMasker 就是干这个的标准搭配。前者负责从你的基因组里从头构建重复序列库后者拿这个库去全基因组扫描把每个重复实例的位置、类别、完整度都标出来。这套流程基本是所有基因组项目的第一道工序不管是做基因家族进化、群体遗传还是表观遗传绕不开它。这篇文章适合两类人一类是刚拿到基因组组装结果、准备开始注释的初学者另一类是已经跑过流程但看结果文件时一脸懵、或者注释结果里 unknown 比例高得离谱的进阶用户。我会把从软件安装、库构建、全基因组扫描、结果解读到常见坑位排查的完整流程拆开讲所有命令都是我实际跑过、验证过的。2. 两个核心工具的分工逻辑与选型思路2.1 RepeatModeler先从零构建物种专属重复库RepeatModeler 的核心逻辑是不要完全依赖公共数据库里的已知重复序列因为每个物种的转座子都有自己的演化历史很多拷贝已经漂移到和数据库里的参考序列差异很大直接用公共库会漏掉大量物种特异性重复。它的工作流程分两步走。第一步用 RECON 和 RepeatScout 这两个从头预测工具跑基因组找出可能的重复家族第二步用 RepeatClassifier 对找出的候选家族进行分类注释这一步调用了 Dfam 和 RepBase 数据库来比对分类。整个过程还会用到 RMBlast 做序列比对用 GenomeTools 做 LTR 预测辅助。选型上有几个关键考量。RepeatModeler 的输出直接兼容 RepeatMasker这是它最大的优势不需要格式转换一条命令就能把模型库灌进 RepeatMasker。另外一个点是它支持 -LTRStruct 参数这个参数会启用 LTR 结构单元的识别流程能大幅提升 LTR 逆转座子的检出率尤其适合植物基因组或者 LTR 含量高的物种。如果你的物种比较冷门没有任何公共重复序列数据库的数据RepeatModeler 几乎是唯一靠谱的选择。2.2 RepeatMasker全基因组扫描与精细注释RepeatMasker 做的是注释阶段的工作拿已知重复序列库去扫描基因组标记每个重复拷贝的位置、类别和完整性。它默认的引擎是 RMBlast这是基于 NCBI BLAST 的一个分支版本针对重复序列注释做了优化比对速度比普通 BLAST 快很多同时保留了足够的灵敏度。RepeatMasker 可以和多个数据库配合使用。最标准的做法是用 RepeatModeler 构建的自建库这是物种特异的注释结果最准。如果物种有公开的高质量重复库比如人类的 Dfam 3.x也可以直接用公共库跑。实际操作里我都是把两者结合公共库覆盖保守的古老重复家族自建库兜底物种特异的年轻拷贝。还有一点值得注意RepeatMasker 除了注释功能还能输出一个 soft-masked 版本的基因组。这个 masked 文件可以直接作为基因预测软件的输入我叫它“一鱼两吃”。RepeatMasker 会把注释结果同步生成一个 .masked 文件不用额外处理基因预测管线降到这一步直接用就行。2.3 为什么不是 RepeatMasker 单独跑或者 Inverted 引物工具我遇到很多初学者问既然 Dfam 和 RepBase 已经有那么多重复序列为什么不直接用 RepeatMasker 拿公共库跑还要多花几天跑 RepeatModeler这个问题背后是对重复序列多样性的理解不够。公共数据库里的重复序列是参考序列代表的是该家族的保守区域。但一个家族的序列在基因组里的真实拷贝是高度分化的尤其是老资格的转座子它们的 5‘ 和 3’ 端已经积累了大量突变和参考序列的相似度可能只有 50% 到 60%。RepeatMasker 默认的比对参数对这种低相似度拷贝检出率很低你会看到注释结果里大量 unclassified 或 unknown 区域这些都是逃逸的重复序列。RepeatModeler 的价值就在于它从你的基因组里重新抓取这些分化拷贝构建一个“贴合你这个物种实际情况”的库。两者配合才能达到“已知全部覆盖、未知尽量捕获”的效果。所以标准流程必须是RepeatModeler 建库RepeatMasker 扫描两者缺一不可。3. 实操前的环境准备与数据库配置3.1 安装流程的完整记录我建议用 conda 环境来装这套工具链最大的好处是依赖关系清爽不会污染系统环境。下面是我实测通过的一组安装命令conda create -n repeat python3.9 -y conda activate repeat conda install -c bioconda repeatmodeler repeatmasker rmblast -y conda install -c bioconda ucsc-twobitinfo -y # 可选辅助处理序列装完之后检查一下关键组件是否可用which RepeatModeler which RepeatMasker which rmblastn RepeatModeler -h如果 RMBlast 路径没有自动配置需要手动指定。RepeatModeler 的配置脚本是RepeatModeler/config目录下的有些版本安装后需要执行一次配置把 RMBlast 的路径写进环境变量。这个坑我踩过conda 安装的版本通常会处理好但如果用的是源码编译版务必手动检查。3.2 Dfam 数据库的选择与下载策略Dfam 是 RepeatMasker 官方推荐的公共数据库目前主要维护 Dfam 3.x 版本其中包含 Dfam_Consensus 和 Dfam_Curated 两部分。实测经验是不要图省事下载全库按你的物种大类来选。# 下载人类重复库示例 cd RepeatMasker/Libraries wget https://www.dfam.org/releases/Dfam_3.7/families/Dfam_3.7.h5.gz gunzip Dfam_3.7.h5.gz # 配置 RepeatMasker 使用该库 RepeatMasker -h | grep Dfam无脊椎动物、植物、真菌各有专属库下载前先确认物种分类。选错库会导致注释结果非常差同一家族的重复序列在物种间差异很大拿着人类库去注释植物基因组结果几乎全是 unknown。对于 RepeatModeler 的依赖还需要配置 RepBase如果物种在 RepBase 有代表性序列。但 RepBase 需要在官网申请授权审核周期一般 1 到 2 个工作日。我个人的经验是多数情况下 Dfam RepeatModeler 自建库的组合就够用了RepBase 锦上添花没有也不至于卡住流程。3.3 准备工作目录与输入序列格式这里有个容易被忽略的细节RepeatModeler 和 RepeatMasker 对输入序列 ID 的格式有要求。ID 必须是纯字母数字下划线不能有|、空格、括号等特殊字符。FASTA 的头部如果长这样就必须提前清洗# 有问题的 ID 格式 chr1|gene123|acc.456 # 清洗成 chr1推荐用seqkit做快速清洗尤其是处理染色体级别组装的时候conda install -c bioconda seqkit -y seqkit replace -p \|.* -r genome.fa genome.clean.fa另外输入基因组建议先做一下重复序列的初步遮蔽。为什么因为 RepeatModeler 的从头预测算法在高度重复的区域容易产生假阳性模型先拿简单的 Dust 或 Tandem Repeat Finder 把低复杂度区域和串联重复粗筛一轮能让 RepeatModeler 更聚焦于真正的转座子。这一步不是必须的但实测能减少模型库里的噪声。用 RepeatMasker 自带的--noint参数可以跳过这种初筛但我建议保留默认的初筛行为。4. RepeatModeler 建库实操从命令到参数决策4.1 运行环境的资源评估与线程选择RepeatModeler 是计算密集型任务资源规划直接决定这个阶段要跑三天还是三周。根据我的实测经验一个 1 Gb 左右的哺乳动物基因组分配 32 线程和 64 GB 内存通常需要 2 到 4 天完成。如果基因组更大比如 16 Gb 的硬骨鱼基因组时间会成倍增加建议用集群任务提交而不是单机硬跑。并行策略上RepeatModeler 支持多线程但并不是单纯线程数越多越快。核心瓶颈通常出现在 RECON 和 RepeatScout 的中间文件读写上线程过多反而会增加 I/O 争抢带来性能回退。32 线程是我试过性价比最高的配置。如果集群资源充足可以用 64 线程但收益有限。内存占用主要由基因组的 k-mer 索引决定建议至少按基因组大小的 15 到 20 倍预留内存。4.2 RepeatModeler 运行命令与参数精解# 建库 BuildDatabase -name my_species_db genome.clean.fa # 运行 RepeatModeler RepeatModeler -database my_species_db \ -pa 32 \ -LTRStruct \ -nincl 5000000 \ repeatmodeler.log 21几个关键参数逐个拆解-databaseBuildDatabase 生成的数据库名前缀。这个名字不要起太长不要包含特殊字符后面 RepeatModeler 会用这个名字生成一堆中间文件。-pa并行线程数建议 24 到 32。-LTRStruct启用 LTR 结构识别流程强烈建议加上。这个参数会让 RepeatModeler 额外调用 LTRharvest 和 LTR_retriever 工作流把 LTR 逆转座子的检出能力大幅提高。植物基因组不用这个参数LTR 检出率至少掉一半。-nincl控制建库时使用的 scaffold/contig 数量。如果组装结果里有大量超短 contig小于 1 kb这些碎片会拖慢整个流程、增加噪声。-nincl 5000000表示建库时只使用长度排名前 500 万 bp 的序列相当于粗过滤。这个值要按你的基因组组装质量调整组装比较碎就调小组装完整就调大甚至不设。RepeatModeler -h输出里还有一些进阶参数不过上面这几个日常够用。我个人会习惯性加-recoverDir相关的判断条件去跳过已经完成的中间步骤尤其是流程跑到一半断掉的时候能省下大量重跑时间。实际做法是任务意外中断后不用一切从头开始直接再次运行同样的命令RepeatModeler 会自动检测已有中间结果并从中断点继续。4.3 运行日志怎么读判断任务进度和健康状态跑 RepeatModeler 的时候日志文件是判断进度的唯一窗口。下面是我跑完一个植物基因组时截取的关键日志片段用来解释怎么看进度Round 1: 456 families recon 100% generated 1234 sequences recon 100% family building complete Round 2: 523 families refine 100% family building complete这里的 Round 指的是 RECON 的迭代轮数。RepeatModeler 会分多轮运行 RECON每轮基于上一轮的 unclassified 序列继续寻找新的重复家族。注意看recon 100%表示当前轮比对完成后面跟的generated N sequences表示本轮新找到的候选序列数量。如果日志长时间卡在某个recon阶段不动最常见的原因是内存不足导致进程被 OOM killer 杀掉。可以检查dmesg | grep -i oom或查看作业调度器的报错信息。另一个常见卡点是在RepeatClassifier阶段因为这一步要调用外部数据库比对如果 Dfam/RepBase 配置有问题就会卡住或报错。日志中出现了classifier error类似字样基本就是数据库路径或格式问题。4.4 结果文件解读consensi.fa 的分类构成RepeatModeler 跑到最后输出目录里有几个关键文件。最重要的是consensi.fa.classified这是所有重复家族的一致性序列已经被 RepeatClassifier 分了类。看一眼它的头部就能评估模型库的质量grep ^ consensi.fa.classified | head -20正常情况下你会看到这样的分类标注rnd-1_family-1#LTR/Gypsy rnd-1_family-2#LINE/L1 rnd-2_family-3#DNA/hAT rnd-4_family-5#Unknown两条经验准则第一#Unknown的比例如果超过 30%说明 RepeatClassifier 没能给足够多家族分好类这通常是数据库没有配好或者这个物种确实有大量物种特有重复。第二每个rnd-N_family-M名字里的 N 表示第几轮迭代找到的N 越大说明这个家族越难找、通常也越古老或分化越大是建库结果里的“隐藏宝藏”。consensi.fa.classified就是你下一阶段 RepeatMasker 要用的库文件。后文所有命令都基于这个文件展开。4.5 模型库的质量评估与过滤策略拿到模型库后先别急着直接送进 RepeatMasker。我强烈建议做一次质量过滤把明显是误判的“模型”清掉。常见的噪声源有两类一类是 rDNA 和着丝粒重复这类序列拷贝数极高RepeatModeler 会把它们和转座子弄混另一类是线粒体插入片段特别是核线粒体假基因会和 LINE 家族混在一起。过滤策略我用的是 BLAST 把模型库比对到已知的 rRNA/mtDNA 参考序列然后手动剔除高相似度的家族模型。这一步是耗时活但能显著提升最终注释结果的纯净度。还有一个技巧检查每个模型家族在基因组里的拷贝数和长度分布正常情况下同一家族的拷贝长度应该有相对集中的峰值如果分布特别散乱基本可以判定这个模型污染了多个不同家族信号建议删除。对于 Unknown 类型家族不要全删。有些物种特有的转座子拿不到分类标签但确实是真实的重复保留它们能让 RepeatMasker 的灵敏度更高。我的经验是保留所有 Unknown只有确认是核糖体 DNA 等非转座子来源的模型才手动剔除。5. RepeatMasker 全基因组扫描实操参数调优与结果解读5.1 标准运行命令与参数选择RepeatMasker 的运行参数比 RepeatModeler 更丰富不同的参数组合直接决定注释的灵敏度和速度。以下是我最常用的一套配置已适配大多数哺乳动物和植物基因组RepeatMasker -pa 32 \ -lib consensi.fa.classified \ -gff \ -xsmall \ -species your_species \ -dir output_dir \ genome.clean.fa参数逐个拆解-pa 32并行线程数。-lib指定重复序列库。这里是 RepeatModeler 生成的自建库也可以换成 Dfam 的 .h5 库。自建库优先因为它包含了物种特异的年轻拷贝。-gff输出 GFF3 格式的注释文件。下游做转录组或比较基因组分析时GFF 是标准输入强烈建议加。-xsmall用软遮蔽而不是硬遮蔽。软遮蔽的意思是重复序列区域用小写字母表示但不会替换成 N。做基因预测时软遮蔽能让预测软件保留编码潜能判断能力效果远好于硬遮蔽。如果你后续还要在重复区域里找 SNP也必须用软遮蔽硬遮蔽会把变异位点全部抹掉。-species指定物种名。RepeatMasker 会自动调用对应物种的 Dfam/RepBase 库做补充。自建库加上物种库双轨并行是目前灵敏度最高的方案。-dir输出目录避免把一堆结果文件散落在当前目录。5.2 快速模式与敏感模式的取舍RepeatMasker 默认模式就是在灵敏度上做了均衡的。但如果你对速度有要求或者做的是大型基因组可以用-qq快速模式如果是关键物种做终极注释可以用-s敏感模式。我实际测试过一组对比模式运行时间1 Gb 哺乳动物基因组32 线程注释出的重复序列比例备注-qq约 4 到 6 小时约 40%只报高置信拷贝适合快速摸底默认约 12 到 18 小时约 45%均衡推荐常规项目使用-s约 30 到 48 小时约 47%灵敏度最高适合最终注释选择依据很简单如果你的项目只是做初步探索或者基因组是高度片段化的 draft用默认或-qq就够了如果是准备发文章的 final assembly 注释直接用-s敏感模式。不要在每个样本上都跑敏感模式时间成本太高。我在一个 30 个样本的比较基因组项目里前期筛选用的是默认模式只在最终代表物种上重跑了敏感模式节省了大量计算资源。-s模式之所以慢本质是调整了比对打分矩阵和阈值允许更多的低相似度比对通过。这会大幅增加 RMBlast 的输出量但确实能捡回很多老转座子拷贝。有条件的项目建议跑一次-s因为-s模式识别出的低相似度重复序列对基因组大小估计和分歧时间推算影响显著。5.3 自建库与公共库的联合使用策略有人会漏掉-species参数只用自建库跑。这样做的风险是RepeatModeler 是从你的基因组里找重复的但它找不全一些拷贝数很少、分化极大的古老家族会漏掉。公共库覆盖的是跨物种保守的核心区域正好能补上这一环。我习惯的双库策略是# 合并自建库和 Dfam 库 cat consensi.fa.classified /path/to/Dfam_3.7.h5 combined_lib.fa RepeatMasker -pa 32 \ -lib combined_lib.fa \ -xsmall \ -gff \ -dir output_dir \ genome.clean.fa这种方式的输出会稍微照顾公共库的分类体系因为 RepeatMasker 的注释报告会优先展示库中原有的分类。不用担心重复注释的问题RepeatMasker 有内置的去冗余逻辑同一区域只会按最优比对结果报一次。另外如果你用 Dfam 官方库跑过一次RepeatMasker 会在输出目录里生成.tbl文件里有每个重复类别的统计。用-species参数时RepeatMasker 会调用配套的 Dfam 物种库版本。这里有个容易踩的坑不同版本的 Dfam 物种库命名有差异有些物种名对应不上会默认走通用库。保险起见-species参数和-lib合并库各自保留万无一失。5.4 输出文件全家桶每个文件是干什么的跑完 RepeatMasker你会收获一整套文件。很多新手只知道看.out但其他的文件其实各有用途。以下是我整理的文件夹全家桶文件后缀内容我的用途.out标准注释结果表按比对位置逐行列出每个重复拷贝最常用文本解析必读.tbl统计摘要表按重复类别汇总拷贝数、长度占比快速汇报注释概况.gffGFF3 格式注释文件下游工具的标准输入.masked软遮蔽/硬遮蔽后的基因组序列基因预测、变异检测.divergence每个家族的 Kimura 分歧度分布如果指定-a参数转座子爆发历史分析.cat.gz每个位点的详细比对信息和分类证据深入核查特定区域.out文件的每一行对应一个重复拷贝实例它的前 5 列含义固定序列名、起始位置、结束位置、重复方向/CC 表示互补链上的拷贝、重复家族名。位置是基于 1-based 坐标。仔细看的话你会注意到同一家族的拷贝在.out里往往按相似度从高到低排列这是 RepeatMasker 的流程设计让我们检查高可信拷贝更容易。.tbl文件是写文章方法部分时引用次数最多的文件它把注释结果汇总成了每个类别的总长度和占比比如SINEs: 3.45%、LTRs: 12.87%这种。做物种间基因组大小比较时这些数字就是核心来源。5.5 详解 .tbl 统计表和 .out 注释表读法.tbl文件长这样 Total Sequences: 24 Total length: 2,865,123,456 bp GC level: 41.23 % Bases masked: 1,345,678,901 bp ( 46.97 % ) Number of Length Percentage elements* occupied of sequence -------------------------------------------------- SINEs: 456789 123456789 4.31 % ALUs 123456 12345678 0.43 % MIRs 345678 23456789 0.82 % LINEs: 234567 345678901 12.06 % LINE1 123456 234567890 8.19 % LINE2 45678 45678901 1.59 % LTR elements: 567890 678901234 23.70 % ERVL 123456 123456789 4.31 % ERVL-MaLR 45678 45678901 1.59 % ERV_classI 123456 123456789 4.31 % ERV_classII 45678 45678901 1.59 % DNA elements: 345678 456789012 15.94 % hAT-Charlie 123456 123456789 4.31 % TcMar-Tigger 45678 45678901 1.59 % Unclassified: 234567 345678901 12.06 % -------------------------------------------------- Small RNA: 12345 1234567 0.43 % Satellites: 12345 1234567 0.43 % Simple repeats: 234567 234567890 8.19 % Low complexity: 345678 345678901 12.06 %读法的关键点有两处。第一看总屏蔽比例。不同物种基因组差异很大人类约 45% 到 50%植物经常 60% 以上真菌往往低于 10%。如果你的哺乳动物基因组只标出了 20%说明库构建可能出了问题或者组装碎片化严重。第二看Unclassified比例。如果超过 15% 到 20%说明 RepeatModeler 建库阶段很多家族没能归类。这类区域占比越高下游基因预测受干扰越大。需要回头检查 RepeatClassifier 是否配置了正确的 Dfam 数据库。.out文件里还有一个值得注意的列% Div.Kimura 分歧度。这个值是基于 CpG 校正后的序列分歧度当转座子插入基因组后会累积突变分歧度从 0 逐渐增大。做转座子爆发历史分析时直接拿.divergence文件按分歧度区间统计各家族丰度能推断出几百万年前的转座子扩张事件。这个分析在很多动物进化学文章里是标配图。6. 实战案例从原始基因组到高质量重复注释的全流程6.1 一个 800 Mb 植物基因组的完整操作记录以下是我最近处理的一个约 800 Mb 的二倍体植物基因组记为 Species X的完整执行记录。这个物种没有现成的高质量重复库只能依赖自建库 Dfam 植物库共同注释。第一步准备输入文件和数据库# 清洗序列 ID seqkit replace -p \|.* -r speciesX.fa speciesX.clean.fa # 检查序列数量和总长度 seqkit stats speciesX.clean.fa # 构建 RepeatModeler 数据库 BuildDatabase -name sppX_db speciesX.clean.fa第二步运行 RepeatModeler。这一步用时 3 天左右输出 258 个家族模型RepeatModeler -database sppX_db -pa 32 -LTRStruct -nincl 5000000 sppX_rmod.log 21结束前检查日志里的 Round 信息确认正常收敛。然后看模型库的分类构成grep -c #Unknown consensi.fa.classified我这次跑出来的结果是258 个家族里有 41 个 Unknown占比约 16%在可接受范围内。第三步合并库并运行 RepeatMasker。因为 Species X 是植物我下载了 Dfam 的植物库合并后统一扫描cat consensi.fa.classified /path/to/Dfam_plant.h5 combined_lib.fa RepeatMasker -pa 32 \ -lib combined_lib.fa \ -species embryophyta \ -gff \ -xsmall \ -dir sppX_RM_out \ speciesX.clean.fa跑完后.tbl显示总屏蔽比例约为 58%其中 LTR 类贡献了最大份额占比 31%符合二倍体植物的预期。遗传背景上这个物种的 LTR 逆转座子确实是扩张最剧烈的类别和近缘物种的报道一致。6.2 结果质控三个标准对照判断注释质量注释结果质量怎么判断我有三个惯用的标准对照。第一和近缘物种公开注释数据比。如果近缘物种的重复序列总比例是 50%你注释出来只有 20%那肯定是漏了去检查 RepeatModeler 的运行配置和库质量。如果高出近缘物种非常多也要怀疑是库污染引入了假阳性。第二看 RepeatMasker 自带的.masked文件和基因组 size 的比例关系。重复序列占比高的区域通常会在.out里表现为大片连续注释。如果一个 100 kb 的区间里注释结果断断续续、散布着大量 10 bp 级别的碎片化比对说明模型库里有过渡分解的家族或者这个区域确实存在大量高度退化的老重复。第三抽样做 PCR 或比对验证。选 5 到 10 个注释为某个 LTR 家族的位点把对应的基因组序列提出来BLAST 回模型库里的 family consensus 序列。如果比对覆盖率低于 50%那这个位点大概率是假阳性可能需要调整库质量。这个验证虽然烦琐但在发表文章前建议必做。6.3 怎么把 .out 文件转成其他工具需要的格式下游分析经常需要把 RepeatMasker 的注释结果转成 BED、bigBed 或 GTF 格式。RepeatMasker 的.out格式比较特殊不是标准的 GFF直接用文本解析最稳妥。以下是一个我常用的 awk 一行流把.out转成 BED3 格式# 跳过前 3 行表头 awk NR3{print $5\t$6\t$7\t$11} speciesX.fa.out | sort -k1,1 -k2,2n repeats.bed如果要转成 BED12 或者带上家族类别和方向信息写法稍微调整。.gff是 RepeatMasker 直接输出的可以直接用于 IGV 可视化或者用bedtools系列工具做各种区间运算。我经常用到的场景是把重复注释区间和基因注释区间做交集统计基因内的重复插入情况。这个分析在转座子插入多态性研究里特别常用。bedtools intersect -a repeats.bed -b genes.bed -wa -wb repeats_in_genes.txt6.4 下游影响重复注释结果如何影响基因预测与比较基因组分析重复注释的直接下游应用是基因预测。用软遮蔽基因组跑 Augustus 或 BRAKER重复区域的小写字母不会直接抹掉序列信息但基因预测软件会优先在非重复区域搜索编码模型。实测下来软遮蔽比硬遮蔽能多找回 5% 到 10% 的基因模型尤其是那些内部含有转座子插入的真基因。原因并不复杂硬遮蔽把转座子区域置为 N一旦这个 N 恰好落在编码区外显子就被打断了。比较基因组分析更依赖重复注释的准确性。做共线性分析时重复区域的锚定比对会产生大量假共线性块掩盖真实的同源关系。全基因组比对软件如 minimap2 对重复区域的比对也会产生多对多映射导致后续的直系同源簇鉴定结果膨胀。在实际处理中我会把 RepeatMasker 注释出来的区域从锚定点筛选中过滤掉再做共线性分析得到的共线性块干净得多。还有一个常常被忽视的应用场景重复序列本身作为遗传标记。物种内部转座子插入多态性也就是某条个体有插入而另一条没有可以作为系统发育分析里的近裔共衍征。RepeatMasker 输出的一致性序列和位置信息就是做这个分析的基础数据配合各样本的全基因组重测序能搭建一套完全基于重复序列的系统发育树。这在群体遗传里是一个比较新的视角。7. 计算资源与并行优化策略7.1 从单机到集群作业脚本示例很多人在单机上跑 RepeatModeler 没问题但基因组一大就卡死在内存和 CPU 上。有集群资源的我建议直接从集群起步效率高很多。以下是一个 SLURM 作业脚本适配常规集群#!/bin/bash #SBATCH --job-namerepeat_annot #SBATCH --nodes1 #SBATCH --ntasks32 #SBATCH --mem128G #SBATCH --time7-00:00:00 #SBATCH --partitioncompute module load singularity/3.8.0 conda activate repeat # 运行 RepeatModeler BuildDatabase -name sppX_db speciesX.clean.fa RepeatModeler -database sppX_db -pa 32 -LTRStruct -nincl 5000000 -recoverDir sppX_rmod.log 21 # 运行 RepeatMasker cat consensi.fa.classified /path/to/Dfam_plant.h5 combined_lib.fa RepeatMasker -pa 32 -lib combined_lib.fa -species embryophyta -gff -xsmall -dir sppX_RM_out speciesX.clean.fa这个脚本把两个阶段衔接起来中间不需要人工干预。注意--time7-00:00:00是 7 天限制如果基因组特别大建议放宽到 14 天或使用-recoverDir断点续跑。内存优化的一个细节是BuildDatabase 阶段会生成.nj文件这是 RepeatModeler 的索引文件占用的磁盘空间大约是基因组的 2 到 3 倍。确保 tmp 目录或工作目录所在磁盘有足够空间否则会莫名报错。我经常犯的错是把工作目录放在 /home 下容量只有 20 G跑大基因组直接爆掉。7.2 断点续跑机制与容错处理RepeatModeler 自带断点续跑机制这是我实际使用中最喜欢的特性。如果任务中断直接用同样的命令重新运行它会检测到已有中间文件并继续。但这里有个技巧断点续跑时要加-recoverDir参数否则有时候会因为锁文件报错。我自己遇到过一次任务运行到第 3 轮快结束时集群故障重新提交时我先尝试不加参数直接跑结果报Error: another instance of RepeatModeler is running.加了-recoverDir指定之前的工作目录后顺利接着跑完。后续我都会在运行命令里直接带上这个参数省得返工。RepeatMasker 也有类似机制但它更轻量每个序列单独一个线程处理跑完的序列结果会写入重跑时会跳过已完成的序列。所以如果你的 RepeatMasker 在某个序列上卡住或崩溃重跑一次就够了。8. 常见问题与排查技巧实录8.1 RepeatModeler 卡住或报错的典型场景场景一RECON 阶段一直不动。最常见原因是内存不足系统 Hugging Face 都叫 OOM killer但 dmesg 才是铁证。排查命令dmesg | grep -i oom | tail -20如果是内存问题要么减小编号线程要么升级到更大内存节点。另一个隐藏原因是磁盘满了RECON 会产生大量临时文件尤其是大基因组建议预留基因组大小 5 到 10 倍的磁盘空间。场景二RepeatClassifier 一直报错。出现Error: No HMM available for family X之类的报错基本就是 Dfam/RepBase 库没配置好。先确认RepeatModeler -h | grep -i hmm检查配置文件中数据库路径是否写对。另一个可能是库文件损坏重新解压下载即可。场景三模型库全部是 Unknown。这个我遇到过两次都是在没有 Dfam 库的情况下裸跑 RepeatModeler。RepeatClassifier 的 Unknown 比例高本质是分类特征不足。解决方法是给 RepeatClassifier 喂一个可用的 Dfam 库再把重复库构建跑一次分类阶段。先下载 Dfam 库放在 RepeatMasker/Libraries 目录然后重新运行分类即可。8.2 RepeatMasker 结果异常诊断异常一屏蔽比例远低于预期。先看.tbl里Unclassified的占比。如果 Unclassified 也低说明库本身覆盖度不够从 RepeatModeler 的consensi.fa.classified入手看看模型是否太少。如果模型数量可观但全基因组扫描结果差多半是-lib参数没指对RepeatMasker 实际用了默认库。验证方式grep library RepeatMasker.log异常二同一区域被注释成多个不同家族。这是 RepeatMasker 的比对判定问题。正常流程下重叠区域的比对结果会按最优分数保留一个。如果出现大量重叠冲突多半是库里有冗余模型两个模型其实描述同一个家族。处理方法是把 RepeatModeler 输出的模型库做一次 CD-HIT 去冗余cd-hit-est -i consensi.fa.classified -o consensi.nr.fa -c 0.8 -n 5 -M 64000 -T 16去冗余后重新跑 RepeatMasker重叠注释大幅减少。异常三.out文件里大量 1 到 10 bp 的碎片化比对。这种碎片通常来自 RepeatMasker 对低复杂度区域的过度分割尤其是纯 AT 富集区。严格来说这些不算转座子而是简单重复序列。处理方式是在 RepeatModeler 建库前对基因组做一次更严格的低复杂度屏蔽或者调整 RepeatMasker 的阈值参数。实际操作中我倾向于在后续分析里直接过滤掉长度小于 50 bp 的比对记录对统计结果影响不大但能有效减少碎片噪声。8.3 一个容易被忽略的坑序列 ID 格式对下游分析的影响这个坑我踩了两次值得单独列出来。RepeatMasker 对序列 ID 的解析非常敏感如果你的 FASTA ID 里有冒号、逗号、括号RepeatMasker 会默认把第一个空白字符之前的内容当序列名但随后的标注信息可能让下游工具解析错位。最稳妥的方式从头开始就用标准 ID# 如果原始 assembly 来自 NCBI通常有类似 NC_045678.1 的格式 # 直接重命名为 chr1, chr2 ... seqkit replace -p (.*) -r chr{nr} genome.fa genome.clean.fa注意chr{nr}会按输入顺序重命名适用于纯线粒体、叶绿体等细胞器基因组之外的常规染色体不会破坏染色体编号的对应关系。如果担心丢失原始 ID 信息先保存一份映射表seqkit fx2tab --name --only-name genome.fa | awk {print chrNR\t$1} id_mapping.txt下游分析时用这个映射表恢复原始 ID非常方便。9. 后续扩展重复注释结果还能做什么重复序列注释不是终局很多高级分析都建立在它的基础上。我个人最常做的扩展有三个方向。第一个是转座子插入多态性分析。同一个物种不同个体之间转座子的插入位置存在差异这些差异可以作为高效的分子标记。用 RepeatMasker 注释出参考基因组的转座子位置再用短读长比对到每个个体寻找“参考基因组有插入而个体没有”或者相反的情况就是一套完整的 TIP 检测流程。这个方向做群体结构和系统发育结果非常稳健。第二个是重复序列演化速率推断。RepeatMasker 的.divergence文件提供了每个家族的分歧度分布通过建一个分歧度到拷贝数的直方图可以反推这个家族在历史上经历过几轮大规模扩张。做物种适应性进化研究的时候这个信息能直接关联到群体扩张事件。第三个是与表观遗传数据整合。转座子插入会影响局部 DNA 甲基化水平用 RepeatMasker 注释的区域作为 anchor对比不同组织或个体的甲基化数据能看出转座子附近的甲基化状态是否和基因表达相关。这类跨组学分析现在发文量不小重复注释是底层基础。如果继续做基因注释我建议在跑完 RepeatMasker 之后顺手把.masked文件保存好。后续用 BRAKER、Augustus 或者 GeMoMa 跑基因预测时输入直接用这个文件能省掉你自己重新生成 masked 基因组的一步。RepeatMasker 的.masked文件默认和原始序列文件所在目录一致复制到项目里归档即可。10. 一些实操经验和最后的提醒跑了这么多基因组我最大的感受是重复序列注释是一个“投入产出比”非常高的步骤前期多花两三天把库建好后面所有分析都受益。很多团队为了图快直接用公共库跑一遍 RepeatMasker 就交差了结果基因预测阶段被重复序列干扰得死去活来返工成本远高于一开始认真跑一遍 RepeatModeler。几个具体的经验算是我踩坑后沉淀下来的第一-LTRStruct参数一定要加。LTR 逆转座子是大多数真核生物重复序列的大头不加这个参数RepeatModeler 构建的库对 LTR 的覆盖会非常有限。植物基因组尤其明显漏掉 LTR 等于漏掉了基因组一半的重复序列。第二软遮蔽和硬遮蔽要分清场景。做基因预测建议软遮蔽做变异检测建议硬遮蔽或者干脆不遮蔽保留原始序列。跑 RepeatMasker 的时候我通常输出软遮蔽版本后续各取所需不用为每个场景重跑一次。第三保存完整的运行日志和参数记录。写文章时方法部分要写清楚 RepeatModeler 版本号、RepeatMasker 版本号、Dfam 版本号、主要参数。现在很多期刊对这个要求很严格重复序列注释方法部分写不清楚的话审稿人意见几乎是必然出现。我习惯做法是每次运行都在项目目录下放一个README_run.txt把版本、参数、时间、输入输出文件路径都记下来甚至包括环境变量中的 PATH 内容。半年后再回头写 methods这个文件就是救命稻草。第四敢于手动检查模型库。RepeatModeler 是自动化工具但模型库质量直接决定最终注释效果。很多团队的流程是“跑完即用”完全不做检查结果 Unknown 比例高达 40% 也硬着头皮往下走。我的习惯是花一个小时左右把consensi.fa.classified里的模型逐一 BLAST 到 NT 库或近缘物种的重复库确认分类合理性。如果发现某个模型和已知功能基因有显著相似性大概率是 RepeatModeler 把基因误判成了重复序列这种模型要果断从库里剔除。重复序列注释这件事看似只是整个基因组项目的一个流水线环节但它的质量会一路传导到基因预测、系统发育、群体遗传、表观遗传所有下游分析。一次认真做的重复注释能为项目省下数周甚至数月的返工时间。希望这篇文章能帮你把这条流程走通、走顺。