企业官网建设流程全解析

昨天有个做组学的朋友发来一份人类参考基因组的注释文件，问我说：“我就想把每个基因的染色体、起止坐标整理成一张表，后面要去提启动子、做富集分析，怎么弄最省事？”我第一反应是，这还不简单，把第三列 feature 等于 gene 的行抽出来，取第 1、4、5 列不就行了。但等文件真正摊开以后，才发现坑全藏在第九列那一长串 attribute 里：gene_id 和 gene_name 谁拿来做主键？基因版本号要不要剥掉？负链基因的 TSS 到底在哪一端？BED 的 0-based 坐标和 GTF 的 1-based 坐标怎么换算？任何一个环节没对齐，后面你的启动子提取、基因集富集、peak 注释，全都白做。

这篇文章就围绕“根据 GTF 格式的基因注释文件得到人类所有基因的染色体坐标”这件事，把我实际操作的完整思路、脚本和踩过的坑一条条讲清楚。适合基因组注释接触得不多、正准备把 GTF 当作普通文本硬parse 的人，也适合已经会跑命令但想确认自己没埋雷的老手。

1. 先看懂GTF这九列，坐标才不会抓错

网上搜 GTF 格式能搜到一大堆定义，但我更建议你把它当成一张“每一行注释一个基因组元件”的大表。人类参考基因组的注释文件动辄 1GB 以上，行数上百万，如果不搞清楚每一列的含义，后面的过滤和提取就是在盲人摸象。

1.1 九列字段的职责划分：真正与坐标有关的只有四五项

标准的 GTF 每一行有九列，用 Tab 分隔，典型的一行长这样：

chr1 gene 11869 14409 . + . gene_id "GENE00000001"; gene_name "SYMBOL_X"; gene_biotype "lncRNA";

注意这里的列顺序有个很容易忽略的点：第二列是 source，不是染色体。染色体在第一列。下面这张表我每次做注释文件处理都会先默念一遍：

先说清楚 start 和 end 的一个底层约定：无论基因在正链还是负链，GTF 里存的 start 永远是较小的坐标，end 永远是较大的坐标。所以一个负链基因，它的转录起始位点 TSS 并不在 start 位置，而是在右边界 end 那一边。很多新手拿到负链基因的坐标直接拿 start 当 TSS，后面做启动子提取时方向就全反了，这一点后面专门展开。

染色体命名同样是个暗坑。有的注释文件里是 chr1、chr2 这种带 chr 前缀的写法，有的则是 1、2 这种裸编号，线粒体有的写作 MT，有的写作 chrM，还有随机对照的 contig、alt 染色体的名字五花八门。如果这份坐标表还要跟别的表做 join，染色体命名不一致时你会得到一堆莫名其妙的空匹配。

1.2 feature 类型的选择：为什么必须卡在 gene 这一层

GTF 文件里每一行都有第三列 feature，取值常见的有 gene、transcript、exon、CDS、start_codon、stop_codon、five_prime_utr、three_prime_utr 等等。它们之间的层级关系可以理解成一棵树：一个基因可以有多个转录本，一个转录本由多个外显子组成，外显子里面才有 CDS、UTR 这些更细的元件。

如果我们要的是“基因级坐标”，即这个基因在染色体上占用的大致区间，那么最直接的来源就是 feature 等于 gene 的那些行。注释数据库在生成 gene 行时，通常会把该基因所有转录本的外显子起点取最小、终点取最大，合并成一个完整的基因跨度。也就是说，gene 行自带的 start 和 end，就已经是我们想要的“全基因占据区间”，不需要再去 transcript、exon 层做手工聚合。

反过来，如果你过滤时没限制 feature，只是把文件里所有行都抽出来，一个基因会伴随几十上百行 exon 记录，看起来同一个基因重复出现几十次，根本无法直接使用。过滤到 transcript 层虽然行数变少，但一个多转录本基因会出现好几行，每行的起止区间还不一样，仍然需要额外处理。

