美文网首页TCGA数据挖掘Rcox
TCGA生存模型的构建以及模型预测和评估

TCGA生存模型的构建以及模型预测和评估

作者: Hayley笔记 | 来源:发表于2021-06-20 14:24 被阅读0次

    生存模型的构建方法包括:
    1. Lasso回归;
    2. Cox多因素回归;
    3. 随机森林;
    4. 支持向量机

    可以把log_rank_test或cox筛选出的基因单独做模型构建,也可以取交集之后再做模型构建。

    1. Lasso回归(机器学习算法)

    目的:从若干个基因中挑选真正对生存有影响的基因。
    Lasso回归可以对这些基因进行统计和打分,从而挑出关键基因。

    1.1 准备输入数据(表达矩阵数据和临床信息数据)
    load("TCGA-KIRC_sur_model.Rdata")
    ls()
    exprSet[1:4,1:4]
    meta[1:4,1:4]
    
    表达矩阵数据 临床信息数据
    1.2 构建lasso回归模型

    输入数据是表达矩阵(仅含tumor样本)和每个病人对应的生死(顺序必须一致)。

    x=t(exprSet) #转换成基因在列
    y=meta$event #结局
    library(glmnet)
    
    • 1.2.1 挑选合适的λ值

    Lambda 是构建模型的重要参数。他的大小关系着模型选择的基因个数

    #调优参数
    set.seed(1006) 
    #⚠️不设置的话,每次运行时后面的结果(选出的基因)会不断变化。这也就是说,如果文献中使用了Lasso回归,是没有办法复现结果的。
    cv_fit <- cv.glmnet(x=x, y=y) #cv.glmnet()就是一个调优参数的过程
    plot(cv_fit)
    

    图的横轴是选择的λ值,图最上方的数字代表选择对应的λ值时,模型所使用的基因数。
    两条虚线分别指示了两个特殊的λ值,一个是lambda.min,一个是lambda.1se,这两个值之间的lambda都认为是可用的。lambda.1se构建的模型最简单,即使用的基因数量少,而lambda.min则准确率更高一点,使用的基因数量更多一点。

    #系数图
    fit <- glmnet(x=x, y=y) #glmnet是构建模型的
    plot(fit,xvar = "lambda")
    

    这张图中的每一条线代表一个基因。横坐标和上面那张图一样是log Lambda,纵坐标是系数。从左往右看,系数一直在0附近的,不管λ是几都不会被选到。
    这两张图就是Lasso回归最经典的两张图。

    • 1.2.2 用这两个λ值重新建模

    两个都构建一下,再比较哪个模型更好。实际操作也可直接选择lambda.min到lambda.1se中的任意一个λ值构建模型,λ值的选择不是固定的。

    model_lasso_min <- glmnet(x=x, y=y,lambda=cv_fit$lambda.min) # 得到的模型是一个列表
    model_lasso_1se <- glmnet(x=x, y=y,lambda=cv_fit$lambda.1se) # 得到的模型是一个列表
    View(model_lasso_1se) #从中选一个查看一下格式
    

    选中的基因与系数存放于模型的子集beta中,beta是一个稀疏矩阵,用到的基因有一个s0值(系数值),没用的基因只记录了“.”(系数是0),所以可以用下面代码挑出用到的基因。

    head(model_lasso_min$beta,20)
    choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0] #as.numeric不等于0的基因就是有系数的也就是被选中了的基因。
    choose_gene_1se=rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta)!=0]
    length(choose_gene_min)
    # [1] 35 ##选到了35个基因
    length(choose_gene_1se)
    # [1] 11 ##选到了11个基因
    save(choose_gene_min,file = "lasso_choose_gene_min.Rdata")
    
    • 1.2.3 模型预测
      使用predict函数进行模型预测(提供模型和数据就可以得到预测结果)
      newx参数是预测对象。输出结果lasso.prob是一个矩阵,第一列是min的预测结果,第二列是1se的预测结果,预测结果是概率,或者说百分比,不是绝对的0和1。
      将每个样本的生死和预测结果放在一起,直接cbind即可。
    lasso.prob <- predict(cv_fit, newx=x, s=c(cv_fit$lambda.min,cv_fit$lambda.1se) ) 
    re=cbind(y ,lasso.prob)
    head(re)
    #                              y         1         2
    # TCGA-A3-3307-01A-01T-0860-13 0 0.1151638 0.2300976
    # TCGA-A3-3308-01A-02R-1324-13 0 0.4020126 0.3591050
    # TCGA-A3-3311-01A-02R-1324-13 1 0.2906333 0.2954458
    # TCGA-A3-3313-01A-02R-1324-13 1 0.3831590 0.3456060
    # TCGA-A3-3316-01A-01T-0860-13 0 0.3284941 0.3113726
    # TCGA-A3-3317-01A-01T-0860-13 0 0.3832764 0.3342198
    ## y是每个病人,1和2是给的lamda值是min和1se的时候,模型给每个病人计算出来的预测值评分是多少
    re=as.data.frame(re)
    colnames(re)=c('event','prob_min','prob_1se')
    #将真实结局event和预测值合并在一起
    re$event=as.factor(re$event)
    
    • 1.2.4 模型评估(绘制ROC曲线)

    计算AUC取值范围在0.5-1之间,越接近于1越好。可以根据预测结果绘制ROC曲线。

    library(ROCR)
    library(caret)
    #min
    pred_min <- prediction(re[,2], re[,1])
    auc_min = performance(pred_min,"auc")@y.values[[1]]
    perf_min <- performance(pred_min,"tpr","fpr")
    plot(perf_min,colorize=FALSE, col="blue") 
    lines(c(0,1),c(0,1),col = "gray", lty = 4 )
    text(0.8,0.2, labels = paste0("AUC = ",round(auc_min,3)))
    #1se
    pred_1se <- prediction(re[,3], re[,1])
    auc_1se = performance(pred_1se,"auc")@y.values[[1]]
    perf_1se <- performance(pred_1se,"tpr","fpr")
    plot(perf_1se,colorize=FALSE, col="red") 
    lines(c(0,1),c(0,1),col = "gray", lty = 4 )
    text(0.8,0.2, labels = paste0("AUC = ",round(auc_1se,3)))
    

    把两个模型的图画一起

    plot(perf_min,colorize=FALSE, col="blue") 
    plot(perf_1se,colorize=FALSE, col="red",add = T) 
    lines(c(0,1),c(0,1),col = "gray", lty = 4 )
    text(0.8,0.3, labels = paste0("AUC = ",round(auc_min,3)),col = "blue")
    text(0.8,0.2, labels = paste0("AUC = ",round(auc_1se,3)),col = "red")
    
    一般画成正方形,可以自己调整

    ggplot绘图(和上面的图含义一样,ggplot2画的图更好看,绘图也更灵活)

    tpr_min = performance(pred_min,"tpr")@y.values[[1]]
    tpr_1se = performance(pred_1se,"tpr")@y.values[[1]]
    dat = data.frame(tpr_min = perf_min@y.values[[1]],
                     fpr_min = perf_min@x.values[[1]],
                     tpr_1se = perf_1se@y.values[[1]],
                     fpr_1se = perf_1se@x.values[[1]])
    library(ggplot2)
    ggplot() + 
      geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue") + 
      geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
      theme_bw()+
      annotate("text",x = .75, y = .25,
               label = paste("AUC of min = ",round(auc_min,2)),color = "blue")+
      scale_x_continuous(name  = "fpr")+
      scale_y_continuous(name = "tpr")
    
    
    ggplot() + 
      geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue") + 
      geom_line(data = dat,aes(x = fpr_1se, y = tpr_1se),color = "red")+
      geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
      theme_bw()+
      annotate("text",x = .75, y = .25,
               label = paste("AUC of min = ",round(auc_min,2)),color = "blue")+
      annotate("text",x = .75, y = .15,label = paste("AUC of 1se = ",round(auc_1se,2)),color = "red")+
      scale_x_continuous(name  = "fpr")+
      scale_y_continuous(name = "tpr")
    
    • 1.2.5 切割数据构建模型并预测

    对数据进行切割,一组做训练集,使用表达矩阵和临床信息用来构建模型。一组做测试集,输入表达矩阵来对结局进行预测,测序构建的模型预测的结果。

    用R包caret切割数据,生成的结果是一组代表列数的数字,用这些数字来给表达矩阵和meta取子集即可。

    library(caret) #机器学习R包,可以比较科学的拆分数据。
    set.seed(12345679)
    sam<- createDataPartition(meta$event, p = .5,list = FALSE) #根据event来切分,一半一半和7比3都可以
    head(sam)
    
    train <- exprSet[,sam] #切分表达矩阵,sam是一组,非sam是另一组。
    test <- exprSet[,-sam]
    train_meta <- meta[sam,] #切分临床信息
    test_meta <- meta[-sam,]
    
    #看一下,不要让临床信息差的太多
    prop.table(table(train_meta$stage))
    prop.table(table(test_meta$stage)) 
    prop.table(table(test_meta$race)) 
    prop.table(table(train_meta$race)) 
    

    切割后的train数据集建模

    #计算lambda
    x = t(train)
    y = train_meta$event
    cv_fit <- cv.glmnet(x=x, y=y)
    plot(cv_fit)
    #构建模型
    model_lasso_min <- glmnet(x=x, y=y,lambda=cv_fit$lambda.min)
    model_lasso_1se <- glmnet(x=x, y=y,lambda=cv_fit$lambda.1se)
    #挑出基因
    head(model_lasso_min$beta)
    choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]
    choose_gene_1se=rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta)!=0]
    length(choose_gene_min)
    # [1] 18
    length(choose_gene_1se)
    # [1] 3
    

    模型预测
    用训练集构建模型,预测测试集的生死,注意newx参数变了。

    lasso.prob <- predict(cv_fit, newx=t(test), s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
    re=cbind(event = test_meta$event ,lasso.prob)
    re=as.data.frame(re)
    colnames(re)=c('event','prob_min','prob_1se')
    re$event=as.factor(re$event)
    head(re)
    

    再画ROC曲线

    library(ROCR)
    library(caret)
    # 训练集模型预测测试集
    #min
    pred_min <- prediction(re[,2], re[,1])
    auc_min = performance(pred_min,"auc")@y.values[[1]]
    perf_min <- performance(pred_min,"tpr","fpr")
    
    #1se
    pred_1se <- prediction(re[,3], re[,1])
    auc_1se = performance(pred_1se,"auc")@y.values[[1]]
    perf_1se <- performance(pred_1se,"tpr","fpr")
    
    tpr_min = performance(pred_min,"tpr")@y.values[[1]]
    tpr_1se = performance(pred_1se,"tpr")@y.values[[1]]
    dat = data.frame(tpr_min = perf_min@y.values[[1]],
                     fpr_min = perf_min@x.values[[1]],
                     tpr_1se = perf_1se@y.values[[1]],
                     fpr_1se = perf_1se@x.values[[1]])
    
    ggplot() + 
      geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue") + 
      geom_line(data = dat,aes(x = fpr_1se, y = tpr_1se),color = "red")+
      geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
      theme_bw()+
      annotate("text",x = .75, y = .25,
               label = paste("AUC of min = ",round(auc_min,2)),color = "blue")+
      annotate("text",x = .75, y = .15,label = paste("AUC of 1se = ",round(auc_1se,2)),color = "red")+
      scale_x_continuous(name  = "fpr")+
      scale_y_continuous(name = "tpr")
    

    2. Cox多因素回归

    如果Lasso回归挑出的基因数目还是太多,就可以通过Cox多因素回归再进行筛选。
    使用Lasso回归挑出的基因作为Cox多因素回归的输入数据,使用逐步回归法去挑选可选范围内最好的模型。(通常来说做数据挖掘的文章构建模型的话最后一步都是多因素Cox)

    2.1 准备输入数据
    if(!require(My.stepwise))install.packages("My.stepwise")
    load("TCGA-KIRC_sur_model.Rdata") #KIRC整理好的生存数据
    load("lasso_choose_gene_min.Rdata") #Lasso回归选出的35个基因(这里没有区分lnc和mRNA)
    
    2.2 构建coxph模型

    将用于建模的基因(例如lasso回归选中的基因)从表达矩阵中取出来,可作为列添加在meta表格的后面,组成的数据框赋值给dat。

    library(stringr)
    e=t(exprSet[choose_gene_min,]) #从矩阵中取出Lasso回归选出的35个基因
    colnames(e)= str_replace_all(colnames(e),"-","_") #因为miRNA的名字中有减号-,而cox的公式中有加号,会造成干扰,因此需要将基因名中的-替换成_。
    dat=cbind(meta,e)
    
    dat$gender=as.numeric(factor(dat$gender)) #把levels转换成数字作为输入
    dat$stage=as.numeric(factor(dat$stage))
    colnames(dat) #dat矩阵行是样本,列是感兴趣的临床信息和挑出的35个基因
    #   [1] "ID"             "event"          "death"          "last_followup"  "race"          
    #  [6] "age"            "gender"         "stage"          "time"           "age_group"     
    # [11] "hsa_mir_101_1"  "hsa_mir_1179"   "hsa_mir_1185_1" "hsa_mir_1248"   "hsa_mir_1269"  
    # [16] "hsa_mir_1277"   "hsa_mir_1305"   "hsa_mir_130a"   "hsa_mir_130b"   "hsa_mir_144"   
    # [21] "hsa_mir_149"    "hsa_mir_181b_2" "hsa_mir_181c"   "hsa_mir_18a"    "hsa_mir_223"   
    # [26] "hsa_mir_2276"   "hsa_mir_27b"    "hsa_mir_28"     "hsa_mir_3149"   "hsa_mir_34c"   
    # [31] "hsa_mir_3613"   "hsa_mir_3614"   "hsa_mir_365_2"  "hsa_mir_3651"   "hsa_mir_3667"  
    # [36] "hsa_mir_3676"   "hsa_mir_376b"   "hsa_mir_3917"   "hsa_mir_548q"   "hsa_mir_612"   
    # [41] "hsa_mir_627"    "hsa_mir_676"    "hsa_mir_9_2"    "hsa_mir_939"    "hsa_mir_95"    
    

    逐步回归法构建最优模型(从若干个因素中挑选出更关键的)

    library(survival)
    library(survminer)
    library(My.stepwise)
    vl <- colnames(dat)[c(6,7,8,11:ncol(dat))] #选择感兴趣的输入因素,这里选了列名中从age开始到最后一列的35个因素作为输入
    My.stepwise.coxph(Time = "time",
                      Status = "event",
                      variable.list = vl, #从这些因素中挑选更关键的因素,输出结果中最后一个模型是最好的
                      data = dat)
    

    使用输出结果里的最后一个模型(逐步回归法经过若干次迭代后选出的最好的模型)

    # 复制输出结果中的最后一个公式
    model = coxph(formula = Surv(time, event) ~ stage + hsa_mir_223 + age + 
        hsa_mir_34c + hsa_mir_3917 + hsa_mir_3651 + hsa_mir_144 + 
        hsa_mir_2276 + hsa_mir_3667 + hsa_mir_3149 + hsa_mir_627 + 
        hsa_mir_101_1, data = dat)
    

    逐步回归法最重要的是选出了~后面的这些因素,作为它认为的最优因素组合来构建模型,~后的因素也可以通过其它方式挑选,公式还是一样的写法。

    2.3 模型可视化-森林图
    ggforest(model,data = dat)
    
    每一行是一个因素,n=516是用来建模的样本数,第三列和右边的图都是表示SR值,以1为界,小于1的是保护因素,大于1的是危险因素。最后一列是p值。图中有一些p值不显著的也被选中了,但这不影响这是逐步回归法选出的一个最优模型。左下角的Events: 158是说有158个人结局是1。Global p-value: 3.0657e-35,这个是全局的p值,显示这个模型还不错。AIC是衡量模型精确度和简洁度的标准,这个值越小越好。逐步回归法选择模型的标准就是AIC。Concordance index和ROC曲线的AUC值是相同的取值范围相同的作用,只是意义不同,它和AUC一样是评价模型好坏的一个指标。模型不好的话Concordance index就更接近于0.5,好的话就更接近于1。

    森林图绘图更多参考:https://www.jianshu.com/p/58c90b2f3910

    2.4 模型预测
    fp <- predict(model,newdata = dat)
    library(Hmisc) #使用Hmisc包中的函数计算C Index
    options(scipen=200)
    with(dat,rcorr.cens(fp,Surv(time, event))) #用1-返回的C Index值是我们所说的C Index,也就是上面森林图中标出的0.81
    #         C Index             Dxy            S.D.               n         missing      uncensored 
    #      0.19234306     -0.61531388      0.03361968    516.00000000      0.00000000    158.00000000 
    #  Relevant Pairs      Concordant       Uncertain 
    #  89174.00000000  17152.00000000 176550.00000000 
    

    C-index用于计算生存分析中的COX模型预测值与真实之间的区分度(discrimination),也称为Harrell's concordanceindex。C-index在0.5-1之间。0.5为完全不一致,说明该模型没有预测作用,1为完全一致,说明该模型预测结果与实际完全一致。

    2.5 切割数据构建模型并预测
    • 2.5.1 切割数据
      用R包caret切割数据,生成的结果是一组代表列数的数字,用这些数字来给表达矩阵和meta取子集即可。
    library(caret)
    set.seed(12345679)
    sam<- createDataPartition(meta$event, p = .5,list = FALSE)
    train <- exprSet[,sam]
    test <- exprSet[,-sam]
    train_meta <- meta[sam,]
    test_meta <- meta[-sam,]
    
    • 2.5.2 切割后的train数据集建模
      和上面的建模方法一样。
    e=t(train[choose_gene_min,])
    colnames(e)= str_replace_all(colnames(e),"-","_")
    dat=cbind(train_meta,e)
    
    dat$gender=as.numeric(factor(dat$gender))
    dat$stage=as.numeric(factor(dat$stage))
    colnames(dat)
    
    #install.packages("My.stepwise")
    # library(My.stepwise)
    # vl <- colnames(dat)[c(6,7,8,11:ncol(dat))]
    # My.stepwise.coxph(Time = "time",
    #                   Status = "event",
    #                   variable.list = vl,
    #                   data = dat)
    
    model = coxph(formula = Surv(time, event) ~ stage + hsa_mir_223 + age + 
        hsa_mir_34c + hsa_mir_181b_2 + hsa_mir_3614, data = dat)
    
    • 2.5.3 模型可视化
    ggforest(model, data =dat)
    
    • 2.5.4 用切割后的数据test数据集验证模型
    e=t(test[choose_gene_min,])
    colnames(e)= str_replace_all(colnames(e),"-","_")
    test_dat=cbind(test_meta,e)
    test_dat$gender=as.numeric(factor(test_dat$gender))
    test_dat$stage=as.numeric(factor(test_dat$stage))
    
    fp <- predict(model,newdata = test_dat)
    library(Hmisc)
    with(test_dat,rcorr.cens(fp,Surv(time, event)))
    #        C Index            Dxy           S.D.              n        missing 
    #     0.22565299    -0.54869403     0.04919064   258.00000000     0.00000000 
    #     uncensored Relevant Pairs     Concordant      Uncertain 
    #    75.00000000 21440.00000000  4838.00000000 44858.00000000 
    

    3 随机森林

    3.1 准备输入数据

    输入数据是肿瘤样本表达矩阵exprSet和临床信息meta

    load("TCGA-KIRC_sur_model.Rdata")
    library(randomForest)
    library(ROCR)
    library(genefilter)
    library(Hmisc)
    
    3.2 构建随机森林模型

    输入数据是表达矩阵(仅含tumor样本)和对应的生死。(和Lasso回归一样)

    x=t(exprSet)
    y=meta$event
    #构建模型,一个叫randomForest的函数,运行时间很长,存Rdata跳过
    tmp_rf="TCGA_KIRC_miRNA_rf_output.Rdata"
    if(!file.exists(tmp_rf)){
      rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )
      save(rf_output,file = tmp_rf)
    }
    load(file = tmp_rf)
    #top30的基因
    varImpPlot(rf_output, type=2, n.var=30, scale=FALSE, 
               main="Variable Importance (Gini) for top 30 predictors",cex =.7)
    

    随机森林算法筛选基因不像Lasso回归那么绝对,不是直接分是用到了还是没用到,随机森林只能帮我们把基因的重要性进行排名 ,要选前多少个基因是我们自己定的。(这里选了前30,如上图)

    rf_importances=importance(rf_output, scale=FALSE)
    head(rf_importances)
    #                   %IncMSE IncNodePurity
    # hsa-let-7a-1 1.852761e-04     0.1787383
    # hsa-let-7a-2 2.167420e-04     0.1916623
    # hsa-let-7a-3 2.218169e-04     0.1858544
    # hsa-let-7b   7.399404e-05     0.1628646
    # hsa-let-7c   7.658155e-05     0.1635053
    # hsa-let-7d   1.974099e-04     0.2382185
    #解释量top30的基因,和图上是一样的,从下往上看。
    choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30)) #选前30个基因,这个数值可变
    
    3.3 模型预测和评估
    rf.prob <- predict(rf_output, x) #rf.prob是预测值,越接近于0就表示模型预测这个人活着,越接近于1就表示模型预测这个人死了。
    re=cbind(y ,rf.prob) #预测值和真实值放在一起查看比较
    head(re)
    #                              y   rf.prob
    # TCGA-A3-3307-01A-01T-0860-13 0 0.1364447
    # TCGA-A3-3308-01A-02R-1324-13 0 0.1793771
    # TCGA-A3-3311-01A-02R-1324-13 1 0.6709712
    # TCGA-A3-3313-01A-02R-1324-13 1 0.7742376
    # TCGA-A3-3316-01A-01T-0860-13 0 0.2035863
    # TCGA-A3-3317-01A-01T-0860-13 0 0.1619938
    

    ROC曲线(随机森林如果不做另外一个数据的验证,它的AUC值是1,所以绘制ROC曲线没有什么意义)

    library(ROCR)
    #library(caret)
    
    pred <- prediction(re[,2], re[,1])
    auc = performance(pred,"auc")@y.values[[1]]
    auc
    # [1] 1
    
    3.4 切割数据构建模型并预测
    • 3.4.1 切割数据

    用R包caret切割数据,生成的结果是一组代表列数的数字,用这些数字来给表达矩阵和meta取子集即可。

    library(caret)
    set.seed(12345679)
    sam<- createDataPartition(meta$event, p = .5,list = FALSE)
    
    train <- exprSet[,sam]
    test <- exprSet[,-sam]
    train_meta <- meta[sam,]
    test_meta <- meta[-sam,]
    
    • 3.4.2 切割后的train数据集建模

    和上面的建模方法一样。

    #建模
    x = t(train)
    y = train_meta$event
    tmp_rf="TCGA_KIRC_miRNA_t_rf_output.Rdata"
    if(!file.exists(tmp_rf)){
      rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )
      save(rf_output,file = tmp_rf)
    }
    load(file = tmp_rf)
    
    choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30))
    head(choose_gene)
    # [1] "hsa-mir-511-1"  "hsa-mir-155"    "hsa-mir-409"    "hsa-mir-1185-1"
    # [5] "hsa-mir-1277"   "hsa-mir-149"
    
    3.5 模型预测

    用训练集构建模型,预测测试集的生死。

    x = t(test)
    y = test_meta$event
    rf.prob <- predict(rf_output, x)
    re=cbind(y ,rf.prob)
    re=as.data.frame(re)
    colnames(re)=c('event','prob')
    re$event=as.factor(re$event)
    

    再看AUC值。

    library(ROCR)
    # 训练集模型预测测试集
    pred <- prediction(re[,2], re[,1])
    auc= performance(pred,"auc")@y.values[[1]]
    auc
    # [1] 0.7121311
    

    模型还可以。

    4 支持向量机(Support Vector Machine, SVM)

    前三种预测方法都是可以把建模时所 使用的基因挑出来,也可以写出公式,SVM比前三种方法要简单很多。

    4.1 准备输入数据
    load("TCGA-KIRC_sur_model.Rdata")
    library(ROCR)
    library(genefilter)
    library(Hmisc)
    library(e1071) #SVM的包名
    
    4.2 构建支持向量机模型
    • 4.2.1 切割数据
      用R包caret切割数据,生成的结果是一组代表列数的数字,用这些数字来给表达矩阵和meta取子集即可。
    library(caret)
    set.seed(12345679)
    sam<- createDataPartition(meta$event, p = .5,list = FALSE)
    
    train <- exprSet[,sam]
    test <- exprSet[,-sam]
    train_meta <- meta[sam,]
    test_meta <- meta[-sam,]
    
    • 4.2.2 train数据集建模
    x=t(train)
    y=as.factor(train_meta$event)
    model = svm(x,y,kernel = "linear")
    summary(model) 
    # 
    # Call:
    # svm.default(x = x, y = y, kernel = "linear")
    # 
    # 
    # Parameters:
    #    SVM-Type:  C-classification 
    #  SVM-Kernel:  linear 
    #        cost:  1 
    # 
    # Number of Support Vectors:  176
    # 
    #  ( 108 68 )
    # 
    # 
    # Number of Classes:  2 
    # 
    # Levels: 
    #  0 1
    
    • 4.2.3 模型预测
      用训练集构建模型,预测测试集的生死。不同于其他模型,这个预测结果是分类变量,直接预测生死,而不是prob。
    x=t(test)
    y=as.factor(test_meta$event)
    pred = predict(model, x)
    table(pred,y)
    #     y
    # pred   0   1
    #    0 142  47
    #    1  41  28
    

    上面的结果列的0和1是预测的存活和死亡,行的0和1是真实的存活和死亡。预测值是0的有142+41个,预测值是1的有47+28个。真实是0的有142+47个,真实是1的有41+28个。41(false negative)和47(false positive)就属于误判。


    代码来自2021生信技能树数据挖掘课

    相关文章

      网友评论

        本文标题:TCGA生存模型的构建以及模型预测和评估

        本文链接:https://www.haomeiwen.com/subject/lwidyltx.html