如何应对MR放射组学中的多中心变异性？验证校正程序

  目标：测试一种实用的重新校准方法，以补偿MR放射组学特征的技术变异性。

     方法：对两台扫描仪获取的水模图像、18名脑瘤患者的FLAIR和增强T1加权（CE-T1w）图像（分别在1.5-T和3-T扫描仪上扫描），以及在两个中心之一扫描的36张T2加权（T2w）前列腺癌患者图像进行研究。使用ComBat程序对放射组学特征进行调和。在调和前后，测试了1.5-T和3-T图像间特征值的统计分布差异。前列腺研究用于确定调和对区分Gleason等级（GGs）的影响。

     结果：在水模数据中，在调和之前，42个放射组学特征值中有40个在两台扫描仪之间显著不同，而在调和后则没有。在白质区域，FLAIR和CE-T1w图像中的42个特征中有37个特征在1.5-T和3-T图像之间的统计分布显著不同（p < 0.05）。调和后，未观察到统计学上显著的差异。在脑肿瘤中，未经调和的1.5-T和3-T图像中，42个特征中有41个（FLAIR）或36个（CE-T1w）特征显著不同，而经过调和后，只有1个（FLAIR）或没有（CE-T1w）特征差异显著。在前列腺研究中，调和后有636个放射组学特征在Gleason分级（GGs）之间显著不同，而调和前为461个。使用放射组学特征区分GGs的能力在调和后增加。

      结论：ComBat调和能有效地消除放射组学特征值中的跨中心技术不一致性，并增加了使用多台扫描仪数据的研究的敏感性。

    本文发表在European Radiology杂志。（可添加微信号19962074063或18983979082获取原文及补充材料，另思影提供免费文献下载服务，如需要也可添加此微信号入群)

关键点：

• 使用不同MR扫描仪或成像协议获得的放射组学特征值可以通过结合现成的图像标准化和特征重新校准程序来调和。

• 调和后的放射组学特征使得可以汇集不同扫描仪和中心的数据，而不会因中心内和中心间的变异性而显著损失统计功效。

• 所提出的重新校准方法适用于来自不同MR序列和肿瘤类型的放射组学特征，且不依赖于任何水模采集。

缩写：

• CE-T1w：增强T1加权（Contrast-enhanced T1-weighted）

• CT：计算机断层扫描（Computed Tomography）

• D1/D2：前列腺癌患者数据库1/2

• hWS：混合白色条纹（Hybrid White Stripe）

• LDA：线性判别分析（Linear Discriminant Analysis）

• MRI：磁共振成像（Magnetic Resonance Imaging）

• PET：正电子发射断层扫描（Positron Emission Tomography）

• ROI：感兴趣区域（Region of Interest）

• T2w：T2加权（T2-weighted）

• VOI：感兴趣体积（Volume of Interest）

• WM：白质（White Matter）

