项目目标

将 RepeatMasker 的 GFF3 注释文件转换为 RepeatMasker 标准输出格式（.out / .cat），以便使用 ProcessRepeats 工具生成 .tbl 汇总统计表。

项目结果

项目无法继续

最终选择：寻求公司使用原始的 RepeatMasker 运行环境和中间文件重新进行统计。

输入与输出

输入文件

目标输出

• .out 格式（cross_match 输出格式，14 列）
• .cat 格式（RepeatMasker 内部格式，14-15 列，含 lineageId）
• .tbl 格式（汇总统计表）

工作历程

第一阶段：GFF → .out 格式转换

编写 gff2out.sh，将 GFF 的 9 列解析并映射到 .out 格式的 14 列。

成功映射的字段（13/14 列）：

• SW score、PercDiv/Del/Ins（分歧率）、seqid、query 起止坐标、strand、repeat 名称、class、consensus 起止坐标、元素 ID

无法映射的字段：

• Col 8（query remaining）：需要 scaffold 全长 → 通过用户提供的 lengths.txt 解决
• Col 12（subjRemaining，共有序列剩余碱基）：需要 repeat 库中共有序列的全长 → 最终未解决

第二阶段：发现 ProcessRepeats 格式兼容性问题

运行 ProcessRepeats output.out 时报错：

The input file was generated using an older RepeatMasker format.

原因：ProcessRepeats v4.x 期望的是 新的 .cat 格式（含 lineageId 字段），而非旧的 .out 格式。

新旧格式的关键区别：

修复：更新脚本，输出新 .cat 格式，添加 lineageId，将 class 嵌入 repeat name，修复正向链省略 strand 等。

第三阶段：添加文件头信息

运行更新后的脚本时，ProcessRepeats 报警告：

Warning: Calculation of statistics hampered by missing sequence sizes

原因：.cat 文件需要包含 ## Total Sequences 和 ## Total Length 头信息行，供 ProcessRepeats 计算百分比。ProcessRepeats 的解析逻辑（ProcessRepeats 第 323-352 行）通过读取这些 ## 开头的头信息行来获取基因组总长和序列数。

修复：在输出文件开头添加头信息：

## Total Sequences: 1
## Total Length: 20000
RepeatMasker version 4.0.6

第四阶段：统计结果偏小问题的根因分析

头信息修复后，.tbl 统计结果仍小于公司提供的参考值。

通过阅读 ProcessRepeats 源代码定位到关键逻辑（ProcessRepeats 第 1758-1767 行）：

# Don't bother if consensus remaining is lower
next if (
    ($currentEquiv->getOrientation() ne "C"
     && $currentEquiv->getSubjRemaining() < 20)
    || ($currentEquiv->getOrientation() eq "C"
        && $currentEquiv->getQueryStart() < 20)
);

这段代码在元素合并（element joining）阶段执行：对于正向链匹配，如果 subjRemaining < 20，该匹配将被跳过，不参与元素合并。由于我们在转换时将所有条目的 subjRemaining 都填充为 0，导致所有正向链匹配都被错误地过滤掉了，从而使统计结果偏小。

第五阶段：尝试绕过 ProcessRepeats 直接生成 .tbl

编写 gff2tbl.sh，从 GFF + scaffold 全长直接计算 .tbl 统计。

结果：bp 总数和百分比、元素个数（number of elements）均偏大。 原因是 GFF 中每个片段被计为一个元素，而 ProcessRepeats 会将属于同一次插入事件的多个片段合并计数，而且会对部分不符合“评分条件”的进行过滤。

第六阶段：寻求 subjRemaining 的正确值

要正确计算 subjRemaining，需要知道每个 repeat 家族共有序列（consensus sequence）的全长：

subjRemaining = consensus_total_length - max(consensus_start, consensus_end)

共有序列全长来源于重复序列库（Repeat Library），这是一个 FASTA 文件，包含各重复家族的代表性序列。

获取途径：

