【原】生信编程16.对有临床信息的表达矩阵批量进行生存分析

健明 2021-07-14

展开全文

有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享：

200个生信工程师面试考题

生存分析简介

在生物医学研究中，生存分析是非常重要和常见的分析方法。生存分析经常用在癌症等疾病的研究中，例如在对某种抗癌药物做临床实验时，会首先筛选一部分癌症患者随机分为两组，一组服用该实验药物，一组用对照药物，服药之后开始统计每个患者从服药一直到死亡的生存时间，通过考察两组之间的病人在生存时间上是否有统计学差异来判断试验药物是否有效。在这种情况下，死亡是整个实验中重点观测的事件，即event。对于每个病人，需要记录它们发生该事件的具体时间。因此，生存分析可以抽象的概述为，研究在不同情况下，特定时间发生的时间的关系是否存在差异。这些具体事件可以是死亡，也可以是肿瘤转移、复发、病人出院、重新入院等任何可以明确识别的事件，而不同条件即为不同的分组依据，可以是年龄、性别、地狱、某个基因表达量的高低、某个突发的携带是否等等。

两个重要的概念：

生存概率：即Survival probability，指的是研究对象从试验开始知道某个时间点仍然存活的概率，可见它是一个对时间t的函数，我们将其定义为S(t)
风险概率：Hazard probability，指的是研究对象从试验开始到某个特定时间t之间的存活，但在t时间点发生观测事件如死亡的概率，它也是对时间t的函数，定义为H(t)。

我觉得这个关于生存分析的教程写得真心不错，像我这样第一次接触生存分析的人也能理解：http://thisis./blog/index.php/2020/04/06/survival-analysis/

首先使用模拟数据来了解下生存分析和风险分析

生存分析

目的，挑选出表达高低与生存差异存在显著相关的基因

#首先生存模拟数据，50个病人，200个基因
dat = matrix(rnorm(10000,mean = 8,sd=4),nrow = 200)
rownames(dat) = paste0('gene_',1:nrow(dat)) #行为基因
colnames(dat) = paste0('sample_', 1:ncol(dat)) #列
os_years = abs(floor(rnorm(ncol(dat),mean = 50,sd=20))) #年龄
os_status = sample(rep(c(0,1),25)) #生存状态,0代表死亡，1代表存活

library(survival)
#创建一个生存对象，
mysurv = Surv(os_years,os_status )
## [1] 40  11  24+  5  52+ 31+ 40  46  47+ 57  41+ 77+ 48+  5+ 47+ 25   5  29+ 73+

#根据公式（Kaplan-Meier）、之前拟合的Cox模型或者先前拟合的加速失效时间模型创建生存曲线
fit.KM = survfit(mysurv~1)

#根据基因的表达情况进行分组，高表达组和低表达组，然后确定每个基因的高表达组和低表达组的病人之间是否存在生存的显著差异，计算P值
log_rank_p <- apply(dat, 1, function(values1){
  group=ifelse(values1> median(values1), 'high', 'low')
  kmfit2 <- survfit(mysurv~group)
  #plot(kmfit2)
  data.survdiff <- survdiff(mysurv~group)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) -1)
})

#挑选出表达高低能够显著区分生存的基因，这些基因可以用来绘制生存曲线图
names(log_rank_p[log_rank_p<0.05])

挑选差异最显著的两个基因进行生存分析绘图

#将p值设置为0.01的时候，只能筛选出两个基因 ，我们使用这两个基因进行分析
diff_genes<-names(log_rank_p[log_rank_p<0.01]) 
new_dat <- as.data.frame(t(dat[diff_genes,]))
new_dat$DSS.time <- os_years*30
new_dat$DSS <- os_status
new_dat$patient <- rownames(new_dat)
library(survminer)

#确定连续表达量的最优切割点
new.cut <- surv_cutpoint(
  new_dat,
  time = 'DSS.time',
  event = 'DSS',
  variables = c("gene_37","gene_91")
)
summary(new.cut)
plot(new.cut,'gene_37',palette = 'npg')
#以切割位点为阈值，将信息分为high，low
new.cut <- surv_categorize(new.cut)#使用两个基因颚高低表达信息来判断是否与病人的生存有关
library(survival)
fit <- survfit(Surv(DSS.time,DSS) ~ gene_37+gene_91, data = new.cut)

ggsurvplot(fit,pval = TRUE, conf.int = FALSE, title = 'METABRIC', 
           xlim = c(0,10000), break.time.by = 2000, risk.table.y.text.col=T, 
           risk.table.y.text = FALSE)

