学徒数据挖掘之谁说生存分析一定要按照表达量中位值或者平均值分组呢？

rm(list=ls())
options(stringsAsFactors = F)
##读取数据
mRNA_HiSeqV2 = read.table('HiSeqV2.gz',header = T,sep = 't',check.names = F,row.names = 1)
dim(mRNA_HiSeqV2)
mRNA_HiSeqV2[1:4,1:4]

#查看NA的数据
dim(mRNA_HiSeqV2)
exp = mRNA_HiSeqV2
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ]
dim(exp)
exp[1:4,1:4]


##临床信息（在这里没有用到）
mRNA_clinical = read.table('LIHC_clinicalMatrix',header = T,sep = 't',row.names = 1)
dim(mRNA_clinical)
mRNA_clinical[1:4,1:4]


##生存相关信息
mRNA_survival = read.table('LIHC_survival.txt.gz',header = T,sep = 't',row.names = 1)
dim(mRNA_survival)
mRNA_survival <- mRNA_survival[, -1] 
mRNA_survival[1:4,1:4]


save(mRNA_survival,mRNA_clinical,exp,group_list,file='LIHC.Rdata')


#01–09是癌症，10–19是正常，20–29是癌旁
根据样本ID的第14-15位，给样本分组（tumor和normal）
library(stringr)
table(str_sub(colnames(exp),14,15))
group_list = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
group_list = factor(group_list,levels = c("normal","tumor"))
table(group_list)


##xena下载的数据是经过log的count值，所以要还原才可以进行差异分析
exprSet = exp
exprSet <- 2^exprSet-1
dim(exprSet)
exprSet[1:4,1:4]

#并且还要取整数，现在就变成count值了
exprSet <- floor(exprSet)
exprSet[1:4,1:4]


save(mRNA_survival,mRNA_clinical,exp,group_list,exprSet,file='LIHC1.Rdata')



rm(list=ls())
options(stringsAsFactors = F)
load('LIHC1.Rdata')
library("stringr")


先看以count数据的结果 OS
tumor_sample <- colnames(exprSet)[substr( colnames(exprSet),14,15) < 10] 
mRNA_survival = mRNA_survival[,-1]
pheno = rownames(mRNA_survival)[substr(rownames(mRNA_survival),14,15) < 10] 
fin_tumor = intersect(tumor_sample,pheno)

# 提取表达矩阵和临床信息 
exprSet = exprSet[,fin_tumor] 
pheno_fina = mRNA_survival[fin_tumor,] 

### 这里又归一化了？
exprSet = log2(exprSet+1)

# 开始生存分析
library(ggrisk) 
library(survival) 
library(survminer)

# 先做无病生存期，首先把DFI不明（NA）的样本给去掉 
DFI =  pheno_fina[,3:6] 
dfi = DFI[!is.na(DFI$DFI),] 
dfi_count = exprSet[,rownames(dfi)]
dd = scale(as.numeric( dfi_count["LHPP",]))

# 绘图 
mySurv<-Surv(dfi$DFI.time, dfi$DFI) 
### 以LHPP的表达量z-score后，<= -1为界，哦z-score是这样用的啊
dfi_group<-ifelse(dd <= -1 ,"low","high")
fi <-survfit(mySurv~dfi_group,data = dfi) 
survminer::ggsurvplot(fi, data = dfi, risk.table = TRUE, 
                      pval = T,  )

# 接下来绘制os的 
os = pheno_fina[,c(1,2,7,8)] 
os = os[!is.na(os$OS),] 
os_count =  exprSet[,rownames(os)]

dd = scale(as.numeric( os_count["LHPP",]))

mySurv<-Surv(os$OS.time, os$OS) 
dfi_group<-ifelse(dd <= -1 ,"low","high")
fi <-survfit(mySurv~dfi_group,data = dfi) 
survminer::ggsurvplot(fi, data = dfi, risk.table = TRUE, 
                      pval = T,  )