测序公司给的分析报告中为什么没有Ensembl GTF注释版本?
导读:
在做RNA-seq等测序分析时,通常会用到Ensembl的GTF基因注释文件,它有不同的版本。我们调研了十几家测序公司的报告,令人惊讶的是,在报告中明确说明GTF版本的公司几乎没有,甚至个别公司连基因组版本都没有明确说明。本文以human为例,比较了不同版本GTF所包含的基因情况,以帮助大家更好地理解基因注释的变化。
GTF注释数据库
在进行RNA-seq定量分析时,常用的定量软件,例如featureCounts,HTSeq等都会要求输入一个GTF基因注释文件。不同的数据库和研究机构(如NCBI、Ensembl、GENCODE等)使用不同的算法和标准进行基因注释,因此即使是同一基因组版本,不同机构提供的注释文件也可能存在差异。
Ensembl数据库GTF不同版本基因的变化
随着科技的进步和研究的深入,对基因组的认识不断更新,新的基因、转录本和其他基因组特征的发现需要反映在注释文件中,因此会定期更新GTF注释文件以包含最新的基因信息。Ensembl GTF约3-6个月更新一次,每次更新都会删除部分基因或者添加部分基因。因此,不同的GTF版本中包含的基因数不一样。
以人类hg38基因组的GTF注释为例,从2014年7月的v76到2024年8月的v113版,前前后后总共有38个版本,基因数(包括protein_coding,lncRNA等)从58051个(v85)到78932个(v113),极差为20881个基因。
数据分析时,GTF版本不匹配带来的问题
由于我们查基因一般都是在NCBI或Ensembl等数据库在线查询,在线版一般默认都是最新版,如果我们的报告所使用的版本跟在线版不匹配,就会存在基因查不到,或者在线查询结果跟测序报告不一致的情况。好在Ensembl有archive版本,如果明确了基因注释版本,就可以直接使用archive版本在线查询。
同样,在对数据进行重新分析时,使用的GTF版本不匹配,会造成部分基因的缺失。例如GEO数据库中很多RNA-seq上传的都是Ensembl基因的原始count矩阵,在将原始count转成FPKM或者TPM时,必需知道基因的长度,如果GTF版本不匹配,就会出现部分基因长度未知,无法计算该基因的FPKM或者TPM了。
例如,在使用rnanorm将count矩阵转成TPM时,会弹出如下警告信息:
UserWarning: X contains 120 genes
that are not in GTF <_io.TextIOWrapper name='Homo_sapiens.GRCh38.104.chr.gtf'
mode='r' encoding='UTF-8'>. This will result in NaN values for missing genes
in the output
最终生成的结果中,有120个基因的tpm是空值。
120个基因影响不大,约占全部基因的2‰左右。然而,参考上述统计结果,最大的差别可以高达2w个基因,是不是瞬间觉得“非常可怕”?
Ensembl GTF版本预测思路
此问题由来已久,这里给出了一个预测Ensembl GTF版本的思路:
1,下载所有版本的GTF文件
2,提取其包含的基因
3,将原始表达矩阵里边的基因与各个版本的GTF基因取交集
4,若提供的表达矩阵里边的基因是某个版本GTF基因的子集,则表达矩阵使用的GTF就有可能是该GTF版本。
注:这里可能会出现预测多个版本的情况,不影响使用。
Python代码如下(相关文件加好友索取):
因此,在出报告时必需明确GTF注释版本,最好在报告中直接给出分析时所使用的GTF文件。这样可以确保结果的可重复性,同时也为进一步信息挖掘提供方便。
小伙伴们赶紧看下你手上的报告,找一找基因组版本和注释版本吧!
微生信助力高分文章,用户220000+,谷歌学术4400+
共有 0 条评论