风险分析

目的：挑选出能够与患病风险有相关性的基因

gender= ifelse(rnorm(ncol(dat))>1, 'male','female')
age = abs(floor(rnorm(ncol(dat), mean = 50, sd=20)))

cox_results <- apply(dat, 1, function(values1){
  group = ifelse(values1 > median(values1), 'high','low')
  survival_dat <- data.frame(group=group, gender=gender, age= age, stringsAsFactors = F)
  #拟合比例风险回归模型，使用Andersen和Gill技术公示将时间因变量，时间因式，每个主题的多个时间以及其他或者的内容都纳入计算
  m = coxph(mysurv~age + gender + group, data = survival_dat)
  
  beta <- coef(m)
  se = sqrt(diag(vcov(m)))
  HR <- exp(beta)
  HRse = HR * se
  
  tmp = round(cbind(coef = beta, se=se, z = beta/se, 
                    p = 1-pchisq((beta/se)^2,1), 
                    HR = HR, HRse= HRse, 
                    HEz = (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[,cox_results[4,]<0.05]

使用CGDS数据库的乳腺癌表达数据进行生存分析

library(cgdsr)
#install.packages("cgdsr")
#DT是DataTables的R接口，使得R的数据对象可以在HTML界面中以表格的形式展示
library(DT)

#cgdsr是一个用来访问MSKCC Cancer Genomics Data Server（CGDS）中的数据的包
#从CGDS的终端URL创建一个连接对象
mycgds = cgdsr::CGDS("http://www./")
#为CGDA函数的调用设置详细的日志级别
setVerbose(mycgds, TRUE)
#从CGDS中获取可用的癌症研究
all_TCGA_studies = getCancerStudies(mycgds)
#从CGDS中查询到的所有可用的数据信息以网页的形式展示
DT::datatable(all_TCGA_studies)
#从某个特定的癌症研究中获得可用的遗传数据
all_gen_pro = getGeneticProfiles(mycgds, 'brca_metabric')
#获得特定癌症研究中可用的病例列表
DT::datatable(all_gen_pro)
all_tables <- getCaseLists(mycgds, 'brca_metabric')

#通过查看all_gen_pro可以看到很多可用的数据，我们挑选mRNA的表达数据
my_dataset <- 'brca_metabric_mrna'
my_table <- 'brca_metabric_mrna'
#检索基因以及遗传特征的基因组数据,第二个参数为目标基因，是必须要求的，这里以
#"BCL11A","MBNL1"两个基因为例
exp<-getProfileData(mycgds,c("BCL11A","MBNL1"), my_dataset,my_table)
exp$patient = rownames(exp)

#获取临床信息
cil_dat = getClinicalData(mycgds, my_table)
head(cil_dat)
DT::datatable(cil_dat,
              extensions = "FixedColumns",
              options = list(
                scrollX=TRUE,
                fixedColumns = TRUE
              ))

#查看临床信息，然后取出重要的两列
colnames(cil_dat)
cil_dat_1 <- cil_dat[,c("OS_MONTHS","VITAL_STATUS")]
cil_dat_1$patient <- rownames(cil_dat)

#将表达数据和临床信息根据病人进行合并
tmp <- merge(exp,cil_dat_1,by='patient')
tmp = na.omit(tmp)
#将OS_MONTHS替换为DSS.time, VITAL_STATUS替换为DSS
colnames(tmp)[4] = 'DSS.time'
colnames(tmp)[5] = 'DSS'
tmp$DSS = ifelse(tmp$DSS=='Died of Disease',1,0)
tmp$DSS.time = tmp$DSS.time*30

library(survminer)
#确定连续表达量的最优切割点
tmp.cut <- surv_cutpoint(
  tmp,
  time = 'DSS.time',
  event = 'DSS',
  variables = c("BCL11A","MBNL1")
)
summary(tmp.cut)
plot(tmp.cut,'BCL11A',palette = 'npg')
#以切割位点为阈值，将信息分为high，low
tmp.cut <- surv_categorize(tmp.cut)#使用两个基因颚高低表达信息来判断是否与病人的生存有关
library(survival)
fit <- survfit(Surv(DSS.time,DSS) ~ BCL11A+MBNL1, data = tmp.cut)

ggsurvplot(fit,pval = TRUE, conf.int = FALSE, title = 'METABRIC', 
           xlim = c(0,10000), break.time.by = 2000, risk.table.y.text.col=T, 
           risk.table.y.text = FALSE)