这里要提两个现实中的例外。第一，部分年代较早或来源特殊的注释文件里没有显式的 gene 行，只有 transcript 和 exon 行，这时只能按 gene_id 聚合所有元件的起止来重建基因跨度，方法我放在第 3 节。第二，feature 的写法在不同版本里有细微差别，绝大多数现代 GTF 是小写 gene，但我见过老 GFF 转出来的文件里写成 Gene，所以正式处理前最好先跑一句看清楚：

cut -f3 annotation.gtf | sort | uniq -c | sort -nr | head

先搞清楚这个文件里 feature 到底有哪些取值、大小写长什么样，再写过滤条件，能避免很多低级错误。

1.3 attribute 列是重灾区：gene_id、gene_name、gene_biotype 不能混着用

第九列 attribute 是一串分号隔开的键值对，也是处理 GTF 时唯一需要“解析”而不是“切列”就能拿到的信息。常见键包括：

• gene_id：唯一的稳定 ID，比如形如 ENSG 开头的标识，这才是我们做外键匹配时的真正主键；
• gene_name：人类可读的基因符号，展示用；
• gene_biotype：基因类型，protein_coding、lncRNA、pseudogene 等，筛子基因或者统计基因类型分布时非常有用。

为什么一定要拿 gene_id 当主键而不是 gene_name？因为 gene_name 理论上不是严格唯一的：不同数据库之间可能存在同名冲突，某些新注释基因可能没有可读符号，直接显示为空或类似 unknown 的占位。gene_id 则是注释系统内部保证唯一的标识，用它做 join 最稳。输出的坐标表里，我会把 gene_id、gene_name、gene_biotype 同时保留，谁也别丢，后面无论你想按符号查还是按 ID 查、想筛 biotype 还是不过滤，都有现成的列可用。

另外一个需要盯防的是 gene_id 的版本号后缀。不少注释文件的 gene_id 长得像GENE00000001.17这种，点后面跟的是版本号。如果你后续要用这个 ID 去和别的注释资源做匹配，而另一份表里存的是不带版本号的裸 ID，两边直接匹配会大量掉数据。稳妥做法是在提取的时候就剥掉版本号，规则见 3.3 节。

2. 解析的三条路线：awk、Python，以及那个我不推荐的 grep

明确了 GTF 的结构以后，接下来就是选工具把目标行抽出来、把 attribute 拆开、输出成你想要的模样。同一个任务可以有三种主流干法，效率、灵活性和可读性各有侧重。

2.1 awk 一行流：快、省内存、适合一次出表

针对“从 GTF 提取基因坐标”这种相对固定的操作，awk 是我最常用的方案。它的优势在于单遍扫描、不把整个文件读进内存，对 GB 级文件非常友好，通常几秒到几十秒就能跑完。下面这段脚本可以直接存成一个文件，也可以一行敲进终端：

awk -F'\t' 'BEGIN{OFS="\t"} $1 !~ /^#/ && $3=="gene" { gid=""; gname=""; biotype=""; n=split($9, arr, ";"); for (i=1; i<=n; i++) { gsub(/^[ ]+/, "", arr[i]); if (arr[i] ~ /^gene_id/) { gid=arr[i]; sub(/^gene_id[ ]+/, "", gid); gsub(/"/, "", gid); } else if (arr[i] ~ /^gene_name/) { gname=arr[i]; sub(/^gene_name[ ]+/, "", gname); gsub(/"/, "", gname); } else if (arr[i] ~ /^gene_biotype/) { biotype=arr[i]; sub(/^gene_biotype[ ]+/, "", biotype); gsub(/"/, "", biotype); } } if (gid != "") { print $1, gid, gname, biotype, $7, $4, $5, $5 - $4 + 1; } }' annotation.gtf > human_gene_coords.tsv head human_gene_coords.tsv

拆解一下关键点。-F'\t'是把分隔符指定成 Tab，避免默认按空白切分时 attribute 内部的空格把列切碎。$1 !~ /^#/是跳过文件头部的注释行，很多 GTF 文件开头有几行以 # 开头的元信息，不跳的话会混进很脏的行。$3=="gene"把范围卡死在基因层。attribute 的解析我用了一个循环：先按分号切成若干段，每一段去掉开头空格，然后判断键名，再把键名、值两侧的引号剥掉。判断时用了^gene_id的前缀匹配，这是为了避免某个键名恰好以别的键开头导致误匹配。

