A radiomics model from joint FDG-PET and MRI texture features for the prediction of lung metastases in soft-tissue sarcomas of the extremities

2.1.1. Patient cohort.

2.1.2. Imaging data.

Contours defining the 3D tumour region for each patient were manually drawn slice-by-slice on T2FS scans by an expert radiation oncologist. For patients with visible edema in the vicinity of the tumours (n = 32), two contours were drawn: one incorporating the visible edema and one excluding it, as shown in figure 2. Contours were propagated to FDG-PET and T1 scans using rigid registration with the commercial software MIM^® (MIM software Inc., Cleveland, OH). The results presented in this work were obtained from texture analysis performed on the volume of interest of each patient as defined by the contour containing no edema.

Zoom In Zoom Out Reset image size

Prior to texture analysis, FDG-PET and MRI DICOM data were transferred into MATLAB^® (The MathWorks Inc., Natick, MA) format using the software CERR (Deasy et al 2003). All subsequent analyses were performed in MATLAB^®. FDG-PET scans were first converted to standard uptake value (SUV) maps, followed by the application of a square-root transform to help in the stabilization of the PET noise in the images. MRI scans were kept in raw data form, and voxels within the tumour region with intensities outside the range μ ± 3σ were rejected and not considered in subsequent texture computations, as suggested by Collewet et al (2004) for making MRI texture measurements more reliable.

The fusion of FDG-PET and MRI volumes first starts with the registration of the scans as described in section 2.2. The 3D discrete wavelet transform (DWT) is then used to combine the spatial and frequency characteristics of the two modalities as follows:

• (i)  
  Downsample the MRI volume (raw data, no pre-processing) to the resolution of the FDG-PET volume (pre-processed, see section 2.3) using cubic interpolation. Normalize the intensity range of FDG-PET and MRI tumour regions between 0 and 255. Invert MRI intensities if needed.
• (ii)  
  Apply the 3D DWT to the FDG-PET and MRI volumes up to one decomposition level using the wavelet basis function symlet8.
• (iii)  
  Apply the μ ± 3σ normalization scheme (see section 2.3) respectively to the wavelet coefficients of the tumour region of the different MRI sub-bands. The rejected MRI wavelet coefficients are then replaced by the spatially corresponding coefficients of the FDG-PET sub-bands.
• (iv)  
  Perform a weighted average of the spatially corresponding wavelet coefficients of all PET and MRI sub-bands to obtain a single set of fused wavelet coefficients. If the weight given to MRI wavelet coefficients is denoted as 'MRI weight' and the weight given to PET wavelet coefficents is denoted as 'PET weight', MRI weight ranges from 0 to 1 and PET weight = 1 − MRI weight.
• (v)  
  Apply the 3D inverse DWT to the set of fused wavelet coefficients using the reconstruction wavelet basis function symlet8 to obtain a fused FDG-PET/MRI tumour volume.

The choice of the wavelet basis function symlet8 is based on our previous work (Vallières 2013), in which that basis function was shown to produce fused textures with best predictive value. This could be explained from the fact that symlets is a family of orthogonal and compactly supported wavelets with the least asymmetry and highest number of vanishing moments for a given support width, which would help in the local preservation of spatial characteristics of images.

The fusion process yields two new types of scans: FDG-PET/T1 and FDG-PET/T2FS. We also tested if the inversion of MRI intensities prior to fusion with FDG-PET could enhance texture characteristics in the fused volumes. Figure 3 shows an example of the fusion of FDG-PET and T2FS scans for a patient of the LungMets group.

Zoom In Zoom Out Reset image size

The methodology used to extract the imaging features from the tumour region of the pre-treatment FDG-PET and MRI scans is described below.

2.5.1. Non-texture features (9).

2.5.2. Texture features (41).

2.5.3. Texture extraction parameters (6).

The influence of the following six extraction parameters on the predictive value of textures was investigated.

Fusion parameters (2).

The following two parameters apply only to fused scans (FDG-PET/T1 and FDG-PET/T2FS scans):

• (i)  
  Inversion of MR imaging intensities in the FDG-PET/MRI fusion process (see section 2.4). This parameter is denoted as 'MRI Inv'. MRI Inv. of 0 and 1 (no inversion/inversion) were tested.
• (ii)  
  Weight applied to MRI wavelet sub-bands in the FDG-PET/MRI fusion process (see section 2.4). This parameter is denoted as 'MRI weight'. MRI weight values of $\frac{1}{4}$, $\frac{1}{3}$, $\frac{1}{2}$, $\frac{2}{3}$ and $\frac{3}{4}$ were tested in this work.

Wavelet band-pass filtering (1).

This parameter is denoted as 'WBPF'. This operation is carried out by applying a different weight to band-pass sub-bands (LHL, LHH, LLH, HLL, HHL and HLH) of the tumour region as compared to low- and high-frequency sub-bands (LLL and HHH) in the wavelet domain. The ratio of the weight applied to band-pass sub-bands to the weight applied to the other sub-bands is defined by R. Ratios of $\frac{1}{2}$, $\frac{2}{3}$, 1, $\frac{3}{2}$ and 2 were tested.

