image.png
image.png
image.png

image.png
event需要是0、1
image.png

生存分析前的数据整理

  1. 表达矩阵只需要tumor数据，不要normal，将其去掉，新表达矩阵数据命名为exprSet；
  2. 临床信息需要进一步整理，成为生存分析需要的格式，新临床信息数据命名为meta。
  3. rm(list=ls())
  4. proj = "TCGA-KIRC"
  5. load(paste0(proj,".Rdata"))
  6. library(stringr)
  7. 1.整理表达矩阵
  8. 1.1去除normal样本
  9. table(Group)
  10. #> Group
  11. #> normal  tumor 
  12. #>     72    535
  13. exprSet = exp[,Group=='tumor']
  14. ncol(exp)
  15. #> [1] 607
  16. ncol(exprSet)
  17. #> [1] 535
  18. 1.2 基因过滤
  19. 再次进行基因过滤.
  21. (1)标准1：至少要在50%的样本里表达量大于0（最低标准）。
  23. exp(600样本)满足“至少在300个样本里表达量>0”,不能等同于 exprSet(500样本)满足“至少在250个样本里表达量>0”
  25. k = apply(exprSet,1, function(x){sum(x>0)>0.5*ncol(exprSet)});table(k)
  26. #> k
  27. #> FALSE  TRUE 
  28. #>   155 31798
  29. exprSet = exprSet[k,]
  30. nrow(exprSet)
  31. #> [1] 31798
  32. (2)标准2：至少在一半以上样本里表达量＞10(其他数字也可，酌情调整)
  34. k = apply(exprSet,1, function(x){sum(x>10)>0.5*ncol(exprSet)});table(k)
  35. #> k
  36. #> FALSE  TRUE 
  37. #> 11718 20080
  38. exprSet = exprSet[k,]
  39. nrow(exprSet)
  40. #> [1] 20080
  41. 1.3 使用logCPM或logTPM数据
  42. exprSet[1:4,1:4]
  43. #>            TCGA-BP-4342-01A TCGA-A3-3387-01A TCGA-BP-4167-01A TCGA-B8-4620-01A
  44. #> WASH7P                   77               69               23               29
  45. #> AL627309.6              282               36               26               88
  46. #> WASH9P                  222              147               74               47
  47. #> AP006222.1               77               23               23               14
  48. exprSet=log2(edgeR::cpm(exprSet)+1)
  49. exprSet[1:4,1:4]
  50. #>            TCGA-BP-4342-01A TCGA-A3-3387-01A TCGA-BP-4167-01A TCGA-B8-4620-01A
  51. #> WASH7P             1.049412        1.0031558        0.5304210        0.5510009
  52. #> AL627309.6         2.297940        0.6078858        0.5871819        1.2698370
  53. #> WASH9P             2.029993        1.6506568        1.2807483        0.8104724
  54. #> AP006222.1         1.049412        0.4166163        0.5304210        0.2922303
  55. 1.4 样本与病人
  56. 有的病人会有两个或两个以上的肿瘤样本，就有重复。两种可行的办法：
  57. （1）以病人为中心，对表达矩阵的列按照病人ID去重复，每个病人只保留一个样本。
  59. exprSet = exprSet[,sort(colnames(exprSet))]
  60. k = !duplicated(str_sub(colnames(exprSet),1,12));table(k)
  61. #> k
  62. #> FALSE  TRUE 
  63. #>     5   530
  64. exprSet = exprSet[,k] 
  65. ncol(exprSet)
  66. #> [1] 530
  67. （2）以样本为中心，如果每个病人有多个样本则全部保留。(删掉上面这一段代码即可)
  69. 2.整理生存信息和临床信息
  70. 2.1数据合并
  71. xena单独整理了生存信息，和临床信息是分开存放的。后续构建模型可能会需要纳入一些临床信息，所以要合并到一起。
  73. library(dplyr)
  74. nrow(surv)
  75. #> [1] 979
  76. nrow(clinical)
  77. #> [1] 985
  78. meta = left_join(surv,clinical,by = c("sample"= "submitter_id.samples"))
  79. nrow(surv)
  80. #> [1] 979
  81. 2.2 样本过滤
  82. 去掉生存信息不全或者生存时间小于30天的样本，样本纳排标准不唯一，且差别很大
  84. k1 = meta$OS.time >= 30;table(k1)
  85. #> k1
  86. #> FALSE  TRUE 
  87. #>    20   959
  88. k2 = !(is.na(meta$OS.time)|is.na(meta$OS));table(k2)
  89. #> k2
  90. #> TRUE 
  91. #>  979
  92. meta = meta[k1&k2,]
  94. 2.3 选列、简化列名
  95. 小技巧：搜索列名
  96. tmp = data.frame(colnames(meta))
  97. #View(tmp)
  98. 选择需要的列，简化列名
  100. meta = meta[,c(
  101.   'sample',
  102.   'OS',
  103.   'OS.time',
  104.   'race.demographic',
  105.   'age_at_initial_pathologic_diagnosis',
  106.   'gender.demographic' ,
  107.   'tumor_stage.diagnoses'
  108. )]
  109. colnames(meta)=c('ID','event','time','race','age','gender','stage')
  110. str(meta)
  111. #> 'data.frame':    959 obs. of  7 variables:
  112. #>  $ ID    : chr  "TCGA-B8-5162-01A" "TCGA-B4-5834-01A" "TCGA-CW-6087-01A" "TCGA-CW-6087-11A" ...
  113. #>  $ event : int  0 0 1 1 1 1 1 1 0 1 ...
  114. #>  $ time  : int  36 38 41 41 42 42 43 43 43 51 ...
  115. #>  $ race  : chr  "white" "white" "white" "white" ...
  116. #>  $ age   : int  62 59 61 61 75 75 65 65 35 72 ...
  117. #>  $ gender: chr  "male" "male" "male" "male" ...
  118. #>  $ stage : chr  "stage ii" "stage i" "stage iv" "stage iv" ...
  119. 2.4 简化、规范内容
  120. 结局事件
  121. 生存分析的输入数据里，要求结局事件必须用0和1表示，1表示阳性结局。
  123. xena的数据是整理好的，其他来源的需要自行检查和整理。
  125. table(meta$event)
  126. #> 
  127. #>   0   1 
  128. #> 622 337
  129. #自行整理的例子：
  130. a = c("Dead", "Dead", "Alive", "Dead", "Alive")
  131. ifelse(a=="Alive",0,1)
  132. #> [1] 1 1 0 1 0
  133. 生存时间
  134. 认清生存时间的单位（通常是月，也可以用天和年）;
  136. range(meta$time)
  137. #> [1]   36 4537
  138. meta$time = meta$time/30
  139. range(meta$time)
  140. #> [1]   1.2000 151.2333
  141. (3)检查各列的信息是否规范，没有冗余信息。缺失的信息用NA代替
  142. table(meta$stage)
  143. #> 
  144. #> not reported      stage i     stage ii    stage iii     stage iv 
  145. #>            4          467          101          226          161
  146. meta$stage = meta$stage %>% 
  147.   str_replace("not reported",as.character(NA)) %>% 
  148.   str_remove("stage ") %>% 
  149.   str_to_upper()
  150. table(meta$stage,useNA = "always")
  151. #> 
  152. #>    I   II  III   IV <NA> 
  153. #>  467  101  226  161    4
  154. table(meta$race)
  155. #> 
  156. #>                     asian black or african american              not reported 
  157. #>                        16                        69                        14 
  158. #>                     white 
  159. #>                       860
  160. meta$race = meta$race %>% 
  161.   str_replace("not reported",as.character(NA))
  162. table(meta$race,useNA = "always")
  163. #> 
  164. #>                     asian black or african american                     white 
  165. #>                        16                        69                       860 
  166. #>                      <NA> 
  167. #>                        14
  168. table(meta$gender,useNA = "always")
  169. #> 
  170. #> female   male   <NA> 
  171. #>    327    632      0
  172. str(meta)
  173. #> 'data.frame':    959 obs. of  7 variables:
  174. #>  $ ID    : chr  "TCGA-B8-5162-01A" "TCGA-B4-5834-01A" "TCGA-CW-6087-01A" "TCGA-CW-6087-11A" ...
  175. #>  $ event : int  0 0 1 1 1 1 1 1 0 1 ...
  176. #>  $ time  : num  1.2 1.27 1.37 1.37 1.4 ...
  177. #>  $ race  : chr  "white" "white" "white" "white" ...
  178. #>  $ age   : int  62 59 61 61 75 75 65 65 35 72 ...
  179. #>  $ gender: chr  "male" "male" "male" "male" ...
  180. #>  $ stage : chr  "II" "I" "IV" "IV" ...
  181. 3.实现表达矩阵与临床信息的匹配
  182. 即：meta的每一行与exprSet每一列一一对应
  184. rownames(meta) <- meta$ID
  185. s = intersect(rownames(meta),colnames(exprSet));length(s)
  186. #> [1] 513
  187. exprSet = exprSet[,s]
  188. meta = meta[s,]
  189. dim(exprSet)
  190. #> [1] 20080   513
  191. dim(meta)
  192. #> [1] 513   7
  193. identical(rownames(meta),colnames(exprSet))
  194. #> [1] TRUE
  195. save(meta,exprSet,proj,file = paste0(proj,"_sur_model.Rdata")) #将_sur_model.Rdata弄好

