最近想重新进行免疫浸润计算,用的是TCGA的数据。然后注意到cibersort算法对于测序数据用的是tpm值。
我以前是直接下载小洁老师存在网盘中的tpm的Rdata,然后我发现我分析的LAML数据总共是151例样本,但是小洁老师不知道为什么漏了一例,只有150例。强迫病犯了,少一例不能忍,所以就自己下载TCGA官网的数据,并进行count到tpm的数据转换。
说明:该安装的包自行安装,99.9999%的代码来自于生信技能树公众号和数据挖掘班上课内容。我不具备生产代码的能力,只有按照自己实际需要改造代码的能力。
第一步,下载TCGA_LAML官网count数据下载
### 1.下载并整理表达矩阵
##共151个sample,60000+基因
rm(list = ls())
library(TCGAbiolinks)
library(SummarizedExperiment)
library(stringr)
query = GDCquery(project = "TCGA-LAML",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "STAR - Counts")
GDCdownload(query)
dat = GDCprepare(query)
laml_count = assay(dat)
laml_count[1:4,1:4]
save(laml_count,file = "Rdata/LAML_count.Rdata")
第2步,获取基因有效长度。
首先需要去NCBI官网下载基因注释文件,网址在这:
GDC Reference Files | NCI Genomic Data Commons
下载这个,放到工作目录下面
###获取基因有效长度
rm(list = ls())
library(parallel)
cl <- makeCluster(0.75*detectCores())
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("gencode.v36.annotation.gtf.gz",
format="gtf")
exons_gene <- exonsBy(txdb, by = "gene")
head(exons_gene)
exons_gene_lens <- parLapply(cl,exons_gene,function(x){sum(width(reduce(x)))})
exons_gene_lens[1:10]
geneid_efflen <- data.frame(geneid=names(exons_gene_lens),
efflen=as.numeric(exons_gene_lens))
dim(geneid_efflen)
save(geneid_efflen,file = "Rdata/gene_efflen.Rdata")
总共是60660行基因,和我们TCGA下载的表达矩阵的的基因是一一对应的,如果数量不对,那后面的所有结果都会有偏差,因为tpm的计算涉及样本测序深度的标准化,需要用到样本所有基因count数求和进行标化。
第三步:count转化为tpm
###count转tpm
rm(list = ls())
load("Rdata/gene_efflen.Rdata")
load("Rdata/LAML_count.Rdata")
identical(geneid_efflen$geneid,rownames(laml_count))
###这里需要注意表达矩阵的id和基因长度的id是否对应,运行结果必须为True!
counts2TPM <- function(count, efflength=geneid_efflen$efflen){
RPK <- count/(efflength/1000) #每千碱基reads (reads per kilobase) 长度标准化
PMSC_rpk <- sum(RPK)/1e6 #RPK的每百万缩放因子 (“per million” scaling factor ) 深度标准化
RPK/PMSC_rpk
}
laml_tpm <- apply(laml_count,2,counts2TPM)
save(laml_tpm,file = "Rdata/LAML_tpm.Rdata")
到这里,count和tpm数据都下载保存好了,后续就是id的转化,以及去除重复值。
还特意去比对了一下和小洁老师计算的tpm有无差异,发现一模一样。
上述内容参照自以下三篇公众号文章:
以Counts转TPM为例, 感受一下生信自学中阳性对照的重要性-果子学生信
Counts FPKM RPKM TPM CPM 的转化-生信技能树
如有不对,请多批评指正。