做网站的大公司有哪些,小程序模板商城,系部网站建设中期检查表,安装wordpress数据库连接时出错在正式分析之前#xff0c;对于数据的处理是至关重要的#xff0c;这种重要性是体现在很多方面#xff0c;其中有一点是要求分析者采用正确的数据类型。
对于芯片数据#xff0c;原始数据进行log2处理之后可以进行很多常见的分析#xff0c;比如差异分析、热图、箱线图、…在正式分析之前对于数据的处理是至关重要的这种重要性是体现在很多方面其中有一点是要求分析者采用正确的数据类型。
对于芯片数据原始数据进行log2处理之后可以进行很多常见的分析比如差异分析、热图、箱线图、PCA分析、生存分析、模型构建聚类分析和相关性分析等。
对于转录组数据在上述的常见分析中只有差异分析时需要采用count值其他的分析是需要采用log2后的cpmtpmfpkmrpkm数据。
除了上述的常规分析在使用不同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,formatgtf)
# 获取每个基因id的外显子数据
exons.list.per.gene - exonsBy(txdb,bygene)
# 对于每个基因将所有外显子减少成一组非重叠外显子计算它们的长度(宽度)并求和
exonic.gene.sizes - sum(width(GenomicRanges::reduce(exons.list.per.gene)))
# 得到geneid和长度数据
gfe - data.frame(gene_idnames(exonic.gene.sizes),lengthexonic.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 或 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.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 -