实操

  1. 1.准备输入数据
  2. rm(list = ls())
  3. proj = "TCGA-KIRC"
  4. load(paste0(proj,"_sur_model.Rdata"))
  5. ls()
  6. #> [1] "exprSet" "meta"    "proj"
  7. exprSet[1:4,1:4]
  8. #>            TCGA-B8-5162-01A TCGA-B4-5834-01A TCGA-CW-6087-01A TCGA-B0-4698-01A
  9. #> WASH7P           0.38330266        1.0823971       0.24141414        0.8200756
  10. #> AL627309.6       0.04121953        0.1330091       0.60612548        1.7747381
  11. #> WASH9P           0.95751073        1.8118992       0.66255436        1.7390920
  12. #> AP006222.1       0.15825052        0.2713707       0.08503973        0.5589702
  13. str(meta)
  14. #> 'data.frame':    513 obs. of  7 variables:
  15. #>  $ ID    : chr  "TCGA-B8-5162-01A" "TCGA-B4-5834-01A" "TCGA-CW-6087-01A" "TCGA-B0-4698-01A" ...
  16. #>  $ event : int  0 0 1 1 1 0 1 0 1 0 ...
  17. #>  $ time  : num  1.2 1.27 1.37 1.4 1.43 ...
  18. #>  $ race  : chr  "white" "white" "white" "white" ...
  19. #>  $ age   : int  62 59 61 75 65 35 72 50 84 76 ...
  20. #>  $ gender: chr  "male" "male" "male" "male" ...
  21. #>  $ stage : chr  "II" "I" "IV" "IV" ...
  23. 2.KM-plot
  24. 简单版本和进阶版本
  25. library(survival)
  26. library(survminer)
  28. sfit <- survfit(Surv(time, event)~gender, data=meta)
  29. ggsurvplot(sfit,pval=TRUE)
  32. ggsurvplot(sfit,
  33.            palette = "jco",
  34.            risk.table =TRUE,
  35.            pval =TRUE,
  36.            conf.int =TRUE)
  39. 连续型信息怎么作KM分析？例如年龄,基因？
  40. 连续型数据的离散化
  42. 年龄
  43. group = ifelse(meta$age>median(meta$age,na.rm = T),"older","younger")
  44. table(group)
  45. #> group
  46. #>   older younger 
  47. #>     251     262
  48. sfit=survfit(Surv(time, event)~group, data=meta)
  49. ggsurvplot(sfit,pval =TRUE, data = meta, risk.table = TRUE)
  52. 基因
  53. g = rownames(exprSet)[1];g
  54. #> [1] "WASH7P"
  55. meta$gene = ifelse(exprSet[g,]> median(exprSet[g,]),'high','low')
  56. sfit=survfit(Surv(time, event)~gene, data=meta)
  57. ggsurvplot(sfit,pval =TRUE, data = meta, risk.table = TRUE)