引言 

      放射组学面临着可重复性不足的关键问题，这阻碍了将放射组学模型发现成功转化为更好的诊断、患者分类或基于放射组学的监测工具。事实上，放射组学特征值显著受到成像设备和协议的技术设置的影响，这已在正电子发射断层扫描（PET）、计算机断层扫描（CT）和磁共振成像（MRI）中得到证明[1,2,3]。为了解决由不同技术设置引起的放射组学特征的变异性，可以设计包含大多数技术设置的各种图像的放射组学模型，或者在设计模型之前对图像或放射组学特征值进行调和。忽略中心效应，正如在许多论文中常见的那样，导致放射组学模型缺乏普遍性[4]。

      在前瞻性研究中，可以在中心之间上游协调成像协议，以最小化成像协议对特征值的影响[5, 6]，尽管在不同代的机器之间协调通常伴随着降低最新扫描仪实现的图像质量[5]。在回顾性研究中，这种方法不可行。几个小组已经提出通过将图像重采样到共同的体素大小或通过过滤图像以匹配空间分辨率来减少变异性[7, 8]。然而，这需要事后访问图像，且过滤程序会降低使用最新设备获取的图像的质量。其他人应用z得分转换[9]到每个特征值，基于每个中心该特征的平均值和标准偏差，但这假设不同中心产生的图像是在类似的患者样本中获得的（例如，同比例的晚期和早期肿瘤），这有时难以实现。

      在基因组学中，研究人员面临着一个类似的问题，称为批次效应，由不同实验室、不同技术人员或不同天处理样本导致，可能会掩盖个体变异。为了解决基因组学中的这个问题，Johnson等人[10]引入了ComBat重新校准方法。该方法将所有数据重新校准到一个单一空间中，其中批次效应被丢弃，而不改变生物学信息。这种方法已经在支持调和重要性的研究中，成功地验证了从PET[11]和CT[12, 13]图像中测量的患者或水模数据的放射组学特征。

     在MR（磁共振成像）中，挑战更加艰巨，因为与PET（正电子发射断层扫描）和CT（计算机断层扫描）不同，其中图像分别以kBq/mL和Hounsfield单位表示，MR并没有标准的灰度强度级别。这意味着即使在相同的MR成像协议、相同身体区域、同一扫描仪获取的图像，以及同一患者的情况下，也缺乏组织特异性的绝对强度数值含义。因此，为了比较基于强度的特征值，对患者间图像强度的标准化是绝对必要的。在脑部MR中，已经提出了一些标准化方法来纠正强度的变异性[3, 14,15,16]。特别是，混合白色条纹法（hWS）在神经退行性疾病和癌症研究中被证明是成功的[17, 18]。ComBat已在MRI中被验证，用于跨扫描仪调和皮层厚度测量[19]。尽管它已在MR放射组学研究中被使用[20,21,22,23]，但在这一高度具有挑战性的背景下，它从未被验证。

     在这里，我们扩展了ComBat方法，提供了一种适用于任何放射组学特征的调和程序。我们展示了通过结合图像标准化（如混合白色条纹法）与ComBat重新校准，MR放射组学特征可以在不受多种变异源负面影响的情况下汇总，从而确保了多中心MR放射组学研究的更高敏感性和特异性。

材料和方法 

     本研究按照世界医学会赫尔辛基宣言进行。对于涉及MRI脑部研究的实验2，Fondation Ophtalmologique A. Rothschild的机构审查委员会批准了本研究（IRB编号1512-016-726），由于研究的回顾性质，放弃了获得书面知情同意的要求。对于使用MRI前列腺研究的实验3，所有患者数据公开可用，并最初在[24]中研究。所有患者数据均已匿名化。所有作者对提交出版的数据和信息有控制权。

实验1：水模研究 

      本研究中使用的水模数据来自RIDER水模MRI研究[25]，在TCIA平台[26]上公开可用。该水模由19个含有钆基对比剂的掺杂凝胶填充管组成。我们使用了1.5-T扫描仪和3-T扫描仪（分别为[25]中的扫描仪B和D，详见补充数据1）获得的T1采集。对于每张图像，我们在每个分析管道上绘制了19个3.5毫升的球形感兴趣体积（VOIs）。我们使用LIFEx软件[27]计算了42个放射组学特征（补充数据2），该软件包括一个符合图像生物标志物标准化倡议指南[28]的开源放射组学协议。放射组学特征的计算使用固定的箱体大小[3, 29]，设置为所有VOIs中测量的信号强度的平均标准偏差，从最小强度到最大强度。这个离散化步骤是必需的，以将具有相似强度的体素设置为相同的值，从而减少噪声的影响。

