参考基因组与基因注释文件：质控的基石
做RNA-seq质控，你绕不开两个东西：参考基因组和基因注释文件。说白了，它们就是你比对和定量的“地图”和“图例”。地图错了，读段就贴不对位置；图例错了，基因表达量就是一堆废数。我这些年踩过的坑，有一半都跟这两个文件有关。
参考基因组：版本不同，结果天差地别
参考基因组本质上是某个物种的完整DNA序列拼装结果。人类基因组现在主流是GRCh38（也叫hg38），但很多人还在用GRCh37（hg19）。我跟你说个真事：有次一个项目比对率只有65%，换了GRCh38后直接跳到92%。为啥？GRCh37里那些gap和未定位的scaffold，会让大量读段变成unmapped。
基因组拼装有不同层级：
| 组装级别 | 描述 | 常用场景 |
|---|---|---|
| Complete Genome | 每条染色体完整无gap | 模型物种（人、小鼠） |
| Chromosome-level | 染色体骨架+部分gap | 大多数脊椎动物 |
| Scaffold-level | contig拼接为scaffold | 新测序物种 |
| Contig-level | 仅连续片段 | 快速组装项目 |
警告：千万别混用不同版本的基因组做验证。比如你用GRCh38建的索引，却用GRCh37的注释文件去定量——基因坐标都对不上，结果全是错的。我曾经因此浪费了两周时间。
GTF/GFF：基因组上的“标签”文件
基因注释文件告诉你：哪个位置是基因、哪个片段是外显子、哪里是UTR。我最常用的就是GTF（Gene Transfer Format），它基于GFF2衍生而来，每行9列，用tab分隔。
# 典型GTF行示例（来自Ensembl）
1 ensembl gene 11869 14409 . + . gene_id "ENSG00000223972"; gene_name "DDX11L1";
1 ensembl transcript 11869 14409 . + . gene_id "ENSG00000223972"; transcript_id "ENST00000456328";
1 ensembl exon 11869 12227 . + . gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; exon_number 1;GFF3是更新版本，属性字段更规范。两者差别不大，但注意很多工具只认GTF。我个人习惯一律用GTF，省得转换。
关键列含义：
- seqid：染色体或scaffold名称，必须和参考基因组的FASTA文件里的名字完全一致。
- source：注释来源（如ensembl, havana）。
- feature：基因、转录本、外显子、CDS等。这里最容易藏问题——有些GTF里feature名称不规范，比如用“exon”写成了“exon_part”，会导致一些工具识别失败。
- start/end：1-based坐标。注意和BED的0-based区别。
小技巧：拿到新GTF后，先跑 grep -c "exon" 看看数量级。小鼠mm10的GTF大概有100万+ exon行，如果少了一个数量级，基本可以确定文件损坏或下载不完整。
基因注释质量如何影响比对率与表达分布
你可能觉得比对率只跟基因组序列有关，跟注释没关系。错了！很多比对器（比如STAR、RSEM）在比对过程中会利用注释信息来优化junction检测。我来给你拆解一下：
- 注释越全，splice-aware比对越准——GTF里提供的已知剪切位点，能帮STAR优先把跨越内含子的读段定位对。
- 基因重叠区域定量的歧义——同一读段落在一个基因的exon和另一个基因的intron重叠区，注释文件决定了你是计数给基因A还是丢弃。
- 基因长度计算决定FPKM分布——FPKM = 每千碱基每百万读段。如果GTF里外显子总长度算错了，表达量分布就会整体偏移。我遇到过GTF里某个长链非编码RNA exon坐标多了一位，导致长度翻倍，FPKM直接砍半。
来看一张流程图，直观感受参考基因组和注释文件在整个质控流程中的位置：
实战：GTF质量检查三步法
我每次拿到新的注释文件，都会跑下面三个检查。
grep -v "^#" gencode.v38.annotation.gtf | cut -f3 | sort | uniq -c | sort -rn正常输出应该看到大量 exon、CDS、transcript、gene。如果 exon 数量突然很少，八成是下载时丢数据了。
# 提取GTF中所有seqid
cut -f1 gencode.v38.annotation.gtf | grep -v "^#" | sort -u
# 提取FASTA中所有染色体名
grep "^>" GRCh38.primary_assembly.genome.fa | sed 's/>//' | sort -u两个列表必须完全匹配。我曾经因为FASTA里是"chr1"而GTF里是"1"，导致STAR报错“no chromosomes in common”。
awk '$3=="gene" && $10~"ENSG00000141510"' gencode.v38.annotation.gtf拿一个已知基因（比如TP53）看看它的起止位置和已知数据库是否一致。如果chr17:7,660,000-7,690,000这种经典区间都差了几Mb，赶紧换来源。
选哪个数据库？
| 数据库 | 特点 | 我的使用建议 |
|---|---|---|
| GENCODE | 最全面，包含人工审核 | 人类和小鼠首选，28版本以后很好用 |
| Ensembl | 持续更新，覆盖广 | 非模式物种只能靠这个 |
| RefSeq | NCBI官方，保守 | 需要和NCBI其他资源联用时用 |
| UCSC | 常用浏览器格式 | 他们提供的GTF有时feature类型偏少 |
我个人的习惯：人类数据全部走GENCODE，其他哺乳动物用Ensembl。如果你是在做跨物种比较，记得保持所有物种的注释版本号一致，不然基因坐标体系不同会导致假阳性。
你在比对率突降或者表达分布异常的时候，优先检查参考基因组和注释文件。先看版本对不对，再看是否完整，最后看染色体命名。我在项目里用这几招解决过不下十次线上bug。嗯，这就是地基牢不牢的问题——地基不稳，上面盖得多漂亮都没用。