image.png

image.png
image.png

  1. 2.log-rank test
  2. KM的p值是log-rank test得出的，可批量操作
  4. source("KM_cox_function.R")
  5. logrankfile = paste0(proj,"_log_rank_p.Rdata")
  6. if(!file.exists(logrankfile)){
  7.   log_rank_p <- apply(exprSet , 1 , geneKM)
  8.   log_rank_p=sort(log_rank_p)
  9.   save(log_rank_p,file = logrankfile)
  10. }
  11. load(logrankfile)
  12. table(log_rank_p<0.01) 
  13. #> 
  14. #> FALSE  TRUE 
  15. #> 13023  7057
  16. table(log_rank_p<0.05) 
  17. #> 
  18. #> FALSE  TRUE 
  19. #> 10232  9848
  22. 3.批量单因素cox
  23. coxfile = paste0(proj,"_cox.Rdata")
  24. if(!file.exists(coxfile)){
  25.   cox_results <-apply(exprSet , 1 , genecox)
  26.   cox_results=as.data.frame(t(cox_results))
  27.   save(cox_results,file = coxfile)
  28. }
  29. load(coxfile)
  30. table(cox_results$p<0.01)
  31. #> 
  32. #> FALSE  TRUE 
  33. #> 11911  8169
  34. table(cox_results$p<0.05)
  35. #> 
  36. #> FALSE  TRUE 
  37. #>  9373 10707
  39. lr = names(log_rank_p)[log_rank_p<0.01];length(lr)
  40. #> [1] 7057
  41. cox = rownames(cox_results)[cox_results$p<0.01];length(cox)
  42. #> [1] 8169
  43. length(intersect(lr,cox))
  44. #> [1] 6165
  45. save(lr,cox,file = paste0(proj,"_logrank_cox_gene.Rdata"))

