入门篇：医学影像组学与机器学习(进入飞桨Aistudio平台，不用配置环境，开箱即用)_影像组学和机器学习的区别

背景

影像组学的概念最早由荷兰学者在2012年提出，其强调的深层次含义是指从影像（CT、MRI、PET等）中高通量地提取大量影像信息，实现肿瘤分割、特征提取与模型建立，凭借对海量影像数据信息进行更深层次的挖掘、预测和分析来辅助医师做出最准确的诊断。

影像组学大致流程是：影像数据的收集->兴趣区ROI的勾画->兴趣区域的分割->提取特征、量化特征->建立模型->分类/预测。

发现特征处理、建立模型、分类/预测这些步骤方法涉及机器学习。

所以一边学习一边做记录，过一下影像组学的流程。所以就有这个项目。可以到哔哔哩哩看李老师的影像组学视频。看看能不能在编程方面减轻医生的研究工作。因为也是接触没有多久，假如有错误，请大佬指正。

#首次要去掉注释运行一下
#解压数据
#只用到训练集
!cd /home/aistudio/
# !unzip -o data/data67772/MICCAI_BraTS_2019_Data_Training.zip -d work/
1
2
3
4
5

#安装医学影像特征提取库pyradiomics，和医学图像处理库SimpleITK
!pip install pyradiomics
!pip install SimpleITK
1
2
3

#导入常用库
import sys 
import pandas as pd
import os
import random
import shutil
import sklearn 
import scipy
import numpy as np
import radiomics  #这个库专门用来提取特征
from  radiomics import featureextractor
import SimpleITK as sitk  #读取nii文件
import matplotlib.pyplot as plt
from sklearn.linear_model import LassoCV#导入Lasso工具包LassoCV
from sklearn.preprocessing import StandardScaler#标准化工具包StandardScaler
%matplotlib inline
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

数据介绍

BraTS 数据集是胶质瘤脑肿瘤分割比赛数据集，HGG是高级别胶质瘤，LGG是低级别胶质瘤。每个病例都有四个MR模态。分别是T1、T2、Flair、T1增强。

1.因为这个braTs数据集是用作分割任务的。现在我用来做影像组学。我只用到里面的T1增强和Seg两个数据。这个项目目的是分类HGG和LGG

2.Seg.nii数据里面分了标签0、1、2、4。0是背景，其他是胶质瘤病区。我要把1、2、4标签都变成1，重新制作Mask文件。

3.因为这个数据集本身是分割数据集，所有分割的ROI不用自己画，如果是自己的医疗数据集，可以通过ITK-SNAP、LabelX这些软件读取医疗数据自己画ROI，然后再做影像组学

#重新制作Mask标签
def createNewMask(file_path,dst_path):
    file_base_name = os.path.basename(file_path)
    sitkImage = sitk.ReadImage(file_path)          #读取nii.gz文件
    npImage = sitk.GetArrayFromImage(sitkImage)    #simpleITK 转换成numpy
    npImage[npImage > 0 ] =1                       #把大于0的标签都变成1，就是所有病区都要
    outImage = sitk.GetImageFromArray(npImage)     #numpy 转换成simpleITK
    outImage.SetSpacing(sitkImage.GetSpacing())    #设置和原来nii.gz文件一样的像素空间
    outImage.SetOrigin(sitkImage.GetOrigin())      #设置和原来nii.gz文件一样的原点位置
    sitk.WriteImage(outImage,os.path.join(dst_path,file_base_name))#保存文件

if not os.path.exists('/home/aistudio/work/MyData'):
    os.makedirs('/home/aistudio/work/MyData')

source_data_path = '/home/aistudio/work/MICCAI_BraTS_2019_Data_Training'
dst_data_path = '/home/aistudio/work/MyData'
kinds = ['HGG','LGG']
index = 0
for kind in kinds:
    kind_path = os.path.join(source_data_path, kind)    #/home/aistuido/work/MICCAI_BraTS_2019_Data_Training/HGG
    for folder in os.listdir(kind_path):
        file_dir_path =  os.path.join(kind_path,folder )#/home/aistuido/work/MICCAI_BraTS_2019_Data_Training/HGG/BraTS19_CBICA_APY_1
        dst_path =  os.path.join(dst_data_path, kind,str(index)) #/home/aistuido/work/MyData/HGG/1
        for file_name in os.listdir(file_dir_path):
            if not os.path.exists(dst_path):
                os.makedirs(dst_path)
            if 't1ce' in file_name:
                dst_file_path = os.path.join(dst_path, file_name)#/home/aistuido/work/MyData/HGG/1/BraTS19_2013_5_1_t1ce.nii.gz
                shutil.copy(os.path.join(file_dir_path,file_name),dst_file_path)
                index += 1
            elif 'seg' in file_name:
                createNewMask(os.path.join(file_dir_path,file_name), dst_path)
