主讲老师第十六课:Count、TPM、FPKM数据转换
在转录组测序分析中,有三个经典的数值,count,FPKM 以及 TPM 值,在 TCGA 数据库中,其提供了count 和 FPKM 两种结果形式,而平时的分析过程中,FPKM 和 TPM 往往是我们比较常用的数据标准化方法。首先,我们来简单看一下三者的基本概念。
count:原始测序得到的 count 数就是比对到某个基因 i上的总数目。不知道大家是否了解测序的简单过程?在测序分析过程中,我们首先是将测得的短 reads 比对到参考基因组上,然后通过软件来计算该片段上比对到 reads 的数量,也就是说呢,count 是一个整数值FPKM:我们把比对到的某个基因 的 Fragment 数目,除以基因的长度,其比值再除以所有基因的总长度,最终得到 FPKM 值。这里需要注意一点,严格来讲,这里的基因长度是指基因外显子的总长度TPM:与 FPKM 不同的地方在于,其基因的比值是再除以(基因的总数目/基因的总长度),因此,其得到的结果是一个相对的比值而且,从定义里面也可以看到,这几个值主要是针对编码基因而言的。比较三者的定义,我们可以发现,FPKM 和 TPM 两种标准化方法的计算公式,其分子是完全相同的,唯一的区别在于对于分母处的处理方式。如果已知 FPKM 的话,那么 TPM 的值就是可以通过 FPKM除以 FPKM 值的总和,再乘以 10 的 6 次方而得到。相对来讲,就标准化程度而言,FPKM值是不如TPM 值的,因此在后续的分析过程中我们一般是推荐使用 TPM 值的下面,我们将介绍一下如何使用 R 语言来进行 count,FPKM 和 TPM 三者之间的转换1 . 基因长度的计算首先,我们来看一下前面提到的基因 Length 的计算方法,这部分基因长度的计算,大家熟悉即可。我已经算好把结果给大家,后面直接用结果即可以了。相信大家必然听说过可变剪切的概念,也正是因为可变剪切的存在,同一个基因会产生不同的转录本。所以呢,又会产生两种不同的分析思路:思路 1:计算基因在染色体的起始和结束之差;思路 2:计算每个基因的最长转录本或所有外显子之和。
接着,如何去获取基因在染色体上的位置信息呢?
如果基因组测序或分析有些了解的小伙伴应该对这几个文件类型有所接触:fasta 文件是保存基因 DNA或 RNA 的测序信息;gtf、gff 以及 gf3 文件则都是保存基因注释的文件。我们来一起看一下基因注释文件里面所包含的内容。
在该文件中,我们可以看到,这是一份人类参考基因组注释文件,也是 TCGA 数据库测序文件的参考注释文件,GRCH38,版本是 22,存在于 GeneCode 网站上。里面包含了基因的染色体位置,起始和终止位置。GRCH38 是人类基因注释文件名字。根据思路 1,一个基因的染色体上位置相减即为基因长度,同时,+代表其为正向,-则代表为反向。如果基因是反向的,其长度需要后面位置减去前面位置,文件往后是基因注释的 id 号,比如 Ensemble号,symbol 号等具体信息。接着,我们来演示一下如何用 R 来进行计算,这里使用的 R 包是 biomaRt 包。
加载完这个 R 包后,通过 listMarts() 函数查看可以连接到的 BioMart 数据库首先来查看一下基因组的参数,在此我们使用的是 Ensemble id 号。
通过 useMart() 函数设置要连接的 BioMart 数据集。
随后,这里以人类基因组为例,获取基因组信息。
注意不同的物种之间是不一样的, getBM() 函数是检索数据用到的主要函数。
首先需要对它的 4 个参数 (filters, attributes, values, mart) 及参数选项进行查看和选择,然后,我们从输入数据中提取基因名,接着,我们对每个基因计算有效长度 eff_length 。
取染色体上位置相减的绝对值作为每个基因的长度,最终得到了第 1 种思路的基因长度结果。
而对于第 2 种思路来讲,其存在不同的转录本,同一个基因就可以存在多个转录本,分别相减取其中最长的一个,或者将这些外显子分别计算,再进行相加从而得到结果接着,导入 GTF 或者 GFF3 文件,ensembl 或者 gencode 网站 GTF 注释皆可。
读取 gtf 文件,并通过 exonBy() 函数提取基因的外显子信息接着,使用 lapply()函数来计算每个外显子的长度,同时,我们进一步通过 reduce() 函数来避免重复计算重叠区并且,对生成的基因长度结果赋予 geneID,即 ensemble 编号,这样一份完整的基因长度文件就准备完成了,此时,得到了两份不同思路计算得到的基因长度文件。在此,推荐大家在后续 的分析过程中,使用第二种方法得到的基因长度文件,而且,因为每个基因的长度是基本固定的,所以大家在后续分析时不需要再自己进行计算,直接用就可以
2. Count 值转换成 FPKM 值
随后,我们来对示例数据中的 count 值进行转换结果显示,在 count 文件中,一共包含了 6 例样本,56830 个不同的基因表达。
接着,将之前准备好的基因长度文件进行读取。
这里,我们选取第二种方法计算得到的基因长度文件。
在正式的转换之前,我们需要将表达文件和基因长度文件做一个预处理。
使用 intersect()函数取共同基因,然后取对应的子集,使得两者的基因名字和排序一致。
最终一共得到了 56830 个共同基因。
根据定义,设定一个用于转换的自定义函数。
结合 apply()函数,对其直接进行转换,最终得到了每个基因对应的 FPKM 值。
3.FPKM 值转 TPM 值
接下来,同样的,设置了 FPKM 值转换成 TPM 值的函数同样的,使用 apply()函数进行转换。
这样,count 值,FPKM 值和 TPM 值之间的转换就完成了。在单基因分析模式中,一般还是推荐使用
TPM 值进行后续的分析计算