1. 从原始 RepeatModeler 输出目录的 consensi.fa.classified 文件提取
2. 从使用的 RepBase 数据库获取已知元素的序列
3. 从基因组重新运行 RepeatModeler 进行 de novo 预测
4. 从 GFF 自身近似推断（取每个 repeat 的最大 consensus_end 作为下限）

最终结论

由于以下因素，项目无法继续：

1. 共识序列全长的缺失：转换时只能将 subjRemaining 填充为 0，触发了 ProcessRepeats 的过滤逻辑
2. 中间文件已释放：公司原始分析中生成的 RepeatModeler 中间产物（含 consensus 全长信息）已不可获取
3. 数据库版本未知：无法确定公司使用的 Dfam/RepBase 数据库的具体版本及从头预测的参数配置
4. 近似方法不可靠：用 GFF 推断的共识序列长度是下限值，可能导致元素合并逻辑仍有偏差

最终选择：寻求公司使用原始的 RepeatMasker 运行环境和中间文件重新进行统计。

涉及的生物信息学知识

1. 重复序列注释概述

真核生物基因组中含有大量重复序列（repetitive elements），主要包括：

• 转座子（Transposons）：可移动的 DNA 序列
• 逆转座子（Retrotransposons）：通过 RNA 中间体转座，包括 LINE、SINE、LTR 等
• DNA 转座子：直接通过 DNA 切割和插入转座
• 简单重复（Simple repeats）：如微卫星序列（microsatellites）
• 低复杂度区域（Low complexity regions）

RepeatMasker 是最广泛使用的重复序列注释工具。

2. GFF3 格式（Generic Feature Format version 3）

生物信息学中描述基因组特征（如基因、重复序列）的通用格式。9 列制表符分隔：

第 9 列属性中与本项目相关的 key：

• Target：格式 repeat_name:start-end，表示比对到 repeat 库共有序列上的坐标范围
• Class：重复序列分类（如 LTR/Gypsy、LINE/L1）
• PercDiv：替换率（% substitution）
• PercDel：缺失率（% deletion）
• PercIns：插入率（% insertion）

3. cross_match 与 .out 格式

cross_match 是 Phil Green 开发的序列比对程序，RepeatMasker 使用它（或替代引擎）将基因组序列比对到重复序列库。

.out 文件是 cross_match 的标准输出格式，14 列空格分隔：

4. .cat 格式及其与 .out 的区别

.cat 文件是 RepeatMasker 的内部格式，作为 ProcessRepeats 的直接输入。

新 .cat 格式（v4.x）的关键特征：

• 正向链省略 strand 列（无 + 标识），反向链保留 C 标识
• 无独立的 class 列；class 信息通过 repeatName#class/subclass 格式嵌入在 repeat name 字段中
• 包含 lineageId 字段：格式为 m_bNsNiM（正义链）或 c_bNsNiM（互补链），其中 N 为替换数估算值，M 为插入数估算值，b 后的数字为批次号
• 反向链的共识坐标顺序为 subjRemaining subjEnd subjStart（注意 start 和 end 的位置与正向链不同）
• 文件开头需包含 ## Total Sequences: 和 ## Total Length: 头信息行
• 文件头信息段以版本行结束（如 RepeatMasker version 4.0.6）

格式检测逻辑（ProcessRepeats 源码第 705-749 行）：通过检查特定列的值（是否有 + 号、是否存在 [mc]_b 模式的 lineageId）来区分 .out、.cat 和原始 cross_match 输出。

5. ProcessRepeats 与 .tbl 格式

ProcessRepeats 是 RepeatMasker 的后处理程序，主要功能：

• 元素合并（element joining）：将属于同一次转座插入事件的多个片段（因插入其他元件而被分割）识别并合并，计数为一个元素
• 分类汇总：按 repeat 类别（SINE、LINE、LTR、DNA 等）统计元素数和占用碱基数
• 生成 .tbl 统计表：输出标准格式的汇总报告

处理流程：