print("完成")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

完成
1

1.从医疗影像学数据提取特征

主要使用pyradiomics 这个库自动提取，当然也可以人工实际情况、实际病种提取特有特征。

pyradiomics 提取的特征具体是什么，可以到它的官网查看了解

#这一步比较耗时，如果不想运行也用已经生成的csv文件
#HGG.csv 和LGG.csv文件

kinds = ['HGG','LGG']
para_path = '/home/aistudio/MR_1mm.yaml'#这个是特征处理配置文件，具体可以参考pyradiomics官网网站
extractor = featureextractor.RadiomicsFeatureExtractor(para_path) 
for kind in kinds:
    print("{}:开始提取特征".format(kind))
    features_dict = dict()
    df = pd.DataFrame()
    path = '/home/aistudio/work/MyData/' + 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)
            else:
                lab_path = os.path.join(path, folder, f)
        features = extractor.execute(ori_path,lab_path)  #抽取特征
        for key, value in features.items():  #输出特征
            if 'diagnostics_Versions' in  key or 'diagnostics_Configuration' in key:#这些都是一些共有的特征，可以去掉
                continue
            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'.format(kind),index=0)
print("完成")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

2.处理提取出来的特征数据

#读取HGG.csv和LGG.csv ，因为是分类任务，新增标签0是LGG，1标签是HGG
#还有各从HGG.csv 和LGG.csv随机抽取10%作为测试集

def split_df(df, ratio):
    #用来分割数据集，保留一定比例的数据集当做最终的测试集
    cut_idx = int(round(0.1 * df.shape[0]))
    data_test, data_train = df.iloc[:cut_idx], df.iloc[cut_idx:]
    return (data_train, data_test)

test_ratio = 0.1
random_state = 2021 #固定随机种子
hgg_data = pd.read_csv('/home/aistudio/HGG.csv')
lgg_data = pd.read_csv('/home/aistudio/LGG.csv')

hgg_data.insert(0,'label', 1) #插入标签
lgg_data.insert(0,'label', 0) #插入标签

hgg_data = hgg_data.sample(frac=1.0, random_state=random_state)  # 全部打乱
lgg_data = lgg_data.sample(frac=1.0, random_state=random_state)  # 全部打乱

#因为有些特征是字符串，直接删掉
cols=[x for i,x in enumerate(hgg_data.columns) if type(hgg_data.iat[1,i]) == str]
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]
lgg_data=lgg_data.drop(cols,axis=1)

hgg_data_train, hgg_data_test = split_df(hgg_data,test_ratio) #返回train 和test数据集
lgg_data_train, lgg_data_test = split_df(lgg_data,test_ratio) #返回train 和test数据集

#保存测试集为cvs 后面最终验证使用
hgg_data_test.to_csv('HGG_test.csv',index=False)
lgg_data_test.to_csv('LGG_test.csv',index=False)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

#把hgg_data_train 和lgg_data_train 并在一起并且打乱。
data = pd.concat([hgg_data_train, lgg_data_train])
data = data.sample(frac=1.0,random_state=random_state)  # 全部打乱
print("一共有{}行特征数据".format(len(data)))
print("一共有{}列不同特征".format(data.shape[1]))
#再把特征值数据和标签数据分开
x = data[data.columns[1:]]
y = data['label']
#取X的5行看看数据
# x.head()

1
2
3
4
5
6
7
8
9
10
11

一共有301行特征数据
一共有106列不同特征
1
2

3.特征筛选方法

T 检验 特征筛选

使用其中的两独立样本t检验，检验两组独立样本的平均数和其分布是否具有显著性差异（显著性差异p<0.05）

假如两组样本中，里面有涉及平均数的特征，可以使用T检验，判断这些平均数的特征是否有显著性差异，有显著性差异就可以留着做分类任务。没有显著性差异就可以去掉这项平均数特征

