本文关键词:GEO基因count转FPKM
说实话,搞生信分析的头一个月,最让人头秃的大概就是拿到GEO数据库那一堆raw数据时的懵逼感了。你下载下来是个压缩包,解压一看,全是count值或者fpm值,但很多时候原始数据根本就没给标准化后的数据,或者给了个奇怪的文件名你根本不知道这是啥。这时候我就在想,既然大家都能从GEO下count值,那能不能自己转成FPKM来用呢?这问题听起来简单,真动手的时候坑多着呢。
首先得搞明白个概念,count值是原始的read count,受测序深度和基因长度的影响极大。FPKM则是把这些影响因素给标准化了。很多人为了省事,看到网上有人写代码能跑,照抄一遍。比如用R语言的edgeR或者DESeq2包。注意啊,DESeq2内部其实用的是自己的标准化方法(median of ratios),虽然它也能输出类似归一化后的值,但那不是严格意义上的FPKM。如果你非要FPKM,还得看你的目的。是为了画热图?还是做后续的差异表达分析?如果是画图,FPKM确实直观,因为它是每百万reads中千碱基长度的转录本丰度,不同样本间可比性稍微好点,至少不用担心某个样本测序深度特别深导致基因表达看起来特别高。
我有个搞科研的朋友,之前做RNA-seq,从GEO下了两个样本的count矩阵,直接拿来跑差异分析,结果发现两个样本的总counts差了近十倍,他还没意识到问题严重性,直接进了后续流程,最后PCA图分得乱七八糟,怎么都解释不通。后来查半天才发现,原来是没做标准化。这时候如果你用FPKM,或者TPM,至少在视觉上和初步分析上能纠正一部分由于测序量不同带来的偏差。
具体怎么转呢?咱们不整那些虚的,直接说逻辑。FPKM的计算公式大概是这样的:Read Count除以(基因长度(kb)乘以 Total Mapped Reads(百万))。在R里,你可以用Deseq2的estimateSizeFactors函数处理一下,但那不是FPKM。要得到纯粹的FPKM,得自己算,或者用特定的包。比如edgeR里面的cpm函数可以算CPM,再手动除以基因长度就能得类似的东西。或者更简单的,用tximport结合Salmon或Kallisto的输出,那些工具本来就能输出TPM和Estimated Counts,要是只要FPKM,稍微改改公式就行。
但这里有个大坑,就是基因长度的取值。不同的基因组注释文件(GTF/GFF),同一基因的“长度”可能不一样。有的算的是外显子总长,有的算的是转录本加权平均长度。你要是随便找个长度矩阵,算出来的FPKM肯定也是歪的。我上次帮一个学生看数据,他就是用的错误的基因组注释长度,导致长基因的表达量被严重低估,短基因被高估。这事儿得格外小心,去UCSC或者Ensembl确认下你用的注释版本,别搞混了。
还有啊,别指望FPKM能解决所有问题。如果你想做样本间的聚类,或者PCA,FPKM对低丰度基因的噪声还是太敏感。有时候用vst转换后的数据效果比FPKM好得多。但我懂大家的心情,导师就要FPKM的热图,你就得给。别争辩,做了再说。
总结一下,GEO下载count转FPKM,核心在于两点:一是确认基因长度数据来源的准确性,二是确认你使用的标准化公式是否真正消除了测序深度的影响。别盲目跟风网上的脚本,得明白每一行代码在干嘛。要是你手里数据实在搞不定,或者懒得自己写代码检查那些繁琐的矩阵乘法,找个人帮你看一眼也好,毕竟踩坑两次比读十遍文档管用。要是你有具体的数据格式搞不定,欢迎随时来问,咱们一起拆解下问题所在。