输出我设计了八列：染色体、gene_id、gene_name、gene_biotype、链方向、start、end、基因长度。第八列长度纯粹是顺手算的，后面统计基因长度分布时会用得上。

2.2 Python 逐行解析：灵活、可扩展，适合后面接复杂管线

如果这次不只是提取坐标，后面还要做更复杂的处理，比如按 biotype 拆分、生成多种输出格式、或者接进已有的分析流程，我会用 Python 写一个低层解析脚本。所谓低层，是指不借助 pandas 的 read_csv 直接硬读，而是用open()配合for line in逐行处理，这样内存占用稳定，不会因为文件太大而把机器拖垮。

#!/usr/bin/env python3 import sys def parse_attributes(attr_text): attrs = {} parts = attr_text.rstrip(";").split(";") for part in parts: part = part.strip() if not part: continue key, sep, value = part.partition(" ") if sep: attrs[key] = value.strip().strip('"') else: attrs[key] = "" return attrs def main(gtf_path, out_path): with open(gtf_path, "r") as fin, open(out_path, "w") as fout: for line in fin: if line.startswith("#") or not line.strip(): continue cols = line.rstrip("\n").split("\t") if len(cols) < 9: continue chrom, _, feature, start, end, _, strand, _, attr = cols[:9] if feature != "gene": continue attrs = parse_attributes(attr) gene_id = attrs.get("gene_id", "") gene_name = attrs.get("gene_name", "") biotype = attrs.get("gene_biotype", "") if not gene_id: continue start = int(start) end = int(end) fout.write( f"{chrom}\t{gene_id}\t{gene_name}\t{biotype}\t{strand}\t{start}\t{end}\t{end-start+1}\n" ) if __name__ == "__main__": main(sys.argv[1], sys.argv[2])

用法是：

python3 extract_gene_coords.py annotation.gtf human_gene_coords.tsv

这段脚本的逻辑和 awk 版本一一对应，只是 attribute 解析用了 Python 的字符串处理，更直白一些。如果你后面还想顺手筛掉非标准染色体、按 biotype 输出多个文件，甚至直接生成 BED 格式，在这个框架上加条件非常方便，这也是我推荐 Python 的原因。

2.3 为什么我不推荐grep gene然后接 cut

很多刚上手的人会写这样一串命令：

grep gene annotation.gtf | cut -f1,4,5,7 > gene_coords.txt

我不推荐这个写法，原因很简单：GTF 里几乎每一行 transcript、exon 的 attribute 里都带gene_id字样，普通grep gene会把所有包含“gene”这个子串的行全捞出来，输出结果里一个基因重复几十次，坐标还都是外显子或转录本的，不是你想要的基因跨度。

有人会说，那我用grep -w gene精确匹配单词总行了吧？确实，正确的 feature 列恰好是独立的单词 gene，grep -w能把这个条件限定住，至少在大多数规范文件里能命中。但随之而来的第二个问题是：你仍然要处理第九列 attribute 才能拿到 gene_id、gene_name，cut -f9切出来的是一整串带引号的键值对，还得再套一层 sed 或 awk 去拆。一圈下来，相当于自己用两个不完整的工具拼了一个本可以一步完成的任务，中间还更容易在各种转义、引号处理上出错。

真想在 shell 里硬写的，我建议用更精确的正则匹配第五列的区域边界：

grep -P '^chr1\t.*?\tgene\t' annotation.gtf | head

但说实话，既然要精确过滤第三列，awk 一个$3=="gene"就解决了，完全没必要退回去用 grep。很多时候看似“简单”的连续管道，反而是出错率最高的地方。

2.4 pandas 能读吗？能，但不推荐当主力

顺便说一句 pandas。pd.read_csv(annotation.gtf, sep="\t", comment="#", header=None)确实能把这个文件读进 DataFrame，注释行也能通过 comment 参数跳过。问题在于两点：一是文件太大，动辄上 GB 的文件塞进 DataFrame 会占用好几个 GB 内存，机器差一点的直接卡死；二是第九列 attribute 的解析在 pandas 里并不省事，你需要对整列做 str 分割、strip、正则提取，写起来不比 Python 纯循环短，而且向量化处理字符串时很容易因为缺失值、引号不统一而抛异常。

所以我的经验是：一次性提取坐标，awk 就是首选；后续要做复杂处理，用 Python 的逐行解析；pandas 留给坐标提取之后的数据分析阶段，比如读入做好的 TSV 做统计画图。这样每一步都用最合适的工具，不为一时的“能用”留下隐性问题。

3. 坐标抓出来后，先绕开这四个高频坑

坐标文件拿到手，不等于任务完成。我在实际项目中至少撞过以下四类问题，每一个都会让下游分析结果不可信，而且往往是隐蔽的，不主动查根本发现不了。

3.1 有些 GTF 没有 gene 行：按 gene_id 聚合出全基因跨度

前面说现代注释文件基本都有 gene 行，但“基本”不等于“一定”。有的旧文件、某些物种的注释文件、或者从其他来源转出来的 GTF，压根没有 feature 等于 gene 的行，只有 transcript 和 exon。这时候你就得自己根据 gene_id 把所有元件聚合起来。

聚合的逻辑很简单：同一个 gene_id 下，start 取所有行里的最小值，end 取所有行里的最大值。用 awk 的关联数组可以很轻松实现：

awk -F'\t' 'BEGIN{OFS="\t"} $1 !~ /^#/ { n=split($9, arr, ";"); gid=""; gname=""; for (i=1; i<=n; i++) { gsub(/^[ ]+/, "", arr[i]); if (arr[i] ~ /^gene_id/) { gid=arr[i]; sub(/^gene_id[ ]+/, "", gid); gsub(/"/, "", gid); } else if (arr[i] ~ /^gene_name/) { gname=arr[i]; sub(/^gene_name[ ]+/, "", gname); gsub(/"/, "", gname); } } if (gid == "") next; if (!(gid in gstart)) { gstart[gid] = $4; gend[gid] = $5; gchrom[gid] = $1; gstrand[gid] = $7; gname_map[gid] = gname; } else { if ($4 < gstart[gid]) gstart[gid] = $4; if ($5 > gend[gid]) gend[gid] = $5; } } END { for (gid in gstart) print gchrom[gid], gid, gname_map[gid], gstrand[gid], gstart[gid], gend[gid], gend[gid] - gstart[gid] + 1; }' annotation.gtf > aggregated_gene_coords.tsv

人类基因总数通常在 6 万上下（包含各类非编码 RNA 和假基因），这个量级放在 awk 的关联数组里没有任何压力，跑起来和过滤 gene 行差不多快。用这段脚本时心里要有数：它会把所有 feature 类型都当作聚合输入，所以如果同一个基因在不同行里出现了跨染色体的矛盾记录，聚合出来的结果就会有严重问题。我建议聚合完以后顺手检查一下每个基因是否只出现在同一条染色体上，方法对着第 4 列的碱基数和第 1 列的染色体名去 set 去重，发现异样再回溯。

3.2 负链基因的 TSS 不是 start：方向感必须在线

我在 1.1 节埋了一个扣，这里正式解开。GTF 里 start 始终是较小的坐标，但对负链基因来说，基因转录方向是从右往左，真正的转录起始位点在右侧边界。用数字举个例子：某个负链基因的记录是 chr1 1000 5000 -，那么它的转录方向是从 5000 往 1000 走，TSS 在 5000 这一端，TES 在 1000 那一端。

如果你要提取启动子区域，或者判断一个变异位点落在基因上游还是下游，就必须做方向感知。基因级坐标拿到手后我是这样定义 TSS 的：

# 正链基因 TSS 在 start，负链基因 TSS 在 end if (strand == "+") { tss = start } else { tss = end }

这里再补一个精度相关的提醒：用基因特征行的边界当 TSS，是一个在“只有基因坐标表”前提下的折中做法。真要精细到单个转录本的启动子，TSS 应该取最靠外的那个外显子起点，而且不同转录本可能还有不同的启动子。所以如果你的下游是单基因单启动子的精细分析，建议回到 transcript 或 exon 层去构建转录本结构；如果只是做全基因组的粗略 promoter 扫描、富集分析，用基因特征行当 TSS 完全够用，市面上很多流程也是这么干的。