KM-cox

  1. geneKM = function(gene){
  2.   meta$group=ifelse(gene>median(gene),'high','low')  
  3.   data.survdiff=survdiff(Surv(time, event)~group,data=meta)
  4.   p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  5.   return(p.val)
  6. }
  7. genecox = function(gene){
  8.   meta$gene = gene
  9.   #可直接使用连续型变量
  10.   m = coxph(Surv(time, event) ~ gene, data =  meta)
  11.   #也可使用二分类变量
  12.   #meta$group=ifelse(gene>median(gene),'high','low') 
  13.   #meta$group = factor(meta$group,levels = c("low","high"))
  14.   #m=coxph(Surv(time, event) ~ group, data =  meta)
  16.   beta <- coef(m)
  17.   se <- sqrt(diag(vcov(m)))
  18.   HR <- exp(beta)
  19.   HRse <- HR * se
  21.   #summary(m)
  22.   tmp <- round(cbind(coef = beta, 
  23.                      se = se, z = beta/se, 
  24.                      p = 1 - pchisq((beta/se)^2, 1),
  25.                      HR = HR, HRse = HRse,
  26.                      HRz = (HR - 1) / HRse, 
  27.                      HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1),
  28.                      HRCILL = exp(beta - qnorm(.975, 0, 1) * se),
  29.                      HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)
  31.   return(tmp['gene',]) 
  32.   #return(tmp['grouphigh',])#二分类变量
  33. }

生存模型

image.png

Lasso回归

作用：①将重要的基因筛选出来，
②可得出公式，得出参与基因的系数，代入后可得出风险值
image.png