Isotropic voxel size (1).

This parameter is denoted as 'Scale'. Prior to the computation of texture features, all volumes were resampled to an isotropic voxel size set to a desired resolution using cubic interpolation. Scale values of 1 mm, 2 mm, 3 mm, 4 mm, 5 mm and initial in-plane resolution (denoted 'in-pR') were tested. For example, if the desired Scale was set to 5 mm, a FDG-PET volume with voxel size of 5.47 × 5.47 × 3.27 mm³ would be isotropically resampled to a voxel size of 5 × 5 × 5 mm³. If the desired Scale was set to in-pR, a MRI volume with voxel size of 0.86 × 0.86 × 5 mm³ would be isotropically resampled to a voxel size of 0.86 × 0.86 × 0.86 mm³. Note that Global texture features are extracted after isotropically resampling to the in-plane resolution without further processing.

Quantization of gray levels (2).

Prior to the computation of texture features, the full intensity range of the tumour region was quantized to a smaller number of gray levels N_g. The quantization process maps the voxel values of a volume to a finite set $\boldsymbol{r}=\left\{{{r}_{k}}\in \mathbb{R}:k=1,2,\ldots ,{{N}_{\text{g}}}\right\}$ of reconstruction levels by defining a set $\boldsymbol{t}=\left\{{{t}_{k}}\in \mathbb{R}:k=1,2,\ldots ,{{N}_{\text{g}}}+1\right\}$ of decision levels. Two extraction parameters are related to the quantization of gray levels in the volumes:

• (i)  
  Quantization algorithm. This parameter is denoted as 'Quant. algo'. Equal-probability and Lloyd-Max quantization algorithms were implemented in this work using the functions histeq and lloyds of MATLAB^®, respectively. Equal-probability quantization attempts to define decision thresholds in the volume such that the number of voxels with reconstructed level r_k is the same in the quantized volume for all k (i.e. for all gray levels), whereas Lloyd-Max quantization attempts to minimize the mean-squared quantization error of the output (Max 1960, Lloyd 1982).
• (ii)  
  Number of gray levels (N_g) in the quantized volume. N_g of 8, 16, 32 and 64 were tested.

Texture extraction summary.

Considering the full set of texture extraction parameters of Global features and higher-order texture features (GLCM, GLRLM, GLSZM and NGTDM), a total of 27 405 and 182 700 scan-texture-parameter combinations were computed in this work for separate and fused scans, respectively. Figure 4 presents a summary of the workflow of extraction of texture features.

Zoom In Zoom Out Reset image size

Univariate association between the whole set of features (9 non-texture features and 210 105 scan-texture-parameter features) and lung metastases development in STSs was investigated using Spearman's rank correlation (r_s). To correct for multiple tests comparisons, the Bonferroni correction method was applied: the significance level was lowered to a value p < α/K, where K is the number of comparisons and α the significance level set to 0.05.

The process of combining features into a multivariable model was achieved using the logistic regression utilities of the software DREES (El Naqa et al 2006). We are interested in finding a linear combination of p variables such that the multivariable model of interest takes the form:

$$\begin{eqnarray}&&g({{\mathbf{x}}_{i}})={{\beta}_{0}}+\underset{j=1}{\overset{p}{\sum}}\,{{\beta}_{j}}{{x}_{ij}},\quad \text{for}\;i=1,2,\ldots ,N,\end{eqnarray} \tag{ 1 }$$

where the vector of input variables (imaging data) of the ith patient is ${{\mathbf{x}}_{i}}=\left\{{{x}_{ij}}\in \mathbb{R}:j=1,2,\ldots ,p\right\}$, for a number N of patients. The set $\boldsymbol{\beta}=\left\{{{\beta}_{j}}\in \mathbb{R}:j=0,1,\ldots ,p\right\}$ is the set of regression coefficients of the model to be determined such that the conditional probability of the set of outcome states {0,1} given the input data x_i is maximized for i = 1, 2, ..., N. This operation is carried out using a logistic regression model (logit transformation) of the form:

$$\begin{eqnarray}&&\pi ({{\mathbf{x}}_{i}})=\text{P}({{y}_{i}}=1\vert {{\mathbf{x}}_{i}})=\frac{\text{exp}\left[g({{\mathbf{x}}_{i}})\right]}{1+\text{exp}\left[g({{\mathbf{x}}_{i}})\right]},\quad \text{for}\;i=1,2,\ldots ,N.\end{eqnarray} \tag{ 2 }$$

