TCGA數(shù)據(jù)挖掘三：針對不同因素和不同方法的生存分析

加載并處理數(shù)據(jù)

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')
load("D:/R/R TCGA/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=t(exprSet)
 rownames(exprSet)<-group_list
 exprSet=t(exprSet)
 exprSet= expr[,colnames(exprSet)=="tumor"]#選出腫瘤樣本,生存分析不針對正常人做
 exprSet=na.omit(exprSet)#刪除不要的值
  head(meta)
  colnames(meta)
  meta[,3][is.na(meta[,3])]=0#把第3列NA變?yōu)镺
  meta[,4][is.na(meta[,4])]=0#把第4列NA變?yōu)镺
  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")#改變?nèi)〕龅膸捉M的名字
  # R里面實現(xiàn)生存分析非常簡單！
  
  # 用my.surv <- surv(OS_MONTHS,OS_STATUS=='DECEASED')構(gòu)建生存曲線。
  # 用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)#把狀態(tài)改為數(shù)字，死亡為1，生存為0
  meta$age=as.numeric(meta$age)#年齡
  library(stringr) 
  meta$stage=str_split(meta$stage,' ',simplify = T)[,2]
  #對字符串進行處理，把腫瘤分級用空格分開，取后面的部分
  table(meta$stage)
  boxplot(meta$age)
  meta$agegroup=ifelse(meta$age>median(meta$age),'older','younger')#把年齡根據(jù)中位數(shù)分為兩組
  table(meta$agegroup)
  meta$time=meta$days/30#把日變成月
  phe=meta
  meta
  head(phe)
  phe$ID=toupper(phe$ID) #變成大寫，因為前面是大寫
  phe=phe[match(substr(colnames(exprSet),1,12),phe$ID),]
  #substr(colnames(exprSet),1,12)取列名的1到12位，match把臨床數(shù)據(jù)種樣本和表達矩陣樣本匹配，把前面的id找到后面位置排序
  head(phe)
  exprSet[1:4,1:4]
  
  save(exprSet,phe,
       file = 'TCGA-KIRC-miRNA-survival_input.Rdata')
  
}
# 上面被關閉的代碼，就是在整理臨床信息和生存分析的表達矩陣。
# 整理好的數(shù)據(jù)，直接加載即可
load(  file = 'TCGA-KIRC-miRNA-survival_input.Rdata')

針對臨床資料某一因素，如年齡，性別等進行生存分析，并畫圖

head(phe)
exprSet[1:4,1:4]
# 利用ggsurvplot快速繪制漂亮的生存曲線圖
sfit <- survfit(Surv(time, event)~age_group, data=phe)#根據(jù)性別畫圖
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
## more complicate figures.
ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
           risk.table =TRUE,pval =TRUE,
           conf.int =TRUE,xlab ="Time in months", 
           ggtheme =theme_light(), 
           ncensor.plot = TRUE)
## 多個 ggsurvplots作圖生存曲線代碼合并 
sfit1=survfit(Surv(time, event)~gender, data=phe)
sfit2=survfit(Surv(time, event)~age_group, data=phe)
splots <- list()
splots[[1]] <- ggsurvplot(sfit1,pval =TRUE, data = phe, risk.table = TRUE)
splots[[2]] <- ggsurvplot(sfit2,pval =TRUE, data = phe, risk.table = TRUE)
# Arrange multiple ggsurvplots and print the output
arrange_ggsurvplots(splots, print = TRUE,  ncol = 2, nrow = 1, risk.table.height = 0.4)
dev.off()
# 可以很明顯看到，腫瘤病人的生存受著診斷癌癥的年齡的影響，卻與性別無關。
# 在相對年長的時候診斷的癌癥患者通常會死的快一點。
Rplot.jpeg

12.png

針對基因的生存分析：方法一：挑選感興趣的基因做生存分析

# 來自于文章：2015-TCGA-ccRCC-5-miRNAs-signatures
# Integrated genomic analysis identifies subclasses and prognosis signatures of kidney cancer
# miR-21,miR-143,miR-10b,miR-192,miR-183
tmp=as.data.frame(rownames(exprSet))
g1='hsa-mir-21' # p value = 0.0059
g2='hsa-mir-143' # p value = 0.0093
g3='hsa-mir-192' # p value = 0.00073
g4='hsa-mir-183' # p value = 0.00092
g5='hsa-mir-10b' # p value < 0.0001
gs=c('hsa-mir-21','hsa-mir-143','hsa-mir-192',
     'hsa-mir-183','hsa-mir-10b') 
splots <- lapply(gs, function(g){
  phe$gene=ifelse(exprSet[g1,]>median(exprSet[g1,]),'high','low')#用基因的中位數(shù)分組
  table(phe$gene)
  sfit1=survfit(Surv(time, event)~gene, data=phe)
  ggsurvplot(sfit1,pval =TRUE, data = phe, risk.table = TRUE)
}) 
arrange_ggsurvplots(splots, print = TRUE,  
                    ncol = 2, nrow = 3, risk.table.height = 0.4)
dev.off()

