Features module tutorial

This tutorial shows how to perform an exploratory analysis of the features data: to visualize features distribution in classes, plot the feature correlation matrix, check Mann-Whitney U-test p-values, plot univariate ROC (and calculate AUC) for each feature, perform volumetric analysis, and save all the scores. We will use PyRadiomics features as variables and binary survival label as an outcome.

Importing modules:

import os,sys
from pmtool.AnalysisBox import AnalysisBox
import pandas as pd

Required data

To run the tutorial, you will need the following files, available in the 'data/features' folder of the precision-medicine-toolbox repository:
- 'extended_clinical_df.xlsx' - this is the clinical data file provided with the Lung1 dataset (Aerts et al., 2014); we added binary variables of 1-, 1.5-, and 2-years survival, based on the presented data,
- 'extracted_features_full.xlsx' - this is a file with the PyRadiomic features values (can be extracted in the imaging module tutorial) for all the patients from Lung1 dataset.

Binary classes functionality

If you have 2 classes in your dataset, the following functionality is availabe for your analysis. We will divide Lung1 dataset by 1 year survival label, therefore having 2 classes.

Set up the parameters to get the data:

parameters = {
    'feature_path': "../data/features/extracted_features_full.xlsx", # path to csv/xls file with features
    'outcome_path': "../data/features/extended_clinical_df.xlsx", #path to csv/xls file with outcome
    'patient_column': 'Patient', # name of column with patient ID
    'patient_in_outcome_column': 'PatientID', # name of column with patient ID in clinical data file
    'outcome_column': '1yearsurvival' # name of outcome column
}

Initialise the feature set (you will see a short summary):

fs = AnalysisBox(**parameters)

Number of observations: 149
Class labels: ['0' '1']
Classes balance: [0.4228187919463087, 0.5771812080536913]

Print some attributes of the feature set - first 10 patient IDs and first 10 feature names:

print ('Patient IDs: ', fs._patient_name[:10])
print ('\nFeature names: ', fs._feature_column[:10])

Patient IDs:  ['LUNG1-001_20180209_CT_2_GTV-1_mask', 'LUNG1-002_000000_GTV-1_mask', 'LUNG1-002_20180526_CT_1_GTV-1_mask', 'LUNG1-003_000000_GTV-1_mask', 'LUNG1-003_20180209_CT_1_GTV-1_mask', 'LUNG1-003_20180209_CT_1_GTV-2_mask', 'LUNG1-003_20180209_CT_1_GTV-3_mask', 'LUNG1-004_000000_GTV-1_mask', 'LUNG1-006_000000_GTV-1_mask', 'LUNG1-008_000000_GTV-1_mask']

Feature names:  ['original_shape_Elongation', 'original_shape_Flatness', 'original_shape_LeastAxisLength', 'original_shape_MajorAxisLength', 'original_shape_Maximum2DDiameterColumn', 'original_shape_Maximum2DDiameterRow', 'original_shape_Maximum2DDiameterSlice', 'original_shape_Maximum3DDiameter', 'original_shape_MeshVolume', 'original_shape_MinorAxisLength']

Exclude patients with the unknown outcomes, if they are represented in the dataset (feature set parameters are re-initialised and short summary is printed again):

fs.handle_nan(axis=0)

Number of observations: 149
Class labels: ['0' '1']
Classes balance: [0.4228187919463087, 0.5771812080536913]

Print the head of the composed dataframe, containing both the variables and the outcome:

fs._feature_outcome_dataframe.head(5)

Visualize feature values distribution in classes for the first 12 features (will pop up in an interactive .html report):

fs.plot_distribution(fs._feature_column[:8])

Shapiro-Wilcoxon test to get a list of normally distributed features:

features_norm_distr = fs.normality_check()
print (len(features_norm_distr))
print (features_norm_distr)

59
['original_shape_Sphericity', 'original_glcm_Idm', 'original_glcm_Id', 'original_glrlm_RunLengthNonUniformityNormalized', 'original_glszm_HighGrayLevelZoneEmphasis', 'original_glszm_SmallAreaHighGrayLevelEmphasis', 'original_gldm_SmallDependenceHighGrayLevelEmphasis', 'log-sigma-1-0-mm-3D_glcm_DifferenceEntropy', 'log-sigma-2-0-mm-3D_firstorder_MeanAbsoluteDeviation', 'log-sigma-3-0-mm-3D_firstorder_InterquartileRange', 'log-sigma-3-0-mm-3D_firstorder_RobustMeanAbsoluteDeviation', 'log-sigma-3-0-mm-3D_glcm_Id', 'log-sigma-4-0-mm-3D_firstorder_InterquartileRange', 'log-sigma-4-0-mm-3D_glcm_Idm', 'log-sigma-4-0-mm-3D_glcm_Id', 'log-sigma-4-0-mm-3D_glszm_SizeZoneNonUniformityNormalized', 'log-sigma-5-0-mm-3D_firstorder_Skewness', 'log-sigma-5-0-mm-3D_glcm_Imc1', 'log-sigma-5-0-mm-3D_glcm_Idm', 'log-sigma-5-0-mm-3D_glcm_Id', 'wavelet-LLH_firstorder_Entropy', 'wavelet-LLH_glcm_SumEntropy', 'wavelet-LHL_glcm_DifferenceAverage', 'wavelet-LHH_firstorder_10Percentile', 'wavelet-LHH_firstorder_MeanAbsoluteDeviation', 'wavelet-LHH_glcm_Correlation', 'wavelet-LHH_glcm_DifferenceAverage', 'wavelet-LHH_glcm_Idm', 'wavelet-LHH_glcm_Id', 'wavelet-LHH_glrlm_RunVariance', 'wavelet-HLL_glcm_DifferenceEntropy', 'wavelet-HLL_glcm_Imc1', 'wavelet-HLL_glcm_Idm', 'wavelet-HLL_glcm_Id', 'wavelet-HLL_glrlm_RunLengthNonUniformityNormalized', 'wavelet-HLL_glrlm_RunPercentage', 'wavelet-HLH_glcm_Correlation', 'wavelet-HLH_glcm_Idm', 'wavelet-HLH_glcm_Id', 'wavelet-HLH_glcm_InverseVariance', 'wavelet-HHL_firstorder_Entropy', 'wavelet-HHL_firstorder_MeanAbsoluteDeviation', 'wavelet-HHL_firstorder_Uniformity', 'wavelet-HHL_glcm_Correlation', 'wavelet-HHL_glcm_DifferenceEntropy', 'wavelet-HHL_glcm_JointEntropy', 'wavelet-HHL_glcm_Imc2', 'wavelet-HHL_glcm_MaximumProbability', 'wavelet-HHL_glcm_SumEntropy', 'wavelet-HHL_glrlm_ShortRunEmphasis', 'wavelet-HHH_firstorder_10Percentile', 'wavelet-HHH_firstorder_90Percentile', 'wavelet-HHH_firstorder_InterquartileRange', 'wavelet-HHH_firstorder_MeanAbsoluteDeviation', 'wavelet-HHH_firstorder_RobustMeanAbsoluteDeviation', 'wavelet-HHH_glcm_Imc2', 'wavelet-HHH_glcm_MaximumProbability', 'wavelet-LLL_glszm_HighGrayLevelZoneEmphasis', 'wavelet-LLL_gldm_SmallDependenceHighGrayLevelEmphasis']