Following the work by Sahiner et al (2008), we adopted the 0.632+ bootstrap method and the area under the receiver operating characteristic curve (AUC) metric to estimate which models learned from our patient cohort would best predict lung metastases on new prospective data from the whole (or true) STS population. Let AUC(S_train, S_test) denote the value of the test AUC obtained when the classifier is trained on set S_train (computing logistic regression coefficients) and tested in set S_test (testing g(x_i)). Also, let the observed sample (imaging data of our patient cohort) be denoted as the matrix X={x_i : i = 1, 2, ..., N}. A bootstrap sample denoted as ${{\mathbf{X}}^{*}}=\left\{\mathbf{x}_{i}^{*}:i=1,2,\ldots ,N\right\}$ is a sample of input variables x_i of N patients randomly drawn with replacement from the available sample X. The set of original data vectors that do not appear in X* is denoted as X*(0). The generation of a large number B of randomly drawn bootstrap samples X*^b for b = 1, 2, ..., B is used to estimate a statistical quantity of interest on the unknown true population distribution. With the 0.632+ bootstrap method, the estimated AUC is then calculated as:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{[\widehat{{\text{AUC}}}\,]}_{0.632+}}=\frac{1}{B}\underset{b=1}{\overset{B} \sum}\,[(1-\alpha (b))\cdot \text{AUC}(\mathbf{X},\mathbf{X})+\alpha (b)\cdot \text{AU}{{\text{C}}^{\prime}}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))], \\ \text{where}\qquad\qquad \text{AU}{{\text{C}}^{\prime}}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))=\text{max}\left\{0.5,\text{AUC}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))\right\},\,\alpha (b)=\frac{0.632}{1-0.368\cdot R(b)}, \\ \text{and}\qquad\quad R(b)=\left\{\begin{array}{*{35}{l}} 1 & \text{if}~\text{AUC}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))\leqslant 0.5, \\ \frac{\text{AUC}(\mathbf{X},\mathbf{X})-\text{AUC}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))}{\text{AUC}(\mathbf{X},\mathbf{X})-0.5} & \text{if}~2>\frac{\text{AUC}(\mathbf{X},\mathbf{X})}{\text{AUC}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))}>1, \\ 0 & \text{otherwise}.\end{array} \right.\\ \end{array}\end{eqnarray} \tag{ 3 }$$

In this work, each time a bootstrap sample X*^b was drawn from X in the multivariable analysis, the probability of choosing a negative instance (NoLungMets patient group class) was made equal to the probability of choosing a positive instance (LungMets patient group class), a procedure hereby denoted as 'imbalance-adjusted bootstrap resampling'.

Prediction models were constructed for three different types of initial feature sets: (i) 9 non-texture features + 9 135 FDG-PET scan-texture-parameter features; (ii) 9 non-texture features + 27 405 separate FDG-PET and MRI scan-texture-parameter features; and (iii) 9 non-texture features + 182 700 fused FDG-PET/MRI scan-texture-parameter features. First, feature set reduction was performed through a stepwise forward feature selection scheme in order to create reduced feature sets containing 25 different scan-texture features from larger initial sets, a procedure carried out using the Gain equation:

$$\begin{eqnarray}\begin{array}{*{20}{l}}{{\widehat{\text{Gain}}}_{j}}&=\gamma \cdot \left\vert {{\widehat{r}}_{\text{s}}}({{\mathbf{x}}_{j}},\mathbf{y})\right\vert\\ &+{{\delta}_{a}}\cdot \left[\underset{k=1}{\overset{f}{\sum}}\,\left(\frac{2(f-k+1)}{f(f+1)}\right)\widehat{\text{PIC}}({{\mathbf{x}}_{k}},{{\mathbf{x}}_{j}})\right]\\ &+{{\delta}_{b}}\cdot \left[\frac{1}{F}\underset{l=1}{\overset{F}{\sum}}\,\widehat{\text{PIC}}({{\mathbf{x}}_{l}},{{\mathbf{x}}_{j}})\right],\\ &\text{where}\quad {{\widehat{r}}_{\text{s}}}({{\mathbf{x}}_{j}},\mathbf{y})=\frac{1}{B}\underset{b=1}{\overset{B} \sum}\,{{r}_{\text{s}}}(\mathbf{x}_{j}^{*b},\mathbf{y}),\\ &\text{and}\quad \widehat{\text{PIC}}({{\mathbf{x}}_{k}},{{\mathbf{x}}_{j}})=\frac{1}{B}\underset{b=1}{\overset{B} \sum}\,\text{PIC}(\mathbf{x}_{k}^{*b},\mathbf{x}_{j}^{*b}).\end{array}\end{eqnarray} \tag{ 4 }$$