3.3 gene_id 带着版本号：做外键匹配前先剥掉

前面提到过，GTF 的 gene_id 往往带版本号后缀，例如GENE00000001.17。这种设计在注释系统内部是合理的，同一个基因在新版本注释里可能出现位置微调，加上版本号可以精确追溯。但是，一旦你要跟外部数据表做匹配，比如把基因坐标表拿去和某个芯片注释平台、某个甲基化探针注释、某个文献里整理好的基因集做 join，两边 ID 格式经常不一致，一方带版本号、一方不带，匹配率断崖式下跌。

我的习惯是在提取坐标这一步就把版本号剥掉，规则是匹配结尾的点加数字：

# awk 里在输出 gene_id 前执行 gid_clean = gid; sub(/\.[0-9]+$/, "", gid_clean);

如果你是用 Python，直接re.sub(r'\.[0-9]+$', '', gene_id)即可。剥版本号的时候只剥“点号加数字结尾”这种模式，不要一见点号就切，否则某些 ID 里本来包含点号的特殊情况会被误伤。剥完之后，后续所有表都用不带版本的 gene_id 做外键，匹配失败的概率会大幅下降。

3.4 染色体命名全局统一：chr前缀、MT/M、scaffold 一个都不能少

染色体命名是另一个隐蔽但致命的问题。不同来源的文件对染色体的写法可能完全不同：

如果只在 GTF 内部提取坐标，命名不统一问题不大，因为全文件都是同一套规则。但一旦这份坐标表要和 peak 文件、变异注释、芯片位点数据做交集或合并，两边染色体名就得先统一。我一般的做法是全部转成带 chr 前缀的标准写法，线粒体统一成 chrM，这样可以跟绝大多数基因组浏览器和 bedtools 兼容：

awk -F'\t' 'BEGIN{OFS="\t"} { chrom = $1; if (chrom ~ /^chr/) { chrom = chrom } else if (chrom == "MT" || chrom == "M") { chrom = "chrM" } else { chrom = "chr" chrom } $1 = chrom; print; }' human_gene_coords.tsv > human_gene_coords_std.tsv

另外要决定一件事：标准的常染色体和性染色体以外的 contig、alt 染色体、随机链要不要保留。很多全基因组分析工具对非标准染色体支持很差，或者你统计基因数时根本不想让那些 scaffold 干扰。我会在生成最终表之前加一步过滤，只保留^chr([0-9]+|X|Y|M)$这类规则命中的染色体，把其余的行直接丢掉。做这个决定前想清楚：这一步是真正的数据取舍，不是单纯的格式转换，丢掉的基因会在后续分析中永远消失，所以最好在注释文件里记录一下过滤掉了多少行、哪些类别，方便追溯。

4. 坐标系换算与输出格式：别让1之差毁掉整张表

坐标提取完成以后，很多人的下一步是拿着这张表去对接各种工具。这时候最大的坑是坐标系混用。GTF 内部用的是 1-based 全闭区间，而 BED 格式、bedtools、UCSC 浏览器用的是 0-based 半开区间，差之 1 毫厘，谬以千里。

4.1 GTF 是 1-based，BED 是 0-based：换算公式就一行

先花 30 秒把概念理清：

• 1-based 全闭区间：染色体上第一个碱基的位置是 1，start 和 end 都是实际包含的碱基位点；
• 0-based 半开区间：染色体上第一个碱基的位置是 0，区间写作 [start, end)，表示包含 start 位置的碱基，但不包含 end 位置的碱基。

因此，把一个 GTF 基因坐标转成 BED 坐标的规则只有一个：

bed_start = gtf_start - 1 bed_end = gtf_end

举例：

为什么 BED 的 end 保持不变？因为半开区间右端不包含在内，BED 的 end 需要比最后一个实际碱基的 0-based 索引大 1。最后一个实际碱基在 1-based 是 5000，0-based 索引是 4999，所以 BED end 写 5000。这就是换算公式里只减 start 不减 end 的原因。

