Htseq Count To Fpkm

我们通过HTseq-count对hisat2比对后的bam文件进行计数后,会得到每个基因上比对上的reads数,也就是通常所说的count数。接着如果需要比较不同样本同个基因上的表达丰度情况,则需要对count数进行标准化,因为落在一个基因区域内的read counts数目一般可以认为取决于length of the gene(基因长度)和sequencing depth(测序深度),所以引申出两种标准化方法:RPKM和FPKM。前者是以每个reads作为一个单位,在单端测序中应用较多;而后者是以fragment作为一个单位,主要应用在双端测序后的分析。

除了上述两种标准化方法外,最近几年开始应用TPM(Transcripts Per Million)较多,Transcripts Per Kilobase of exonmodel per Million mapped reads (个人觉得从翻译上不太好理解,总而言之是先对reads进行基因长度矫正,然后再除以所有矫正后的reads的总和),认为优化的RPKM计算方法,可以用于同一物种不同组织的比较。

但下面我主要还是以FPKM为例,如何从Htseq-count的count数计算FPKM(其实TPM也是类似的)

FPKM (Fragments Per Kilobase Million)的定义:Fragment Per Kilobase of exon model per Million mapped reads(每1百万个map上的reads中map到外显子的每1K个碱基上的Fragments个数)。如果一对paired reads都比对上了,那么这对reads当做一个fragments;如果paired reads中一个比对上,另外一个没比对上,那么比对上的那个reads当做一个fragments。

由于FPKM与RPKM的唯一差别在于前者在reads map上的情况下只计数1,而后者会计数2;所以两者的公式其实是一样的:

FPKM=  read counts / (mapped reads (Millions) * exon length(KB))

在转录分析时,如果是有参分析,一般使用htseq-count计数后是没有FPKM值的,需要我们通过公式来计算;如果是无参分析的话,使用RSEM定量会给出FPKM值和TPM值。

从公式上可看出,我们可以从htseq-count结果文件中获取每个基因上的read counts数;但后面两者参数mapped reads和exon length就需要自行提供了。

我在网上搜过不少信息,对于这两个参数的解读真是版本众多:

对于mapped reads这个参数而言,一些人认为是clean reads,也就是比对前的clean reads总数;但大多数人还是定义为有效的reads,即mapped reads。那么对于htseq-count结果而言则就是所有基因上计数后的reads总数

对于exon length这个参数而言,一些人粗暴地理解为gene length(是因为好获取,容易计算吗?);但一般人还是理解为所有exon的长度总和。那么怎么获取所有外显子的长度和呢,一般可以从gtf文件入手,但是gtf文件中的exon的区域很多都是重合的(不同转录本的可变剪切)

针对上述问题,可以用GenomicFeatures包来解决;当然也可以自己写脚本对gtf处理,虽然会有点麻烦,但是值得一试

  1. 载入GTF文件

    library(GenomicFeatures)
    txdb <- makeTxDbFromGFF("hg38.gtf",format="gtf")
    
  2. 通过exonsBy获取每个gene上的所有外显子的起始位点和终止位点,然后用reduce去除掉重叠冗余的部分,最后计算长度

    exons_gene <- exonsBy(txdb, by = "gene")
    exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
    

获得每个基因下所有外显子的总长度后,就可以利用上述公式计算FPKM了。。

其实这篇我主要想介绍GenomicFeatures包的,使用学习包之前需要了解TxDb对象,因为这个包主要是围绕着TxDb对象来操作的。

其不仅可以提取上面所说的外显子的位置信息,当然还可以转录本的位置信息,启动子信息等等

其不仅可以读取gtf文件,还可以根据需求使用UCSC Genome Bioinformatics数据

总而言之其主要是用来retrieves and manages transcript-related features

具体可以参看
http://bioconductor.org/packages/GenomicFeatures
http://bioconductor.org/packages/release/bioc/vignettes/GenomicFeatures/inst/doc/GenomicFeatures.pdf

本文出自于http://www.bioinfo-scrounger.com转载请注明出处