训练集和测试集不同

  1. 1.准备输入数据
  2. rm(list = ls())
  3. proj = "TCGA-KIRC"
  4. load(paste0(proj,"_sur_model.Rdata"))
  5. ls()
  6. #> [1] "exprSet" "meta"    "proj"
  7. exprSet[1:4,1:4]
  8. #>            TCGA-B8-5162-01A TCGA-B4-5834-01A TCGA-CW-6087-01A TCGA-B0-4698-01A
  9. #> WASH7P           0.38330266        1.0823971       0.24141414        0.8200756
  10. #> AL627309.6       0.04121953        0.1330091       0.60612548        1.7747381
  11. #> WASH9P           0.95751073        1.8118992       0.66255436        1.7390920
  12. #> AP006222.1       0.15825052        0.2713707       0.08503973        0.5589702
  13. str(meta)
  14. #> 'data.frame':    513 obs. of  7 variables:
  15. #>  $ ID    : chr  "TCGA-B8-5162-01A" "TCGA-B4-5834-01A" "TCGA-CW-6087-01A" "TCGA-B0-4698-01A" ...
  16. #>  $ event : int  0 0 1 1 1 0 1 0 1 0 ...
  17. #>  $ time  : num  1.2 1.27 1.37 1.4 1.43 ...
  18. #>  $ race  : chr  "white" "white" "white" "white" ...
  19. #>  $ age   : int  62 59 61 75 65 35 72 50 84 76 ...
  20. #>  $ gender: chr  "male" "male" "male" "male" ...
  21. #>  $ stage : chr  "II" "I" "IV" "IV" ...
  22. load(paste0(proj,"_logrank_cox_gene.Rdata"))
  23. g = cox[1:1000] #这里是示例，取了1000个单因素cox筛选的基因，活学活用
  24. exprSet = exprSet[g,]
  27. 2.构建lasso回归模型
  28. 输入数据是表达矩阵(仅含tumor样本)和每个病人对应的生死（顺序必须一致）。
  30. x=t(exprSet)
  31. y=meta$event
  32. library(glmnet)
  33. #> Warning: package 'glmnet' was built under R version 4.1.3
  34. #> Warning: package 'Matrix' was built under R version 4.1.3
  36. 2.1挑选合适的λ值
  37. Lambda 是构建模型的重要参数。他的大小关系着模型选择的基因个数
  38. #调优参数
  39. set.seed(1006)
  40. cv_fit <- cv.glmnet(x=x, y=y)
  41. plot(cv_fit)
  43. #系数图
  44. fit <- glmnet(x=x, y=y)
  45. plot(fit,xvar = "lambda")
  48. 两条虚线分别指示了两个特殊的λ值,一个是lambda.min,一个是lambda.1se,这两个值之间的lambda都认为是合适的。
  49. lambda.1se构建的模型最简单，即使用的基因数量少，而lambda.min则准确率更高一点，使用的基因数量更多一点。
  51. 2.2 用这两个λ值重新建模
  52. model_lasso_min <- glmnet(x=x, y=y,lambda=cv_fit$lambda.min)
  53. model_lasso_1se <- glmnet(x=x, y=y,lambda=cv_fit$lambda.1se)
  54. 选中的基因与系数存放于模型的子集beta中，用到的基因有一个s0值，没用的基因只记录了“.”，所以可以用下面代码挑出用到的基因。
  56. head(model_lasso_min$beta,20)
  57. choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]
  58. choose_gene_1se=rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta)!=0]
  59. length(choose_gene_min)
  60. #> [1] 41
  61. length(choose_gene_1se)
  62. #> [1] 33
  63. save(choose_gene_min,file = paste0(proj,"_lasso_choose_gene_min.Rdata"))
  64. save(choose_gene_1se,file = paste0(proj,"_lasso_choose_gene_1se.Rdata"))
  66. 3.模型预测和评估
  67. newx参数是预测对象。输出结果lasso.prob是一个矩阵，第一列是min的预测结果，第二列是1se的预测结果，预测结果是概率，或者说百分比，不是绝对的0和1。
  69. 将每个样本的生死和预测结果放在一起，直接cbind即可。
  71. lasso.prob <- predict(cv_fit, newx=x , s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
  72. head(lasso.prob)
  73. #>                           s1          s2
  74. #> TCGA-B8-5162-01A  0.35441206  0.37673630
  75. #> TCGA-B4-5834-01A  0.06975962  0.07579627
  76. #> TCGA-CW-6087-01A  0.72387660  0.64176153
  77. #> TCGA-B0-4698-01A  0.63402199  0.54280316
  78. #> TCGA-B0-4690-01A  0.83592571  0.75407756
  79. #> TCGA-AS-3778-01A -0.13649093 -0.02911398
  81. ROC曲线
  82. library(pROC)
  83. library(ggplot2)
  84. m <- roc(meta$event,lasso.prob[,1])
  85. g <- ggroc(m,legacy.axes = T,size = 1,color = "#2fa1dd")
  86. auc(m)
  87. #> Area under the curve: 0.8687
  89. g + theme_bw() +
  90.   geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), 
  91.                colour = "grey", linetype = "dashed")+
  92.   annotate("text",x = .75, y = .25,
  93.            label = paste("AUC of min = ",format(round(as.numeric(auc(m)),2),nsmall = 2)),color = "#2fa1dd")
  96. 计算AUC取值范围在0.5-1之间，越接近于1越好。可以根据预测结果绘制ROC曲线。
  98. 两个模型的曲线画在一起
  100. m2 <- roc(meta$event, lasso.prob[,2])
  101. auc(m2)
  102. #> Area under the curve: 0.8507
  103. g <- ggroc(list(min = m,se = m2),legacy.axes = T,size = 1)
  105. g + theme_bw() +
  106.   scale_color_manual(values = c("#2fa1dd", "#f87669"))+
  107.   geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), 
  108.                colour = "grey", linetype = "dashed")+
  109.   annotate("text",x = .75, y = .25,
  110.            label = paste("AUC of min = ",format(round(as.numeric(auc(m)),2),nsmall = 2)),color = "#2fa1dd")+
  111.   annotate("text",x = .75, y = .15,
  112.            label = paste("AUC of 1se = ",format(round(as.numeric(auc(m2)),2),nsmall = 2)),color = "#f87669")
  115. 5.切割数据构建模型并预测
  116. 5.1 切割数据（数据>500）
  117. 用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。
  119. library(caret)
  120. #> Warning: package 'caret' was built under R version 4.1.3
  121. set.seed(12345679)
  122. sam<- createDataPartition(meta$event, p = .5,list = FALSE)
  123. head(sam)
  124. #>      Resample1
  125. #> [1,]         5
  126. #> [2,]         9
  127. #> [3,]        13
  128. #> [4,]        17
  129. #> [5,]        19
  130. #> [6,]        21
  131. 可查看两组一些临床参数切割比例
  133. train <- exprSet[,sam]
  134. test <- exprSet[,-sam]
  135. train_meta <- meta[sam,]
  136. test_meta <- meta[-sam,]
  138. prop.table(table(train_meta$stage))
  139. #> 
  140. #>         I        II       III        IV 
  141. #> 0.5058824 0.1019608 0.2117647 0.1803922
  142. prop.table(table(test_meta$stage)) 
  143. #> 
  144. #>         I        II       III        IV 
  145. #> 0.4941176 0.1176471 0.2470588 0.1411765
  146. prop.table(table(test_meta$race)) 
  147. #> 
  148. #>                     asian black or african american                     white 
  149. #>                0.01587302                0.09920635                0.88492063
  150. prop.table(table(train_meta$race)) 
  151. #> 
  152. #>                     asian black or african american                     white 
  153. #>                0.01574803                0.10236220                0.88188976
  155. 5.2 切割后的train数据集建模
  156. 和上面的建模方法一样。
  158. #计算lambda
  159. x = t(train)
  160. y = train_meta$event
  161. cv_fit <- cv.glmnet(x=x, y=y)
  162. plot(cv_fit)
  165. #构建模型
  166. model_lasso_min <- glmnet(x=x, y=y,lambda=cv_fit$lambda.min)
  167. model_lasso_1se <- glmnet(x=x, y=y,lambda=cv_fit$lambda.1se)
  168. #挑出基因
  169. head(model_lasso_min$beta)
  170. #> 6 x 1 sparse Matrix of class "dgCMatrix"
  171. #>            s0
  172. #> WASH7P      .
  173. #> AL627309.6  .
  174. #> WASH9P      .
  175. #> MTATP6P1    .
  176. #> LINC01409   .
  177. #> LINC00115   .
  178. choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]
  179. choose_gene_1se=rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta)!=0]
  180. length(choose_gene_min)
  181. #> [1] 31
  182. length(choose_gene_1se)
  183. #> [1] 21
  185. 4.模型预测
  186. 用训练集构建模型，预测测试集的生死，注意newx参数变了。
  188. lasso.prob <- predict(cv_fit, newx=t(test), s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
  189. head(lasso.prob)
  190. #>                           s1         s2
  191. #> TCGA-B8-5162-01A  0.39653904 0.43398261
  192. #> TCGA-B4-5834-01A  0.10082314 0.16288716
  193. #> TCGA-CW-6087-01A  0.69858720 0.63793147
  194. #> TCGA-B0-4698-01A  0.52026571 0.54582571
  195. #> TCGA-AS-3778-01A -0.02908933 0.07121892
  196. #> TCGA-B2-4098-01A  0.38169229 0.35982935
  198. 再画ROC曲线
  199. library(pROC)
  200. library(ggplot2)
  201. m <- roc(test_meta$event, lasso.prob[,1])
  202. g <- ggroc(m,legacy.axes = T,size = 1,color = "#2fa1dd")
  203. auc(m)
  204. #> Area under the curve: 0.7696
  206. g + theme_bw() +
  207.   geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), 
  208.                colour = "grey", linetype = "dashed")+
  209.   annotate("text",x = .75, y = .25,
  210.            label = paste("AUC of min = ",format(round(as.numeric(auc(m)),2),nsmall = 2)),color = "#2fa1dd")
  213. 计算AUC取值范围在0.5-1之间，越接近于1越好。可以根据预测结果绘制ROC曲线。
  215. 两个模型的曲线画在一起
  217. m2 <- roc(test_meta$event, lasso.prob[,2])
  218. auc(m2)
  219. #> Area under the curve: 0.7569
  220. g <- ggroc(list(min = m,se = m2),legacy.axes = T,size = 1)
  222. g + theme_bw() +
  223.   scale_color_manual(values = c("#2fa1dd", "#f87669"))+
  224.   geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), 
  225.                colour = "grey", linetype = "dashed")+
  226.   annotate("text",x = .75, y = .25,
  227.            label = paste("AUC of min = ",format(round(as.numeric(auc(m)),2),nsmall = 2)),color = "#2fa1dd")+
  228.   annotate("text",x = .75, y = .15,
  229.            label = paste("AUC of 1se = ",format(round(as.numeric(auc(m2)),2),nsmall = 2)),color = "#f87669")

