TCGA数据挖掘十：风险因子联动图
2019-07-15 本文已影响129人
mayoneday
一.加载表达数据,并分组
这里没有按照全部处理，直接加载了处理好的数据，前面处理数据的代码详解看见前面《TCGA数据挖掘一》
# https://mp.weixin.qq.com/s/J-vaFq1Vv-zR1LC0Wop1zw
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5051943/figure/djw200-F2/
# https://www.ncbi.nlm.nih.gov/pubmed/24893932
# Sixteen of the 111 most significantly altered miRNAs were associated with OS across different clinical subclasses of the TCGA-derived LUAD cohort.
# A linear prognostic model of eight miRNAs (miR-31, miR-196b, miR-766, miR-519a-1, miR-375, miR-187, miR-331 and miR-101-1) was constructed and weighted by the importance scores from the supervised principal component method to divide patients into high- and low-risk groups.
# Patients assigned to the high-risk group exhibited poor OS compared with patients in the low-risk group (hazard ratio [HR]=1.99, P <0.001).
# The eight-miRNA signature is an independent prognostic marker of OS of LUAD patients and demonstrates good performance for predicting 5-year OS (Area Under the respective ROC Curves [AUC] = 0.626, P = 0.003), especially for non-smokers (AUC = 0.686, P = 0.023).
rm(list=ls())
options(stringsAsFactors = F)
Rdata_dir='Rdata/'
Figure_dir='figures/'
library(survival)
library(survminer)
# 这里举例的文章不一样，所以不再使用前面步骤的数据。
# 而是对TCGA-LUAD-miRNA重新处理，正好跟前面的步骤对应学习。
if(F){
library(RTCGA.miRNASeq)
??RTCGA.miRNASeq
s=rownames(LUAD.miRNASeq)[seq(1,nrow(LUAD.miRNASeq),by=3)]
expr <- expressionsTCGA(LUAD.miRNASeq)
dim(expr)
expr[1:40,1:4]
expr=as.data.frame(expr[seq(1,nrow(expr),by=3),3:ncol(expr)])
mi=colnames(expr)
expr=apply(expr,1,as.numeric)
colnames(expr)=s
rownames(expr)=mi
expr[1:4,1:4]
expr=na.omit(expr)
dim(expr)
expr=expr[apply(expr, 1,function(x){sum(x>1)>10}),]
dim(expr)
library(RTCGA.clinical)
??RTCGA.clinical
meta <- LUAD.clinical
tmp=as.data.frame(colnames(meta))
meta[(grepl('patient.bcr_patient_barcode',colnames(meta)))]
meta[(grepl('patient.days_to_last_followup',colnames(meta)))]
meta[(grepl('patient.days_to_death',colnames(meta)))]
meta[(grepl('patient.vital_status',colnames(meta)))]
## patient.race # patient.age_at_initial_pathologic_diagnosis # patient.gender
# patient.stage_event.clinical_stage
meta=as.data.frame(meta[c('patient.bcr_patient_barcode','patient.vital_status',
'patient.days_to_death','patient.days_to_last_followup',
'patient.race',
'patient.age_at_initial_pathologic_diagnosis',
'patient.gender' ,
'patient.stage_event.pathologic_stage')])
#meta[(grepl('patient.stage_event.pathologic_stage',colnames(meta)))]
save(expr,meta,file = file.path(Rdata_dir,'TCGA-LUAD-miRNA-example.Rdata'))
}
真正加载表达数据的代码：
load(file = file.path(Rdata_dir,'TCGA-LUAD-miRNA-example.Rdata'))#直接加载处理好的数据
group_list=ifelse(as.numeric(substr(colnames(expr),14,15)) < 10,'tumor','normal')#分组
table(group_list)
二.加载生存数据：
这里没有按照全部处理，直接加载了处理好的数据，前面处理数据的代码详解看见前面《TCGA数据挖掘一》
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")
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,na.rm = T),'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-LUAD-survival_input.Rdata'))
}
真正加载生存数据的代码：
load(file=file.path(Rdata_dir,'TCGA-LUAD-survival_input.Rdata'))
head(phe)
exprSet[1:4,1:4]
dim(exprSet)
# 这个时候是515个病人的673个miRNA表达矩阵。
三.挑选感兴趣的基因构建coxph模型 ,顺便画了个森林图
## 挑选感兴趣的基因构建coxph模型
#由前面那些算法所得
# miR-31, miR-196b, miR-766, miR-519a-1, miR-375, miR-187, miR-331 and miR-101-1
# hsa-mir-31,hsa-mir-196b,hsa-mir-766,hsa-mir-519a-1,hsa-mir-375,hsa-mir-187,hsa-mir-331,hsa-mir-101-1
e=t(exprSet[c('hsa-mir-31','hsa-mir-196b','hsa-mir-766','hsa-mir-519a-1','hsa-mir-375','hsa-mir-187','hsa-mir-331','hsa-mir-101-1'),])#挑选基因并进行行列转换
e=log2(e+1)#log基因数值
colnames(e)=c('miR31','miR196b','miR766','miR519a1','miR375','miR187','miR331','miR101')#加上列名为基因名
挑选且处理过的表达矩阵
dat=cbind(phe,e)#按照列合并了phe和e
得到的dat包括临床信息及表达数据列
dat$gender=factor(dat$gender)#把gender变成因子，才能进行下面的分析
dat$stage=factor(dat$stage)#把stage变成因子，才能进行下面的分析
colnames(dat) #查看dat的行
s=Surv(time, event) ~ miR31+miR196b+miR766+miR519a1+miR375+miR187+miR331+miR101#time是生存时间，event是生存状态，~后接生存模型需要分析哪些东西
model <- coxph(s, data = dat )#进行cox回归
summary(model,data=dat)
options(scipen=1)
ggforest(model, data =dat,
main = "Hazard ratio",
cpositions = c(0.10, 0.22, 0.4),
fontsize = 1.0,
refLabel = "1", noDigits = 4)#画森林图
dev.new()
image.png
四.用C-index评价模型预测情况
# 对于生存数据预后模型评价很多采用C-index ，但c-index展示没有roc曲线图来的直观
new_dat=dat
fp <- predict(model,new_dat,type="risk");boxplot(fp)
#predict,放形函数，放进去什么就预测什么，此处的model是coxh
fp <- predict(model,new_dat,type="expected") ;boxplot(fp)#expected预测时间相关
plot(fp,phe$days)
fp <- predict(model,new_dat) ;boxplot(fp)
#查看predict.coxph说明书可以知道：预测值的类型（意思就是预测的什么东西）：选择包括线性预测因子(“lp”)、风险评分exp(“risk”)、给定协变量和随访时间的预期事件数(“expected”)和线性预测因子的术语(“terms”)。一个对象的生存概率等于exp(-expected)。
basehaz(model)
library(Hmisc)
options(scipen=200)
with(new_dat,rcorr.cens(fp,Surv(time, event) ))#我理解是比较原始数据和预测数据判断准确性得到C-index
# http://dni-institute.in/blogs/cox-regression-interpret-result-and-predict/
生存数据预后模型评价C-index值
五.画风险因子联动图
library(cowplot)
library(pheatmap)
# https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html
fp_dat=data.frame(s=1:length(fp),v=as.numeric(sort(fp )))#把fp根据大小排序，S是它的顺序，这个顺序很重要，因为三个图都要根据它来画
fp_dat
sur_dat=data.frame(s=1:length(fp),
t=phe[names(sort(fp )),'time'] ,
e=phe[names(sort(fp )),'event'] )
#sort(fp )吧fp排序，fp从小到大排列好了，原来的排列中每一个fp对应的行名名字的那一行的的event被取出来，这个取出来的方式的有先后的，就是fp的顺序，所有最后状态就按照了fp的对应fp的大小的顺序排序了
#简单说就是把生存时间和生存状态按照样本fp大小顺序排序，并生成数据框，这就是把一列按照另一列排序的方法！！
QQ截图20190715213037.jpg
sur_dat$e=ifelse(sur_dat$e==0,'alive','death')#把生死数据进行处理
QQ截图20190715213212.jpg
exp_dat=new_dat[names(sort(fp )),10:17]#挑出表达矩阵
new_dat
可以看到new_dat是这样的，所有取表达矩阵要取10:17列才是表达数据，且names(sort(fp ))就是刚介绍的按照fp进行排序了的，所以最后得到的exp_dat就是按照fp排序了的表达矩阵
exp_dat
plot.point=ggplot(fp_dat,aes(x=s,y=v))+geom_point();print(plot.point)
第一个风险因子图
plot.sur=ggplot(sur_dat,aes(x=s,y=t))+geom_point(aes(col=e));print(plot.sur)
生存时间图.png
mycolors <- colorRampPalette(c("black", "green", "red"), bias = 1.2)(100)
tmp=t(scale(exp_dat))
tmp[tmp > 1] = 1
tmp[tmp < -1] = -1
plot.h=pheatmap(tmp,col= mycolors,show_colnames = F,cluster_cols = T)
热图一列聚类了
plot.h=pheatmap(tmp,col= mycolors,show_colnames = F,cluster_cols = F)
热图一列没有聚类
plot_grid(plot.point, plot.sur, plot.h$gtable,
labels = c("A", "B","C"),
align = 'v',ncol = 1)
风险因子联动图就结合起来了
