我们通过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处理，虽然会有点麻烦，但是值得一试
载入GTF文件
library(GenomicFeatures) txdb <- makeTxDbFromGFF("hg38.gtf",format="gtf")通过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转载请注明出处