### Create: Jianming Zeng
### Date: 2019-04-02 21:59:01
### Email: jmzeng1314@163.com

### https://github.com/jmzeng1314/GEO/blob/master/GSE11121/step5-surivival.R

rm(list=ls())
options(stringsAsFactors = F)

Rdata_dir='../Rdata/'
Figure_dir='../figures/'
# 加載上一步從RTCGA.miRNASeq包里面提取miRNA表達矩陣和對應的樣本臨床信息产阱。
load( file = 
        file.path(Rdata_dir,'TCGA-KIRC-miRNA-example.Rdata')
)
dim(expr)
dim(meta)
# 可以看到是 537個病人朋譬，但是有593個樣本，每個樣本有 552個miRNA信息。
# 當然绷雏，這個數(shù)據(jù)集可以下載原始測序數(shù)據(jù)進行重新比對窗轩，可以拿到更多的miRNA信息

# 這里需要解析TCGA數(shù)據(jù)庫的ID規(guī)律桩撮，來判斷樣本歸類問題雷绢。
group_list=ifelse(as.numeric(substr(colnames(expr),14,15)) < 10,'tumor','normal')

table(group_list)
exprSet=na.omit(expr)

library(survival)
library(survminer)

# 這里做生存分析，已經(jīng)不需要正常樣本的表達矩陣了课梳，所以需要過濾距辆。
# 而且臨床信息余佃，有需要進行整理。
### survival analysis only for patients with tumor.
if(F){
  exprSet=na.omit(expr)
  exprSet=exprSet[,group_list=='tumor']
  
  head(meta)
  colnames(meta)
  meta[,3][is.na(meta[,3])]=0
  meta[,4][is.na(meta[,4])]=0
  meta$days=as.numeric(meta[,3])+as.numeric(meta[,4])
  meta=meta[,c(1:2,5:9)]
  colnames(meta)
  colnames(meta)=c('ID','event','race','age','gender','stage',"days")
  # R里面實現(xiàn)生存分析非常簡單跨算！
  
  # 用my.surv <- surv(OS_MONTHS,OS_STATUS=='DECEASED')構建生存曲線咙冗。
  # 用kmfit2 <- survfit(my.surv~TUMOR_STAGE_2009)來做某一個因子的KM生存曲線。
  # 用 survdiff(my.surv~type, data=dat)來看看這個因子的不同水平是否有顯著差異漂彤，其中默認用是的logrank test 方法。
  # 用coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung) 來檢測自己感興趣的因子是否受其它因子(age,gender等等)的影響灾搏。
  
  library(survival)
  library(survminer)
  meta$event=ifelse(meta$event=='alive',0,1)
  meta$age=as.numeric(meta$age)
  library(stringr) 
  meta$stage=str_split(meta$stage,' ',simplify = T)[,2]
  table(  meta$stage)
  boxplot(meta$age)
  meta$age_group=ifelse(meta$age>median(meta$age),'older','younger')
  table(meta$race)
  meta$time=meta$days/30
  phe=meta
  
  head(phe)
  phe$ID=toupper(phe$ID) 
  phe=phe[match(substr(colnames(exprSet),1,12),phe$ID),]
  head(phe)
  exprSet[1:4,1:4]
  
  save(exprSet,phe,
       file = 
         file.path(Rdata_dir,'TCGA-KIRC-miRNA-survival_input.Rdata')
  )
}
# 上面被關閉的代碼挫望，就是在整理臨床信息和生存分析的表達矩陣。
# 整理好的數(shù)據(jù)狂窑，直接加載即可
load(  file = 
         file.path(Rdata_dir,'TCGA-KIRC-miRNA-survival_input.Rdata')
)
head(phe)
exprSet[1:4,1:4]


library(survival)
library(survminer)

批量生存分析 使用 coxph 回歸方法

