FPKM，TPM，count基本计算

三种获得基因长度的方法(qq.com)

说明：虽然下面有两种计算基因长度的方法，但是目前看来都没有上面这个好

转换公式

countToTpm <- function(counts, effLen)
{
  rate <- log(counts) - log(effLen)
  denom <- log(sum(exp(rate)))
  exp(rate - denom + log(1e6))
}

countToFpkm <- function(counts, effLen)
{
  N <- sum(counts)
  exp( log(counts) + log(1e9) - log(effLen) - log(N) )
}

fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}

countToEffCounts <- function(counts, len, effLen)
{
  counts * (len / effLen)
}

这里的efflen就是基因长度或者说总外显子长度

如果你要计算的话，切记用apply对每一列，也就是每个样本进行计算,不要直接把矩阵丢进去

apply(dataexp,2,function)

获取外显子的长度方法1

10.计算FPKM和RPKM (qq.com)

library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb<-TxDb.Hsapiens.UCSC.hg38.knownGene
## 定义基因的长度为非冗余exon长度之和
if (F) {
  exon_txdb=exons(txdb) #取出人类基因组中的全部外显子
  genes_txdb=genes(txdb)
  o=findOverlaps(exon_txdb,genes_txdb) # 找出外显子和基因之间的重叠部分
  o
  t1=exon_txdb[queryHits(o)]# 将与gene重叠的exon提取出来
  t2=genes_txdb[subjectHits(o)]
  t1=as.data.frame(t1)
  t1$geneid=mcols(t2)[,1]
  #lapply:历遍列表向量中的每个元素，并且使用指定的函数来对其元素进行处理。返回向量列表。
  # 函数split()可以按照分组因子，把向量，矩阵和数据框进行适当的分组；返回值为列表
  g_l=lapply(split(t1,t1$geneid), function(x){
    head(x)
    tmp=apply(x, 1,function(y){
      y[2]:y[3]
    })
    length(unique(unlist(tmp)))
  })
  head(g_l)
  g_l=data.frame(gene_id=names(g_l), length=as.numeric(g_l))
  save(g_l,file ='hg38_g_l.RData')
}
load('hg38_g_l.RData')
## 下面是定义基因长度为最长转录本长度
if(F) {
  t_l=transcriptLengths(txdb)
  head(t_l)
  t_l=na.omit(t_l)
  t_l=t_l[order(t_l$gene_id, t_l$tx_len, decreasing = T),]
  str(t_l)
  t_l=t_l[!duplicated(t_l$gene_id),]
  head(t_l)
  g_l=t_l[,c(3,5)]
}
 head(g_l)
 library(org.Hs.eg.db)
 s2g=toTable(org.Hs.egSYMBOL)
 head(s2g)
 gl=merge(g_l,s2g,by='gene_id') # merge函数可以实现对两个数据框根据共同的列名来进行连接。

获取外显子的长度方法2

Htseq Count To Fpkm | KeepNotes blog (bioinfo-scrounger.com)

library(GenomicFeatures)
txdb <- makeTxDbFromGFF("hg38.gtf",format="gtf")
exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})

如何优雅的统计基因外显子长度 - 云+社区 - 腾讯云 (tencent.com)

能不能反过来算

很多的时候我们从数据库中得到的不是原始数据count，而是fpkm或者是log2(fpkm+1)，那么我们就想知道能不能反过来进行计算？

这里说一下结论，之前简单推了一下（不知道有没有错），发现是没有办法反过去算的，比如你只知道FPKM，想逆过来算count，按照上面的公式，你会发现，你得先知道Count的内容，比如总的count的数量，才能逆过来计算Count,否则没法计算。这就是死循环。