#通过T检验从106个特征筛选出85个特征
from scipy.stats import levene, ttest_ind
counts = 0
columns_index =[]
for column_name in hgg_data_train.columns[1:]:
    if levene(hgg_data_train[column_name], lgg_data_train[column_name])[1] > 0.05:
        if ttest_ind(hgg_data_train[column_name],lgg_data_train[column_name],equal_var=True)[1] < 0.05:
            columns_index.append(column_name)
    else:
        if ttest_ind(hgg_data_train[column_name],lgg_data_train[column_name],equal_var=False)[1] < 0.05:
            columns_index.append(column_name)

print("筛选后剩下的特征数：{}个".format(len(columns_index)))
1
2
3
4
5
6
7
8
9
10
11
12
13

筛选后剩下的特征数：85个


/opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages/scipy/stats/morestats.py:2352: RuntimeWarning: invalid value encountered in double_scalars
  W = numer / denom
1
2
3
4
5

Lasso 特征筛选

Lasso回归是最小二乘法回归的基础上加上L1正则表达式，L1正则表达式可以防止模型过拟合。

这里的用途是把不重要的特征变为0，进行特征筛选(降维，从一堆特征中，挑出一些有用的特征做分类任务)

#数据只保留从T检验筛选出的特征数据，重新组合成data
if  not 'label' in columns_index:
    columns_index = ['label'] + columns_index
hgg_train = hgg_data_train[columns_index]  
lgg_train = lgg_data_train[columns_index]  

data = pd.concat([hgg_train, lgg_train])
data = data.sample(frac=1.0,random_state=random_state)  # 全部打乱

#再把特征值数据和标签数据分开
x = data[data.columns[1:]]
y = data['label']
#先保存X的列名
columnNames = x.columns
x = x.astype(np.float32)#把x数据转换成np.float格式
x = StandardScaler().fit_transform(x)#对x进行均值-标准差归一化
x = pd.DataFrame(x,columns=columnNames)#转 DataFrame 格式
# 自己建立Lasso进行alpha选择的范围
# 形成10为底的指数范围 10**（-10）~ 10**（1）,200个数值
alpha_range = np.linspace(-1,1,10)
#alpha_range在这个参数范围里挑出aplpha进行训练，cv是把数据集分5分，进行交叉验证，max_iter是训练10000轮
lassoCV_model = LassoCV(alphas=alpha_range,cv=5,max_iter=10000)
#进行训练
lassoCV_model.fit(x,y)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

LassoCV(alphas=array([-1.        , -0.77777778, -0.55555556, -0.33333333, -0.11111111,
        0.11111111,  0.33333333,  0.55555556,  0.77777778,  1.        ]),
        copy_X=True, cv=5, eps=0.001, fit_intercept=True, max_iter=10000,
        n_alphas=100, n_jobs=None, normalize=False, positive=False,
        precompute='auto', random_state=None, selection='cyclic', tol=0.0001,
        verbose=False)
1
2
3
4
5
6

#打印训练找出来的入值
print(lassoCV_model.alpha_)
# print("Coefficient of the model:{}".format(lassoCV_model.coef_) )
# print("intercept of the model:{}".format(lassoCV_model.intercept_))
coef = pd.Series(lassoCV_model.coef_, index=columnNames)
print("从原来{}个特征，筛选剩下{}个".format(len(x),sum(coef !=0)))
print("分别是以下特征")
print(coef[coef !=0])
index = coef[coef !=0].index
x = x[index]
# x.head()
1
2
3
4
5
6
7
8
9
10
11

0.11111111111111116
从原来301个特征，筛选剩下3个
分别是以下特征
original_firstorder_RootMeanSquared    0.043582
original_firstorder_Skewness           0.046152
original_glszm_SmallAreaEmphasis       0.094545
dtype: float32
1
2
3
4
5
6
7

#绘制特征相关系数热力图
import seaborn as sns
f, ax= plt.subplots(figsize = (6, 5))
sns.heatmap(x.corr(),annot=True,annot_kws={'size':10,'weight':'bold', },ax=ax)#绘制混淆矩阵
ax.set_xticklabels(ax.get_xticklabels(), rotation=45,va="top",ha="right")
ax.set_yticklabels(ax.get_yticklabels(), rotation=45)
plt.show()
1
2
3
4
5
6
7

#画一个特征系数的柱状图
weight = coef[coef !=0].to_dict()
#根据值大小排列一下
weight = dict(sorted(weight.items(),key=lambda x:x[1],reverse=False))
plt.figure(figsize=(8,6))#设置画布的尺寸
plt.title('characters classification weight',fontsize=15)#标题，并设定字号大小
plt.xlabel(u'weighted value',fontsize=14)#设置x轴，并设定字号大小
plt.ylabel(u'feature')
plt.barh(range(len(weight.values())), list(weight.values()),tick_label = list(weight.keys()),alpha=0.6, facecolor = 'blue', edgecolor = 'black', label='feature weight')
plt.legend(loc=4)#图例展示位置，数字代表第几象限
plt.show()
1
2
3
4
5
6
7
8
9
10
11