In (4), r_s(x_j, y) is the Spearman's rank correlation computed between feature j defined as ${{\mathbf{x}}_{j}}=\left\{{{x}_{ij}}\in \mathbb{R}:i=1,2,\ldots ,N\right\}$ and the outcome vector y = {y_i ∈ {0: NoLungMets, 1: LungMets} : i = 1, 2, ..., N}. PIC(x_k, x_j) is the potential information coefficient defined as PIC(x_k, x_j) = 1 − MIC(x_k, x_j), where MIC(x_k, x_j) is the maximal information coefficient between feature j and k as defined by Reshef et al (2011). The sum over k is a sum over all f features that have already been chosen to be part of the reduced feature set (employed in forward selection schemes), whereas the sum over l is a sum over all F features that have not yet been removed from a larger initial set (employed in backward selection schemes). The sum over the k features is always done in order of appearance of the different features in the reduced set in order to favour the features from the larger initial set with the least dependence with the features chosen first in the reduced set. In this work, γ was set to 0.5, δ_a to 0.5 and δ_b to 0. Every time a new feature had to be chosen in the reduced set from a larger initial set, a new Gain was calculated for all remaining features in the larger initial set using imbalance-adjusted bootstrap resampling (1000 samples). Note that (4) allows to rank specific scan-texture-parameter features, as part 1 of the Gain equation uses Spearman's rank correlations varying over the whole set of texture extraction parameters. However, to speed up calculations, average scan-texture features over all texture extraction parameters were used in part 2 (and 3 if needed) of the Gain equation.

From the reduced feature sets, stepwise forward feature selection was then carried out by maximizing the 0.632+ bootstrap AUC. For a given model order and a given reduced feature set, the feature selection step was divided into 25 experiments. In each of these experiments, a different feature from the reduced set was used as a different 'starter'. For a given starter, 1000 logistic regression models g(x_i) or order 2 were first created using imbalance-adjusted bootstrap resampling (1000 samples) for each of the remaining features in the reduced feature set. Then, the single remaining feature that maximized the 0.632+ bootstrap AUC of the models was chosen, and the process was repeated up to model order 10. Finally, for each model order, the experiment that yielded the highest 0.632+ bootstrap AUC was identified, and combination of features were chosen for model orders of 1–10 (only features of the models were selected, but logistic regression coefficients were not yet computed).

The feature reduction and feature selection processes were repeated for all possible combinations of texture extraction parameter types being allowed to vary or not (i.e. degrees of freedom): 2⁴ combinations in the case of the initial feature set containing textures extracted from FDG-PET scans and the initial feature set containing textures extracted from separate scans, and 2⁶ combinations in the case of the initial feature set containing textures extracted from fused scans. For example, in an experiment where only MRI weight and Quant. algo. texture extraction parameters were allowed to vary, the FDG-PET initial feature set contained 9 non-texture features + 79 scan-texture-parameter features, the separate scan initial feature set contained 9 non-texture features + 237 scan-texture-parameter features, and the fused scan initial feature set contained 9 non-texture features + 790 scan-texture-parameter features. Then, separately for each model order of 1–10 for the three feature sets, the degree of freedom on texture extraction parameters yielding the multivariable models with the highest 0.632+ bootstrap AUC was found. For the experiments in which specific extraction parameters were not allowed to vary, baseline parameters had to be defined. Table 4 presents the baseline texture extraction parameters used for the five different types of scans.

Table 4. Baseline texture extraction parameters.

MRI Inv.	MRI weight	WBPF	Scale	Quant. algo.	Ng
No	$\frac{1}{2}$	R = 1	in-pR	Lloyd-Max	32

Note: MRI Inv.: MRI Inversion, MRI weight: weight in FDG-PET/MR fusion process, WBPF: wavelet band-pass filtering, Scale: isotropic voxel size, Quant. algo.: quantization algorithm, N_g: Number of gray-levels, R: ratio of the weight applied to band-pass sub-bands to the weight applied to low- and high-frequency sub-bands. in-pR: in-plane resolution.

Once optimal combination of features were identified for model orders of 1–10 for the three different types of feature sets, imbalance-adjusted bootstrap resampling (1000 samples) was again performed for all models. Prediction performance was then estimated using the 0.632+ bootstrap method in terms of AUC as defined in (3), and in terms of sensitivity and specificity metrics as defined in (5). Using the prediction estimates for the three initial feature sets, a single combination of features possessing the best parsimonious properties was then determined.

$$\begin{eqnarray}\begin{array}{*{20}{l}}&{{[\hat{\text{S}}]}_{0.632+}}=\frac{1}{B}\underset{b=1}{\overset{B} \sum}\,[(1-\alpha (b))\cdot \text{S}(\mathbf{X},\mathbf{X})+\alpha (b)\cdot \text{S}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))],\\ &\text{where}\qquad\quad \alpha (b)=\frac{0.632}{1-0.368\cdot R(b)},\\ &\text{and}\qquad\quad R(b)=\left\{\begin{array}{c} \frac{\text{S}(\mathbf{X},\mathbf{X})-\text{S}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))}{\text{S}(\mathbf{X},\mathbf{X})} & \ \text{if}\ \frac{\text{S}(\mathbf{X},\mathbf{X})}{\text{S}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))}>1, \\ 0 & \text{otherwise}, \\ \end{array} \right.\\ &\text{for}\qquad\quad \text{S}:\text{Sensitivity},\text{Specificity}.\end{array}\end{eqnarray} \tag{ 5 }$$