实验2：MRI脑部研究 

       对于实验2，我们回顾性地选择了18名患者（13男；平均年龄50 ± 18岁；年龄范围26-85岁；表1），这些患者患有III级和IV级胶质瘤，时间为2017年1月至2018年5月，来自于一个机构数据库。所有患者均接受了两次使用相同协议的MRI扫描：一次在1.5-T扫描仪（飞利浦Achieva，飞利浦医疗系统）上，另一次在3-T扫描仪（飞利浦Ingenia）上。两次扫描之间的中位间隔时间为30天（范围4-93天），期间没有进行化疗、手术、放疗，且两次扫描间肿瘤和肿瘤异质性没有可见变化。获得了两种MR序列（详见补充数据1）：3D FLAIR扫描（17名患者）和3D增强T1加权（CE-T1w）扫描（14名患者）。

表1 病人特征

     对于每位患者和每种序列，3-T图像使用FSL-FLIRT中的刚性变换被配准到1.5-T图像上[30]。使用N4算法[31]校正了磁场不均匀性，该算法通过公开可用的ANTs软件实现，并采用了超参数的标准设置。

      对于每种序列，肿瘤病变基于两位放射科医师（A.L.和L.D.，分别拥有9年和2年的经验）的共识在1.5-T图像上手动分割，得到的区域被复制到3-T图像上。在每个肿瘤中选择三个切片（顶部、中部、底部）以获得每个肿瘤三个2D感兴趣区域（ROI），共计FLAIR图像的54个肿瘤ROI（=3×18个肿瘤；一名患者有两个不同的病变）和CE-T1w图像的51个肿瘤ROI（=3×17个肿瘤；一名患者有两个不同的病变，另一名有三个）。此外，在每位患者中，分别在白质（WM）中绘制了6个0.5毫升的区域，共计FLAIR图像的102个WM-VOI和CE-T1w图像的84个WM-VOI，这些被复制到3-T图像上。

      使用先前描述的混合白色条纹法（hWS）[17]，每位患者的图像体积无视其他患者进行了标准化[3]。hWS方法对基于正常出现的WM强度分布的脑体素值应用z得分转换。

      对于基于原始和hWS标准化图像（重采样为1×1×1 mm³）的每个ROI和VOI，我们使用LIFEx计算了42个放射组学特征。放射组学特征的计算使用固定的箱体大小[3, 29]，设置为WM信号强度的平均标准偏差，介于所有WM和肿瘤VOI中测量的最小和最大强度之间，每个序列分别计算（详见补充数据1）。

实验3：MRI前列腺研究 

      使用了两个公开可用放射组学特征的前列腺癌患者数据库（D1和D2；表1）[24]。这两个数据库最初由一个独立团队完整研究，以探讨从MR图像和数字化组织图像计算的特征之间的关系，以区分前列腺癌等级，但没有考虑到MRI扫描是在两个不同中心获取的[24]。在这里，我们准确地研究了考虑中心效应实际上如何改变每个MRI特征区分肿瘤等级的能力。实验方案获得克利夫兰大学医院机构审查委员会IRB方案#02-13-42C的批准。患者在根治性前列腺切除术前接受了T2加权（T2w）MRI。在D1中，宾夕法尼亚大学的23名患者在2009年至2011年间进行了扫描（3 T Verios, Siemens Healthcare; 回波时间107-127毫秒；重复时间3690-7090毫秒）。在D2中，圣文森特医院的13名患者在2012年至2014年间进行了扫描（11名患者[3 T, Philips Medical Systems; 回波时间67-100毫秒；重复时间2525-3567毫秒]和2名患者[1.5 T, Siemens Healthcare; 回波时间119毫秒；重复时间3760毫秒]）。手术后，病理学家分析了切除的前列腺腺体，以确定Gleason评分，将其分类为低风险（3 + 3分）或中/高风险（3 + 4分，4 + 3分，4 + 4分或更高）。D1包括21个低风险区域和44个中/高风险区域，而D2包括26个低风险区域和14个中/高风险区域（表1）。基于与组织学图像的配准，相应的肿瘤区域由放射科医师手动分割。MR图像根据D1中患者的前列腺内像素强度的中位数标准化为模板分布[32]。对于每个区域，使用自制软件计算了2379个放射组学特征（特征计算的详细信息在[24]中描述），我们为分析选择了所有患者均可用的2326个特征。