Visualize mutual feature correlation coefficient (Spearman's) matrix for the first 12 features (in .html report):

fs.plot_correlation_matrix(fs._feature_column[:8])

Visualize Mann-Whitney (Bonferroni corrected) p-values for binary classes test (in .html report):

fs.plot_MW_p(fs._feature_column[:8])

Vizualize univariate ROC curves:

fs.plot_univariate_roc(fs._feature_column[:8], auc_threshold=0.70)

Calculate the basic statistics for each feature (save the values in data/features/extracted_features_full_basic_stats.xlsx): number of NaN, mean, std, min, max; if applicable: MW-p, univariate ROC AUC, volume correlation:

fs.calculate_basic_stats(volume_feature='original_shape_VoxelVolume')

Check the excel table:

print('Basic statistics for each feature')
pd.read_excel('../data/features/extracted_features_full_basic_stats.xlsx')

Basic statistics for each feature

Volume analysis will show you if your features have a high correlation to volume and if volume itself is a predictive feature in separation of 2 classes. You need to have a volume feature in your dataset and send it as a function parameter (in our case it is 'original_shape_VoxelVolume').

fs.volume_analysis(volume_feature='original_shape_VoxelVolume')

Multi-class

If you have more than 2 classes in your dataset, the following functionality is availabe for your analysis. We will use a disease stage as a class label.

Set up the parameters to get the data:

parameters = {
    'feature_path': "../data/features/extracted_features_full.xlsx", # path to csv/xls file with features
    'outcome_path': "../data/features/extended_clinical_df.xlsx", # path to csv/xls file with outcome
    'patient_column': 'Patient', # name of column with patient ID
    'patient_in_outcome_column': 'PatientID', # name of column with patient ID in clinical data file
    'outcome_column': 'Overall.Stage' # name of outcome column
}

Initialise the feature set:

fs = AnalysisBox(**parameters)

Number of observations: 149
Class labels: ['I' 'II' 'IIIa' 'IIIb' 'nan']
Classes balance: [0.24161073825503357, 0.09395973154362416, 0.2348993288590604, 0.4228187919463087, 0.006711409395973154]

Check the patients with the absent outcomes (because there is a NaN value among class labels in the data summary):

fs._feature_outcome_dataframe[fs._feature_outcome_dataframe['Overall.Stage'].isnull()]

Exclude the patients with the unknown outcomes:

fs.handle_nan(axis=0)

Number of observations: 148
Class labels: ['I' 'II' 'IIIa' 'IIIb']
Classes balance: [0.24324324324324326, 0.0945945945945946, 0.23648648648648649, 0.42567567567567566]

Visualize feature values distribution in classes for the first 12 features (in .html report):

fs.plot_distribution(fs._feature_column[:8])

Visualize feature values distribution in classes for the first 12 features only for the selected classes (in .html report):

fs.plot_distribution(fs._feature_column[:8], ['I', 'IIIb'])

Visualize mutual feature correlation coefficient (Spearman's) matrix for the first 12 features (in .html report):

fs.plot_correlation_matrix(fs._feature_column[:12])

Even though we have more than 2 classes and cannot perform Mann-Whitney test for all the classes, we can perform it for the pairs of the classes (result in .html report):

fs.plot_MW_p(fs._feature_column[:100], ['I', 'IIIb'])

The same situation is with univariate analysis - we can perform it for the selected pairs of the classes:

fs.plot_univariate_roc(fs._feature_column[:100], ['I', 'IIIb'])

Calculate the basic statistics for each feature (save the values in data/features/features_basic_stats.xlsx); for multi-class analysis the following are available: number of NaN, mean, std, min, max

fs.calculate_basic_stats(volume_feature='original_shape_VoxelVolume')

Check the table:

print('Basic statistics for each feature')
pd.read_excel('../data/features/extracted_features_full_basic_stats.xlsx')

Basic statistics for each feature

Volume analysis (result in .html report):

fs.volume_analysis(volume_feature='original_shape_VoxelVolume')