3. 重复序列注释评估：覆盖率、分类准确性与转座子识别
重复序列注释，说白了就是给基因组里的“暗物质”做标记。你想想看，一个人类基因组里超过50%的序列都是重复的，要是这步搞砸了，后面的基因预测、变异检测全得跟着跑偏。
我个人习惯在评估重复注释质量时，死磕三个核心指标：重复覆盖率、重复分类准确性，以及转座子识别。今天咱们就把这三个硬骨头啃下来。
3.1 重复覆盖率：注释到底“够不够”
重复覆盖率，衡量的是你注释出来的重复序列占整个基因组的比例。我在项目里遇到过不少新手，上来就问我：“覆盖率是不是越高越好？” 嗯，答案没那么简单。
如果一个基因组真实重复比例是45%，你注释出来只有20%，那肯定是漏了很多。但如果注释出来80%，多半是误把非重复区也标成了重复——这叫过注释，比漏注释更麻烦。
那怎么量化注释的完整性呢？我个人习惯用下面这组命令：
# 计算基因组总长度
samtools faidx genome.fa
awk '{sum+=$2} END {print sum}' genome.fa.fai > genome_length.txt
# 用RepeatMasker的输出统计
# .out 文件第6列是“重复覆盖的碱基数”
awk 'NR>3 {sum+=$6} END {print sum}' genome.fa.out > repeat_bases.txt
# 覆盖率 = 重复碱基数 / 基因组总长度
coverate=$(awk 'NR==FNR{a=$1;next}{print $1/a*100}' genome_length.txt repeat_bases.txt)
echo "重复覆盖率: $coverate%"不过这里有个大坑。RepeatMasker的.out文件里，有些低复杂度区域（比如微卫星）也算进去了。如果你的研究重点是转座子，建议单独统计。我曾经帮一个水稻项目做评估，一上来覆盖率看着挺高（52%），结果去掉低复杂度序列之后，真实转座子覆盖率直接掉到38%。
3.2 重复分类准确性：别把LINE当成SINE
分类准确性比覆盖率更难搞。你注释出来的覆盖率再高，如果大部分重复序列都被标成了“Unknown”或者“Other”，那这个注释的质量就大打折扣了。
我记得有次处理一个植物基因组，结果发现超过60%的重复序列被标成了“Unknown”，一下子就把注释质量打回了原形。我花了整整两周重建重复库，才把Unknown比例降到15%以下。
那怎么定量评估分类准确性呢？我一般用RepeatMasker的.tbl文件来查：
# 查看分类统计头
head -50 genome.fa.tbl
# 输出示例
==================================================
file name: genome.fa
sequences: 5000
total length: 250000000 bp
GC level: 42.00 %
bases masked: 130000000 bp (52.00 %)
==================================================
type count length percentage
SINE 12000 3500000 1.40 %
LINE 45000 15000000 6.00 %
LTR 40000 32000000 12.80 %
DNA_TE 15000 7500000 3.00 %
Unknown 8000 4200000 1.68 %
Low_complexity 2000 600000 0.24 %
Satellite 500 250000 0.10 %
Simple_repeat 3000 300000 0.12 %
Other 5000 2500000 1.00 %我会重点关注Unknown类别的占比。如果这个值超过20%，注释的质量就悬了。
为了方便对比不同方法的结果，我画过一张分类准确性评估表，你直接拿去用：
| 评估维度 | 高质量标准 | 可接受标准 | 需改进标准 |
|---|---|---|---|
| Unknown占比 | < 5% | 5% - 15% | > 15% |
| 主要超家族覆盖率 | 全部覆盖率 > 90% | 主要超家族覆盖率 > 80% | 多个超家族 < 60% |
| 与参考注释的一致性 | 分类标签一致率 > 85% | 70% - 85% | < 70% |
| 转座子家族数量 | 覆盖已知所有主要家族 | 缺少1-2个次要家族 | 缺少3个以上家族 |
你还可以用RepeatClassifier（RepeatModeler自带的工具）做再分类，看看之前被标成Unknown的那些序列能不能进一步归类：
# 将Unknown序列提取出来
grep "Unknown" genome.fa.out | awk '{print $5}' > unknown_headers.txt
seqtk subseq genome.fa unknown_headers.txt > unknown_seqs.fa
# 用RepeatClassifier分类
RepeatClassifier -consensi unknown_seqs.fa -engine ncbi
# 输出新的分类结果
cat unknown_seqs.fa.classified > reclassified.lib3.3 转座子识别：把“跳跃基因”抓出来
转座子（TE）是重复序列里的主角。我评估转座子识别质量时，关注三件事：全基因组TE覆盖率、TE家族完整性、以及TE插入时间分布。
先说TE覆盖率。你需要在RepeatMasker输出中单独统计各超家族的占比：
# 从.out文件统计各超家族覆盖度
awk 'NR>3 && $11=="DNA" {sum_dna+=$6}
NR>3 && $11=="LINE" {sum_line+=$6}
NR>3 && $11=="LTR" {sum_ltr+=$6}
END {print "DNA:",sum_dna,"LINE:",sum_line,"LTR:",sum_ltr}' genome.fa.out再来说TE家族完整性。一个完整TE家族应该有：完整的两端重复序列（LTR或TIR）、目标位点重复（TSD），以及编码区没有提前终止密码子。
我习惯用LTR_Finder和LTR_retriever来鉴定完整LTR转座子：
# 运行LTR_Finder
ltr_finder -D 15000 -d 1000 -L 7000 -l 100 -p 20 -M 0.9 genome.fa > genome.ltr.finder.scn
# 用LTR_retriever评估完整性
LTR_retriever -genome genome.fa -infinder genome.ltr.finder.scn -threads 20
# 查看结果摘要
cat genome.fa.LTR_retriever.pass.list | wc -l # 完整LTR数量
cat genome.fa.LTR_retriever.fa | grep ">" | head
# 输出示例
>LTR_R1_chr1_123456_134567
>LTR_R2_chr2_234567_245678最后是TE插入时间。LTR转座子的两个LTR序列的差异可以推算插入时间，这是评估TE活性的重要指标：
# 用LTR_retriever估算插入时间
# 查看LTR序列差异分布
awk '{if(NR>0) print $5}' genome.fa.LTR_retriever.pass.list | sort -n > ltr_divergence.txt
# 用R或Python画个直方图
# 假设突变率为1.3e-8 per site per year
# 时间 = 差异度 / (2 * 突变率)我当年做玉米基因组注释的时候，发现B73自交系有接近75%的LTR转座子是近100万年内插入的。这个信息直接告诉我们该物种的基因组还处于“活跃扩张期”。
我把整个重复序列评估的逻辑画了一张图，你看一眼就清楚了：
嗯，这章就讲到这里。这三个指标——覆盖率、分类准确性、转座子识别——都是实打实能通过命令和数据验证的。我建议你在项目前期就把这些检测点埋进去，避免最后出问题倒查时，才发现是重复注释这步埋了雷。