1. 解析 .cat 文件头信息（序列数、总长等）
2. 按类别层级（class/subclass）将注释分配给相应统计桶
3. 多轮处理（cycle 1、cycle 2 等）：去除边缘效应、合并碎片、重命名等
4. 生成 .tbl 输出

元素合并（element joining）逻辑：ProcessRepeats 判断两个注释片段是否来自同一次插入事件的依据包括：

• 比对到同一 repeat 家族的相近区域
• 在查询序列中位置相邻
• 在共识序列中坐标相邻或重叠
• subjRemaining 值（共识序列剩余碱基）——如果值过小（< 20），表示匹配已到达共识序列末端，不作为合并候选

6. 重复序列库（Repeat Library）与共有序列（Consensus Sequence）

重复序列库是包含各重复家族代表性序列（共有序列）的 FASTA 文件。

来源：

• Dfam（开源）：重复序列家族数据库，包含 curated 的共有序列和 HMM 模型
• RepBase（商业）：最全面的真核生物重复序列数据库
• RepeatModeler：基于基因组 de novo 预测重复序列的工具，产出包含 RepeatModeler_rnd-*_family-* 等命名的共有序列

共有序列全长的重要性：Target 属性中的坐标（consensus_start-consensus_end）是该注释在共有序列上比对到的区域。但仅凭起止坐标无法知道注释覆盖了共有序列的多少比例——subjRemaining = 全长 - max(坐标) 代表匹配终点之后还有多少碱基未被覆盖。这个值影响：

• 是否将注释视为"全长匹配"
• 元素合并时是否将其作为合并候选

7. 链方向（Strand）约定

• +（正义链 / Watson 链）：匹配方向与共有序列方向一致
• C（互补链 / Crick 链）：匹配方向与共有序列方向相反

对共识坐标的影响：

• + 链：匹配从较低坐标走向较高坐标（subjStart < subjEnd）
• C 链：在 top-strand 编号下匹配从较高坐标走向较低坐标（subjStart > subjEnd）

GFF 中 Target 的坐标始终是 start < end（与链方向无关），链信息在 GFF 第 7 列单独提供。

8. Kimura 双参数距离与分歧率

在重复序列注释中，"分歧率"（divergence）用于衡量注释片段与共有序列之间的演化距离。常用的度量包括：

• PercDiv：替换率（substitutions），即碱基替换的百分比
• PercDel：缺失率（deletions），query 序列中相对于共有序列的缺失
• PercIns：插入率（insertions），共有序列中相对于 query 序列的缺失

这些值与Kimura 双参数距离（Kimura 2-parameter distance）相关，后者是分子演化中用于估计序列分歧时间的标准方法，考虑了转换（transitions）和颠换（transversions）的不同速率。

9. 重复序列分类体系

标准的重复序列分类（.tbl 中使用的层级结构）：

10. RepeatModeler 与 de novo 重复序列预测

RepeatModeler 是 de novo（从头）重复序列识别工具，工作流程：

1. 构建基因组索引：BuildDatabase -name db_name genome.fa
2. 运行预测：使用 RECON 和 RepeatScout 两种算法分别识别重复序列
3. 聚类与分类：将相似的序列聚类为家族，生成共有序列
4. 输出：consensi.fa.classified 包含分类后的共有序列 FASTA 文件

缺点：只能识别在基因组中出现多次且相似度较高的重复家族；无法为已知名称的元件（如 Alu、LINE1）提供准确的学术命名。

11. lineageId 格式

lineageId 是 RepeatMasker v4.x 引入的批次标识符，格式为：

m_b<批次号>s<替换数>i<插入数>    # 正义链
c_b<批次号>s<替换数>i<插入数>    # 互补链

示例：m_b7s551i6 表示正义链、批次 7、约 551 个替换、约 6 个插入。 m_ 前缀表示 main strand（正义链），c_ 表示 complement strand（互补链）。

它在 ProcessRepeats 中用于判断链方向（PRSearchResult.pm 第 244-252 行），并且在格式检测中作为新版本文件的标志（必须匹配 [mc]_b\d+s\d+i\d+ 模式）。