很多人会在这里精神内耗：我到底是减 start 还是减 end？我建议你记住下面这句话：把 GTF 坐标往 BED 转，永远是把 left 边界减 1，right 边界不动；把 BED 坐标往 GTF 转，则反过来是 left 边界加 1。不要凭感觉改，每一步都对照上面的公式写。

4.2 四种常见输出：TSV、BED、GFF、区间表怎么选

坐标表最终以什么形式落地，取决于下游工具要吃什么。

如果你只是自己在 R 或 Python 里做后续分析，一个整洁的 TSV 最合适，列名起清楚，读进来直接当 data.frame 或 DataFrame 用。我常用的列顺序是：

chrom gene_id gene_name biotype strand start end length

如果你要对接 bedtools、UCSC 浏览器、或者做区间交集，那就转成 BED。BED 至少要有前三列（chrom、start、end），通常建议给上第 4 列 name（这里可以放 gene_name）和第 6 列 strand。从上面那张 TSV 转 BED 的命令非常直接：

# 输入列: chrom gene_id gene_name biotype strand start end length # 输出 : chrom start-1 end gene_name 0 strand awk -F'\t' 'BEGIN{OFS="\t"} {print $1, $6 - 1, $7, $3, 0, $5}' human_gene_coords.tsv > human_gene_coords.bed

如果你习惯用 UCSC 风格做可视化，直接把这个 BED 文件拖进浏览器就能看。至于 GFF 格式，一般只需要在极少数要求“必须提供注释原点”的工具里才会用到，日常任务基本用不上，这里不展开。

还有一类情况是工具本身要求 1-based 坐标，比如 R 里的某些 Bioconductor 对象，那就直接用 GTF 风格的 start/end 列，不要再转。判断标准很简单：读工具的文档，确认它内部用的是 0-based 还是 1-based，再决定用哪个版本的表。

4.3 实操：一张坐标表直接产出启动子区域 BED

这里给一个我经常用的完整场景：从基因坐标表生成全基因组启动子区域 BED，然后抓序列。所谓启动子，我这里定义成“TSS 上游 2000bp 到 TSS 下游 500bp”。

方向感知很重要，负链基因的 TSS 在 end，而且上游在坐标更大的方向。于是生成启动子的 awk 脚本要分正负链处理：

# 输入列: chrom gene_id gene_name biotype strand start end length # 输出 : 启动子区域的 BED，TSS 上游 2000bp，下游 500bp awk -F'\t' 'BEGIN{OFS="\t"} { chrom = $1; gene = $3; strand = $5; start = $6; end = $7; if (strand == "+") { tss = start } else { tss = end } if (strand == "+") { gstart = tss - 2000; gend = tss + 500 } else { gstart = tss - 500; gend = tss + 2000 } if (gstart < 1) gstart = 1; print chrom, gstart - 1, gend, gene, 0, strand; }' human_gene_coords.tsv > human_gene_promoters.bed

解释一下负链分支的写法。负链基因的 TSS 在 end，它上游 2000bp 是坐标更大的那一边，下游 500bp 是坐标更小的那一边。所以 gstart 取 tss - 500，gend 取 tss + 2000，这样区间才真正对应“上游-下游”的生物学方向。如果直接把正链的 gstart = tss - 2000、gend = tss + 500 套到负链上，得到的区间实际上是下游 2000bp 加 上游 500bp，方向完全反了。

生成 BED 之后，抓启动子序列就一行 bedtools：

bedtools getfasta -fi reference.fa -bed human_gene_promoters.bed -s -name > human_gene_promoters.fa

这里-s表示按链方向取序列，-name用 BED 第 4 列的基因名作为 fasta 的序列名。如果你要做的是区间交集，比如“手头有一批变异位点，看看落在哪些基因的启动子区域”，直接拿 my_variants.bed 和这个 promoters.bed 做bedtools intersect -wa -wb就行。

5. 交付前的质检：坐标表也要有“验收标准”

经常有人把坐标表发给我，说“跑完了，你帮我看看对不对”。我从来不会直接拿去做分析，而是先在心里过一遍验收清单。坐标表这种东西，错一个碱基都能让下游分析出现离群值，所以宁愿多花五分钟查一遍，也不要后面花五小时重新跑流程。