#绘制误差棒图
MSEs = lassoCV_model.mse_path_
mse = list()
std = list()
for m in MSEs:
    mse.append(np.mean(m))
    std.append(np.std(m))

plt.figure(figsize=(8,6))
plt.errorbar(lassoCV_model.alphas_, mse, std,fmt='o:',ecolor='lightblue',
			elinewidth=3,ms=5,mfc='wheat',mec='salmon',capsize=3)
plt.axvline(lassoCV_model.alpha_, color='red', ls='--')
plt.title('Errorbar')
plt.xlabel('Lambda')
plt.ylabel('MSE')
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

#这个图显示随着lambda 的变化，系数的变化走势
x = data[data.columns[1:]]
y = data['label']
#先保存X的列名
columnNames = x.columns
x = x.astype(np.float32)#把x数据转换成np.float格式
x = StandardScaler().fit_transform(x)#对x进行均值-标准差归一化
x = pd.DataFrame(x,columns=columnNames)#转 DataFrame 格式
coefs = lassoCV_model.path(x,y, alphas=alpha_range, cv=5,max_iter=10000)[1].T
plt.plot(lassoCV_model.alphas,coefs,'-')
plt.axvline(lassoCV_model.alpha_, color='red', ls='--')
plt.xlabel('Lambda')
plt.ylabel('coef')
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14

4.使用随机森林进行分类

from sklearn.model_selection import train_test_split #分割训练集和验证集
from sklearn.ensemble import RandomForestClassifier #导入随机森林分类器
from sklearn.externals import joblib #用来保存 sklearn 训练好的模型
#把数据分成训练集和验证集，7：3比例
index = coef[coef !=0].index
x = x[index]
x_train,x_test, y_train, y_test = train_test_split(x,y,test_size=0.3)
model_forest = RandomForestClassifier(n_estimators=30,random_state=random_state).fit(x_train,y_train)
score = model_forest.score(x_test, y_test)
print("在验证集上的准确率：{}".format(score))
#把随机森林的模型保存下来
joblib.dump(model_forest, 'model_forest1.model')
1
2
3
4
5
6
7
8
9
10
11
12

在验证集上的准确率：0.9120879120879121





['model_forest1.model']
1
2
3
4
5
6
7

5.在测试集上进行验证

from sklearn.externals import joblib #用来保存 sklearn 训练好的模型
hgg_test = pd.read_csv('/home/aistudio/HGG_test.csv')
lgg_test = pd.read_csv('/home/aistudio/LGG_test.csv')

#再把特征值数据和标签数据分开
x_test_data = data_test[data_test.columns[1:]]
columnNames = x_test_data.columns
x_test_data = x_test_data.astype(np.float32)
x_test_data = StandardScaler().fit_transform(x_test_data)
x_test_data = pd.DataFrame(x_test_data,columns=columnNames)
y_test_data = data_test['label']
#只提取之前Lasso 筛选后的特征
index = coef[coef !=0].index
x_test_data = x_test_data[index]
print("测试集一共有{}行特征数据，{}列不同特征,包含HGG:{}例，LGG:{}例".format(len(x_test_data),x_test_data.shape[1],len(hgg_data_test),len(lgg_data_test)))
#加载保存后的模型，然后进行预测
model_forest = joblib.load('/home/aistudio/my_model_forest.model') #这是自己训练模型，记得替换自己的。
score = model_forest.score(x_test_data, y_test_data)
print("在测试集上的准确率：{}".format(score))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

测试集一共有34行特征数据，3列不同特征,包含HGG:26例，LGG:8例
在测试集上的准确率：0.9117647058823529
1
2

6.绘制混淆矩阵和ROC图表

有些概念比较混淆，例如准确率、精准率等。

准确率(accurary)：（正确预测0的个数+正确预测1的个数）/ 所有样本数

精确率(precision)=查准率：预测为1中真实为1的样本数 / 预测为1的样本数

真阳性率(True Positive Rate, TPR)=灵敏度(sensitivity)=敏感度=召回率(recall)=标签1的查全率：正确预测1的个数 / 实际1的样本数

假阳性率(False Positive Rate, FPR)=1-真阴率 ：预测为0中真实为1的样本数 / 实际0的样本数

