实现影像组学全流程

admin管理员组

文章数量:1794759

实现影像组学全流程

对一篇影像组学的的论文（《Development and validation of an MRI-based radiomics nomogram for distinguishing Warthin’s tumour from pleomorphic adenomas of the parotid gland》）中方法进行复现。完整地跑通影像组学全流程，对临床+影像组学特征进行建模并绘制Lasso回归图和ROC曲线、诺模图、校准曲线、决策曲线等。 需要完整代码到这里aistudio.baidu/aistudio/projectdetail/3270880 fork项目之后下载glioma.7z压缩包，然后用RStudio的RMarkdown进行运行

1.项目介绍

2012年荷兰学者提出影像组学(radiomics）概念：借助计算机，从医学影像中挖掘海量定量影像特征，使用统计学/机器学习方法，筛选出最有价值的影像学特征，用来解析临床信。 可以说影像组学是人工智能在医学领域的一种特定研究方式。

做影像组学一般会经过以下步骤1.收集数据，2.标注数据，3.特征提取，4.特征工程，5.模型设计，6.结果分析和绘图。

1.1 项目目的

而对于大多数医学专业的朋友来说，毕竟学医的对写代码不一定擅长，或者第一次接触影像组学。因此想打算写一个影像组学相关的项目，对某篇组学论文方法进行复现。通用这个项目的代码希望可以稍微帮助到他们，可以让他们把时间更多花在疾病的分析上，而不是解决某个bug花大量的时间。

1.2 项目大概情况

这个项目主要是对《Development and validation of an MRI-based radiomics nomogram for distinguishing Warthin’s tumour from pleomorphic adenomas of the parotid gland》这篇影像组学论文中的方法进行复现，完整的跑一个影像组学流程。包括：

1.对临床特征进行建模

2.提取影像组学特征，通过LassoCv进行特征筛选再建模

3.结合影像组学特征+age临床特征进行建模

4.绘制三个模型的ROC曲线进行对比

5.对组合模型绘制诺模图+决策曲线+校准曲线

1.3运行环境

由于有些图例，例如诺模图(列线图)需要用到R语言来绘制，但是提取特征用到python语言的Pyradiomics库，所以这个项目是一个混合编程语言的项目。刚好RMarkdown可以同时运行Python和R语言，所以请有需要运行项目的请离线安装RMarkdown来运行。

虽然不能通过在线BML运行项目，但是可以展示复现论文方法的代码，主要可以容易复制粘贴代码~~~ 。

注：对于如何安装python和安装R语言和如何使用RMarkdown不是这篇项目的范围，自己通过百度自学学习。

注：这个项目也是我在学习影像组学中做下的笔记，如有错误请纠正和谅解

2.关于数据

虽然是论文复现，但是论文中涉及的腮腺MR数据是不公开的，无法和论文实验数据一致，不过这个项目主要是为了实现论文中大部分方法。因此采用了一个公开分割比赛的数据集作为这个项目的数据。

在数据上用到的是一个公开的胶质瘤数据集–【BraTS2019】。是一个两分类的数据集，Hgg是高级别胶质瘤，Lgg是低级别的胶质瘤，数量上分别是240和76。每个病例有四个模态的磁共振图像。Hgg患者是带有age的年龄数据，但是lgg的患者是没有的。为了不影响临床+影像组学的建模。在（25到55）之间随机生成一个随机数赋予个Lgg患者作为他的年龄。所以Lgg患者的年龄是虚拟数据。提取影像组学的特征只用到T1增强的模态。

因为原始的数据集Hgg和Lgg的比例差不多达到3比1，所以为了搭到数据均衡一点，把Hgg删掉一部分，剩下101例。Lgg为76例。

因此在实际运用到自己的数据的时候，收集病例尽量做到数据均衡一点。不然有点自己坑自己的感觉。

3.复现的论文 《Development and validation of an MRI-based radiomics nomogram for distinguishing Warthin’s tumour from pleomorphic adenomas of the parotid gland》 3.1论文地址

论文地址

3.2论文大概内容

这篇论文主要对腮腺肿瘤（Warthin和多形性腺瘤）的磁共振图像使用计算机提取大量高通量的特征，通过对特征进行筛选和建立模型来鉴别MR图像是Warthin还是多形性腺瘤。

两种腺瘤的治疗方法和预后不同，建立一个良好的分类模型对临床来说至关重要。论文通过建立三种不同特征的模型，分别是1.单纯临床特征模型，2.单纯影像组学特征模型，3.影像组学特征+年龄模型。经过实验得到模型3的auc最好。并根据模型3构建影像组学列线图（如下图），使预测模型的结果更具有可读性，让模型更加方便运用实际临床决策之中。

3.2论文中的方法流程：

4.开始影像组学 4.1目录结构和数据准备

重要的一点就是临床特征的index列需要和每一例病例文件夹要一一对应

4.2 特征提取

特征提取是采用pyradiomics库进行提取，可以通过配置文件来设置提取那些特征，例如设置提取特征前是否重采样，设置滤波核的大小等等。如下图

可以到pyradiomics官网查看地址

""" 经过一两个小时提取后会生成HGG.csv和LGG.csv文件, 生成的csv文件每一行都有接近一千个特征，数量会根据不同yaml文件设置不同而不同 """ #导入相关的库 import sys import pandas as pd import os import random import shutil import numpy as np import radiomics from radiomics import featureextractor import SimpleITK as sitk kinds = ['HGG','LGG'] #这个是特征处理配置文件，具体可以参考pyradiomics官网 para_path = 'yaml/MR_1mm.yaml' extractor = featureextractor.RadiomicsFeatureExtractor(para_path) dir = 'data/MyData/' for kind in kinds: print("{}:开始提取特征".format(kind)) features_dict = dict() df = pd.DataFrame() path = dir + kind # 使用配置文件初始化特征抽取器 for index, folder in enumerate( os.listdir(path)): for f in os.listdir(os.path.join(path, folder)): if 't1ce' in f: ori_path = os.path.join(path,folder, f) break lab_path = ori_path.replace('t1ce','seg') features = extractor.execute(ori_path,lab_path) #抽取特征 #新增一列用来保存病例文件夹名字 features_dict['index'] = folder for key, value in features.items(): #输出特征 features_dict[key] = value df = df.append(pd.DataFrame.from_dict(features_dict.values()).T,ignore_index=True) print(index) df.columns = features_dict.keys() df.to_csv('csv/' +'{}.csv'.format(kind),index=0) print('Done') print("完成") """ 再对HGG.csv和LGG.csv文件进行处理,去掉字符串特征，插入label标签。 HGG标签为1，LGG标签为0 """ import matplotlib.pyplot as plt import seaborn as sns hgg_data = pd.read_csv('csv/HGG.csv') lgg_data = pd.read_csv('csv/LGG.csv') hgg_data.insert(1,'label', 1) #插入标签 lgg_data.insert(1,'label', 0) #插入标签 #因为有些特征是字符串，直接删掉 cols=[x for i,x in enumerate(hgg_data.columns) if type(hgg_data.iat[1,i]) == str] cols.remove('index') hgg_data=hgg_data.drop(cols,axis=1) cols=[x for i,x in enumerate(lgg_data.columns) if type(lgg_data.iat[1,i]) == str] cols.remove('index') lgg_data=lgg_data.drop(cols,axis=1) #再合并成一个新的csv文件。 total_data = pd.concat([hgg_data, lgg_data]) total_data.to_csv('csv/TotalOMICS.csv',index=False) #简单查看数据的分布 fig, ax = plt.subplots() sns.set() ax = sns.countplot(x='label',hue='label',data=total_data) plt.show() print(total_data['label'].value_counts())

4.3划分数据

对影像组学数据集进行划分训练集合测试集。比例为8:2。同样把临床特征也进行同样的划分（顺序和影响组学特征是一致的）。

#导入常用R包 library(glmnet) library(rms) library(foreign) library(ggplot2) library(pROC) #设置种子为了保证每次结果都一样 set.seed(888) data <- read.csv("csv/TotalOMICS.csv") nn=0.8 print(paste('总样本数:',length(data[,1]))) sub<-sample(1:nrow(data),round(nrow(data)*nn)) trainOmics<-data[sub,]#取0.8的数据做训练集 testOmics<-data[-sub,]#取0.2的数据做测试集 print(paste('训练集样本数:',length(trainOmics[,1]))) print(paste('测试集样本数:',length(testOmics[,1]))) write.csv(trainOmics,"csv/trainOmics.csv",row.names = FALSE ) write.csv(testOmics,"csv/testOmics.csv",row.names = FALSE ) """ R语言输出结果 [1] "总样本数: 177" [1] "训练集样本数: 142" [1] "测试集样本数: 35" """ """ 根据上面对影像组学TotalOMICS.csv的数据划分，对TotalClinic.csv同样的顺序划分 """ import pandas as pd def compare(data1, data2,filename): # 读取两个表 dt1 = pd.read_csv(data1,encoding='utf-8') dt2 = pd.read_csv(data2,encoding='gb18030') dt2.head() df = pd.DataFrame() dt1_name = dt1['index'].values.tolist() dt2_name = dt2['index'].values.tolist() for i in dt1_name: if i in dt2_name: dt2_row = dt2.loc[dt2['index'] == i] df = df.append(dt2_row) df.to_csv('./csv/'+filename+'.csv',header=True,index=False,encoding="utf_8_sig") data_train= "./csv/trainOmics.csv" data_test = "./csv/testOmics.csv" data_clinic= "./csv/TotalClinic.csv" compare(data_train,data_clinic,"trainClinic") compare(data_test,data_clinic,"testClinic") 4.4对单纯临床特征进行建模

因为这个数据集只有age这个临床特征，所以没有做进一步的特征分析和筛选，直接用age进行建模

4.4.1处理临床数据 #读取临床数据trainClinic.csv trainClinic<-read.csv("./csv/trainClinic.csv",fileEncoding = "UTF-8-BOM") trainClinic<- data.frame(trainClinic) #把age变量转换成数值型。 trainClinic$Age <- as.numeric(trainClinic$Age) #打印临床特征名 names(trainClinic) """ R语言输出结果 [1] "index" "Age" "Label" """ 4.4.2 使用逻辑回归对临床特征进行建模 #建模，使用逻辑回归 model.Clinic<-glm(Label~Age,data = trainClinic,family=binomial(link='logit')) summary(model.Clinic) """ R语言输出结果 Call: glm(formula = Label ~ Age, family = binomial(link = "logit"), data = trainClinic) Deviance Residuals: Min 1Q Median 3Q Max -1.7106 -0.6673 0.1870 0.6294 2.6505 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -6.6037 1.1252 -5.869 4.39e-09 *** Age 0.1419 0.0233 6.090 1.13e-09 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 194.57 on 141 degrees of freedom Residual deviance: 121.48 on 140 degrees of freedom AIC: 125.48 Number of Fisher Scoring iterations: 5 """ 4.4.2查看模型在训练集的验证情况 probClinicTrain<-predict.glm(object =model.Clinic,newdata=trainClinic,type = "response") predClinicTrain<-ifelse(probClinicTrain>=0.5,1,0) #计算模型精度 error=predClinicTrain-trainClinic$Label #accuracy:判断正确的数量占总数的比例 accuracy=(nrow(trainClinic)-sum(abs(error)))/nrow(trainClinic) #precision:真实值预测值全为1 / 预测值全为1 --- 提取出的正确信条数/提取出的信条数 precision=sum(trainClinic$Label & predClinicTrain)/sum(predClinicTrain) #recall:真实值预测值全为1 / 真实值全为1 --- 提取出的正确信条数 /样本中的信条数 recall=sum(predClinicTrain & trainClinic$Label)/sum(trainClinic$Label) #P和R指标有时候会出现的矛盾的情况，这样就需要综合考虑他们，最常见的方法就是F-Measure（又称为F-Score） #F-Measure是Precision和Recall加权调和平均，是一个综合评价指标 F_measure=2*precision*recall/(precision+recall) #输出以上各结果 #精确率和召回率和F_measure是预测Hgg的值 #可以模型预测HGG的能力比较强，但是预测Lgg的能力比较弱 print(paste('准确率accuracy:',accuracy)) print(paste('精确率precision:',precision)) print(paste('召回率recall:',recall)) print(paste('F_measure:',F_measure)) table(trainClinic$Label,predClinicTrain) """ R语言输出结果 [1] "准确率accuracy: 0.78169014084507" [1] "精确率precision: 0.795180722891566" [1] "召回率recall: 0.825" [1] "F_measure: 0.809815950920245" predClinicTrain 0 1 0 45 17 1 14 66 """ 4.4.3 查看模型在测试集的验证情况 testClinic<-read.csv("./csv/testClinic.csv",fileEncoding = "UTF-8-BOM") testClinic$Age <- as.numeric(testClinic$Age) probClinicTest<-predict.glm(object =model.Clinic,newdata=testClinic,type = "response") predClinicTest<-ifelse(probClinicTest>=0.5,1,0) error=predClinicTest-testClinic$Label accuracy=(nrow(testClinic)-sum(abs(error)))/nrow(testClinic) precision=sum(testClinic$Label & predClinicTest)/sum(predClinicTest) recall=sum(predClinicTest & testClinic$Label)/sum(testClinic$Label) F_measure=2*precision*recall/(precision+recall) print(paste('准确率accuracy:',accuracy)) print(paste('精确率precision:',precision)) print(paste('召回率recall:',recall)) print(paste('F_measure:',F_measure)) table(testClinic$Label,predClinicTest) """ R语言输出结果 [1] "准确率accuracy: 0.828571428571429" [1] "精确率precision: 0.857142857142857" [1] "召回率recall: 0.857142857142857" [1] "F_measure: 0.857142857142857" predClinicTest 0 1 0 11 3 1 3 18 """ 4.5 对单纯影像组学进行建模

先对影像组学特征进行T检验筛选一部分特征，然后用lasso回归进一步筛选特征，最后剩下5个特征用来建立逻辑回归方程

4.5.1先做T检验筛选一部分特征 """ 用python做T检验,筛选后剩下的特征数：562个 """ from scipy.stats import levene, ttest_ind tData = pd.read_csv('./csv/trainOmics.csv') classinformation = tData["label"].unique() for temp_classinformation in classinformation: temp_data = tData[tData['label'].isin([temp_classinformation])] exec("df%s=temp_data"%temp_classinformation) counts = 0 columns_index =[] for column_name in tData.columns[2:]: if levene(df1[column_name], df0[column_name])[1] > 0.05: if ttest_ind(df1[column_name],df0[column_name],equal_var=True)[1] < 0.05: columns_index.append(column_name) else: if ttest_ind(df1[column_name],df0[column_name],equal_var=False)[1] < 0.05: columns_index.append(column_name) print("筛选后剩下的特征数：{}个".format(len(columns_index))) #print(columns_index) #数据只保留从T检验筛选出的特征数据，重新组合成data if not 'label' in columns_index: columns_index = ['label'] + columns_index if not 'index' in columns_index: columns_index = ['index'] + columns_index df1 = df1[columns_index] df0 = df0[columns_index] tData = pd.concat([df1, df0]) tData.to_csv('./csv/tData_train.csv',header=True,index=False,encoding="utf_8_sig") 4.5.2再用lasso回归进一步筛选特征 tData_train <- read.csv("csv/tData_train.csv",fileEncoding = "UTF-8-BOM") dim(tData_train) Y <-as.data.frame(tData_train$label) #[,-1]是为了去掉截距 Y <- model.matrix(~.,data=Y)[,-1] #除去因变量，提取自变量 yavars<-names(tData_train) %in% c("label","index") X <- as.data.frame(tData_train[!yavars]) X <- model.matrix(~.,data=X)[,-1] #Lasso回归 fit <- glmnet(X,Y, alpha=1, family = "binomial") plot(fit, xvar="lambda", label=TRUE)

这个图显示随着lambda变化，特征数量的变化

#采用k折交叉验证进行lasso回归 cv.fit <- cv.glmnet(X,Y, alpha=1,nfolds = 10,family="binomial") plot(cv.fit) abline(v=log(c(cv.fit$lambda.min, cv.fit$lambda.lse)), lty=2) #把lambda取1倍标准误时的值绘制到图中 plot(cv.fit$glmnet.fit,xvar="lambda") abline(v=log(cv.fit$lambda.1se), lty=2,)

这个图显示随着lambda增大，MSE的变化，右边的垂直虚线是1倍标准误时lambda的取值。

4.5.3经过lasso回归筛选抽出5个特征

分别是

[2] “wavelet.LHL_glcm_Imc2” [3] “wavelet.LHL_glszm_ZoneEntropy” [4] “wavelet.HLL_glszm_ZoneEntropy” [5] “wavelet.LLL_firstorder_90Percentile” [6] “wavelet.LLL_firstorder_RootMeanSquared”

#如果取1倍标准误时,获取筛选后的特征 lambda = cv.fit$lambda.1se Coefficients <- coef(fit, s = lambda) Active.Index <- which(Coefficients != 0) Active.Coefficients <- Coefficients[Active.Index] Active.Index Active.Coefficients row.names(Coefficients)[Active.Index] """ R语言输出结果 [1] 1 273 300 396 507 518 [1] -1.135967e+01 5.551516e+00 5.053583e-01 4.978627e-01 4.710207e-04 [6] 1.236226e-03 [1] "(Intercept)" [2] "wavelet.LHL_glcm_Imc2" [3] "wavelet.LHL_glszm_ZoneEntropy" [4] "wavelet.HLL_glszm_ZoneEntropy" [5] "wavelet.LLL_firstorder_90Percentile" [6] "wavelet.LLL_firstorder_RootMeanSquared" """ 4.6 对5个影像组学特征建立逻辑回归 #建立公式 formulalse <-as.formula(label ~wavelet.LHL_glcm_Imc2+wavelet.LHL_glszm_ZoneEntropy+wavelet.HLL_glszm_ZoneEntropy+wavelet.LLL_firstorder_90Percentile+wavelet.LLL_firstorder_RootMeanSquared) #逻辑回归 model.Omics <- glm(formula=formulalse,data=tData_train,family=binomial(link="logit")) #查看结果 summary(model.Omics) """ R语言输出结果 Call: glm(formula = formulalse, family = binomial(link = "logit"), data = tData_train) Deviance Residuals: Min 1Q Median 3Q Max -3.8353 -0.3845 0.1660 0.3471 2.3193 Coefficients: Estimate Std. Error z value (Intercept) -31.638928 7.898076 -4.006 wavelet.LHL_glcm_Imc2 14.026526 5.980053 2.346 wavelet.LHL_glszm_ZoneEntropy 0.428890 1.381598 0.310 wavelet.HLL_glszm_ZoneEntropy 2.143145 1.333109 1.608 wavelet.LLL_firstorder_90Percentile 0.002704 0.003126 0.865 wavelet.LLL_firstorder_RootMeanSquared 0.004004 0.004434 0.903 Pr(>|z|) (Intercept) 6.18e-05 *** wavelet.LHL_glcm_Imc2 0.019 * wavelet.LHL_glszm_ZoneEntropy 0.756 wavelet.HLL_glszm_ZoneEntropy 0.108 wavelet.LLL_firstorder_90Percentile 0.387 wavelet.LLL_firstorder_RootMeanSquared 0.367 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 194.566 on 141 degrees of freedom Residual deviance: 84.983 on 136 degrees of freedom AIC: 96.983 Number of Fisher Scoring iterations: 6 """ 4.6.1 查看在训练集的情况 probOmicsTrain<-predict.glm(object =model.Omics,newdata=tData_train,type = "response") predOmicsTrain<-ifelse(probOmicsTrain>=0.5,1,0) error=predOmicsTrain-tData_train$label accuracy=(nrow(tData_train)-sum(abs(error)))/nrow(tData_train) precision=sum(tData_train$label & predOmicsTrain)/sum(predOmicsTrain) recall=sum(predOmicsTrain & tData_train$label)/sum(tData_train$label) F_measure=2*precision*recall/(precision+recall) print(paste('准确率accuracy:',accuracy)) print(paste('精确率precision:',precision)) print(paste('召回率recall:',recall)) print(paste('F_measure:',F_measure)) table(tData_train$label,predOmicsTrain) """ R语言输出结果 [1] "准确率accuracy: 0.908450704225352" [1] "精确率precision: 0.935064935064935" [1] "召回率recall: 0.9" [1] "F_measure: 0.917197452229299" predOmicsTrain 0 1 0 57 5 1 8 72 """ 4.6.2 查看模型在测试集的结果 tData_test <- read.csv("csv/testOmics.csv",fileEncoding = "UTF-8-BOM") probOmicsTest<-predict.glm(object =model.Omics,newdata=tData_test,type = "response") predOmicsTest<-ifelse(probOmicsTest>=0.5,1,0) error=predOmicsTest-tData_test$label accuracy=(nrow(tData_test)-sum(abs(error)))/nrow(tData_test) precision=sum(tData_test$label & predOmicsTest)/sum(predOmicsTest) recall=sum(predOmicsTest & tData_test$label)/sum(tData_test$label) F_measure=2*precision*recall/(precision+recall) print(paste('准确率accuracy:',accuracy)) print(paste('精确率precision:',precision)) print(paste('召回率recall:',recall)) print(paste('F_measure:',F_measure)) table(tData_test$label,predOmicsTest) """ R语言输出结果 [1] "准确率accuracy: 0.857142857142857" [1] "精确率precision: 0.863636363636364" [1] "召回率recall: 0.904761904761905" [1] "F_measure: 0.883720930232558" predOmicsTest 0 1 0 11 3 1 2 19 """ 4.7结合临床特征和影像组学特征进行建模

先处理数据，把临床特征的age赋值到影像组学中新增age一列。然后根据lasso模型计算每个病例的得分。把这个得分和age作为新的特征建立逻辑回归模型。

4.7.1处理数据 #读取经过T检验筛选后的影像组学训练集文件 tData_train2 <- read.csv("csv/tData_train.csv",fileEncoding = "UTF-8-BOM") #把临床特征的Age放到影像组学中作为新的一列 #放的时候是根据病例名一一对应的 tData_train2$Age <- trainClinic$Age for (i in trainClinic$index){ if(tData_train2[tData_train2$index == i,]$index == i){ Age <- trainClinic[trainClinic$index == i,]$Age tData_train2[tData_train2$index == i,]$Age <- Age } } #对测试集也这样处理 tData_test2 <- read.csv("csv/testOmics.csv",fileEncoding = "UTF-8-BOM") testClinic2<-read.csv("./csv/testClinic.csv",fileEncoding = "UTF-8-BOM") testClinic2<- data.frame(testClinic2) testClinic2$Age <- as.numeric(testClinic2$Age) columns <- colnames(tData_train2)[1:564]#经过T检验筛选后的562个特征+index+label tData_test2 <- as.data.frame(tData_test2[,columns]) tData_test2$Age <- testClinic2$Age for (i in testClinic2$index){ if(tData_test2[tData_test2$index == i,]$index == i){ Age <- testClinic2[testClinic2$index == i,]$Age tData_test2[tData_test2$index == i,]$Age <- Age } } 4.7.2根据lasso模型计算得分 #分别计算影像组学得分（RS） y_vad <-as.data.frame(tData_test2$label) y_vad <- model.matrix(~.,data=y_vad)[,-1] #除去因变量，提取自变量 yavars<-names(tData_test2) %in% c("label") x_vad <- as.data.frame(tData_test2[!yavars]) columns <- colnames(tData_train2)[3:564]#T检验筛选后的564个特征不要index+label x_vad <- as.data.frame(tData_test2[,columns]) x_vad <- model.matrix(~.,data=x_vad)[,-1] x_dev <- X y_dev <- Y #fit是lassoCV模型 tData_train2$RS <-predict(fit,type="link", newx=x_dev,newy=y_dev,s=cv.fit$lambda.1se)#训练集的RS tData_test2$RS<-predict(fit,type="link", newx=x_vad,newy=y_vad,cv.fit$lambda.1se)#测试集的RS 4.7.3建立逻辑回归模型 #因为后续要绘制诺模图，所以用rms这包建立逻辑回归 #通过下列例子，发现两种的逻辑回归结果是一样的 tData_train2$RS <- as.numeric(tData_train2$RS) model.and1 = glm(label ~ RS+Age,data=tData_train2,binomial(link='logit')) print(model.and1$coef) print("#####") model.and2 <- lrm(label ~ RS+Age,data=tData_train2,x=TRUE,y=TRUE) print(model.and2$coef) """ R语言输出结果 (Intercept) RS Age -5.8730358 2.4543423 0.1192463 [1] "#####" Intercept RS Age -5.8730331 2.4543415 0.1192462 """ 4.7.4查看模型在临床结合影像组学的训练集的情况 probOmicsTrainAnd<-predict.glm(object =model.and1,newdata=tData_train2,type = "response") predOmicsTrainAnd<-ifelse(probOmicsTrainAnd>=0.5,1,0) error=predOmicsTrainAnd-tData_train2$label accuracy=(nrow(tData_train2)-sum(abs(error)))/nrow(tData_train2) precision=sum(tData_train2$label & predOmicsTrainAnd)/sum(predOmicsTrainAnd) recall=sum(predOmicsTrainAnd & tData_train2$label)/sum(tData_train2$label) F_measure=2*precision*recall/(precision+recall) print(paste('准确率accuracy:',accuracy)) print(paste('精确率precision:',precision)) print(paste('召回率recall:',recall)) print(paste('F_measure:',F_measure)) table(tData_train2$label,predOmicsTrainAnd) """ R语言输出结果 [1] "准确率accuracy: 0.929577464788732" [1] "精确率precision: 0.948717948717949" [1] "召回率recall: 0.925" [1] "F_measure: 0.936708860759494" predOmicsTrainAnd 0 1 0 58 4 1 6 74 """ 4.7.5 查看模型在临床结合影像组学的测试集的情况 tData_test2$RS <- as.numeric(tData_test2$RS) probOmicsTestAnd<-predict.glm(object =model.and1,newdata=tData_test2,type = "response") predOmicsTestAnd<-ifelse(probOmicsTestAnd>=0.5,1,0) error=predOmicsTestAnd-tData_test2$label accuracy=(nrow(tData_test2)-sum(abs(error)))/nrow(tData_test2) precision=sum(tData_test2$label & predOmicsTestAnd)/sum(predOmicsTestAnd) recall=sum(predOmicsTestAnd & tData_test2$label)/sum(tData_test2$label) F_measure=2*precision*recall/(precision+recall) print(paste('准确率accuracy:',accuracy)) print(paste('精确率precision:',precision)) print(paste('召回率recall:',recall)) print(paste('F_measure:',F_measure)) table(tData_test2$label,predOmicsTestAnd) """ R语言输出结果 [1] "准确率accuracy: 0.857142857142857" [1] "精确率precision: 0.9" [1] "召回率recall: 0.857142857142857" [1] "F_measure: 0.878048780487805" predOmicsTestAnd 0 1 0 12 2 1 3 18 """ 4.8绘制ROC

把三个模型都绘制到同一个ROC曲线里面，进行比较

4.8.1绘制训练集的ROC曲线 #计算roc和auc rocClinic <- roc(trainClinic$Label, probClinicTrain) rocOmics <- roc(tData_train$label, probOmicsTrain) rocClinicAndOmics <-roc(tData_train2$label, probOmicsTrainAnd) # 绘图 plot(rocClinic, print.auc=TRUE, # 图像上输出AUC的值 print.auc.x=0.4, print.auc.y=0.6, # 设置AUC值坐标为（x，y） auc.polygon=TRUE, # 将ROC曲线下面积转化为多边形 auc.polygon.col="#fff7f7", # 设置ROC曲线下填充色 grid=FALSE, print.thres=FALSE, # 图像上输出最佳截断值 main=" Train ROC curves", # 添加图形标题 col="red", # 设置ROC曲线颜色 legacy.axes=TRUE,# 使x轴从0到1，表示为1-特异度 xlim=c(1,0), mgp=c(1.5, 1, 0), lty=3) # 再添加1条ROC曲线 plot.roc(rocOmics, add=TRUE, # 增加曲线, col="green", # 设置ROC曲线颜色 print.thres=FALSE, # 图像上输出最佳截断值 print.auc=TRUE, # 图像上输出AUC print.auc.x=0.4,print.auc.y=0.5,# AUC的坐标为（x，y） lty=2) # 再添加1条ROC曲线 plot.roc(rocClinicAndOmics, add=TRUE, # 增加曲线, col="blue", # 设置ROC曲线颜色 print.thres=FALSE, # 图像上输出最佳截断值 print.auc=TRUE, # 图像上输出AUC print.auc.x=0.4,print.auc.y=0.4,# AUC的坐标为（x，y） lty=1) # 添加图例 legend(0.45, 0.30, # 图例位置x，y bty = "n", # 图例样式 legend=c("Clinical model","Radiomics signature","Clinical and Radiomics"), # 添加分组 col=c("red","green","blue"), # 颜色跟前面一致 lwd=1, lty=c(3,2,1)) # 线条粗细

4.8.1绘制测试集的ROC曲线 rocClinicTest <- roc(testClinic$Label, probClinicTest) rocOmicsTest <- roc(tData_test$label, probOmicsTest) rocClinicAndOmicsTest <-roc(tData_test2$label, probOmicsTestAnd) plot(rocClinicTest, print.auc=TRUE, # 图像上输出AUC的值 print.auc.x=0.4, print.auc.y=0.6, # 设置AUC值坐标为（x，y） auc.polygon=TRUE, # 将ROC曲线下面积转化为多边形 auc.polygon.col="#fff7f7", # 设置ROC曲线下填充色 grid=FALSE, print.thres=FALSE, # 图像上输出最佳截断值 main=" Test ROC curves", # 添加图形标题 col="red", # 设置ROC曲线颜色 legacy.axes=TRUE,# 使x轴从0到1，表示为1-特异度 xlim=c(1,0), mgp=c(1.5, 1, 0), lty=3) # 再添加1条ROC曲线 plot.roc(rocOmicsTest, add=TRUE, # 增加曲线, col="green", # 设置ROC曲线颜色 print.thres=FALSE, # 图像上输出最佳截断值 print.auc=TRUE, # 图像上输出AUC print.auc.x=0.4,print.auc.y=0.5,# AUC的坐标为（x，y） lty=2) # 再添加1条ROC曲线 plot.roc(rocClinicAndOmicsTest, add=TRUE, # 增加曲线, col="blue", # 设置ROC曲线颜色 print.thres=FALSE, # 图像上输出最佳截断值 print.auc=TRUE, # 图像上输出AUC print.auc.x=0.4,print.auc.y=0.4,# AUC的坐标为（x，y） lty=1) # 添加图例 legend(0.45, 0.30, # 图例位置x，y bty = "n", # 图例样式 legend=c("Clinical model","Radiomics signature","Clinical and Radiomics"), # 添加分组 col=c("red","green","blue"), # 颜色跟前面一致 lwd=1, lty=c(3,2,1)) # 线条粗细

4.8绘制诺模图

临床结合影像组学建立的逻辑回归模型整体来说比较良好。 把这个模型绘制成诺模图，使模型可视化。

library(rms) formula <- as.formula(label ~ Age+RS) #数据打包 dd = datadist(tData_train2) options(datadist="dd") fitnomogram <- lrm(formula,data=tData_train2, x=TRUE, y=TRUE) nom1 <- nomogram(fitnomogram,#模型对象 fun=function(x)1/(1+exp(-x)),#保持不变，logistic固定 lp=F,#是否显示线性概率 fun.at=c(0.1,0.2,0.5,0.7,0.9),#预测结果坐标轴 funlabel="Risk")#坐标轴名称 #可以使用Cairo导出pdf #library(Cairo) #CairoPDF("nomogram.pdf",width=6,height=6) plot(nom1)

诺模图怎样看？现在知道某个患者的age师多少岁，找到对应顶部Points的分数，同样知道患者的RS得分，也找到对应Points的得分，把两个等分加起来。在Total Points找到底部Risk得分，这个得分就是判断患者是HGG的概率。

4.9绘制校准曲线 cal1 <- calibrate(fitnomogram, cmethod="hare", method="boot", B=1000) plot(cal1, xlim=c(0,1.0), ylim=c(0,1.0), xlab = "Predicted Probability", ylab = "Observed Probability", legend = FALSE,subtitles = FALSE) #abline对角线 abline(0,1,col = "black",lty = 2,lwd = 2) #再画一条模型预测的实际曲线 lines(cal1[,c("predy","calibrated.orig")], type = "l",lwd = 2,col="red",pch =16) #再画一条模型Bias-corrected是校准曲线 lines(cal1[,c("predy","calibrated.corrected")], type = "l",lwd = 2,col="green",pch =16) legend(0.55,0.35, c("Ideal","Apparent","Bias-corrected"), lty = c(2,1,1), lwd = c(2,1,1), col = c("black","red","green"), bty = "n") # "o"为加边框

彩色的曲线越接近点斜线，模型的性能就比较好

4.10 绘制决策曲线

把三个模型都绘制在同一个图中，进行比较

library(rmda) formulClinic<-as.formula(Label~Age) formulaOmics <- as.formula(label ~wavelet.LHL_glcm_Imc2+wavelet.LHL_glszm_ZoneEntropy+wavelet.HLL_glszm_ZoneEntropy+wavelet.LLL_firstorder_90Percentile+wavelet.LLL_firstorder_RootMeanSquared) formulaAnd <- as.formula(label~Age+RS) model_Clinic <- decision_curve(formulClinic, data=trainClinic, family=binomial(link='logit'), thresholds=seq(0,1,by=0.01), confidence.intervals=0.95, study.design = 'case-control', population.prevalence=0.3) model_Omics <- decision_curve(formulaOmics, data=tData_train, family=binomial(link='logit'), thresholds=seq(0,1,by=0.01), confidence.intervals=0.95, study.design = 'case-control', population.prevalence=0.3) model_And <- decision_curve(formulaAnd, data=tData_train2, family=binomial(link='logit'), thresholds=seq(0,1,by=0.01), confidence.intervals=0.95, study.design = 'case-control', population.prevalence=0.3) model_all <- list(model_Clinic,model_Omics,model_And) plot_decision_curve(model_all,curve.names=c('Clinical model','Radiomics signature','Clinical and Radiomics'), xlim=c(0,0.8), cost.benefit.axis=F,col=c('green','red','blue'), confidence.intervals = F, standardize = F, ain2, family=binomial(link='logit'), thresholds=seq(0,1,by=0.01), confidence.intervals=0.95, study.design = 'case-control', population.prevalence=0.3) model_all <- list(model_Clinic,model_Omics,model_And) plot_decision_curve(model_all,curve.names=c('Clinical model','Radiomics signature','Clinical and Radiomics'), xlim=c(0,0.8), cost.benefit.axis=F,col=c('green','red','blue'), confidence.intervals = F, standardize = F, legend.position="bottomleft")

一般曲线越靠左上角的，模型的性能就比较好，可以看到蓝色曲线（临床结合影像组学)的模型是在最外面的，可见这个模型性能比较好。

4.11 总结

到这一步基本结束这个项目，基本的方法和图表都已经绘制出来，代码都已经全部给出了。如果想在本地运行整个项目可以到这个项目地址（aistudio.baidu/aistudio/projectdetail/3270880）中下载【glioma.7z】压缩包，解压后用Rstudio软件打开就可以运行（前提已经配置好环境)

本文标签： 影像流程组学全