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 python=3.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 2>&1几个关键参数逐个拆解:
-database:BuildDatabase 生成的数据库名前缀。这个名字不要起太长,不要包含特殊字符,后面 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 | 统计摘要表,按重复类别汇总拷贝数、长度占比 | 快速汇报注释概况 |
.gff | GFF3 格式注释文件 | 下游工具的标准输入 |
.masked | 软遮蔽/硬遮蔽后的基因组序列 | 基因预测、变异检测 |
.divergence | 每个家族的 Kimura 分歧度分布(如果指定-a参数) | 转座子爆发历史分析 |
.cat.gz | 每个位点的详细比对信息和分类证据 | 深入核查特定区域 |
.out文件的每一行对应一个重复拷贝实例,它的前 5 列含义固定:序列名、起始位置、结束位置、重复方向(+/C,C 表示互补链上的拷贝)、重复家族名。位置是基于 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 2>&1结束前检查日志里的 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 'NR>3{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-name=repeat_annot #SBATCH --nodes=1 #SBATCH --ntasks=32 #SBATCH --mem=128G #SBATCH --time=7-00:00:00 #SBATCH --partition=compute 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 2>&1 # 运行 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这个脚本把两个阶段衔接起来,中间不需要人工干预。注意--time=7-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 "chr"NR"\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 把基因误判成了重复序列,这种模型要果断从库里剔除。
重复序列注释这件事,看似只是整个基因组项目的一个流水线环节,但它的质量会一路传导到基因预测、系统发育、群体遗传、表观遗传所有下游分析。一次认真做的重复注释,能为项目省下数周甚至数月的返工时间。希望这篇文章能帮你把这条流程走通、走顺。