TCGA新版数据count的下载及转换为tpm

最近想重新进行免疫浸润计算，用的是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 的转化-生信技能树

基因长度之多少-生信技能树

如有不对，请多批评指正。