The last step in the construction of the final prediction model was to compute the coefficients of the optimal combination of features using imbalance-adjusted bootstrap resampling (1000 samples). Let the logistic regression coefficient of feature j computed in a bootstrap sample X*^b and modeling an outcome vector y be denoted as β_j(X*^b, y) for j = 0, 1, ..., p, where p is the multivariable model order and j = 0 refers to the offset of the model g(x_i). The computation of the different coefficient estimates of the final prediction model was then performed as follows:

$$\begin{eqnarray}&&\widehat{{{\beta}_{j}}}=\frac{1}{B}\underset{b=1}{\overset{B} \sum}\,{{\beta}_{j}}\left({{\mathbf{X}}^{*b}},\mathbf{y}\right),\quad \text{for}\;j=0,1,\ldots ,p.\end{eqnarray} \tag{ 6 }$$

Figure 5 summarizes the workflow of multivariable analysis.

Zoom In Zoom Out Reset image size

Table 5 presents the Spearman's rank correlation (r_s) between the nine non-texture features and lung metastases development in STSs. Table 6 presents the Spearman's rank correlation between the 205 different scan-texture features and lung metastases development in STSs. For each entry in table 6, two values appear: the r_s in the case of texture features extracted using baseline parameters as defined in table 4 (left), and the maximal r_s in the case of texture features extracted using the optimal set of extraction parameters when all parameters are allowed to vary, i.e. with full degrees of freedom (right). In tables 5 and 6, the values in italic font indicate features for which p < α/K < 0.05/(9 non-texture features + 5 scans × 41 textures × 2 extraction parameter degrees of freedom) ≈ 0.000 12 according to the correction for multiple testing comparison.

Table 6. Spearmans rank correlation between texture features and lung metastases development in STSs using baseline | optimal texture extraction parameters.

Type	Texture	FDG-PET	T1	T2FS	FDG-PET/T1	FDG-PET/T2FS
Global	Variance	0.12 |   0.13	0.21 |   0.21	0.03 |   0.05	− 0.06 | −0.31	− 0.44 | −0.49
	Skewness	0.23 |   0.25	− 0.36 | −0.36	0.16 |   0.17	0.13 |   0.28	0.28 |   0.39
	Kurtosis	− 0.02 |   0.06	− 0.28 | −0.31	− 0.15 | −0.15	− 0.11 | −0.24	0.33 |   0.42
GLCM	Energy	− 0.01 |   0.49	− 0.22 | −0.44	0.14 |   0.23	− 0.04 |  $-{\it 0.51}$	0.37 |    $\it 0.53$
	Contrast	− 0.14 | −0.44	− 0.08 | −0.33	− 0.15 |   0.33	− 0.19 |  $-{\it 0.52}$	− 0.36 |  $-{\it 0.51}$
	Entropy	0.01 | −0.43	0.12 |   0.41	− 0.15 |   0.22	− 0.02 | −0.46	− 0.38 |  $-{\it 0.55}$
	Homogeneity	0.28 |   0.49	0.04 |   0.26	0.18 |   0.26	0.23 |    $\it 0.53$	0.42 |    $\it 0.51$
	Correlation	0.28 |   0.42	0.18 |   0.32	0.15 |   0.25	0.26 |   0.42	0.10 |   0.38
	Sum Average	− 0.35 |  $-{\it 0.52}$	0.28 |    $\it 0.51$	− 0.11 | −0.29	− 0.18 |  $-{\it 0.52}$	− 0.27 |  $-{\it 0.52}$
	Variance	0.31 | −0.43	0.32 |   0.49	− 0.07 | −0.26	− 0.20 | −0.50	− 0.32 |  $-{\it 0.59}$
GLRLM	SRE	− 0.33 |  $-{\it 0.53}$	− 0.04 | −0.31	− 0.18 | −0.34	− 0.25 |  $-{\it 0.54}$	− 0.41 |  $-{\it 0.53}$
	LRE	0.34 |    $\it 0.53$	0.02 |   0.34	0.20 |   0.33	0.25 |    $\it 0.53$	0.36 |    $\it 0.53$
	GLN	0.15 |   0.32	0.25 |   0.29	0.25 |   0.30	0.12 |   0.33	0.18 |   0.33
	RLN	0.13 |   0.32	0.26 |   0.30	0.22 |   0.30	0.14 |   0.33	0.14 |   0.33
	RP	− 0.34 |  $-{\it 0.54}$	− 0.02 | −0.33	− 0.19 | −0.34	− 0.25 |  $-{\it 0.54}$	− 0.39 |  $-{\it 0.52}$
	LGRE	0.29 | −0.48	0.06 | −0.37	− 0.08 | −0.33	0.11 | −0.50	− 0.16 |  $-{\it 0.58}$
	HGRE	− 0.23 | −0.32	0.23 | −0.47	− 0.12 | −0.28	− 0.12 |   0.48	− 0.21 |   0.43
	SRLGE	0.26 | −0.49	0.08 | −0.44	− 0.10 | −0.38	0.09 | −0.48	− 0.17 |  $-{\it 0.58}$
	SRHGE	− 0.25 | −0.35	0.10 | −0.44	− 0.17 | −0.31	− 0.14 | −0.40	− 0.25 | −0.44
	LRLGE	0.34 |    $\it 0.51$	0.09 |   0.37	0.08 |   0.40	0.10 |   0.50	− 0.11 |    $\it 0.54$
	LRHGE	− 0.21 |    $\it 0.55$	0.27 |   0.50	0.16 |   0.36	− 0.11 |    $\it 0.52$	− 0.12 |    $\it 0.55$
	GLV	− 0.33 |  $-{\it 0.52}$	− 0.10 | −0.41	− 0.28 | −0.39	− 0.25 |  $-{\it 0.53}$	0.26 |  $-{\it 0.55}$
	RLV	− 0.31 |  $-{\it 0.57}$	− 0.18 | −0.38	− 0.27 | −0.40	− 0.27 | $-{\it 0.56}$	− 0.35 |  $-{\it 0.55}$
