count转fpkm,tpm完整版

参考资料:
1.如何把从TCGA下载的HTseq count转换成FPKM - 简书 (jianshu.com)
2.RNAseq数据分析中count、FPKM和TPM之间的转换-腾讯云开发者社区-腾讯云
3.基因有效长度计算
4.简单小需求：如何将FPKM转换成TPM？

tinyarray包是生信技能树讲师小洁老师写的，里面的很多函数都非常好用，让我们的代码更加简洁易读，让我们感谢小洁老师！

1.处理count矩阵

1.1读取

dat = read.delim("GSE104140_raw_counts_GRCh38.p13_NCBI.tsv.gz",row.names=1)
dat[1:3,1:3]
##           GSM2790523 GSM2790524 GSM2790525
## 100287102         22         29         18
## 653635           599        574        530
## 102466751         32         29         27

1.2行名id转化

行名是Gene ID，把它转化成symbol

dat1 =trans_entrezexp(dat)  #为了防止后面不能查看最原始的数据，还是再建一个变量吧。

1.3基因过滤

表达量为0的基因有点多，我选择去除表达量为0占一半以上样本的基因

nrow(dat1)
## [1] 37734
dat1 = dat1[apply(dat1, 1, function(x) sum(x > 0) > 0.5*ncol(dat1)), ]
nrow(dat1)
## [1] 27615

2.count转化为fpkm

2.1获取基因有效长度

这里的有效长度指的是非冗余外显子长度之和。
（环境内存不够怎么办？有效数据存起来，另开一个session）

if(T){
  if(!require(GenomicFeatures))BiocManager::install("GenomicFeatures")
  library(GenomicFeatures)
  txdb <- makeTxDbFromGFF("gencode.v45.annotation.gtf.gz",format="gtf") 
  exons_gene <- exonsBy(txdb,by="gene")
  exons_gene_lens<-lapply(exons_gene,function(x){sum(width(reduce(x)))})
  length(exons_gene_lens) #查看基因数量
  exons_gene_lens1<-as.data.frame(exons_gene_lens)
  #转换格式：列表转为数据框
  exons_gene_lens1 <- t(exons_gene_lens1)#转置一下
  save(exons_gene_lens1,file = "efflength.Rdata")
}

2.2行名id转化

load("efflength.Rdata")
head(exons_gene_lens1)
##                    [,1]
## ENSG00000000003.16 4530
## ENSG00000000005.6  1476
## ENSG00000000419.14 9276
## ENSG00000000457.14 6883
## ENSG00000000460.17 5970
## ENSG00000000938.13 3382

发现问题：exon是一个只有一列的矩阵，比较特殊，不符合行名转化的示例格式。只好再添加一列，最大程度的模仿示例格式（要不然使用函数会报错）

exons_gene_lens1 = as.data.frame(exons_gene_lens1)
exons_gene_lens1$rnorm <- rnorm(n=nrow(exons_gene_lens1))
exons1 = trans_exp_new(exons_gene_lens1)

3.合并count矩阵和exon长度表

3.1行名对应相等

colnames(exons1)[1] <- "length" #改一下列名
p = identical(rownames(dat1),rownames(exons1));p #行名不一致
## [1] FALSE

s = intersect(rownames(dat1),rownames(exons1)) #取交集
dat1 = dat1[s,]
exons1 = exons1[s,]
table(rownames(dat1) == rownames(exons1)) #检查一下行名是不是对应相等
## 
##  TRUE 
## 21367

3.2合并

合并完成，exon length是新count表的最后一列啦！

count_with_length <- cbind(dat1,length = exons1[,1])

4.count转fpkm

方法一：使用函数（常用推荐）
不合并也可以，但是要行名对应相等。

countToFpkm <- function(counts, effLen)
{
    N <- sum(counts)
    exp( log(counts) + log(1e9) - log(effLen) - log(N) )
}
fpkms <- apply(dat1,2,countToFpkm,effLen = exons1[,1] )
head(fpkms[1:3,1:3])

##           GSM2790523 GSM2790524 GSM2790525
## DDX11L1     1.225466   1.498593   1.073968
## WASH7P     15.280702  13.584247  14.482132
## MIR6859-1  16.566724  13.928099  14.972371

用fpkm转tpm，列加和相等验证一下。

fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
TPM <- apply(fpkms,2,fpkmToTpm)
table(colSums(TPM)) 
## 999999.999999999            1e+06 
##               14               18

方法二：需要先合并

kb <- count_with_length$length/1000
countdata <- count_with_length[,1:32]
rpk <- countdata/kb
fpkms2 <- t(t(rpk)/colSums(countdata)*10^6)
head(fpkms2[1:3,1:3])

##           GSM2790523 GSM2790524 GSM2790525
## DDX11L1     1.225466   1.498593   1.073968
## WASH7P     15.280702  13.584247  14.482132
## MIR6859-1  16.566724  13.928099  14.972371

终于大功告成啦!

启发

不要无脑复制粘贴别人的代码，理解最重要。
毕竟别人的代码不一定能适应你的数据。自己的数据自己去处理，灵活一点。
最靠谱的学习来源：官网和帮助文档。
网上的各种各样的代码最终还是由官方提供的演化来的，那为什么不从源头开始学习呢？
报错怎么处理？
根据R的提示，一步一步排查问题。不符合函数要求的格式？参数没设置对？