Count数据转换为TPM数据方法整理-常规方法、DGEobj.utils和IOBR包
原创Count数据转换为TPM数据方法整理-常规方法、DGEobj.utils和IOBR包
原创
在正式分析之前，对于数据的处理是至关重要的，这种重要性是体现在很多方面，其中有一点是要求分析者采用正确的数据类型。
对于芯片数据，原始数据进行log2处理之后可以进行很多常见的分析，比如差异分析、热图、箱线图、PCA分析、生存分析、模型构建，聚类分析和相关性分析等。
对于转录组数据，在上述的常见分析中只有差异分析时需要采用count值，其他的分析是需要采用log2后的cpm，tpm，fpkm，rpkm数据。
除了上述的常规分析，在使用不同R包进行分析之前务必浏览一下输入数据的要求。
那么芯片数据还好说，毕竟后续进行log2处理后就可以做很多分析。但是转录组数据的可选项就比较多了。
但目前常用的只有tpm和cpm
count数据转换为cpm数据非常简单
# exprSet是count表达矩阵
# 一句代码搞定
exprSet = log2(edgeR::cpm(exprSet)+1)比较难的就是count数据转换为tpm数据，因此搬运了常规的流程和R包的方法，做个对比
首先要去获取基因长度文件，因为后续需要用这个数据去矫正基因长度。
网址：https://gdc.cancer.gov/about-data/gdc-data-processing/gdc-reference-files
方法一：常规方法-来自生信技能树
1、获取gencode.v36.annotation并处理
if(!require(GenomicFeatures))BiocManager::install("GenomicFeatures")
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("gencode.v36.annotation.gtf.gz",format="gtf")
# 获取每个基因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)[1:5,1:2]
# gene_id length
# ENSG00000000003.15 ENSG00000000003.15 4536
# ENSG00000000005.6 ENSG00000000005.6 1476
# ENSG00000000419.13 ENSG00000000419.13 1207
# ENSG00000000457.14 ENSG00000000457.14 6883
# ENSG00000000460.17 ENSG00000000460.17 5970
save(gfe,file = "gfe.Rdata")2、使用TCGA-HNSC数据
load("hnsc_exp.Rdata")
head(hnsc)[1:3,1:3]
# TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# ENSG00000000003.15 2090 1680 5433
# ENSG00000000005.6 0 0 0
# ENSG00000000419.13 2098 3872 2240
identical(rownames(hnsc),rownames(gfe))
#[1] TRUE 行名是能够对上的使用了TCGA-HNSC中的count数据，检测了一下count数据和下载的基因信息的顺序是一致的。
3、曾老师写的代码进行count/tpm转化
load("gfe.Rdata")
#提取基因长度列
effLen = gfe$length
#转化
Counts2TPM <- function(counts, effLen){
rate <- log(counts) - log(effLen)
denom <- log(sum(exp(rate)))
exp(rate - denom + log(1e6))
}
hnsc_tpm_raw <- apply(hnsc, 2, Counts2TPM, effLen = effLen)
head(hnsc_tpm_raw)[1:3,1:3]
# TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# ENSG00000000003.15 31.13713 15.74248 72.37614
# ENSG00000000005.6 0.00000 0.00000 0.00000
# ENSG00000000419.13 117.46366 136.35307 112.142314、加载既往处理好的数据进行对比
load("hnsc_tpm.Rdata")
head(hnsc_tpm)[1:3,1:3]
# TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# ENSG00000000003.15 31.13713 15.74248 72.37614
# ENSG00000000005.6 0.00000 0.00000 0.00000
# ENSG00000000419.13 117.46366 136.35307 112.14231结果是一致的
方法二：R包 DGEobj.utils
load("gfe.Rdata")
#提取上面处理好的基因长度列
effLen = gfe$length
library(DGEobj.utils)
CC_res <- convertCounts(
exp,
unit = "TPM", # CPM、FPKM、FPK 或 TPM
geneLength = effLen, #这里还是需要下在基因长度数据
log = FALSE, #默认为False，设为TRUE时返回Log2值
normalize = "none", #默认为none，调用edgeR的calcNormFactors()进行标准化
prior.count = NULL #为避免取0的对数，对每个观测值添加平均count。仅当 log = TRUE 时使用，
)
head(CC_res)[1:3,1:3]
# TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
#ENSG00000000003.15 31.13713 15.74248 72.37614
#ENSG00000000005.6 0.00000 0.00000 0.00000
#ENSG00000000419.13 117.46366 136.35307 112.14231与前面得到的分析结果是一致的。
需要下载基因长度数据，但是前期处理完之后后面可以很方便的转化为各种想要的格式(CPM、FPKM、FPK 或 TPM)。
方法三：R包 IOBR
load("hnsc_exp.Rdata")
exp <- hnsc
library(IOBR)
exp_tpm <- count2tpm(
exp, #表达矩阵， 行为基因
idType = "Ensembl", #"Ensembl" "ENTREZ","SYMBOL"
)
head(exp_tpm)[1:3,1:3]
# TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# TSPAN6 31.89988 16.16771 73.85059
# TNMD 0.00000 0.00000 0.00000
# DPM1 15.65888 18.22162 14.88931
# 可以看到这个函数还自动把gene_id给转换成了symbol
# 为了更好的对比，我们也把hnsc_tpm_raw中的gene_id转换成symbol
# 使用小洁老师的trans_exp_new函数
library(tinyarray)
a <- trans_exp_new(hnsc_tpm_raw)
head(a)[1:3,1:3]
# TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# DDX11L1 0.000000 0.0000000 0.000000
# WASH7P 1.150477 0.8180029 2.191641
# MIR6859-1 0.993794 0.0000000 4.443138
# 排个序再看一下
p <- identical(rownames(a),rownames(exp_tpm));p
if(!p) {
s = intersect(rownames(a),rownames(exp_tpm))
a = a[s,]
exp_tpm = exp_tpm[s,]
}
head(a)[1:3,1:3]
# TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# DDX11L1 0.000000 0.00000000 0.000000
# WASH7P 1.150477 0.81800290 2.191641
# FAM138A 0.000000 0.03486849 0.000000
head(exp_tpm)[1:3,1:3]
# TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# DDX11L1 0.00000 0.00000000 0.00000
# WASH7P 22.63327 5.95811827 15.88102
# FAM138A 0.00000 0.03581036 0.00000IOBR包中的count2tpm函数只能进行tpm转化(github上搜了这个R包内容似乎没有转化为其他数据格式的函数了)。用默认的方式进行运算得到的结果存在一定的偏差，而且我个人觉得这个偏差有点大... 但是我暂时不知道是什么原因？是内置的基因长度顺序有问题？还是我某个参数设的不对？求高手指点~
综合上述分析，暂时还是选择常规方法和DGEobj.utils R包中的convertCounts函数吧~
参考资料：
1、https://hbctraining.github.io/Training-modules/planning_successful_rnaseq/lessons/sample_level_QC.html
2、https://rdrr.io/cran/DGEobj.utils/man/convertCounts.html
3、https://rdrr.io/github/IOBR/IOBR/man/count2tpm.html
4、生信技能树推文：https://mp.weixin.qq.com/s/IUV9dSbRBK1nvetixKOCRw
致谢：感谢曾老师，小洁老师以及生信技能树团队全体成员。
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟
- END -
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
