批量生存分析 提取生存正相关基因 生存负相关基因 survival_genes_split 生存分析得到的显著性基因 高低表达与生存期的关系 0=高表达生存预期差；1=高表达生存预期好s_正相关和负相关 wgcna 0和1

TCGA学习笔记13-对所有差异表达基因进行生存分析
 
https://uteric.github.io/TCGA/TCGA13Survival/
 
 
 
load("G:/r/duqiang_IPF/surval_analysis_3_independent_dataset_IPF_pval_filtered/surval_genes_filtered_log_rank_p.RData")
 
getwd()
dir.create("G:/r/duqiang_IPF/survival_genes_high_low_association")
setwd("G:/r/duqiang_IPF/survival_genes_high_low_association")
 
load("G:/r/duqiang_IPF/surval_analysis_3_independent_dataset_IPF/combined_data_for_surval.RDdata")
table(complete.cases(rownames(expr.freibrug.IPF)))
 
 
 
## 批量生存分析 使用  logrank test 方法  批量生存分析 并获得正负相关基因 高低表达情况
getwd()
## 批量生存分析 使用  logrank test 方法
phe=phe_final_3
head(phe)[1:3,1:19]
mySurv=with(phe,Surv(time, event))
head(mySurv)
log_rank_p=list()
 
s <- vector() # 定义变量s，放surv，即生存预期
s
 
for (eachgene in colnames(phe)[5:ncol(phe)] ) {
  
  data.survdiff=survdiff(mySurv~phe[,eachgene],data=phe)
 
   p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  log_rank_p[[eachgene]]=p.val
  
  s[eachgene]<-ifelse(data.survdiff$obs[1]/data.survdiff$n[[1]] > 
                        data.survdiff$obs[2]/data.survdiff$n[[2]],
                      0,1) # 0=高表达生存预期差；1=高表达生存预期好，；
  
}
1 >
  3
head(s)
table(s)
head(s[s==0])
str(s)
#table(names(s))
head(names(s))
names(s)
 
length(log_rank_p)
head(names(log_rank_p))
table(log_rank_p<0.01)
head(log_rank_p)
log_rank_p=unlist(log_rank_p)
log_rank_p=sort(log_rank_p)
head(log_rank_p)
log_rank_p[1000]
table(log_rank_p<0.005)
boxplot(log_rank_p) 
phe$H6PD=ifelse(exprSet['H6PD',]>median(exprSet['H6PD',]),'high','low')
table(phe$H6PD)
ggsurvplot(survfit(Surv(time, event)~MRVI1, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~ACER3, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~MMP12, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~SPP1, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~GPNMB, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~CSF1, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~PDGFRA, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~BMP6, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~IBSP, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~S100A14, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~STAB1, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~LOC284751, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~ANKRD22, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~AMICA1, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~FLJ43903, data=phe), conf.int=F, pval=TRUE)
ggsurvplot(survfit(Surv(time, event)~RANBP3L, data=phe), conf.int=F, pval=TRUE)
 
 
 
intersect(names(   s[s==1]   ),
          names(    log_rank_p[log_rank_p<0.01]      ))
 
 
getwd()
#save(s,log_rank_p,file ="G:/r/duqiang_IPF/survival_genes_high_low_association/survival_genes_split.RData" )
 
load("G:/r/duqiang_IPF/survival_genes_high_low_association/survival_genes_split.RData")