用自己的数据进行逻辑回归- lasso分析
导入数据格式
image.png

  1. rm(list = ls()) 
  2. library(glmnet)
  3. library(readxl)
  4. library(plyr)
  5. library(caret)
  6. library(corrplot)
  7. library(ggplot2)
  8. library(Hmisc)
  9. library(openxlsx)
  11. data <- read.xlsx("/Users/mac/Desktop/r/data1.xlsx")
  12. x <- as.matrix(data[,-c(1:2)])
  13. y <- as.double(data$Group)
  14. fit <- glmnet(x,y,family = "binomial",nlambda = 1000,alpha = 1)
  15. print(fit)
  16. plot(fit,xvar = "lambda")
  17. #筛选最合适的变量
  18. lasso_fit <- cv.glmnet(x,y,family="binomial",alpha = 1,type.measure = "auc",nlambda = 1000)
  19. plot(lasso_fit)
  20. print(lasso_fit)
  22. #使用lambda.lse 获取的模型较好
  23. lasso_best <- glmnet(x=x,y=y,alpha = 1,lambda = lasso_fit$lambda.1se)
  24. coef(lasso_best)
  26. coefficient <- coef(lasso_best,s=lasso_best$lambda.lse)
  27. coe <- coefficient@x
  28. coe <- as.data.frame(coe)
  29. Active_Index <- which(as.numeric(coefficient)!=0)
  30. active_coefficients <- as.numeric(coefficient)[Active_Index]
  31. variable <- rownames(coefficient)[Active_Index]
  32. variable <- as.data.frame(variable)
  33. variable <-cbind(variable,coe)
  35. #使用lambda.min 适用于样本质量较差，纳入模型中的变量较多
  36. lasso_best <- glmnet(x=x,y=y,alpha = 1,lambda = lasso_fit$lambda.min)
  37. coef(lasso_best)
  39. coefficient <- coef(lasso_best,s=lasso_best$lambda.min)
  40. coe <- coefficient@x
  41. coe <- as.data.frame(coe)
  42. Active_Index <- which(as.numeric(coefficient)!=0)
  43. active_coefficients <- as.numeric(coefficient)[Active_Index]
  44. variable <- rownames(coefficient)[Active_Index]
  45. variable <- as.data.frame(variable)
  46. variable <-cbind(variable,coe)

