TCGA数据挖掘三：针对不同因素和不同方法的生存分析
2019-07-04 本文已影响15人
mayoneday
加载并处理数据
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信息。
# 当然，这个数据集可以下载原始测序数据进行重新比对，可以拿到更多的miRNA信息
# 这里需要解析TCGA数据库的ID规律，来判断样本归类问题。
group_list=ifelse(as.numeric(substr(colnames(expr),14,15)) < 10,'tumor','normal')
#得到分组信息
table(group_list)
exprSet=na.omit(expr)#删除不要的值
library(survival)
library(survminer)
# 这里做生存分析，已经不需要正常样本的表达矩阵了，所以需要过滤。
# 而且临床信息，有需要进行整理。
### 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变为O
meta[,4][is.na(meta[,4])]=0#把第4列NA变为O
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里面实现生存分析非常简单！
# 用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)#把状态改为数字，死亡为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')#把年龄根据中位数分为两组
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把临床数据种样本和表达矩阵样本匹配，把前面的id找到后面位置排序
head(phe)
exprSet[1:4,1:4]
save(exprSet,phe,
file = 'TCGA-KIRC-miRNA-survival_input.Rdata')
}
# 上面被关闭的代码，就是在整理临床信息和生存分析的表达矩阵。
# 整理好的数据，直接加载即可
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)#根据性别画图
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')#用基因的中位数分组
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])
把分析出来的生存结果可视化：利用选出来的生存差异基因做图
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()函数去中心化和标准化，热图必备
#对每个探针的表达量进行去中心化和标准化
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包包装成一个简单的函数，复杂的原理后面讲解。
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信息。
# 当然，这个数据集可以下载原始测序数据进行重新比对，可以拿到更多的miRNA信息
# 这里需要解析TCGA数据库的ID规律，来判断样本归类问题。
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))#组合生存状态和时间
# 对五百多个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,]#选出P<0.05的基因
根据调出来的基因画图
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()函数去中心化和标准化
#对每个探针的表达量进行去中心化和标准化
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包，视频就不继续没完没了的讲解了。
library("FactoMineR")
library("factoextra")
## 这里的PCA分析，被该R包包装成一个简单的函数，复杂的原理后面讲解。
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"
)