5.1 数量级先对一遍：基因总数和 protein_coding 占比

第一个要核对的是行数。不同注释版本之间基因数有波动，但量级是稳定的。人类注释文件里 gene 行总数大概在 6 万上下，其中 protein_coding 在 2 万左右。如果跑完发现输出只有几千行，或者有几百万行，那基本说明过滤条件写错了，比如把 transcript 当成了 gene，或者把 gene 以外的行也混了进来。

核对命令很简单，直接在原始 GTF 上数：

# 基因总数 awk -F'\t' '$3=="gene"' annotation.gtf | wc -l # 蛋白质编码基因数 awk -F'\t' '$3=="gene" && $9 ~ /gene_biotype "protein_coding"/' annotation.gtf | wc -l

我还会顺带统计一下 biotype 分布的 summary，避免某些冷门 biotype 的数量异常到能说明文件版本不对。做完这些数量核对，再往下做逐行校验。

5.2 逐行校验：start < end、strand 合法、无空 gene_id

数量对得上只能说明大致结构没坏，还要逐行检查数据本身的合理性。我通常一次性跑几组检查命令，全部通过才算合格：

# 1. 找出所有 start >= end 的行 awk -F'\t' '$3=="gene" && $4 >= $5 {print}' annotation.gtf | head # 2. 找出链方向不是 + 或 - 的行 awk -F'\t' '$3=="gene" && $7 != "+" && $7 != "-" {print}' annotation.gtf | head # 3. 找出没有 gene_id 的行 awk -F'\t' '$3=="gene" && $9 !~ /gene_id/ {print}' annotation.gtf | head

以上每一组如果返回空，说明原始的 gene 行在这三个基本规则上是干净的。但如果返回了行，就去对照该文件的原始注释看是不是有特殊记录。虽然理论上 GTF 规范不允许这些情况，但真实文件里偶尔会出现来源不明的 malformed 行，宁肯早期发现也不要留到下游。

对生成好的坐标表，我也会检查重复键和可疑的基因名：

# 检查 gene_id 是否有重复 cut -f2 human_gene_coords.tsv | sort | uniq -d | head # 检查基因名是否有空值 awk -F'\t' '$3 == "" {print}' human_gene_coords.tsv | head

这一步看似平凡，但在实际项目中救过我很多次。

5.3 我平时一定会顺手做的三件小事

最后分享几个我自己的固定习惯，很难写进任何正式文档，但确实帮我避过不少雷。

第一，保存注释文件的来源版本和校验值。GTF 坐标必须绑定具体的参考基因组版本才有意义，同一个基因在旧版本和新版本的坐标可能差出几百上千个碱基。我会随手用md5sum annotation.gtf算个校验值存进说明文件，同时记录文件名称、下载日期和版本信息。以后任何人问我“这张基因坐标表是哪来的”，我都能在五分钟内给出准确回答。

第二，做一次“回环验证”。拿一个自己信得过的基因集合（比如从公共资源挑几个坐标非常明确的基因），把它们的坐标从表里抽出来，再去参考基因组注释里比对一下是否一致。或者更直接一点，把生成好的基因 BED 和原始 GTF 的 exon 行做一次 intersect，检查是不是每个基因都能覆盖到自己的外显子。如果某个基因的 BED 区间里完全碰不到外显子，那说明这个基因的坐标可能聚合错了。

第三，挑两个基因名，去基因组浏览器上肉眼对一下。别笑，这个方法虽然原始，但非常有效。我会随机挑一个正链基因、一个负链基因，把表里的 chr、start、end、strand 输进去看染色体位置，然后跟文件里的注释比对。肉眼确认过一次以后，对整张表的信心是完全不一样的。

我现在做任何基因坐标表，都会在脚本里顺手输出一版带 biotype 统计的 summary 文件，再附上来源版本、过滤条件和校验值。这个习惯曾经帮我抓出过一次因为旧版注释文件和新版参考基因组混用导致的系统性坐标偏移，当时只差一点点就拿着错误坐标跑完了全流程。坐标提取本身不难，难的是每一步都清楚自己在干什么、为什么这么干。把这篇文章里的细节都过一遍，你产出的那张基因坐标表，就可以放心交出去了。