timeROC

image.png
image.png

cox-forest

  1. 1.准备输入数据
  2. rm(list = ls())
  3. proj = "TCGA-KIRC"
  4. if(!require(My.stepwise))install.packages("My.stepwise")
  5. #> Warning: package 'My.stepwise' was built under R version 4.1.3
  6. load(paste0(proj,"_sur_model.Rdata"))
  7. load(paste0(proj,"_lasso_choose_gene_1se.Rdata"))
  8. g = choose_gene_1se
  10. 2.构建coxph模型
  11. 将用于建模的基因（例如lasso回归选中的基因）从表达矩阵中取出来，，可作为列添加在meta表噶的后面,组成的数据框赋值给dat。
  13. library(stringr)
  14. e=t(exprSet[g,])
  15. colnames(e)= str_replace_all(colnames(e),"-","_")
  16. dat=cbind(meta,e)
  18. dat$gender=as.numeric(factor(dat$gender))
  19. dat$stage=as.numeric(factor(dat$stage))
  20. colnames(dat)
  21. #>  [1] "ID"         "event"      "time"       "race"       "age"       
  22. #>  [6] "gender"     "stage"      "AL391244.2" "AJAP1"      "AL357140.2"
  23. #> [11] "CPLANE2"    "RCC2"       "AL137127.1" "NIPAL3"     "KDF1"      
  24. #> [16] "LINC01389"  "ANGPTL3"    "GADD45A"    "LINC01725"  "LINC02609" 
  25. #> [21] "CHI3L2"     "LINC01719"  "KCNN3"      "MSTO1"      "PMF1_BGLAP"
  26. #> [26] "CRABP2"     "CD1E"       "FCER1A"     "NUF2"       "DUTP6"     
  27. #> [31] "IGFN1"      "NEK2"       "PROX1"      "DUSP5P1"    "HS1BP3_IT1"
  28. #> [36] "OTOF"       "CENPA"      "SLC5A6"     "PRKCE"      "SPTBN1"
  29. 逐步回归法构建最优模型
  30. 输出结果行数太多，所以我注释掉了
  32. library(survival)
  33. #> Warning: package 'survival' was built under R version 4.1.3
  34. library(survminer)
  35. # 不能允许缺失值
  36. dat2 = na.omit(dat)
  37. library(My.stepwise)
  38. vl <- colnames(dat)[c(5:ncol(dat))]
  39. # My.stepwise.coxph(Time = "time",
  40. #                   Status = "event",
  41. #                   variable.list = vl,
  42. #                   data = dat2)
  43. 使用输出结果里的最后一个模型
  45. model = coxph(formula = Surv(time, event) ~ stage + AL357140.2 + AL391244.2 + 
  46.     age + OTOF + IGFN1 + CRABP2 + PMF1_BGLAP + SLC5A6 + PROX1 + 
  47.     LINC01719 + DUTP6, data = dat2)
  49.  3.模型可视化-森林图
  50. ggforest(model,data = dat2)
  52. 4.模型预测
  53. fp <- predict(model,newdata = dat2)
  54. library(Hmisc)
  55. #> Warning: package 'Hmisc' was built under R version 4.1.3
  56. options(scipen=200)
  57. with(dat2,rcorr.cens(fp,Surv(time, event))) #with函数下可将列名当作变量用
  58. #>         C Index             Dxy            S.D.               n         missing 
  59. #>      0.17850596（实为1- C Index 的值，故须换算）     -0.64298808      0.02857663    503.00000000      0.00000000 
  60. #>      uncensored  Relevant Pairs      Concordant       Uncertain 
  61. #>    166.00000000 105834.00000000  18892.00000000 146654.00000000
  62. C-index用于计算生存分析中的COX模型预测值与真实之间的区分度（discrimination），也称为Harrell’s concordanceindex。C-index在0.5-1之间。0.5为完全不一致,说明该模型没有预测作用,1为完全一致,说明该模型预测结果与实际完全一致。
  65. 5.切割数据构建模型并预测
  66. 5.1 切割数据
  67. 用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。
  69. library(caret)
  70. #> Warning: package 'caret' was built under R version 4.1.3
  71. set.seed(12345679)
  72. sam<- createDataPartition(meta$event, p = .5,list = FALSE)
  73. train <- exprSet[,sam]
  74. test <- exprSet[,-sam]
  75. train_meta <- meta[sam,]
  76. test_meta <- meta[-sam,]
  77. 5.2 切割后的train数据集建模
  78. 和上面的建模方法一样。
  80. e=t(train[g,])
  81. colnames(e)= str_replace_all(colnames(e),"-","_")
  82. dat=cbind(train_meta,e)
  84. dat$gender=as.numeric(factor(dat$gender))
  85. dat$stage=as.numeric(factor(dat$stage))
  86. colnames(dat)
  87. #>  [1] "ID"         "event"      "time"       "race"       "age"       
  88. #>  [6] "gender"     "stage"      "AL391244.2" "AJAP1"      "AL357140.2"
  89. #> [11] "CPLANE2"    "RCC2"       "AL137127.1" "NIPAL3"     "KDF1"      
  90. #> [16] "LINC01389"  "ANGPTL3"    "GADD45A"    "LINC01725"  "LINC02609" 
  91. #> [21] "CHI3L2"     "LINC01719"  "KCNN3"      "MSTO1"      "PMF1_BGLAP"
  92. #> [26] "CRABP2"     "CD1E"       "FCER1A"     "NUF2"       "DUTP6"     
  93. #> [31] "IGFN1"      "NEK2"       "PROX1"      "DUSP5P1"    "HS1BP3_IT1"
  94. #> [36] "OTOF"       "CENPA"      "SLC5A6"     "PRKCE"      "SPTBN1"
  95. library(My.stepwise)
  96. dat2 = na.omit(dat)
  97. vl <- colnames(dat2)[c(5:ncol(dat2))]
  98. # My.stepwise.coxph(Time = "time",
  99. #                   Status = "event",
  100. #                   variable.list = vl,
  101. #                   data = dat2)
  102. model = coxph(formula = Surv(time, event) ~ stage + AL357140.2 + IGFN1 + 
  103.     AJAP1 + ANGPTL3 + SLC5A6 + HS1BP3_IT1 + AL137127.1 + LINC01719 + 
  104.     CRABP2, data = dat2)
  106. 5.3 模型可视化
  107. ggforest(model, data =dat2)
  109. 5.4 用切割后的数据test数据集验证模型
  110. e=t(test[g,])
  111. colnames(e)= str_replace_all(colnames(e),"-","_")
  112. test_dat=cbind(test_meta,e)
  113. test_dat$gender=as.numeric(factor(test_dat$gender))
  114. test_dat$stage=as.numeric(factor(test_dat$stage))
  115. fp <- predict(model,newdata = test_dat)
  116. library(Hmisc)
  117. with(test_dat,rcorr.cens(fp,Surv(time, event)))
  118. #>        C Index            Dxy           S.D.              n        missing 
  119. #>     0.22187105    -0.55625790     0.04736798   255.00000000     1.00000000 
  120. #>     uncensored Relevant Pairs     Concordant      Uncertain 
  121. #>    86.00000000 28476.00000000  6318.00000000 36294.00000000