参考资料:
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的提示，一步一步排查问题。不符合函数要求的格式？参数没设置对？