真阴性率(True Negative Rate,TNR)=特异度(specificity)=标签0的查全率：正确预测0的个数 / 实际0的样本数

#绘制混淆矩阵
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
import seaborn as sns

predict_label = model_forest.predict(x_test_data) #预测的标签
label = y_test_data.to_list()  #真实标签
confusion = confusion_matrix(label, predict_label)#计算混淆矩阵

plt.figure(figsize=(6,5))
sns.heatmap(confusion,cmap='Blues_r',annot=True,annot_kws={'size':20,'weight':'bold', })#绘制混淆矩阵
plt.xlabel('Predict')
plt.ylabel('True')
plt.show()
print("混淆矩阵为：\n{}".format(confusion))
#计算灵敏度(相对标签1，HGG的召回率)
#计算特异度(相对标签0，LGG的召回率)
print("\n计算各项指标：")
print(classification_report(label, predict_label))
print("根据上面显示的报告，得知：")
print("准确率=0.91，灵敏度=1.00，特异度=0.62")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

混淆矩阵为：
[[ 5  3]
 [ 0 26]]

计算各项指标：
              precision    recall  f1-score   support

           0       1.00      0.62      0.77         8
           1       0.90      1.00      0.95        26

    accuracy                           0.91        34
   macro avg       0.95      0.81      0.86        34
weighted avg       0.92      0.91      0.90        34

根据上面显示的报告，得知：
准确率=0.91，灵敏度=1.00，特异度=0.62
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

#绘制ROC曲线
from sklearn.metrics import roc_curve, roc_auc_score
kind = {'HGG':1,"LGG":0}
model_forest = joblib.load('/home/aistudio/my_model_forest.model')#这是自己训练模型，记得替换自己的
label = y_test_data.to_list()  #真实标签
y_predict = model_forest.predict_proba(x_test_data)#得到标签0和1对应的概率
fpr , tpr ,threshold = roc_curve(label, y_predict[:,kind['LGG']], pos_label=kind['LGG'])
fpr1 , tpr1 ,threshold = roc_curve(label, y_predict[:,kind['HGG']], pos_label=kind['HGG'])
plt.plot(fpr, tpr,marker='o', markersize=5,label='LGG')
plt.plot(fpr1, tpr1,marker='o', markersize=5,label='HGG')
plt.title("HGG AUC:{:.2f}".format(roc_auc_score(label,  y_predict[:,1]),fontsize=15))
plt.xlabel('FPR')
plt.ylabel('TPR')
plt.legend(loc=4)
plt.show()
GG'])
plt.plot(fpr, tpr,marker='o', markersize=5,label='LGG')
plt.plot(fpr1, tpr1,marker='o', markersize=5,label='HGG')
plt.title("HGG AUC:{:.2f}".format(roc_auc_score(label,  y_predict[:,1]),fontsize=15))
plt.xlabel('FPR')
plt.ylabel('TPR')
plt.legend(loc=4)
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

7.小结

基于MRI增强T1WI的影像组学特征值，探讨机器学习模型 随机森林(random forest，RF)对胶质瘤分化程度鉴别价值。

【材料与方法】
胶质瘤MRI增强T1WI影像资料335例(高级胶质瘤(HGG)259例，低级胶质瘤(LGG)76例)，以6:3:1的比例分为训练集和验证集和测试集对增强T1WI数据集进行影像组学特征值提取，并对随机森林模型进行训练和测试，采用 ROC曲线评价预测效能。

【结果】
基于增强T1WI图像建立的RF模型在测试集的AUC为0.96，准确度、敏感度、特异度分别为91%、 100%、 62%， 模型中重要性排名前3位的特征分别为平方强度值平均平方根(RMS)、偏度(Skewness)、灰度大小区域矩阵(GLSZM,SAE)。

8.总结

整个项目的目的就是跑一下整个影像组学的流程。也给医生提供了一些代码模板，从怎样获取特征到特征筛选、模型训练、绘制图等等。
运行代码请点击:
项目地址（开箱即用https://aistudio.baidu.com/aistudio/projectdetail/1690500）
进入地址后点击fork

【个人介绍】 广州某医院的放射科的一名放射技师。

只是一位编程爱好者

只是想把自己的爱好融入工作中

只是想让自己通过努力获取成就和快乐

欢迎更多志同道合的朋友一起玩耍~~~

我在AI Studio上获得钻石等级，点亮10个徽章，来互关呀~ https://aistudio.baidu.com/aistudio/personalcenter/thirdview/181096