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.14231

4、加载既往处理好的数据进行对比

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 或 TPMgeneLength = 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.000000head(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.00000

IOBR包中的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 -