針對基因的生存分析：方法二：批量生存分析 使用 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)
})#得出每一個基因生存分析的P值
log_rank_p=sort(log_rank_p)#取出每一個基因生存分析的P值，形成表
head(log_rank_p)
boxplot(log_rank_p)  
table(log_rank_p<0.01)#哪些是P小于0，001的
log_rank_p[log_rank_p<0.01]#選列出那些P<0.001的基因

# 可以看到，文章里面挑選出來的生存分析相關的miRNA基因，在我們的分析里面都是顯著的。

c('hsa-mir-21','hsa-mir-143','hsa-mir-192',
  'hsa-mir-183','hsa-mir-10b')  %in% names(log_rank_p[log_rank_p<0.01])

把分析出來的生存結(jié)果可視化：利用選出來的生存差異基因做圖


library(pheatmap)
choose_gene=names(log_rank_p[log_rank_p<0.01])
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 = 'logRank_genes.heatmap.png')

library(ggfortify)
df=as.data.frame(t(choose_matrix))
df$group=group_list
png('logRank_genes.pca.png',res=120)
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')+theme_bw()
dev.off()
library("FactoMineR")
library("factoextra")  
## 這里的PCA分析，被該R包包裝成一個簡單的函數(shù)，復雜的原理后面講解。
dat.pca <- PCA(t(choose_matrix), graph = FALSE) #'-'表示“非”
fviz_pca_ind(dat.pca,repel =T,
             geom.ind = "point", # show points only (nbut not "text")只顯示點不顯示文本
             col.ind =  group_list, # color by groups 顏色組
             # palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses 集中成橢圓
             legend.title = "Groups"
)
image.png

針對基因的生存分析：方法二：批量生存分析 使用 coxh

把其他因素對于生存的影響也考慮進去了

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

#Rdata_dir='../Rdata/'
#Figure_dir='../figures/'
# 加載上一步從RTCGA.miRNASeq包里面提取miRNA表達矩陣和對應的樣本臨床信息。
load("D:/R/R TCGA/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)
load("D:/R/R TCGA/survival_input.Rdata")
library(survival)
library(survminer)

## 批量生存分析 使用 coxph 回歸方法
# http://www.sthda.com/english/wiki/cox-proportional-hazards-model
colnames(phe)
mySurv=with(phe,Surv(time, event))#組合生存狀態(tài)和時間
# 對五百多個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)#構(gòu)建一個分組和多個因素的生存分析表
  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',])#返回最后一行，也就是關于基因的生存分析結(jié)果
  
})#循環(huán)這個函數(shù)，每個基因都進行一次運算，最后輸出所有基因的運算結(jié)果
cox_results=t(cox_results)
table(cox_results[,4]<0.05)
cox_results[cox_results[,4]<0.05,]#選出P<0.05的基因

根據(jù)調(diào)出來的基因畫圖

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()

## 也可以嘗試其它主成分分析的R包，視頻就不繼續(xù)沒完沒了的講解了。


library("FactoMineR")
library("factoextra")  
## 這里的PCA分析，被該R包包裝成一個簡單的函數(shù)，復雜的原理后面講解。
dat.pca <- PCA(t(choose_matrix), graph = FALSE) #'-'表示“非”
fviz_pca_ind(dat.pca,repel =T,
             geom.ind = "point", # show points only (nbut not "text")只顯示點不顯示文本
             col.ind =  group_list, # color by groups 顏色組
             # palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses 集中成橢圓
             legend.title = "Groups"
)

最后

感謝jimmy的生信技能樹團隊！

感謝導師岑洪老師！

感謝健明、孫小潔，慧美等生信技能樹團隊的老師一路以來的指導和鼓勵！

文中代碼來自生信技能樹jimmy老師!

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內(nèi)容

  • 隨著公共數(shù)據(jù)庫的建立和開放，越來越多的研究者可以接觸到測序數(shù)據(jù)，非常適合想我們這種“三無”研究者（無課題，無經(jīng)費，...
    墨墨如閱讀 34,178評論 5贊 49
  • 書名貪婪的大腦：為何人類會無止境地尋求意義作者（英）丹尼爾·博爾（Daniel Bor）譯者林旭文豆瓣http:/...
    xuwensheng閱讀 15,912評論 9贊 54
  • 天涯詩雨染流年詩詞征集 在這到處彌漫著金錢氣息的時代 520 這個陽光燦爛的日子 太陽 只想曬曬 人世間到底還有多...
    破天船長閱讀 551評論 0贊 0
  • 因一個偶然的寫作靈感，決定以姑媽為中心寫一些家事。第一篇從最幼稚的吃貨和玩精開始寫起，只因喜歡姑媽做的好飯，更喜歡...
    凡客俠行閱讀 178評論 0贊 1
  • 在漢朝的時候，有一個名醫(yī)叫淳于意。他醫(yī)術(shù)高超，每天看病的人很多，他的理想就是治病救人，因此也很少收病人的錢。弄...
    Aaa艾米閱讀 1,826評論 23贊 5

友情鏈接更多精彩內(nèi)容