R语言操作——TCGA数据处理

获取表达矩阵，处理TCGA的count数据，1表示为行。

exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ]
dim(exp)
exp[1:4,1:4]
exp = as.matrix(exp)

导入数据

library(rio)
count <- import("GSE168184_GEO_submission_FPKM.txt",
                format = "\t")

加 ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）

library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人类
dim(deg)
library(dplyr)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
dim(deg)

转化空格为NA

clinical[clinical==""] = NA

用花花的专属TCGA包，ID进行转换

rm(list=ls())
load("TCGA-CHOL_gdc.Rdata")
library(tinyarray)
k = tinyarray::trans_exp(exp) 

把空着的值改为NA

#简化meta的列名
colnames(meta)=c('ID','event','death','last_followup','race','age','gender','stage')
#空着的值改为NA
meta[meta==""]=NA

以病人为中心，表达矩阵按病人ID去重复

# 以病人为中心，表达矩阵按病人ID去重复
k = !duplicated(str_sub(colnames(exprSet),1,12));table(k)
exprSet = exprSet[,k]

去除重复

an = gtf_gene[,c("gene_name","gene_id","gene_type")]
exp = exp[rownames(exp) %in% an$gene_id,]
an = an[match(rownames(exp),an$gene_id),]
identical(an$gene_id,rownames(exp))
k = !duplicated(an$gene_name);table(k)
an = an[k,]
exp = exp[k,]

TPM数据做单个基因的生存分析file:///C:/Users/denghuan/Desktop/The%20learning%20of%20R%20software/Practice/%E7%94%9F%E5%AD%98%E5%88%86%E6%9E%90%20survival%20analysis/6.Survival.html

#单个基因
g = "TRIP13"
meta$gene = ifelse(exprSet[g,]>median(exprSet[g,]),
                   'high',
                   'low')
sfit1=survfit(Surv(time, event)~gene, data=meta)
ggsurvplot(sfit1,pval =TRUE, data = meta, 
           risk.table = TRUE)
ggsurvplot(sfit1,palette = c("#E7B800", "#2E9FDF"),
           risk.table =TRUE,pval =TRUE,
           conf.int =TRUE,xlab ="Time in months", 
           ggtheme =theme_light(), 
           ncensor.plot = TRUE)
##方法2:生信自学网
diff=survdiff(Surv(time, event)~gene, data=meta)
pValue=1-pchisq(diff$chisq,df=1)
if(pValue<0.001){
              pValue=paste0("=",round(pValue,4))
          }else{
              pValue=paste0("=",round(pValue,3))
          }
surPlot=ggsurvplot(sfit1, 
                   data=meta,
                   conf.int=TRUE,
                   pval=paste0("p",pValue),
                   pval.size=6,
                   risk.table=T,
                   legend.labs=c("high","low"),
                   legend.title=paste0(g," level"),
                   xlab="Time(years)",
                   break.time.by = 1,
                   risk.table.title="",
                   palette=c("red", "blue"),
                   risk.table.height=.25) 
print(surPlot)

stringr::str_replace_all()
str_detect(colnames(exp),"TCGA-W5-AA2R")