CN107480675A - 一种基于放射组学的肝静脉压力梯度计算模型的构建方法 - Google Patents

Abstract

本发明提供了一种基于放射组学的肝静脉压力梯度(radiomics‑based hepatic venous pressure gradient,rHVPG)计算模型的构建方法，以此为基础可以构建一种更具优势的肝静脉压力梯度计算模型，为计算门脉高压患者的早期无创指标提供一种新的途径。本发明还提供了一种基于放射组学的肝静脉压力梯度的计算系统。本发明能够克服现有无创测量方法的诸多局限。本发明提出的rHVPG计算系统切实可行，有望为门脉压力的无创测量提供新的思路，为改善门脉高压患者的生活质量、减轻家庭和社会的疾病负担起到积极推动作用。

Description

一种基于放射组学的肝静脉压力梯度计算模型的构建方法

技术领域

本发明涉及一种基于放射组学的肝静脉压力梯度计算模型的构建方法。

背景技术

在全球范围内，肝脏疾病(病毒性感染、酒精性肝病、非酒精性脂肪肝及相关的肝硬化、肝细胞肝癌)是导致死亡的主要原因之一，仅在中国肝脏疾病就威胁了近3亿人的健康和生命，严重加剧了全球疾病的负担。以乙型肝炎病毒感染引起的慢性肝病为例，据世界卫生组织报道，全球约20亿人曾感染乙肝病毒，其中2.4亿人为慢性乙肝病毒感染者，每年约有65万人死于乙肝病毒感染所致的肝硬化、肝细胞肝癌等。在全球肝硬化患者中，由乙肝病毒感染引起的为30％。而在我国，由于感染乙肝病毒而导致的肝硬化比例高达60％。因此，肝硬化的高发病率和较差的临床结局已造成了国内外严重的社会公共健康问题。

门脉高压是肝硬化的重要表现之一，临床上门脉高压的发现和诊断常由于出现晚期严重的并发症，如：急性静脉曲张破裂出血、顽固性腹水、肝性脑病、门脉高压性胃肠病、肝肾综合症、肝肺综合症和继发感染等。因此，门脉高压及其并发症严重影响了肝硬化患者的生活质量和长期预后，对于门脉高压(尤其是早期无症状性门脉高压)的诊断和监测是肝硬化等终末期肝病治疗链中最重要的环节之一。