# http://www.sthda.com/english/wiki/cox-proportional-hazards-model
colnames(phe)
mySurv=with(phe,Surv(time, event))
# 對五百多個miRNA基因進行批量運行cox媳板，需要一點點時間。
# 如果是mRNA-seq的表達矩陣泉哈， 通常耗時更長蛉幸。
# 注意，如果是某些基因表達量為恒定丛晦，比如在所有樣本為0奕纫，這個代碼會爆倉
# 需要去除這樣的基因，沒有分析的必要性烫沙。

cox_results <-apply(exprSet , 1 , function(gene){
  # gene= exprSet[1,]
  group=ifelse(gene>median(gene),'high','low') 
  survival_dat <- data.frame(group=group,stage=phe$stage,age=phe$age,
                             gender=phe$gender,
                             stringsAsFactors = F)
  m=coxph(mySurv ~ gender + age + stage+ group, data =  survival_dat)
  
  beta <- coef(m)
  se <- sqrt(diag(vcov(m)))
  HR <- exp(beta)
  HRse <- HR * se
  
  #summary(m)
  tmp <- round(cbind(coef = beta, se = se, z = beta/se, p = 1 - pchisq((beta/se)^2, 1),
                     HR = HR, HRse = HRse,
                     HRz = (HR - 1) / HRse, HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1),
                     HRCILL = exp(beta - qnorm(.975, 0, 1) * se),
                     HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)
  return(tmp['grouplow',])
  
})
cox_results=t(cox_results)
table(cox_results[,4]<0.05)
cox_results[cox_results[,4]<0.05,]

cox_results

批量生存分析 使用 logrank test 方法

mySurv=with(phe,Surv(time, event))
log_rank_p <- apply(exprSet , 1 , function(gene){
  # gene=exprSet[1,]
  phe$group=ifelse(gene>median(gene),'high','low')  
  data.survdiff=survdiff(mySurv~group,data=phe)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  return(p.val)
})

log_rank_p

require("VennDiagram")
VENN.LIST=list(cox=rownames(cox_results[cox_results[,4]<0.05,]),
             log=names(log_rank_p[log_rank_p<0.05]))
venn.plot <- venn.diagram(VENN.LIST , NULL, 
                          fill=c("darkmagenta", "darkblue"), 
                          alpha=c(0.5,0.5), cex = 2, 
                          cat.fontface=4,  
                          main="overlap of coxph and log-rank test")
# 可以看到兩種生存分析挑出來的基因重合度尚可匹层。
grid.draw(venn.plot) 

save(log_rank_p,cox_results ,
     file = 
       file.path(Rdata_dir,'TCGA-KIRC-miRNA-survival_results.Rdata')
)

venn.plot

library(pheatmap)
choose_gene=rownames(cox_results[cox_results[,4]<0.05,])
choose_matrix=expr[choose_gene,]
choose_matrix[1:4,1:4] 
n=t(scale(t(log2(choose_matrix+1))))  #scale()函數(shù)去中心化和標準化
#對每個探針的表達量進行去中心化和標準化
n[n>2]=2 #矩陣n中歸一化后，大于2的項锌蓄，賦值使之等于2（相當于設置了一個上限）
n[n< -2]= -2 #小于-2的項升筏，賦值使之等于-2（相當于設置了一個下限）
n[1:4,1:4]

## http://www.bio-info-trainee.com/1980.html
annotation_col = data.frame( group_list=group_list  )
rownames(annotation_col)=colnames(expr)

pheatmap(n,show_colnames = F,annotation_col = annotation_col,
         filename = 'cox_genes.heatmap.png')

library(ggfortify)
df=as.data.frame(t(choose_matrix))
df$group=group_list
png('cox_genes.pca.png',res=120)
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')+theme_bw()
dev.off()

cox_genes.heatmap

cox_genes.pca

## 也可以嘗試其它主成分分析的R包
library("FactoMineR")
library("factoextra")  

參考來源：生信技能樹

友情鏈接：

課程分享
生信技能樹全球公益巡講
（https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g）
B站公益74小時生信工程師教學視頻合輯
（https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw）
招學徒：
（https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw）

歡迎關注公眾號：青島生信菜鳥團