GLSZM	SZE	− 0.18 | −0.42	− 0.05 |  $-{\it 0.51}$	− 0.05 | −0.35	− 0.41 |  $-{\it 0.63}$	− 0.41 |  $-{\it 0.60}$
	LZE	0.19 |   0.50	0.23 |   0.33	0.21 |   0.36	0.19 |    $\it 0.53$	0.30 |    $\it 0.52$
	GLN	0.15 |   0.33	0.20 |   0.31	0.21 |   0.29	0.10 |   0.32	0.13 |   0.33
	ZSN	0.15 |   0.32	0.22 |   0.31	0.18 |   0.29	0.06 |   0.31	0.04 |   0.34
	ZP	− 0.20 |  $-{\it 0.51}$	− 0.17 | −0.41	− 0.15 | −0.34	− 0.25 |  $-{\it 0.58}$	− 0.40 |  $-{\it 0.56}$
	LGZE	− 0.09 | −0.48	0.31 |   0.36	− 0.01 |   0.36	0.12 | −0.46	− 0.16 |  $-{\it 0.55}$
	HGZE	− 0.12 |   0.38	− 0.19 | −0.40	− 0.04 | −0.29	− 0.17 |   0.50	− 0.06 |    $\it 0.52$
	SZLGE	− 0.22 | −0.46	0.28 | −0.39	− 0.01 |   0.34	0.10 | −0.47	− 0.19 |  $-{\it 0.58}$
	SZHGE	− 0.07 |   0.27	− 0.17 |   0.47	− 0.04 |   0.29	− 0.31 |    $\it 0.52$	− 0.21 |   0.50
	LZLGE	0.32 |    $\it 0.52$	0.20 |   0.34	0.21 |   0.35	0.25 |    $\it 0.55$	0.35 |    $\it 0.55$
	LZHGE	0.01 |   0.47	0.27 |   0.44	0.20 |   0.34	0.12 |   0.47	0.22 |   0.46
	GLV	− 0.18 | −0.45	− 0.25 | −0.34	− 0.26 | −0.33	− 0.21 |  $-{\it 0.53}$	− 0.24 |  $-{\it 0.53}$
	ZSV	− 0.21 | −0.48	− 0.07 | −0.43	− 0.05 | −0.33	− 0.09 |    $\it 0.53$	− 0.01 | −0.48
NGTDM	Coarseness	− 0.06 | −0.26	− 0.22 | −0.28	− 0.21 | −0.30	− 0.08 | −0.29	− 0.13 | −0.34
	Contrast	− 0.14 | −0.39	0.16 |   0.36	0.02 |   0.33	− 0.11 | −0.46	− 0.33 |  $-{\it 0.51}$
	Busyness	0.23 |   0.39	0.22 |   0.28	0.22 |   0.30	0.20 |   0.39	0.17 |   0.39
	Complexity	0.21 |  $-{\it 0.55}$	− 0.16 |   0.39	− 0.13 |   0.40	0.14 |    $\it 0.52$	0.22 | −0.48
	Strength	0.04 | −0.25	− 0.24 | −0.29	− 0.21 | −0.29	− 0.09 | −0.37	− 0.07 | −0.33
Average absolute values		0.20 |   0.42	0.18 |   0.37	0.15 |   0.31	0.16 |   0.46	0.24 |   0.49

Note: The values in italic font indicate features for which p < 0.05/419 ≈ 0.000 12.

3.1.1. Non-texture features.

3.1.2. Texture features.

3.1.3. Texture extraction parameter effect.

The Wilcoxon rank sum test was performed between the set of absolute Spearman's rank correlation coefficients obtained from each of the 41 different texture features extracted using baseline extraction parameters, and the sets of maximal absolute Spearman's rank correlation coefficients obtained from optimal texture-parameters for each of the 41 different textures when one extraction parameter was allowed to vary, and the others set to baseline. The same process was repeated for all parameters of all scans, and also for a Wilcoxon rank sum test comparing baseline extraction parameters to full degrees of freedom on extraction parameters (ALL PARAMs.). Table 7 presents the p-value of the Wilcoxon rank sum tests, with multiple testing corrections (α = 0.05, K = 29). The results point out that the optimization of the Scale texture extraction parameter has the highest impact on the predictive value for lung metastases development in STSs. In general, each extraction parameter seems to positively impact the predictive value of textures.