重新校准方法 

      为了纠正扫描仪效应，使用了ComBat重新校准方法[10]。在放射组学的背景下，ComBat已经为PET[11]和CT特征[12, 13]得到验证。该方法直接应用于放射组学特征值，并通过匹配在每个扫描仪i上测量的感兴趣体积（VOI）j的特征值的统计分布来估计扫描仪效应。

      其中，α是特征y ij的平均值，X是感兴趣协变量的设计矩阵，β是与每个协变量对应的回归系数向量，γ i是添加性扫描仪效应，δ i是受误差项（ε ij）影响的乘法协议效应。模型参数α、β、γ i和δ i使用基于两台扫描仪在实验1和实验2中可用观测值的集合，以及基于两个患者数据库在实验3中，采用最大似然方法估计。校正后的值使用以下方法获得：

       其中

分别是 α, β, γi 和 δi 的估计值。

       该模型的非参数形式被使用，对特征遵循的统计规律没有假设，并且为每个特征单独确定了一个转换。对于实验1和实验2，没有使用生物学协变量（即，X = 0），因为数据来自于在1.5-T和3-T机器上扫描的同一患者或水模，我们分别对患者数据中的WM和肿瘤区域计算的特征值进行重新校准。对于实验3，由于两个数据库中低风险与中/高风险区域的比例差异很大（D1中低风险VOI占32%，D2中占68%；表1），我们将Gleason等级作为二元协变量引入。

       为了便于医学成像专业人员访问ComBat方法，我们提供了一个免费的在线应用程序，名为ComBaTool，其中包含示例输入文件（补充数据3和4）和逐步教程（补充数据5）。这个应用程序嵌入了一个名为ComBat[19]的函数，基于R软件。

统计分析 

统计分析使用R软件（版本3.6.1）进行。

      在实验1中，我们在ComBat重新校准前后对两次水模扫描进行了双侧Friedman检验。在实验2中，我们对每个放射组学特征使用了双侧Friedman检验，以测试从1.5-T和3-T扫描得出的值在以下三种配置中是否在白质（WM）和肿瘤区域显著不同：（C1）原始图像未经ComBat重新校准，（C2）hWS标准化图像未重新校准，（C3）hWS标准化图像经过重新校准。使用Benjamini-Hochberg程序来控制假发现率[33]。p值小于0.05被解释为统计学上显著。绘制Bland-Altman图以展示从1.5-T和3-T扫描计算的特征值的差异。

      在实验3中，我们使用Benjamini-Hochberg程序对所有放射组学特征进行了Wilcoxon检验，以区分在汇总D1和D2患者时的低风险和中/高风险组，分别为未经ComBat重新校准、经过重新校准以及包括Gleason等级作为感兴趣协变量的重新校准。为了证明ComBat不会产生假阳性结果，我们在随机分配标签给每个VOI后重复了这些测试，得到53个虚假低风险VOI和52个虚假中/高风险VOI。为了识别风险组，我们通过使用D1数据集作为训练集，并只包括单变量Wilcoxon检验p值小于5%的特征，通过线性判别分析（LDA）构建了一个多变量标志。我们通过计算Youden指数（=敏感性 + 特异性 - 1）来测试在D2数据上的分类性能。我们在三种配置中重复了这个过程：未经ComBat重新校准、经过重新校准以及包括Gleason等级作为感兴趣协变量的重新校准。

结果 

患者特征见表1。

实验1 

      在水模数据中，42个Friedman检验的p值中有40个在未重新校准的情况下低于5%。两个扫描之间只有两个p值（粗糙度和灰度级非均匀性）大于0.05。ComBat之后，所有Friedman检验的p值都大于0.05，显示协议效应不再可检测。

