1. GFFcompare 到底在解决什么问题做过 RNA-seq 转录本组装的人大概都经历过这样一个阶段用 StringTie、Cufflinks 或者自己写的流程跑出一堆 GTF 文件打开一看里面全是MSTRG.1.1、TCONS_00000001这样的编号完全不知道它们跟已知基因是什么关系。GFFcompare 就是在这个环节里干活的工具它拿你组装出来的转录本去跟参考注释做比对告诉你哪些转录本是真的、哪些是新的、哪些是碎片。核心关键词只有一个给转录本组装结果打分并归类。它的官方定位是 compare, combine and annotate——比较、合并、注释。也就是说它不只是一把尺子还是一台搅拌机。你可以把多个样本的 GTF 一起丢进去它会做跨样本的位点追踪tracking输出一个统一的、带参考注释的分类结果。对做差异可变剪接、做新转录本发掘、做注释质量评估的人来说这是绕不过去的一环。适合谁看已经跑过上游组装流程、手上有一个或多个 GTF 文件的人做基因组注释质检的人做多组学整合时需要对注释做标准化的人。小白也不用怕本文会从分类编码这种最基础的概念讲起一直讲到多文件合并和参数调优。读完你至少能做到三件事看懂.stats里那一堆百分比是什么意思知道-C -M -A这几个参数什么时候该加什么时候不该加遇到结果不对劲的时候知道从哪几个方向去查。GFFcompare 的前身是 Cufflinks 套件里的 Cuffcompare。当年 Cufflinks 那一套工具链Cufflinks、Cuffcompare、Cuffmerge、Cuffdiff在 2010 年前后几乎是 RNA-seq 分析的标准配置后来 StringTie Ballgown 的组合逐渐取代了它。Cuffcompare 的代码被独立出来重写就成了现在的 GFFcompare作者还是原来那批人Geo Pertea 等。相比 CuffcompareGFFcompare 在 GTF 解析的容错性、内存占用、对注释属性的处理上都做了不少改进比如支持 GFF3 的gene_id/transcript_id解析、支持-i列出输入文件、支持把结果重新标注回 GTF。很多人第一次用的时候会误以为它就是一个简单的比对工具跑完拿到.stats就完事了。其实它输出的五六个文件各有用途.tmap是做新转录本筛选的关键.tracking是跨样本一致性的核心.annotated.gtf是能直接拿去下游分析的成品。搞清楚每个文件干什么比背参数重要得多。2. 分类编码体系读懂那些字母的含义2.1 十一种 class code 的完整含义GFFcompare 最有价值的产出之一就是给每个输入转录本打一个 class code。这个字母告诉你它和参考注释之间的关系。我把常用的几种整理成表方便对照。代码名称含义说明完全匹配内含子链与某个参考转录本完全一致就是同一个转录本c包含输入转录本被某个参考转录本完全包含k反向包含参考转录本被输入转录本完全包含方向相反m部分链匹配多外显子至少一个内含子与参考匹配但整体不完全一致n内含子保留与参考相比存在内含子保留j多外显子匹配多外显子转录本至少有一个剪接位点与参考一致e单外显子单外显子转录本与参考外显子有部分重叠o其他重叠同链外显子有重叠但不满足上面任何一种s内含子反义内含子区域与参考内含子重叠且链方向相反x外显子反义外显子区域与参考外显子重叠但链方向相反i内含子内整个转录本落在参考的内含子区间里y含参考输入转录本的内含子里包含了参考转录本p可能聚合位于参考转录本下游 2kb 内的可能聚合酶延伸片段r重复序列至少 50% 的碱基被重复序列覆盖u未知基因间区与参考没有任何交集和j是公认的好成绩。意味着你组装出来的这个转录本跟已知的某个转录本一模一样属于被重现j意味着剪接模式部分对上了可能是同一个基因的不同异构体也可能是组装断点不够准。c、k、m、n这几个是部分匹配的范畴通常出现在组装不完整或者参考注释本身就有冗余的情况下。u是新转录本最直接的候选但要注意u里混杂了大量基因间区的噪声和未注释的非编码转录本。2.2 为什么和j是最关键的两个指标在评估一个转录本组装流程好不好用的时候我一般先看两个数字的占比和j的占比。高说明这个样本里主要表达的都是已知转录本组装流程的重现能力没问题j高但低说明剪接点找得对但外显子边界或者转录本起止位置还有偏差这时候要回去检查测序深度和组装参数。u的数量要辩证看。如果你做的是肿瘤样本或者某个发育阶段特有的组织出现大量u是正常的可能真有大量未注释转录本但如果你做的是常规细胞系u占比超过三成那大概率是组装流程里有噪声没滤干净比如低表达转录本被强行拼出来了。x和s是反义转录本的信号。这两个类别在链特异性建库的数据里出现频率应该很低非链特异建库里会明显偏高因为链信息本身就不可信。如果你用非链特异数据跑出来一堆x基本可以忽略当作链归属错误处理就行。p类特别容易被忽略。它指的是位于参考转录本下游 2kb 范围内的转录本很可能是聚合酶没有正常终止导致的通读产物。做新转录本筛选时这类通常要剔除否则你会得到一堆假阳性的新基因。可以通过-d参数调整这个距离阈值。2.3 分类判定背后的比对逻辑GFFcompare 判定的核心是内含子链intron chain。它先把每个转录本转换成一组内含子区间然后拿两组内含子区间去比对。如果内含子链完全一致就是如果只有部分一致就往j、m这些类别上靠。这个设计有个隐含前提它认为内含子位置比外显子边界更可靠。这个假设在大多数情况下成立因为剪接位点有明确的序列特征GT-AG 规则而外显子边界在组装时容易因为覆盖度不足而漂移。所以 GFFcompare 用内含子链做主判据外显子边界只作为辅助。理解这一点对调参很重要。如果你发现组装结果里外显子边界总是差几个碱基那不是 GFFcompare 的问题是上游组装的问题。GFFcompare 只负责告诉你差在哪不负责帮你修。另外判定时会考虑一个距离容忍度由-d参数控制默认 100bp。也就是说如果两个转录本的剪接位点相差在 100bp 以内可能被判为匹配。这个默认值在多数物种上够用但对剪接位点特别密集的基因比如一些神经系统的基因可能会误判可以适当调小。3. 环境准备与最小可用命令3.1 安装方式与版本选择最省事的装法是用 condaconda install -c bioconda gffcompare这样会自动把依赖处理干净。如果你所在的环境不允许联网或者需要指定版本可以从源码编译git clone https://github.com/gpertea/gffcompare.git cd gffcompare make release编译需要 gclib 这个依赖库make release的时候会自动去拉取。注意编译产物是一个单独的二进制文件gffcompare可以直接拷到别的机器上用不依赖运行时库这一点很方便。版本选择上我建议用 0.12 之后的版本。0.12 之前的版本在处理 GFF3 格式的属性字段时有一些已知问题特别是gene_id和transcript_id的解析。可以这样确认版本gffcompare --version如果输出里显示的版本号低于 0.12建议升级。0.11.x 和 0.12.x 在.tmap文件的列结构上也有一点差异如果你要跟别人的结果做对比版本最好保持一致。3.2 输入文件格式要求GFFcompare 接受的输入是 GTF 或 GFF3。这两种格式虽然都是注释文件但字段结构不一样处理时要注意。GTF 格式每行是seqname source feature start end score strand frame attributes属性部分是key value;的形式。GFF3 的属性部分是keyvalue;的形式。GFFcompare 能自动识别但有个坑如果 GTF 里同一行的属性字段格式不规范比如引号缺失、分号缺失解析会静默失败那个转录本就不会出现在结果里。所以跑之前最好做一次格式校验。参考注释文件用-r指定必须是同源基因组的注释。如果你的比对用的是 GRCh38参考注释也必须是 GRCh38 的版本不能混用。参考注释可以是 GTF 也可以是 GFF3但要注意参考注释里的链信息必须完整否则分类会不准。输入的 query 文件数量不限一个两个都行。多个文件会一起处理GFFcompare 会做跨文件的位点追踪。输入文件里可以没有gene_idGFFcompare 会自动根据位点位置生成XLOC_编号。但如果有gene_id它会尽量保留。# 检查文件是否有明显的格式问题 grep -c $\t sample.gtf awk -F\t NF!9 sample.gtf | head第一条命令数一下有多少行第二条命令找出列数不是 9 的行。如果第二条有输出说明格式有问题得先修。3.3 最常用的命令模板一个最基础的命令长这样gffcompare -r reference.gtf -o cmp_out sample1.gtf-r指定参考-o指定输出前缀最后是 query 文件。输出会在当前目录生成cmp_out.stats、cmp_out.tracking、cmp_out.tmap、cmp_out.annotated.gtf、cmp_out.loci这些文件。如果只有一个 query 文件其实是在做注释评估如果有多个就是在做注释合并加评估。多个文件的写法gffcompare -r reference.gtf -o cmp_out sample1.gtf sample2.gtf sample3.gtf文件多了以后命令行会很长可以用-i指定一个列表文件# list.txt 里每行一个 GTF 路径 gffcompare -r reference.gtf -o cmp_out -i list.txt实际工作中我用得最多的组合是gffcompare -r reference.gtf -C -M -o cmp_out -i list.txt-C丢弃被参考转录本包含的输入转录本-M只保留多外显子转录本。这两个一起用能把大量碎片化的单外显子噪声去掉让输出结果干净很多。后面会详细讲这两个参数。提示-o指定的前缀如果包含目录路径那个目录必须提前存在GFFcompare 不会自动创建。跑批的时候经常在这里翻车。4. 输出文件逐个拆解4.1.stats文件敏感度和精确度怎么读.stats是大部分人第一个打开的文件。它的开头是命令行记录和数据集摘要中间是一张大的指标表格最后是漏检和新增的统计。摘要部分会告诉你 query 有多少个转录本、多少个位点、其中多少是多外显子以及参考里对应的数字。核心表格分六行Base level、Exon level、Intron level、Intron chain level、Transcript level、Locus level。每一行都是敏感度 | 精确度的格式。这六行的意义是逐级放宽的。Base level 是按碱基算最宽松通常都在 90% 以上。Transcript level 最严格要求整个转录本的外显子结构完全一致才算对一般在 30% 到 70% 之间浮动。Intron chain level 是判断剪接模式对不对通常比 transcript level 高一些。敏感度Sensitivity 参考里被找回的比例衡量的是漏没漏。精确度Precision 预测里正确的比例衡量的是错没错。这两个指标天然存在权衡把-C加上去会让精确度升高但敏感度下降因为你在丢掉一部分预测结果。我个人的经验是同一套数据在两个不同组装流程之间做选择时看 Intron chain level 的敏感度就够了。这个指标最能反映组装的质量而且对参数不敏感。Transcript level 的数字波动太大受参考注释完整度影响明显不适合做横向比较。注意如果 reference 里有很多单外显子转录本而你的 query 里几乎没有那么 Base level 的数字会很难看但这不代表组装差只是覆盖的转录本类型不一样。看指标前先搞清楚两边的转录本构成。4.2.tracking文件跨样本追踪表.tracking是做多样本分析时最有用的文件。每一行对应一个位点locus后面跟着这个位点在各样本里的转录本编号以及它和参考的关系。文件的结构是第一列是统一编号TCONS_开头第二列是位点编号XLOC_开头第三列是参考基因和转录本编号第四列是 class code后面每两列对应一个输入样本。判断一个转录本在多个样本里是否一致表达就看它对应的列里有没有值。如果第 5 列和第 7 列都有值说明它在样本 1 和样本 2 里都被组装出来了。这个信息在做可变剪接分析时特别关键因为只有跨样本存在的转录本才值得拿去做差异分析。我自己常用的一个操作是把.tracking转成矩阵形式行是位点列是样本值是有无。这样能快速筛出所有样本都出现的高置信度转录本。可以用 awk 处理awk BEGIN{OFS\t} $4 || $4j {print $2, $3, $4} cmp_out.tracking confident.txt这条命令只保留 class code 是或j的位点也就是跟参考有明确对应关系的。做注释精修的时候这部分是最可靠的种子。4.3.tmap与.refmap一对一的映射关系.tmap是 query 到 reference 的映射表。每一行是一个 query 转录本如果它在参考里有对应的转录本就会在ref_id列填上那个参考转录本的编号同时给出 class code。如果没找到对应ref_id是-class code 是u。.tmap的列比较多常用的几列是ref_gene_id、ref_id、class_code、qry_gene_id、qry_id、num_exons、TPM、cov、len。做新转录本筛选时我就直接筛 class_code 是u并且 TPM 或 cov 高于阈值的行。.refmap是反向的从参考转录本出发列出哪些 query 转录本匹配到了它。这个文件在做参考注释的覆盖度评估时有用比如你想知道参考里有哪些基因在某个样本里根本没被组装出来看.refmap里哪些行是空的就知道了。4.4.annotated.gtf与.loci.annotated.gtf是把分类信息写进注释属性的 GTF 文件每一行的属性里会多出class_code、ref_gene_id、cmp_ref这些字段。这个文件可以直接拿去用 IGV 或者别的工具可视化一眼就能看出每个转录本是什么类别。.loci文件记录的是位点级别的信息包括位点在基因组上的位置、包含多少个转录本、每个转录本的编号。做位点级别的差异表达分析时会用到。顺便说一下.annotated.gtf里的转录本编号是统一后的TCONS_编号跟原始输入文件里的编号不一样。如果你需要把结果对应回原始的 GTF就用.tracking文件做桥梁。这一点新手经常搞混拿到.annotated.gtf发现编号全变了以为是文件出错了。5. 参数调优实战从原始输出到可用注释5.1-C -M -A三件套的取舍这三个参数是过滤用的用途不同不能无脑全加上。-C丢弃被参考转录本包含的输入转录本。所谓被包含指的是输入转录本的整个范围落在某个参考转录本内部而且内含子链也对得上。这类转录本很可能是参考注释里那个转录本的 3 端截断版本保留下来价值不大。加上-C之后精确度会明显上升但敏感度会掉一点因为你也丢掉了一些真实存在的短异构体。-M只保留多外显子转录本。这个参数影响很大。如果研究的是 lncRNA 或者某些以单外显子为主的基因家族加-M会把它们全滤掉。所以我一般只在做编码基因注释评估的时候加-M做全转录组分析时会先用.tmap里的num_exons列自己筛保留控制权。-A比-C更严格它把跟参考完全没有重叠的转录本也丢掉包括i内含子内、u基因间、p聚合延伸这些类别。如果你只想评估已知基因的组装情况加-A最省事。但如果你要发掘新转录本-A是绝对不能加的它会把你要找的东西全滤掉。实际使用时我的建议是分两步走第一步不加任何过滤跑一遍拿到完整的.tmap第二步基于.tmap用脚本做自定义筛选。这样每次筛选条件变了都不用重跑 GFFcompare效率高很多。5.2-e -d -s的比对尺度控制这几个参数控制的是怎么算匹配。-e是在没有参考注释时使用的它让 GFFcompare 只做外显子级别的比对不做转录本级别的。多用于纯 query 之间的比较比如合并多个样本的组装结果时。-d n控制剪接位点的容忍距离默认 100。调小到 50 会让匹配更严格适合剪接位点已知很精确的物种调大到 200 会让匹配更宽松适合组装质量一般的样本。调整这个参数对 Intron chain level 的指标影响比较明显做横向对比时要保证两边用的同一个值。-s 文件用于指定链特异性比较。传入一个文件里面列出哪些序列需要考虑链方向。非链特异建库的数据加这个参数会让分类更准因为x和s这些反义类别会被正确处理。参数使用上要注意这个文件需要自己准备格式是每行一个序列名。还有一个容易被忽视的参数是-R。当你的 query 本身就是参考注释的一个子集时比如你只取了某条染色体的注释去做测试加-R会让 GFFcompare 知道这一点避免把正常的包含关系误判成c。这个参数在做大规模注释评估的时候很有用。5.3 多文件合并与-i列表模式多文件合并的时候GFFcompare 会自动把所有 query 转录本按位置聚成位点然后统一编号。这时候位点级别的追踪信息就体现在.tracking里。列表模式的写法有个细节-i后面的文件里每行是一个路径路径里不能有额外的空格或者制表符空行也不要有。我一般这样生成列表ls -1 ./gtf/*.gtf list.txt gffcompare -r reference.gtf -o cmp_out -i list.txt列表里文件的顺序决定了.tracking里样本列的顺序这个顺序一定要自己记清楚否则后面解读结果时容易搞混哪个列对应哪个样本。还有个参数-L可以给每个输入文件打标签标签会出现在.tracking的列名里。这样结果文件自解释比记住顺序靠谱gffcompare -r reference.gtf -o cmp_out -i list.txt -L sample这样输出的列名会带上sample前缀加文件名。标签只影响列名显示不影响数据本身。提示合并多个样本时如果某些样本的组装质量差异很大结果会被差样本拉低。我一般会先单独跑每个样本的.stats把质量明显不合格的样本剔掉再合并。合并前的筛选比合并后补救要省事得多。6. 常见问题与排查实录6.1 常见报错与异常速查跑 GFFcompare 遇到问题时大部分情况出在输入文件上。我把踩过的坑整理成一张表。现象可能原因处理方式输出里转录本数量远少于输入属性字段格式不规范解析失败检查 GTF 每行属性是否有引号和分号全部转录本都是u参考注释和比对基因组不匹配确认-r的版本和 FASTA 同源大量x/s类别非链特异建库或链信息缺失确认建库方式必要时加-s参数运行报错找不到序列序列名大小写不一致统一 chr 命名方式chr1vs1.tmap里ref_id全为空参考注释链方向字段缺失检查 reference 的 strand 列是否完整内存占用飙升输入文件过多或注释极度冗余分批处理或用-M先过滤序列名大小写这个问题特别常见。有的注释文件用chr1有的用1有的用NC_000001.11只要对不上所有转录本都会掉到u里。跑之前用这两条命令核对一下cut -f1 reference.gtf | sort -u | head cut -f1 sample.gtf | sort -u | head两边输出对比一下名字体系必须完全一致。6.2 结果指标异常时的排查思路拿到.stats之后如果 Transcript level 的敏感度特别低比如低于 20%先别急着改参数按下面的顺序排查。先看参考注释是不是过于冗余。有些版本的参考注释包含大量短转录本和预测转录本这些在你的组装结果里本来就不该出现会拉低分母。可以统计一下参考里多外显子转录本和单外显子转录本的比例如果单外显子占比很高指标会被稀释。再看 query 的转录本构成。如果 query 里绝大部分是单外显子而参考里以多外显子为主那 Intron chain level 的指标会很低这是结构差异导致的不是组装错误。这时候去看 Base level 和 Exon level 更有意义。最后看是不是浮动阈值设置的问题。GFFcompare 内部对重复序列有个判定r类如果你的物种重复序列注释不完整一部分转录本会被误判为重复而影响统计。这个可以通过调整输入预处理来规避比如在组装前就把重复区域 masked 掉。6.3 几个容易踩的坑第一个坑是混淆了-C和-A。-C只丢弃被包含的转录本-A还会丢弃不重叠的转录本。做新转录本发掘时错加了-A结果一个u都看不到这是新手最常犯的错误。第二个坑是拿.annotated.gtf当原始 GTF 用。这个文件里的转录本编号是重新分配过的如果你后面还有流程依赖原始的编号直接用它会导致对应关系断裂。要么用.tracking做映射要么就用原始 GTF 加.tmap做筛选。第三个坑是参考注释用了未软屏蔽的版本。某些参考注释里包含了重复序列区域会干扰分类。建议用软屏蔽soft-masked版本的注释。第四个坑是在不同版本之间对比结果。前面提过 0.11 和 0.12 的.tmap列结构有差异。如果你跟合作方各跑一半数据然后合并一定要先确认版本一致。我见过因为版本差异导致列错位、结果完全对不上的情况。第五个坑也是我个人吃过亏的多文件合并时没检查单文件质量直接合并。有一次我把一个深度只有 5x 的样本跟两个 50x 的样本合并跑结果 .tracking 里那个低深度样本贡献了大量噪声位点最后筛出来的跨样本一致转录本里有一大半是假的。后来我改成分别跑.stats把 Intron chain 敏感度低于某个阈值的样本先剔掉结果质量立刻上来了。关于参数调优我再补充一点个人体会-d这个参数看着不起眼但对 Intron chain level 的影响比我预想的大。在剪接位点密集的物种比如某些昆虫上默认的 100 会让不少本不该匹配的转录本被判成匹配。我现在的习惯是先用默认值跑一遍看看j类别的数量如果j占比异常高就把-d调到 50 再跑一次做对比。两次结果差异大的话说明这个数据集对距离参数敏感后续分析里要格外小心。最后再分享一个小技巧。如果你需要把分类结果转成更易读的表可以用tmap文件做透视awk BEGIN{OFS\t} $3!- {print $3, $4, $5} cmp_out.tmap \ | sort | uniq -c | sort -rn | head -20这条命令统计每种 class code 出现的次数按频次排序。跑完一眼就能看出你的结果构成比翻.stats快得多。做批量任务的时候我会把所有样本的这个统计结果汇总成一张表作为质检的第一道关卡。