count转TPM/FPKM实战（GSE229904）

接下来是对学员的答疑部分，学员提了一个问题，他想知道怎么将我们的count值进行标准化转为tpm和fpkm值。我们技能树对这个转换已经介绍过非常多次啦：

• Counts FPKM RPKM TPM CPM 的转化

下面就是用GEO的转录组测序数据进行一下实战，正好GEO数据库对二代高通量测序数据统一进行了自己的处理，每一个数据都有count值，tpm值以及fpkm值。

我们本次实战选的数据为：GSE229904，.cgi?acc=GSE229904

数据下载

这个数据做的是前列腺癌患者原发肿瘤的mRNA和miRNA表达谱，从这里点进去：

将这个页面的 /?acc=GSE229904 下列数据下载下来：

得到基因长度

这里数据库直接提供了配套定量的基因长度文件，不同的参考基因组版本计算出来的长度会有区别，所以我们这里就用Human.GRCh38.p13.annot.tsv.gz

代码语言：javascript代码运行次数：0运行复制

# NCBI的注释文件
anno <- fread("GSE229904/Human.GRCh38.p13.annot.tsv.gz",data.table = F)
colnames(anno)
anno <- anno[,c("GeneID","Symbol","EnsemblGeneID","Length","GeneType")]
save(anno,file = "GSE229904/GRCh38.p13.annot.Rdata")
head(anno)

image-20250328155407494

image-20250328155407494

count转TPM

TPM的标准化方式如下：

先读取count：

代码语言：javascript代码运行次数：0运行复制

############################### count2TPM
# 1.读取count
counts <- fread("GSE229904/GSE229904_raw_counts_GRCh38.p13_NCBI.tsv.gz",data.table = F)
counts[1:4,1:4]
rownames(counts) <- counts[,1]
counts <- counts[,-1]
counts[1:4,1:4]

#            GSM7179964 GSM7179965 GSM7179966 GSM7179967
# 100287102         14          8          4         10
# 653635           986        597        746       1398
# 102466751         38         24         22         57
# 107985730          0          1          0          4

# 2.提取基因长度列
load("GSE229904/GRCh38.p13.annot.Rdata")
head(anno)
rownames(anno) <- anno$GeneID
effLen <- anno[row.names(counts), "Length"]
range(effLen)

接下来定义两个函数，这里有两个不同的写法：

来自：What the FPKM? A review of RNA-Seq expression units | The farrago (wordpress)

代码语言：javascript代码运行次数：0运行复制

# 转化
Counts2TPM <- function(counts, effLen){
  rate <- log(counts) - log(effLen)
  denom <- log(sum(exp(rate)))
  exp(rate - denom + log(1e6))
}

tpm_raw <- apply(counts, 2, Counts2TPM, effLen = effLen)
tpm_raw[1:5,1:5]

# 与官方的tpm比较
# 3.与原文件相比
tpm <- fread("GSE229904/GSE229904_norm_counts_TPM_GRCh38.p13_NCBI.tsv.gz",data.table = F)
tpm[1:5,1:5]

结果一致：

还有一个根据公式写的版本：

代码语言：javascript代码运行次数：0运行复制

#TPM (Transcripts Per Kilobase Million)  每千个碱基的转录每百万映射读取的Transcripts 
counts2TPM2 <- function(count=count, efflength=efflen){   
  RPK <- count/(efflength/1000)   #每千碱基reads (reads per kilobase) 长度标准化   
  PMSC_rpk <- sum(RPK)/1e6        #RPK的每百万缩放因子 (“per million” scaling factor ) 深度标准化   
  RPK/PMSC_rpk                       
}

tpm_raw2 <- apply(counts, 2, counts2TPM2, efflength = effLen)
tpm_raw2[1:5,1:5]

# 3.与原文件相比
tpm[1:5,1:5]

结果一致：

count转FPKM

FPKM的公式如下：

接下来定义两个函数，这里有两个不同的写法：

来自：What the FPKM? A review of RNA-Seq expression units | The farrago (wordpress)

代码语言：javascript代码运行次数：0运行复制

# 转化
countToFpkm <- function(counts, effLen) { 
  N <- sum(counts) 
  exp( log(counts) + log(1e9) - log(effLen) - log(N) ) 
} 
fpkm_raw <- apply(counts, 2, countToFpkm, effLen = effLen)
fpkm_raw[1:5,1:5]


# 3.与原文件相比
fpkm <- fread("GSE229904/GSE229904_norm_counts_FPKM_GRCh38.p13_NCBI.tsv.gz",data.table = F)
rownames(fpkm) <- fpkm[,1]
fpkm <- fpkm[, -1]
fpkm[1:5,1:5]

结果一致：

还有一个根据公式写的版本：

代码语言：javascript代码运行次数：0运行复制

##### 函数
# FPKM/RPKM (Fragments/Reads Per Kilobase Million )  每千个碱基的转录每百万映射读取的Fragments/reads 
# RPKM与FPKM分别针对单端与双端测序而言，计算公式是一样的 
counts2FPKM <- function(count=count, efflength=efflen){    
  PMSC_counts <- sum(count)/1e6   #counts的每百万缩放因子 (“per million” scaling factor) 深度标准化   
  FPM <- count/PMSC_counts        #每百万reads/Fragments (Reads/Fragments Per Million) 长度标准化   
  FPM/(efflength/1000)                                       
} 

fpkm_raw2 <- apply(counts, 2, counts2FPKM, efflength = effLen)
fpkm_raw2[1:5,1:5]

# 3.与原文件相比
fpkm[1:5,1:5]

结果一致：

Note：根据gtf获取基因长度

这个例子当中，GEO数据库提供了现成的基因长度，如果是自己的项目呢，一般 featureCounts 软件的结果会提供一个基因长度，在结果的第六列，就可以完成上面的操作了：

除此之外，我们还可以使用定量的时候配套的gtf文件获取基因长度，这里需要注意一定是同一个版本（不同的版本算出来的基因长度会有点不一样）

gtf文件下载：.html

代码语言：javascript代码运行次数：0运行复制

rm(list=ls())
#BiocManager::install("GenomicFeatures")
library(GenomicFeatures)

# 
# gtf <- "/nas1/zhangj/database/genome/Human/release-113/Homo_sapiens.GRCh38.113.gtf"
#txdb <- makeTxDbFromGFF(gtf,format="gtf")
# 上面的图片定量使用的版本是GRCh38 release 104
txdb <- makeTxDbFromGFF("GSE229904/Homo_sapiens.GRCh38.104.chr.gtf.gz",format="gtf")
txdb

# 获取每个基因id的外显子数据
exons.list.per.gene <- exonsBy(txdb, by="gene")

# 对于每个基因，将所有外显子减少成一组非重叠外显子，计算它们的长度(宽度)并求和
exonic.gene.sizes <- sum(width(GenomicRanges::reduce(exons.list.per.gene)))

# 得到geneid和长度数据
gfe <- data.frame(gene_id=names(exonic.gene.sizes), length=exonic.gene.sizes)
head(gfe)

长度如下，与上面的图片经过了验证，得到的长度与 featureCounts 软件的结果一致：

代码语言：javascript代码运行次数：0运行复制

                        gene_id length
ENSG00000000003 ENSG00000000003   4536
ENSG00000000005 ENSG00000000005   1476
ENSG00000000419 ENSG00000000419   9276
ENSG00000000457 ENSG00000000457   6883
ENSG00000000460 ENSG00000000460   5970
ENSG00000000938 ENSG00000000938   3382

你学会了吗？

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-03-28，如有侵权请联系 cloudcommunity@tencent 删除数据库count函数软件数据