实验2 

       未经hWS标准化和ComBat重新校准的1.5-T和3-T原生FLAIR脑部图像中，从WM-VOI计算的42个放射组学特征中有37个（88%），从肿瘤病变中计算的有41个（98%），Friedman检验的p值小于0.05（表2；补充数据6）。使用MR图像的hWS标准化后，WM区域的29/42（69%）p值和肿瘤病变的25/42（60%）p值小于0.05。结合hWS标准化和ComBat特征分布重新校准，只有一个p值（长区域强调）小于0.05（p = 0.017），表明扫描仪效应对绝大多数放射组学特征不再可检测。图1显示了从灰度共生矩阵计算的相关放射组学特征分布的演变。在原生FLAIR图像上，图表显示分布的偏移，3-T扫描的WM-VOI和肿瘤病变的值大于1.5-T扫描。经过hWS标准化和重新校准后，两个扫描仪之间的分布更好地重叠。为了阐明hWS和ComBat各自的作用，图2显示了基于3-T扫描和1.5-T扫描的FLAIR图像中WM-VOI的平均值的Bland-Altman图。hWS在每个患者内的标准化重新调整了值，使它们在两次扫描之间变得相似。重新校准减少了两者之间的系统差异。

表2 显示了对脑部扫描进行的显著Friedman检验（p < 0.05），包括未经混合白色条纹（hWS）标准化和/或ComBat重新校准的情况。

图1 实验2。18名患有脑部病变的患者在1.5-T和3-T扫描仪上进行了扫描。基于原生或混合白色条纹（hWS）标准化图像，在肿瘤区域和白质区域计算了42个放射组学特征。作为一个例子，下面绘制了基于FLAIR图像的灰度级共生矩阵计算出的相关放射组学特征的概率密度函数（%），展示了未经ComBat重新校准和经过ComBat重新校准（ComBaTool分别应用于两种组织类型：白质和肿瘤）的情况，对于1.5-T MRI（橙色）和3-T MRI（蓝色）。p值是针对两种MRI设备间每种组织的Friedman检验结果。

图2 实验2。基于1.5-T和3-T扫描的FLAIR原生图像（a）、混合白色条纹（hWS）标准化图像（b）和hWS标准化图像经ComBat重新校准（c）的白质区域中计算出的平均值的Bland-Altman图。

相同的趋势也在CE-T1w图像上观察到（表2；补充数据7）。

       实验3 在[24]执行的标准化后，T2w前列腺图像中的2326个放射组学特征中有461个特征的Wilcoxon检验p值小于0.05，用于区分低风险和中/高风险，当汇总两个患者队列（D1 + D2）时。在未使用任何协变量的ComBat之后，2326个p值中有460个小于0.05。在ComBat中使用Gleason等级协变量时，2326个p值中有636个小于0.05。图3展示了在使用ComBat并考虑每个中心招募特异性的协变量后，从两组患者中分别提取的低风险VOI和中/高风险VOI的放射组学特征值更好地对齐。

图3.实验3。特征#20（在[24]中称为Gabor:cos:theta=0:lambda=2:Standard Deviation）的箱线图，用于低风险VOI和中/高风险VOI，在ComBat重新校准前（a, d），ComBat重新校准后无协变量（b, e），以及ComBat重新校准后有协变量（c, f），分别针对前列腺患者队列D1和D2（a–c）或一起（d–f）。p值来自Wilcoxon检验。

      当随机分配风险（低或中/高）给每个VOI时，无论在ComBat重新校准之前还是之后，无论是否包含代表Gleason等级的协变量，都没有p值小于0.05。

      使用LDA在D1数据上识别的多变量放射组学模型，用于区分低风险和中/高风险，在应用于D2患者时，得出的Youden指数为0.12（敏感性=19%，特异性=93%），这是在ComBat之前的结果。在ComBat之后，Youden指数提高到0.20（敏感性=27%，特异性=93%），并在使用Gleason等级作为ComBat中的协变量时提高到0.43（敏感性=58%，特异性=86%）。

