【论文复现】轻松实现影像组学全流程

  • A+
所属分类:医信大学

来源:飞浆

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官网查看地址

In [ ]

"""
经过一两个小时提取后会生成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("完成")

In [ ]


"""
再对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。同样把临床特征也进行同样的划分(顺序和影响组学特征是一致的)。

In [ ]


#导入常用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"
"""

In [ ]

"""
根据上面对影像组学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处理临床数据

In [ ]

#读取临床数据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 使用逻辑回归对临床特征进行建模

In [ ]

#建模,使用逻辑回归
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查看模型在训练集的验证情况

In [ ]

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 查看模型在测试集的验证情况

In [ ]

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检验筛选一部分特征

In [ ]

"""
用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回归进一步筛选特征

In [ ]

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变化,特征数量的变化 

In [ ]

#采用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"

In [ ]

#如果取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个影像组学特征建立逻辑回归

In [ ]

#建立公式
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 查看在训练集的情况

In [ ]

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 查看模型在测试集的结果

In [ ]

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处理数据

In [ ]

#读取经过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模型计算得分

In [ ]

#分别计算影像组学得分(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建立逻辑回归模型

In [ ]

#因为后续要绘制诺模图,所以用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查看模型在临床结合影像组学的训练集的情况

In [ ]

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 查看模型在临床结合影像组学的测试集的情况

In [ ]

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曲线

In [ ]

#计算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曲线

In [ ]

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绘制诺模图

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

In [ ]

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绘制校准曲线

In [ ]

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 绘制决策曲线

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

In [ ]

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,
                    legend.position="bottomleft")

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

4.11 总结

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

  • HealthIT.CN小程序
  • 产学研医创新平台
  • weinxin
  • 健康数据
  • HealthIT公众号
  • weinxin
admin

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: