admin管理员组文章数量:1033969
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
# 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
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函数软件数据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
# 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
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函数软件数据本文标签： count转TPMFPKM实战（GSE229904）
版权声明：本文标题：count转TPMFPKM实战（GSE229904） 内容由热心网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://it.en369.cn/jiaocheng/1748095296a2251694.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论