在最新的2015版Baveno VI共识(门脉高压的风险分层和个体化管理)中，有创性肝静脉压力梯度(hepatic venous pressure gradient,HVPG)的测量继续被推荐为诊断临床显著性门脉高压的金标准，即HVPG≥10mmHg可确诊(de Franchis R#*,Baveno VIFaculty.Expanding consensus in portal hypertension:Report of the Baveno VIConsensus Workshop:Stratifying risk and individualizing care for portalhypertension.J Hepatol.2015；63(3):743-752.)。该侵袭性方法是通过颈内静脉穿刺置管，依次经过颈内静脉、上腔静脉、右心房、下腔静脉进入肝静脉，分别测量肝静脉自由压(free hepatic venous pressure,FHVP)和肝静脉楔入压(wedged hepatic venouspressure,WHVP)，计算两者差值得到HVPG[Bloom S#,Kemp W,Lubel J*.Portalhypertension:pathophysiology,diagnosis and management.Intern Med J.2015；45(1):16-26.]。此外，HVPG对于治疗应答、并发症风险和长期预后也可提供重要信息，其数值的变化被推荐为研究预后的替代指标。此外，HVPG测量还被鼓励用于研究门脉高压创新诊断和治疗方法的临床试验中。然而，测量HVPG也存在一些问题，首先其检测手段为有创性，在门脉高压早期尚无严重并发症时很难被患者所接受；其次，测量过程需由经过专业化培训的介入医生在肝静脉造影的辅助下进行，这不但增加了受试者的射线暴露和造影剂过敏可能，还存在有一定的技术操作风险；加之其高昂的诊断费用，HVPG测定在我国受到很大限制，目前仅有部分三级甲等医院开展。除此之外，门脉高压的有创性诊断技术还包括：超声引导下门脉穿刺测压和开腹手术中直接门脉测压[de Franchis R#*,Dell'EraA.Invasive and noninvasive methods to diagnose portal hypertension andesophageal varices.Clin Liver Dis.2014；18(2):293-302.]。然而，直接测量门脉压力的风险更大、对操作者的技术水平要求也更高，在临床实践中很难被医生和患者所采纳，目前仍多用于动物实验性研究。

门脉压力的无创性评估是目前该领域研究的重要内容，主要包括以下三个方面：(1)检测肝内阻力：肝内阻力增加是导致门脉高压形成的重要因素。内源性缩血管物质的释放和肝组织结构的紊乱可促使肝内阻力的显著升高。有研究提示，血清内皮素-1表达水平与HVPG存在相关性，但该结论尚缺乏大样本量临床试验的支持。另一方面，肝内机械结构的改变也会导致肝内阻力增加，临床上多采用肝瞬时弹性成像(FibroScan)技术，通过对肝脏硬度值的测定来评估肝硬化程度，该方法操作简便、重复性好。但也有报道指出，FibroScan的诊断准确性受到肥胖、肋间隙狭窄、炎症等因素干扰，对门脉高压并发症(如高风险食管静脉曲张)的预测价值不高，需联合其他无创指标进行综合判断。(2)测定门脉循环血量：门脉高压的直接因素为循环血量的积聚升高超过了机体的代偿，当肝硬化伴发重度门脉高压时，心输出量也会明显升高且与HVPG存在高度相关性。有研究利用彩色多普勒超声(colorDoppler ultrasound,CDUS)对肝硬化患者的门脉直径和血流速度进行测定，在此基础上计算门脉系统的循环血量。然而，CDUS测量技术的主观性较强且易受到肥胖等因素干扰，尚无准确的临床证据表明经血流速度和门脉直径计算的门脉压力与HVPG具有较高的诊断一致性。此外，CT血管成像(CT angiography,CTA)和磁共振静脉显像也可以评估门脉系统的循环血量。然而，其测定原理仍是基于门脉直径和侧枝循环来间接评价门脉压力。(3)仿真计算门脉压力：虚拟仿真技术是计算机图形学、医学图像处理、软件工程计算等多学科交叉融合的结晶。以计算机仿真辅助医学影像获得的虚拟人体数据为疾病的介入和治疗提供了全新的思路。目前已有多项国际前瞻性临床试验表明，基于CTA和流体力学计算获得的血流储备分数对于冠状动脉的功能性狭窄具有高度应用价值，相关仿真计算软件也已获得美国食品药品监督管理局的上市批准。然而，国内外该领域研究的主要问题还表现在：评估范围局限于门脉主干，未考虑肝静脉-门静脉系统的整体影响；仅考虑到门脉血液动力学因素，忽略了肝硬化患者的血液流变学改变，如血液黏度的变化和红细胞变形能力的降低都会促使门脉阻力的升高。

门脉压力的各种测定技术目前仍受到各自干扰因素的影响，或者有创风险大、操作难度高，或者干扰因素多、数值变异大。目前国内外相关报道中，门脉高压的诊断方法，特别是早期无创性诊断技术，仍有待进一步研究完善。

发明内容

本发明的目的在于提供一种基于放射组学的肝静脉压力梯度(radiomics-basedhepatic venous pressure gradient,rHVPG)计算模型的构建方法，以此为基础可以构建一种更具优势的肝静脉压力梯度计算模型，为计算门脉高压患者的早期无创指标提供一种新的途径。本发明还提供了一种基于放射组学的肝静脉压力梯度的计算系统。

为实现上述目的，本发明所采取的技术方案：一种基于放射组学的肝静脉压力梯度计算模型的构建方法，包括以下步骤：

A.从标本的肘正中静脉注射造影剂，进行CT血管成像(CTA)，获取包括肝静脉期在内的CTA图层序列，导出图层序列，格式为dcm，层厚1.25mm，图像分辨率512×512像素；

B.将所获取的CTA图层序列导入医学影像控制软件ITK-SNAP，选择门静脉期图层序列，确认图像信息，ITK-SNAP软件自动识别图像序列，生成门静脉期CTA图像序列的冠状位、矢状位和水平位图像；

C.通过Expand this view to occupy the entire window功能将水平位图像放大，使用Zoom Inspector工具，选择zoom to fit调整图像至合适的大小；

D.以寻找水平位图像中的第一肝门作为目标，利用ITK-SNAP软件的PainburshInspector工具以包含目标的CT值、排除目标周围软组织的CT值为原则，提取目标；再反复利用Painbursh Inspector工具，剔除剩余的非目标结构；

E.选择Save segmentation Image将分割出的目标以NiFTI格式(.nii)导出，形成掩模图像的NiFTI格式；

F.打开dcm2niigui软件，选择Compressed FSL[4D NIfTInii]选项，将筛选出的含第一肝门层面的dcm文件序列放到dcm2niigui工作空间中，利用dcm2niigui软件将dcm文件序列转换为NiFTI格式；然后将含第一肝门层面的dcm文件序列、新生成的NiFTI格式文件和步骤D中生成的NiFTI格式文件放到同一文件夹中；

G.打开MATLAB R2016b软件，进入工作界面；在菜单栏处选择SetPath按钮，将所需要的feature_extra.m、readDICOMdir.m、str2matrix.m、NonTextureFeatures、TextureToolbox、NIfTI、Utilities、MultivariableModeling、STUDIES工具箱所在的绝对路径加载到MATLAB搜索路径中；

H.打开feature_extra.m文件，利用NIfTI工具箱中load_nii函数将步骤F中两个NiFTI文件导入MATLAB中，并且以名为data和data1的细胞数组储存；然后将两个细胞数组中的图像提取出来存储到名为mask和image变量中；

I.利用Utilities工具箱中computeBoundingBox函数提取出图像中感兴趣区域(ROI)的边界坐标提取出来储存到变量maskBox中；然后通过矩阵的基本操作将mask和image变量的边界坐标以内的图像信息提取出来放到变量maskBox和image_originalBox中；完成上述操作后，通过矩阵maskBox和image_originalBox的点乘运算得到图像中感兴趣区域，然后保存在ROI变量中，然后通过imshow与contour两个函数将感兴趣区域展示出来；

J.利用readDICOMdir.m函数将含第一肝门层面的dcm文件序列读入到MATLAB中，通过变量索引，将像素大小、层厚提取出来并保存到变量pixelW与sliceS中；

K.利用NonTextureFeatures工具箱中getEccentricity.m、getSize.m、getSolidity.m和getVolume.m函数，将变量ROI、pixelW、sliceS传输到相应函数中；经计算得到eccentricity、sizeROI、solidity、volumeROI 4个非纹理特征；

L.为了进行纹理特征提取，首先对参数进行设置；在TextureFeatures工具箱中一共有四个可变参数，分别为R、Scale、Algo和Ng，每个可变参数分别对应不同的变化空间，其中R对应为[1/2,2/3,1,3/2,2],Scale对应为[‘pixelW’,1,2,3,4,5],Algo对应为[‘Equal’,‘Lloyd’],Ng对应为[8,16,32,64]，而参数的排列组合使得纹理特征有显著变化；利用for循环将不同参数的取值进行排列组合，得到5*6*2*4组设定；通过参数传递向TextureFeatures工具箱传递设定，而每一组设定将得到43种纹理特征，分别存放于名为glcmTextures、glrlmTextures、glszmTextures、ngtdmTextures和globalTextures的结构体中，最后总共5*6*2*4*43＝10320个纹理特征；

M.对于不同病人，在将第一肝门层面感兴趣区域的非纹理特征与纹理特征提取出来后，一共得到了(4+5*6*2*4*43)＝10324个特征，然后将非纹理特征和纹理特征整合成大小为病人个数*特征个数的特征矩阵；

N.重复步骤A到步骤M完成对脾门层面感兴趣区域的非纹理特征与纹理特征提取，最后整合成大小为病人个数*特征个数的特征矩阵；

O.利用str2matrix.m函数，将步骤M与步骤N中的两个特征矩阵横向拼接，组成大小为病人个数*(特征个数*2)的特征矩阵；以各个病人有无显著性门脉高压为标准制作名为outcome的mat文件，其中有显著性门脉高压用1表示，反之用0表示；接着将特征矩阵和outcome的mat文件横向拼接起来，组成一个大小为病人个数*(特征个数*2+1)的矩阵，并且用csvwrite.m函数将已整合的第一肝门层面、脾门层面的图像特征矩阵保存为csv格式；

P.打开R软件，利用read.table函数将cvs格式的特征矩阵导入到R软件的工作空间中，并保存在变量train中；然后利用as.matrix函数将训练集中病人的图像特征矩阵提取出来，保存在train_x中作为建模的输入变量；同样地，利用as.matrix函数将训练集中病人的outcome标签提取出来，保存在train_y中作为建模的输出标签；

Q.利用library函数加载glmnet工具包，然后使用glmnet函数中的lasso降维方法对(4+5*6*2*4*43)*2＝20648个特征进行筛选，并使用dev.new函数和plot函数将降维过程显示出来，图中每个lambda对应着一个线性模型且每个模型中使用的特征也不尽相同；

R.通过设置参数nfold、type.measure和family并利用cv.glmnet函数对训练集做nfold折交叉验证、建立线性的二项分布数学模型，然后通过dev.new函数和plot函数显示每个线性模型的误差曲线；

S.通过cv$lambda.min选择AUC最高的lambda所对应的最优线性模型，然后利用coef函数和coefficients函数输出该最优模型的非零系数和筛选出来的特征；

T.利用library函数加载boot工具包和Hmisc工具包，然后通过boot函数和boot.ci函数对训练集做1000次bootstraping并计算出训练集的C-index指标；

接着加载pROC工具包，利用roc函数计算出训练集的ROC(Receiver OperatingCharacteristic)曲线下面积(AUC)；

U.为了进行模型验证，以未进行特征分析的人群制作测试集，同样经过上诉步骤A到步骤O将测试人群的第一肝门层面和脾门层面的20648维图像特征提取出来保存为csv格式；利用read.table函数将cvs格式的特征矩阵导入到R软件的工作空间中，并保存在变量test中；然后利用as.matrix函数将测试集中病人的图像特征矩阵提取出来，保存在test_x中作为已建好的模型的输入变量；利用glmnet工具包中的predict函数来计算测试集rHVPG的数值，测试结果与临床应用中金标准HVPG相比较，并利用roc函数计算出测试集的AUC。

肝静脉压力梯度(HVPG)是当前世界公认评估临床门脉高压(CSPH)的唯一有效的“金标准”，即它是个中介----能代表是否存在门脉高压。本发明提出的模型旨在与HVPG建立一个对应关系，如果该关系对应得好便意味着本发明也能够指示是否存在门脉高压。所以本方法需要与HVPG相比较，以表明本发明的模型与HVPG有着很好的对应关系。本发明提供了一种基于放射组学的肝静脉压力梯度的计算系统，包括：

数据输入模块，用于将患者第一肝门层面图像特征和脾门层面图像特征的检测结果输入模型计算模块；

模型计算模块，包括基于放射组学的肝静脉压力梯度计算模型，用于根据患者第一肝门层面图像特征和脾门层面图像特征的检测结果以及基于放射组学的肝静脉压力梯度计算模型计算肝静脉压力梯度数值；所述基于放射组学的肝静脉压力梯度计算模型包括基于放射组学的肝静脉压力梯度计算公式，所述基于放射组学的肝静脉压力梯度计算公式：

rHVPG＝-22.76522+14.94622×Complexity_0.5R_pS_LQ_64N_L

+8934.012×ZP_0.5R_1S_LQ_32N_L+0.4419001×Skewness_0.5R_3S_LQ_64N_L

-14.63828×LZE_0.5R_4S_EQ_64N_L+687.5762×ZP_0.5R_5S_EQ_16N_L

-0.00167672×HGZE_0.5R_5S_LQ_32N_L+19.31699×GLN_1R_2S_EQ_16N_L

-0.02605657 HGRE_1R_3S_LQ_8N_L-0.08363389 HGRE_1.5R_3S_LQ_8N_L

-3.724522e-06 Busyness_2R_1S_EQ_32N_L+545.8976 ZP_2R_1SLQ_32N_L

+267.1424 GLN_2R_3S_EQ_8N_L+1361.243 ZP_2R_3S_EQ_16N_L

-1289.11 Strength_0.5R_pS_EQ_8N_S+432.8978 Strength_1.5R_4S_LQ_64N_S；

结果输出模块，用于根据基于放射组学的肝静脉压力梯度计算结果判定患者肝静脉压力；当rHVPG≥0.71时，患者存在临床显著性门脉高压；当rHVPG<0.71时,患者不存在临床显著性门脉高压。

本发明的有益效果在于：

如前所述，门脉高压作为肝硬化的重要表现，其临床明确诊断往往滞后于晚期严重的并发症发生。然而，目前国内外的门脉高压诊断技术或者有创风险大、操作难度高，或者干扰因素多、数值变异大。本发明首先提出了rHVPG的概念，并在其方法学构建和模型验证中取得重要数据，有利于探索并建立一种安全无创、准确量化的门脉压力计算新技术。本发明优化并完善腹部CT分割感兴趣区域(目标为第一肝门层面肝脏和脾门层面脾脏)、特征的提取与筛选和rHVPG计算模型的建立与应用，构建并验证一种更具优势的rHVPG计算新系统，为门脉高压患者的肝静脉压力梯度计算提供一种新的途径。

附图说明

图1为从CT血管成像(CTA)导出包括门静脉期在内的CTA图层序列。

图2为选择导入的门静脉期CTA图层序列，层厚1.25mm。

图3为确认图像信息。

图4为自动生成的冠状位、矢状位和水平位图像。

图5为Expand this view to occupy the entire window功能放大的水平位图像。

图6为Zoom Inspector工具调整大小后的图像。

图7为分割出的感兴趣区域(目标为第一肝门层面的肝脏)。

图8为保存分割图像的界面。

图9为dcm2niigiu软件的工作界面并选择Compressed FSL[4D NIfTInii]选项。

图10为将筛选出的含肝脏图像的dcm文件放到dcm2niigui工作空间。

图11为利用dcm2niigui软件将dcm文件序列转换为NiFTI格式。

图12为将新生成的NiFTI文件和步骤D中生成的NiFTI文件放到同一文件中。

图13为MATLAB R2016b的工作界面。

图14为MATLAB R2016b设置路径界面。

图15为将所需要的工具箱所在的绝对路径加载到MATLAB搜索路径中。

图16为prepare_ROI.m文件的部分代码。

图17为load_nii函数将两个NiFTI文件导入MATLAB。

图18为Utilities工具箱中computeBoundingBox函数部分代码。

图19为提取图像中感兴趣区域并存储变量的部分代码。

图20为部分图像中感兴趣区域(ROI)。

图21为NonTextureFeatures工具箱中部分函数。

图22为利用NonTextureFeatures工具箱提取出的4个非纹理特征。

图23为调用TextureFeatures工具箱中纹理特征的部分代码。

图24为利用TextureFeatures工具箱提取出的43种纹理特征。

图25为outcome的mat文件。

图26为R语言中特征数据的提取与存储。

图27为降维过程的特征筛选图。

图28为每个线性模型的误差曲线。

图29为1000次bootstraping并计算出训练集的C-index指标代码。

图30为predict函数与AUC计算代码。

图31为实施例中的降维过程的特征筛选图。

图32为实施例中的每个线性模型的误差曲线。

具体实施方式

下面结合附图对本发明进行详细说明。

本发明所述的一种基于放射组学的肝静脉压力梯度的计算方法包括：腹部CT分割感兴趣区域(目标为第一肝门层面肝脏和脾门层面脾脏)、特征的提取与筛选和rHVPG计算模型的建立与验证。

1.从标本的肘正中静脉注射造影剂，进行CT血管成像(CTA)，获取包括门静脉期在内的CTA图层序列，导出图层序列，格式为dicom，层厚1.25mm，图像分辨率512×512像素(如图1所示)；

2.将所获取的CTA图层序列导入医学影像控制软件ITK-SNAP，选择门静脉期图层序列(如图2所示)，确认图像信息(如图3所示)，ITK-SNAP软件自动识别图像序列，生成门静脉期CTA图像序列的冠状位、矢状位和水平位图像(如图4所示)；

3.通过Expand this view to occupy the entire window功能将水平位图像放大(如图5所示)，使用Zoom Inspector工具，选择zoom to fit调整图像至合适的大小(如图6所示)；

4.以寻找水平位图像中的第一肝门的肝脏(目标)，利用ITK-SNAP软件的Painbursh Inspector工具以包含目标的CT值、排除目标周围软组织的CT值为原则，提取目标；再反复利用Painbursh Inspector工具，剔除剩余的非目标结构(如图7所示)；

5.将选择Save segmentation Image将分割出的目标以NiFTI格式(.nii)导出，形成掩模图像的NiFTI格式(如图8所示)；

6.打开dcm2niigui软件，选择Compressed FSL[4D NIfTInii]选项(如图9所示)，将筛选出的含第一肝门层面的dcm文件序列放到dcm2niigui工作空间中(如图10所示)，利用dcm2niigui软件将dcm文件序列转换为NiFTI格式(如图11所示)；然后将含第一肝门层面的dcm文件序列、新生成的NiFTI格式文件和步骤D中生成的NiFTI格式文件放到同一文件夹中(如图12所示)；

7.打开MATLAB R2016b，进入工作界面(如图13所示)；在菜单栏处选择SetPath按钮，进入设置路径界面(如图14所示)，将所需要的prepare_ROI、NonTextureFeatures、TextureFeature、NIfTI、Utilities、MultivaribleModeling工具箱所在的绝对路径加载到MATLAB搜索路径中(如图15所示)；

8.打开feature_extra.m文件(如图16所示)，利用NIfTI工具箱中load_nii函数将步骤F中两个NiFTI文件导入MATLAB中(如图17所示)，并且以名为data和data1的细胞数组储存；然后将两个细胞数组中的图像提取出来存储到名为mask和image变量中；

9.利用Utilities工具箱中computeBoundingBox函数提取出图像中感兴趣区域(ROI)的边界坐标提取出来储存到变量maskBox中(如图18所示)；然后通过矩阵的基本操作将mask和image变量的边界坐标以内的图像信息提取出来放到变量maskBox和image_originalBox中(如图19所示)；完成上述操作后，通过矩阵maskBox和image_originalBox的点乘运算得到出图像中感兴趣区域，然后保存在ROI变量中，然后通过imshow与contour两个函数将感兴趣区域展示出来(如图20所示)；

10.利用readDICOMdir.m函数将含第一肝门层面的dcm文件序列读入到MATLAB中，通过变量索引，将像素大小、层厚提取出来并保存到变量pixelW与sliceS中(如图21所示)；

11.利用NonTextureFeatures工具箱中getEccentricity.m、getSize.m、getSolidity.m和getVolume.m函数，将变量ROI、pixelW、sliceS传输到相应函数中；经计算得到eccentricity、sizeROI、solidity、volumeROI 4个非纹理特征(如图22所示)；

12.为了进行纹理特征提取，首先对参数进行设置；在TextureFeatures工具箱中一共有四个可变参数，分别为R、Scale、Algo和Ng，每个可变参数分别对应不同的变化空间，其中R对应为[1/2,2/3,1,3/2,2],Scale对应为[‘pixelW’,1,2,3,4,5],Algo对应为[‘Equal’,‘Lloyd’],Ng对应为[8,16,32,64]，而参数的排列组合使得纹理特征有显著变化；利用for循环将不同参数的取值进行排列组合，得到5*6*2*4组设定；通过参数传递向TextureFeatures工具箱传递设定(如图23所示)，而每一组设定将得到43种纹理特征，分别存放于名为glcmTextures、glrlmTextures、glszmTextures、ngtdmTextures和globalTextures的结构体中，最后总共5*6*2*4*43＝10320个纹理特征(如图24所示)；

13.对于不同病人，在将第一肝门层面的非纹理特征与纹理特征提取出来后，将非纹理特征和纹理特征整合成大小为病人个数*特征个数的特征矩阵

14.重复步骤A到步骤M完成对脾门层面感兴趣区域的非纹理特征与纹理特征提取，最后整合成大小为病人个数*特征个数的特征矩阵

15.利用str2matrix.m函数，将步骤M与步骤N中的两个特征矩阵横向拼接，组成大小为病人个数*(特征个数*2)的特征矩阵；以各个病人有无显著性门脉高压为标准制作名为outcome的mat文件，其中有显著性门脉高压用1表示，反之用0表示(如图25所示)；接着将特征矩阵和outcome的mat文件横向拼接起来，组成一个大小为病人个数*(特征个数*2+1)的矩阵，并且用csvwrite.m函数将已整合的第一肝门层面、脾门层面的图像特征矩阵保存为csv格式；

16.打开R软件，利用read.table函数将cvs格式的特征矩阵导入到R软件的工作空间中，并保存在变量train中；然后利用as.matrix函数将训练集中病人的图像特征矩阵提取出来，保存在train_x中作为建模的输入变量；同样地，利用as.matrix函数将训练集中病人的outcome标签提取出来，保存在train_y中作为建模的输出标签(如图26所示)；

17.利用library函数加载glmnet工具包，然后使用glmnet函数中的lasso降维方法对(4+5*6*2*4*43)＝20648个特征进行筛选，并使用dev.new函数和plot函数将降维过程显示出来(如图27所示)，图中每个lambda对应着一个线性模型且每个模型中使用的特征也不尽相同；

18.通过设置参数nfold、type.measure和family并利用cv.glmnet函数对训练集做nfold折交叉验证、建立线性的二项分布数学模型，然后通过dev.new函数和plot函数显示每个线性模型的误差曲线(如图28所示)；

19.通过cv$lambda.min选择AUC最高的lambda所对应的最优线性模型，然后利用coef函数和coefficients函数输出该最优模型的非零系数和筛选出来的特征；

20.利用library函数加载boot工具包和Hmisc工具包，然后通过boot函数和boot.ci函数对训练集做1000次bootstraping并计算出训练集的C-index指标(如图29所示)；接着加载pROC工具包，利用roc函数计算出训练集的ROC(Receiver OperatingCharacteristic)曲线下面积(AUC)；

21.为了进行模型验证，以未进行特征分析的人群制作测试集，同样经过上诉步骤A到步骤O将测试人群的第一肝门层面和脾门层面的20648维图像特征提取出来保存为csv格式；利用read.table函数将cvs格式的特征矩阵导入到R软件的工作空间中，并保存在变量test中；然后利用as.matrix函数将测试集中病人的图像特征矩阵提取出来，保存在test_x中作为已建好的模型的输入变量；利用glmnet工具包中的predict函数来计算测试集rHVPG的数值，测试结果与临床应用中金标准HVPG相比较，并利用roc函数计算出测试集的AUC(如图30所示)。

实施例1：R语言glmnet工具包进行特征筛选、建模与计算rHVPG

R语言编程界面中，首先用library(glmnet)导入glmnet工具包，并通过read.table("D:/train_list.csv",header＝TRUE,sep＝',')将包含有189例训练集的图像特征的csv数据导入变量train中。通过as.matrix(train[,1:20648])与as.matrix(train[,20649])将train中的特征矩阵与outcome的数据分别存放于变量train_x和train_y中。设置glmnet函数的参数：x＝train_x、y＝train_y、family＝‘binomial’、nfold＝10，整个glmnet函数调用的过程使用了lasso降维，并绘制特征筛选图(如图31所示)与误差曲线(如图32所示)，然后通过cv.glmnet函数进行10折交叉验证，选择参数lambda.min从而选定了最优的线性模型，此时对应着选取了15维特征所建立的线性函数，接着利用boot(test_all，cindex，R＝1000)对训练数据进行了1000次自助抽样，从而计算出该模型的C-index值与95％的置信区间为0.833(0.745-0.921)。接着通过coef函数将模型的系数和变量提取出来。用已经筛选出的特征建立以下rHVPG线性模型：

rHVPG＝-22.76522+14.94622×Complexity_0.5R_pS_LQ_64N_L

+8934.012×ZP_0.5R_1S_LQ_32N_L+0.4419001×Skewness_0.5R_3S_LQ_64N_L

-14.63828×LZE_0.5R_4S_EQ_64N_L+687.5762×ZP_0.5R_5S_EQ_16N_L

-0.00167672×HGZE_0.5R_5S_LQ_32N_L+19.31699×GLN_1R_2S_EQ_16N_L

-0.02605657 HGRE_1R_3S_LQ_8N_L-0.08363389 HGRE_1.5R_3S_LQ_8N_L

-3.724522e-06 Busyness_2R_1S_EQ_32N_L+545.8976 ZP_2R_1S_LQ_32N_L

+267.1424 GLN_2R_3S_EQ_8N_L+1361.243 ZP_2R_3S_EQ_16N_L

-1289.11 Strength_0.5R_pS_EQ_8N_S+432.8978Strength_1.5R_4S_LQ_64N_S；

公式中各个特征的意义如下：

1.常数：模型的修正系数；

2.Complexity：表示纹理特征的复杂度；

3.ZP：表示纹理特征的区域百分比；

4.Skewness：表示纹理特征的偏态；

5.LZE：表示纹理特征的大区域重点；

6.HGZE：表示纹理特征的高灰阶行程重点；

7.GLN：表示纹理特征的灰度不均衡性；

8.HGRE：表示纹理特征的高灰阶区域重点；

9.Busyness：表示纹理特征的冗余度；

Strength：表示纹理特征的强度。

命名规则：特征名_(数字)R_(数字/字母)S_(字符)Q_(数字)N_L/S。

其中，R为带通滤波器的参数、S为层厚、Q为灰阶量化的算法、N为指定量化的灰阶大小、L为肝脏的特征、S为脾脏的特征---具体数值选择参见步骤L或步骤12。

最后将包含有76例测试集的20648维特征的csv数据导入变量test_x中，通过predict(fit，newx＝test_x,s＝best,type＝'response')和最佳Lambda值对测试集test_x进行测试，结果如下：0.87714904、0.83785385、0.80030452、0.90876208、0.91561780、0.63210996、0.85066266、0.65695926、0.71528874、0.86608072、0.84943332、0.87832022、0.84873060、0.73834663、0.85322711、0.92394762、0.88492043、0.83537417、0.80363906、0.79202477、0.80414997、0.85716399、0.78553426、0.54361100、0.92121468、0.89490606、0.91388886、0.87067349、0.67488919、0.72509808、0.40844381、0.78257496、0.86903816、0.74041718、0.85528537、0.43124224、0.78625240、0.41466839、0.60465136、0.70270344、0.82158725、0.88783996、0.81765130、0.54013720、0.73824989、0.57662493、0.62293054、0.77785490、0.90732723、0.76901854、0.88960968、0.72091699、0.84619279、0.88151185、0.83710546、0.32607287、0.05237339、0.81813336、0.85601072、0.89002421、0.88068450、0.86283232、0.84371540、0.54572344、0.72121639、0.89244319、0.56050864、0.71049905、0.90255725、0.74028185、0.50714588、0.68410789、0.69254680、0.54036742、0.61098802、0.83925689选择合适阈值(roc工具包自动完成),此处为0.71，利用阈值对临床诊断有无显著门脉高压进行AUC计算，结果为AUC＝0.799(0.668-0.931)。AUC(ROC曲线下面积)是各个领域公认的评价模型预测/计算性能好坏的一个指标，范围在(0.5,1)之间，越大说明模型的预测/计算性能越好。

比较结果显示，本发明rHVPG对有创HVPG有较高的预测价值。鉴于有创HVPG已经纳入国际指南推荐的金标准，而本发明能够克服现有无创测量方法的诸多局限。因此，本发明提出的rHVPG无创新技术切实可行，有望为门脉压力的无创测量提供新的思路，为改善门脉高压患者的生活质量、减轻家庭和社会的疾病负担起到积极推动作用。

LASSO是一种实现稀疏数据的降维方法。针对线性模型，设该模型的函数表达如下：

其中，w^T表示该线性模型的系数组成的向量，x表示该模型的自变量。不同模型w^T与x各不相同。在临床应用中，往往得到的数据(病例)个数不多，但是每个病例所对应的自变量(特征)一般远大于数据个数的，这样直接用大量的特征进行线性模型的直接构造就会造成过拟合现象。设每个模使用均方误差评价性能好坏：

其中，y_i为期望输出(临床金标准)，n为数据的个数，f(x_i)为该线性模型的输出。当不给w任何限制时，往往导致w不为零的个数很多，造成过拟合现象。Lasso加入第一范数惩罚项限制w的取值，表达如下：

为了使模型性能最佳，等价于求解：

根据凸优化理论，绝大将的w会降维0，从而实现数据的稀疏与特征选择，为模型的构造从性能和复杂程度上讲都极大地优化了。

最后所应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。

Claims (2)

1.一种基于放射组学的肝静脉压力梯度计算模型的构建方法，其特征在于，包括以下步骤：

A.从标本的肘正中静脉注射造影剂，进行CT血管成像(CTA)，获取包括肝静脉期在内的CTA图层序列，导出图层序列，格式为dcm，层厚1.25mm，图像分辨率512×512像素；

B.将所获取的CTA图层序列导入医学影像控制软件ITK-SNAP，选择门静脉期图层序列，确认图像信息，ITK-SNAP软件自动识别图像序列，生成门静脉期CTA图像序列的冠状位、矢状位和水平位图像；

C.通过Expand this view to occupy the entire window功能将水平位图像放大，使用Zoom Inspector工具，选择zoom to fit调整图像至合适的大小；

D.以寻找水平位图像中的第一肝门作为目标，利用ITK-SNAP软件的PainburshInspector工具以包含目标的CT值、排除目标周围软组织的CT值为原则，提取目标；再反复利用Painbursh Inspector工具，剔除剩余的非目标结构；

E.选择Save segmentation Image将分割出的目标以NiFTI格式(.nii)导出，形成掩模图像的NiFTI格式；

F.打开dcm2niigui软件，选择Compressed FSL[4D NIfTInii]选项，将筛选出的含第一肝门层面的dcm文件序列放到dcm2niigui工作空间中，利用dcm2niigui软件将dcm文件序列转换为NiFTI格式；然后将含第一肝门层面的dcm文件序列、新生成的NiFTI格式文件和步骤D中生成的NiFTI格式文件放到同一文件夹中；

G.打开MATLAB R2016b软件，进入工作界面；在菜单栏处选择SetPath按钮，将所需要的feature_extra.m、readDICOMdir.m、str2matrix.m、NonTextureFeatures、TextureToolbox、NIfTI、Utilities、MultivariableModeling、STUDIES工具箱所在的绝对路径加载到MATLAB搜索路径中；

H.打开feature_extra.m文件，利用NIfTI工具箱中load_nii函数将步骤F中两个NiFTI文件导入MATLAB中，并且以名为data和data1的细胞数组储存；然后将两个细胞数组中的图像提取出来存储到名为mask和image变量中；

I.利用Utilities工具箱中computeBoundingBox函数提取出图像中感兴趣区域(ROI)的边界坐标提取出来储存到变量maskBox中；然后通过矩阵的基本操作将mask和image变量的边界坐标以内的图像信息提取出来放到变量maskBox和image_originalBox中；完成上述操作后，通过矩阵maskBox和image_originalBox的点乘运算得到图像中感兴趣区域，然后保存在ROI变量中，然后通过imshow与contour两个函数将感兴趣区域展示出来；

J.利用readDICOMdir.m函数将含第一肝门层面的dcm文件序列读入到MATLAB中，通过变量索引，将像素大小、层厚提取出来并保存到变量pixelW与sliceS中；

K.利用NonTextureFeatures工具箱中getEccentricity.m、getSize.m、getSolidity.m和getVolume.m函数，将变量ROI、pixelW、sliceS传输到相应函数中；经计算得到eccentricity、sizeROI、solidity、volumeROI 4个非纹理特征；

L.为了进行纹理特征提取，首先对参数进行设置；在TextureFeatures工具箱中一共有四个可变参数，分别为R、Scale、Algo和Ng，每个可变参数分别对应不同的变化空间，其中R对应为[1/2,2/3,1,3/2,2],Scale对应为[‘pixelW’,1,2,3,4,5],Algo对应为[‘Equal’,‘Lloyd’],Ng对应为[8,16,32,64]，而参数的排列组合使得纹理特征有显著变化；利用for循环将不同参数的取值进行排列组合，得到5*6*2*4组设定；通过参数传递向TextureFeatures工具箱传递设定，而每一组设定将得到43种纹理特征，分别存放于名为glcmTextures、glrlmTextures、glszmTextures、ngtdmTextures和globalTextures的结构体中，最后总共5*6*2*4*43＝10320个纹理特征；

M.对于不同病人，在将第一肝门层面感兴趣区域的非纹理特征与纹理特征提取出来后，一共得到了(4+5*6*2*4*43)＝10324个特征，然后将非纹理特征和纹理特征整合成大小为病人个数*特征个数的特征矩阵；

N.重复步骤A到步骤M完成对脾门层面感兴趣区域的非纹理特征与纹理特征提取，最后整合成大小为病人个数*特征个数的特征矩阵；

O.利用str2matrix.m函数，将步骤M与步骤N中的两个特征矩阵横向拼接，组成大小为病人个数*(特征个数*2)的特征矩阵；以各个病人有无显著性门脉高压为标准制作名为outcome的mat文件，其中有显著性门脉高压用1表示，反之用0表示；接着将特征矩阵和outcome的mat文件横向拼接起来，组成一个大小为病人个数*(特征个数*2+1)的矩阵，并且用csvwrite.m函数将已整合的第一肝门层面、脾门层面的图像特征矩阵保存为csv格式；

P.打开R软件，利用read.table函数将cvs格式的特征矩阵导入到R软件的工作空间中，并保存在变量train中；然后利用as.matrix函数将训练集中病人的图像特征矩阵提取出来，保存在train_x中作为建模的输入变量；同样地，利用as.matrix函数将训练集中病人的outcome标签提取出来，保存在train_y中作为建模的输出标签；

Q.利用library函数加载glmnet工具包，然后使用glmnet函数中的lasso降维方法对(4+5*6*2*4*43)*2＝20648个特征进行筛选，并使用dev.new函数和plot函数将降维过程显示出来，图中每个lambda对应着一个线性模型且每个模型中使用的特征也不尽相同；

R.通过设置参数nfold、type.measure和family并利用cv.glmnet函数对训练集做nfold折交叉验证、建立线性的二项分布数学模型，然后通过dev.new函数和plot函数显示每个线性模型的误差曲线；

S.通过cv$lambda.min选择AUC最高的lambda所对应的最优线性模型，然后利用coef函数和coefficients函数输出该最优模型的非零系数和筛选出来的特征；

T.利用library函数加载boot工具包和Hmisc工具包，然后通过boot函数和boot.ci函数对训练集做1000次bootstraping并计算出训练集的C-index指标；接着加载pROC工具包，利用roc函数计算出训练集的ROC(Receiver Operating Characteristic)曲线下面积(AUC)；

U.为了进行模型验证，以未进行特征分析的人群制作测试集，同样经过上诉步骤A到步骤O将测试人群的第一肝门层面和脾门层面的20648维图像特征提取出来保存为csv格式；利用read.table函数将cvs格式的特征矩阵导入到R软件的工作空间中，并保存在变量test中；然后利用as.matrix函数将测试集中病人的图像特征矩阵提取出来，保存在test_x中作为已建好的模型的输入变量；利用glmnet工具包中的predict函数来计算测试集rHVPG的数值，测试结果与临床应用中金标准HVPG相比较，并利用roc函数计算出测试集的AUC。

2.一种基于放射组学的肝静脉压力梯度的计算系统，其特征在于，包括：

数据输入模块，用于将患者第一肝门层面图像特征和脾门层面图像特征的检测结果输入模型计算模块；模型计算模块，包括基于放射组学的肝静脉压力梯度计算模型，用于根据患者第一肝门层面图像特征和脾门层面图像特征的检测结果以及基于放射组学的肝静脉压力梯度计算模型计算肝静脉压力梯度数值；所述基于放射组学的肝静脉压力梯度计算模型包括基于放射组学的肝静脉压力梯度计算公式，所述基于放射组学的肝静脉压力梯度计算公式：

rHVPG＝-22.76522+14.94622×Complexity_0.5R_pS_LQ_64N_L

+8934.012×ZP_0.5R_1S_LQ_32N_L+0.4419001×Skewness_0.5R_3S_LQ_64N_L

-14.63828×LZE_0.5R_4S_EQ_64N_L+687.5762×ZP_0.5R_5S_EQ_16N_L

-0.00167672×HGZE_0.5R_5S_LQ_32N_L+19.31699×GLN_1R_2S_EQ_16N_L

-0.02605657HGRE_1R_3S_LQ_8N_L-0.08363389HGRE_1.5R_3S_LQ_8N_L

-3.724522e-06Busyness_2R_1S_EQ_32N_L+545.8976ZP_2R_1S_LQ_32N_L

+267.1424GLN_2R_3S_EQ_8N_L+1361.243ZP_2R_3S_EQ_16N_L

-1289.11Strength_0.5R_pS_EQ_8N_S+432.8978Strength_1.5R_4S_LQ_64N_S；

结果输出模块，用于根据基于放射组学的肝静脉压力梯度计算结果判定患者肝静脉压力；当rHVPG≥0.71时，患者存在临床显著性门脉高压；当rHVPG＜0.71时，患者不存在临床显著性门脉高压。