讨论

     扫描仪效应影响从MR图像中提取的放射组学特征值，为多中心或多协议研究引入重大的混杂因素。在这里，我们验证了一种结合ComBat重新校准和MR图像标准化的调和程序，以共同分析来自不同扫描仪的MR放射组学特征。使用水模数据和为同一患者（在两次扫描之间未视觉检测到任何肿瘤演变）获取的1.5-T和3-T扫描仪的脑部扫描，我们展示了这种调和程序如何重新校准放射组学特征分布并消除T1、FLAIR和CE-T1w图像的扫描仪效应。我们的目标不是测试我们重现3-T MR图像中测量的特征值的能力，因为我们预期两个设备会产生不同的信号，3-T图像会有更多细节（参见图1）。然而，在放射组学的背景下，经常需要汇总使用不同设备以及不同获取和重建协议获得的图像，以增加队列的规模。在这种情况下，我们展示了ComBat可以重新校准特征值，使所有数据都可以一起分析，即使图像是在不同磁场强度下获取的。重要的是要强调，因为成像协议对每种组织的影响不同，所以对每个序列和每种组织类型独立估计不同的ComBat转换。使用在两个中心的不同患者身上获取的前列腺扫描，我们证实了T2w图像的调和有效性，并证明调和没有改变特征所传达的区分信息。这个实验还表明，汇总校正了扫描仪效应的数据可以增加统计功效，识别更多能够区分前列腺病变中低风险和中/高风险区域的放射组学特征，并产生一个更具区分性的多变量模型。重要的是，我们展示了当不期望两组之间存在差异时，这里是在虚假低风险和中/高风险VOI之间，ComBat并没有引入任何假阳性差异。

     ComBat重新校准方法快速且易于使用，直接作用于放射组学特征值（不需要训练集，不需要水模采集，也不需要访问图像）。它适用于从不同MR序列提取的放射组学特征，在图像标准化的第一步之后，如之前所述[3]。我们还展示了在重新校准过程中协变量的附加价值，当患者特征在中心之间不同时（这里是Gleason等级）对单变量和多变量分析。为了处理中心效应，其他作者报告了生成对抗网络将图像从一个成像协议（或域）转换到另一个的潜力[34]。尽管文献中报告了有前景的结果[35, 36]，但这些技术需要访问图像，与ComBat不同。ComBat重新校准方法之前已在MR放射组学研究中使用过[20,21,22,23]，但没有进行任何明确的验证或调查图像标准化的相应作用以及扫描仪/协议效应补偿，如本研究所研究的（图1和图2）。在[20]中，作者报告了在ComBat之后，从表观扩散系数MR图像提取的熵对宫颈癌局部控制的预测准确性增加，与我们的发现完全一致。

       我们的研究有一些局限性。在实验2中，我们只能纳入18名患者，因为患者在肿瘤没有视觉上演变的时间间隔内接受1.5-T和3-T扫描仪的MR检查是非常罕见的。尽管如此，这个小样本还是让我们能够确认使用水模数据获得的结果。此外，如此小的样本量还允许我们证明，即使在病例数量有限的情况下，ComBat也表现良好，从而确认了在基因组学应用中发表的结果[10]。另一个局限是，我们的发现应该还需要针对其他癌症类型、MR序列和设备进行验证。

     总之，我们证明了ComBat重新校准方法结合患者内部图像标准化可以有效地消除扫描仪/协议效应，同时保留水模、脑部和前列腺MR扫描中的个体变异。这种方法使大型MR多中心研究能够调查放射组学分析在患者管理中的附加价值。为了促进大型多中心/多协议放射组学研究，我们提供了ComBat方法作为在线ComBaTool应用程序供研究者使用。