Table 7. p-value of the Wilcoxon rank sum tests asserting the significance of the effects of texture extraction parameters on the correlation of texture features with lung metastases development in STSs.

Scan	MRI Inv.	MRI weight	WBPF	Scale	Quant. algo	Ng	ALL PARAMs.
FDG-PET	—	—	0.0727	$< \it 0.0010$	0.0017	0.0260	$< \it 0.0010$
T1	—	—	0.2406	$< \it 0.0010$	0.1159	0.0369	$\it <0.0010$
T2FS	—	—	0.2105	0.0013	0.0066	0.0352	$< \it 0.0010$
FDG-PET/T1	0.1085	0.0019	0.0034	$< \it 0.0010$	0.0024	0.0143	$< \it 0.0010$
FDG-PET/T2FS	0.5937	0.1259	0.3076	0.0024	0.4143	0.3909	$< \it 0.0010$

Note: The values in italic font indicate features for which p < 0.05/29 ≈ 0.0017. MRI Inv.: MRI Inversion, MRI weight: weight in FDG-PET/MRI fusion process, WBPF: wavelet band-pass filtering, Scale: isotropic voxel size, Quant. algo.: quantization algorithm, N_g: Number of gray-levels, ALL PARAMs.: all texture extraction parameters allowed to vary.

We compared the prediction performance estimation of multivariable models constructed using three different types of initial feature sets: (i) FDG-PET textures + non-texture features; (ii) separate FDG-PET and MRI textures + non-texture features; and (iii) fused FDG-PET/MRI textures + non-texture features. We performed experiments for all degrees of freedom on texture extraction parameters. Figure 6 presents the prediction performance estimation of multivariable models with optimal degrees of freedom on texture extraction parameters, obtained separately for each model order of each initial feature set. Supplementary appendix C (stacks.iop.org/PMB/60/145471/mmedia) also provides detailed comparison between prediction estimates obtained in the experiments using baseline, full and optimal degrees of freedom on texture extraction parameters. Results show that multivariable models constructed with texture features extracted from separate scans provide no significant prediction estimation improvements as compared to multivariables models constructed with texture features extracted from FDG-PET scans only. On the other hand, multivariable models constructed from fused scans significantly improve the prediction performance estimation compared to FDG-PET scans alone.

Zoom In Zoom Out Reset image size

By inspecting the curves in figure 6, we determined that the simplest multivariable model with best predictive properties (best parsimonious model) is obtained by linearly combining 4 texture features extracted from fused FDG-PET/MRI scans. The associated prediction performance estimation of this optimal combination of features using 1000 bootstrap samples yielded an AUC of 0.984 ± 0.002, a sensitivity of 0.955 ± 0.006 and a specificity of 0.926 ± 0.004. These last results were obtained using the 0.632+ bootstrap method, and as a comparison, the same model reached an AUC of 0.976 ± 0.002, a sensitivity of 0.938 ± 0.008 and a specificity of 0.892 ± 0.006 using the ordinary bootstrap method $\left(\widehat{\text{AUC}}=\frac{1}{B}\sum\nolimits_{b=1}^{B}\text{AUC}({{\mathbf{X}}^{*b}},{{\mathbf{X}}^{*b}}(0))\right)$. Next, the logistic regression coefficients of the final prediction model were computed using 1000 bootstrap samples. We hence propose the following complete multivariable model response g(x_i) to be computed from fused FDG-PET/MRI scans at the time of diagnosis of STSs for the prediction of future lung metastases development:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}g({\bf x}_i) = -256 \times \scriptsize \mathsf{FDG-PET/T2FS(\tiny{\it MRI~Inv.} = No Inv., {\it MRI~weight} = 1/2, {\it R} = 3/2, {\it Scale} = 3\,mm, {\it Quant.~algo.} = Lloyd-Max, {\it Ng} = 64\scriptsize) - - ~GLSZM/SZE}\\ \quad +\\ \quad 5360 \times \scriptsize \mathsf{FDG-PET/T1(\tiny{\it MRI~Inv.} = Inv., {\it MRI~weight} = 1/2, {\it R} = 1/2, {\it Scale} = {\it in-pR}, {\it Quant.~algo.} = Lloyd-Max, {\it Ng} = 16\scriptsize) - -~ GLSZM/ZSV}\\ \quad +\\ \quad 1.75 \times \scriptsize \mathsf{FDG-PET/T1(\tiny{\it MRI~Inv.} = Inv., {\it MRI~weight} = 3/4, {\it R} = 1, {\it Scale} = 2\,mm, {\it Quant.~algo.} = Lloyd-Max, {\it Ng} = 8\scriptsize) - - ~ GLSZM/HGZE}\\ \quad + 3.16 \times \scriptsize \mathsf{FDG-PET/T2FS(\tiny{\it MRI~Inv.} = Inv., {\it MRI~weight} = 3/4, {\it R} = 2, {\it Scale} = 1\,mm, {\it Quant.~algo.} = Equal, {\it Ng} = 8\scriptsize) - -~ GLRLM/HGRE}\\ \quad + 26.7\end{array}\end{eqnarray} \tag{ 7 }$$

In order to evaluate the precision of the proposed model, we calculated how its response changes using texture features extracted from tumour contours that include surrounding edema. Supplementary appendix D.1 (stacks.iop.org/PMB/60/145471/mmedia) details the calculations. Overall, an absolute value of ± 4.89 was estimated as the uncertainty of the model due to contouring variations. This uncertainty is constant across all values of g(x_i). Then, to summarize how the model can separate the instances of the two patient classes (LungMets versus NoLungMets), the vector $\mathbf{g}=\left\{g\left({{\mathbf{x}}_{i}}\right)\in \mathbb{R}:i=1,2,\ldots ,N\right\}$ was computed for all patients using the multivariable model response of (7) and was transformed into the posterior probability π(x_i) of observing outcome y_i = 1 (i.e. developing lung metastases) given the input x_i by using the logit transform of (2). Figure 7 displays the plot of π(x_i) versus g(x_i), along with the associated 95% confidence intervals (CIs) on g(x_i) for i = 1, 2,..., N. For each bootstrap sample b used to calculate the final logistic regression coefficients of (7), a new value of $g(\mathbf{x}_{i}^{*b})$ was calculated for i = 1, 2,..., N from the new coefficients computed on $\mathbf{x}_{i}^{*b}$. Then, the lower and upper CI bounds were estimated for each point i by calculating the 2.5 and 97.5 percentiles from the bootstrap distribution of $g(\mathbf{x}_{i}^{*b})$ for b = 1, 2, ..., B. In figure 7, the dots represent patients who eventually developed lung metastases, and the crosses those who did not develop lung metastases. The uncertainty due to contouring variations around the classification threshold g(x_i) = 0 is also shown, and supplementary appendix D.2 (stacks.iop.org/PMB/60/145471/mmedia) provides the detailed data (lung mets status, g(x_i) and CIs) used to construct the figure. It can be seen that the multivariable model of (7) can clearly separate the patients of the two risk groups. Note that the Spearman's rank correlation between the model response vector g and the outcome vector y reached r_s = 0.84, p < 0.001.

Zoom In Zoom Out Reset image size

Finally, further validation of the proposed model was performed using permutation tests (Fisher 1935). This operation was carried out by randomly shuffling the real outcome vector y of our patient cohort (i.e. keeping the same proportion of 0's and 1's). In order to have a direct comparison with the model of (7), only the prediction performance of multivariable models of order 4 were analyzed. For each of the 1000 permutation tests, a different multivariable model was constructed and its prediction performance was assessed. Table 8 presents a summary of the results of the permutation tests, where the p-values were estimated using the Monte Carlo sampling approach described by Ernst (2004). Supplementary appendix E (stacks.iop.org/PMB/60/145471/mmedia) also provides a display of the permutation distributions. Overall, results in table 8 show that the null hypothesis can be rejected. Very few tests yielded prediction estimates higher than the true observed values. The results give strong evidence that the effect observed on the sample data is most likely present in the general STS population.

Table 8. Summary of permutation tests comparing the single value of the performance metrics estimates of the multivariable model proposed in this work found using the real outcome vector $(\widehat{\text{TRUE}})$ to the distribution of the performance metrics estimates of models of order 4 found in 1000 permutation tests using shuffled outcome vectors $(\widehat{\text{PERM}})$. SD: Standard deviation.

Metric	$\text{Valu}{{\text{e}}_{\widehat{\text{TRUE}}}}$	$\text{Mea}{{\text{n}}_{\widehat{\text{PERM}}}}$	$\text{S}{{\text{D}}_{\widehat{\text{PERM}}}}$	$\text{Rang}{{\text{e}}_{\widehat{\text{PERM}}}}$	$\widehat{\text{PERM}}>\text{Valu}{{\text{e}}_{\widehat{\text{TRUE}}}}$	$\widehat{p}$
AUC	0.984	0.895	0.04	[0.745, 0.988]	3 out of 1000	0.004
Sensitivity	0.955	0.798	0.05	[0.634, 0.938]	0 out of 1000	0.001
Specificity	0.926	0.812	0.05	[0.666, 0.948]	6 out of 1000	0.007

Clinical information and imaging data analyzed in this work are available on The Cancer Imaging Archive (TCIA) website under the following DOI: http://dx.doi.org/10.7937/K9/TCIA.2015.7GO2GSKS. All software code implemented in this work is freely shared under the GNU General Public License at: https://github.com/mvallieres/radiomics.