Many methods have been proposed to classify gliomas using MRI images through deep learning methods based on radiological characteristics (Mohan and Subashini, 2018).

Some other methods have also been proposed to automatically classify gliomas based on histological features by deep learning methods, which will greatly improve diagnostic efficiency and improve patient outcomes (Jothi and Rajam, 2017; Mobadersany et al., 2018; Jin et al., 2020).

The MICCAI 2020 CPM-RadPath Challenge provides publicly available H&E stained digital histopathology images and matched multi-sequence radiology images, including three subtypes of gliomas: astrocytoma, glioblastoma, and oligodendroglioma.

Each image in the WSI dataset is very large (e.g., 95,200 × 87,000 pixels) and cannot be directly fed into our network. We crop these WSIs into patches with a size of 1,024 × 1,024 pixels. Then, the OTSU method is adopted to remove non-tissue regions (Otsu, 1979). Following the current studies (Ma and Jia, 2019; Pei et al., 2019; Pei et al., 2020; Pei et al., 2021), we exclude meaningless tissues using a simple but effective threshold technique. Specifically, we first calculate the mean value and standard deviation of each patch in RGB space and maintain patches with a mean value between 100 and 220 and standard deviations above 20 (Pei et al., 2019; Pei et al., 2020; Pei et al., 2021). Then, we convert each patch to the hue saturation value (HSV) space and exclude patches with the mean value below 50 in the H channel (Ma and Jia, 2019).

The volume of the original MRI image is 240 × 240 × 155 pixels. The beginning and end images (slices) in the scan are removed due to their limited brain tissue. As a result, the number of MRI images per sequence is reduced to 128. We also remove the black background at the edges. These MRI images in the four sequences are finally cropped into small images of size 192 × 192 × 128 pixels, which facilitates computational efficiency.

Figure 2 illustrates the overall architecture of our proposed glioma classification system, which is composed of a 2D WSI model (Figure 3) and a 3D MRI model (Figure 4). Both models adopt our two-stage strategy. The first stage classifies MRI images or histopathological images as glioblastoma and non-glioblastoma. Based on the obtained non-glioblastoma data, the second stage aims to distinguish astrocytoma and oligodendroglioma. The reason for using the two-stage strategy is the uneven data distribution. The differences between astrocytoma and oligodendroglioma are so subtle that it can be difficult to distinguish between the two. In contrast, glioblastoma is easier to identify. As confirmed by our experimental results, our two-stage strategy helps to improve the overall accuracy.

Multiple related tasks can help learn from each other by potentially sharing representations, thus improving generalization ability. Therefore, for the 2D WSI model, we develop several multi-task-based convolutional neural network models, where the backbones include EfficientNet-B2, EfficientNet-B3 (Tan and Le, 2019), and SEResNeXt101 (Hu et al., 2018). Specifically, EfficientNet is a benchmark network that achieves performance gains by scaling network width, network depth, and resolution, which greatly reduces the number of parameters and computation complexity of the model. EfficientNet-B1 to B7 are obtained by synthetically optimizing the width, depth, and resolution of the EfficientNet. SEResNeXt101 is derived by embedding squeeze-and-excitation (SE) blocks into the ResNeXt model. ResNeXt replaces the original residual learning block (He et al., 2016) with parallel blocks of the same topology, which improves the model performance without significantly increasing the parameters. The SE block obtains the weight of each feature channel and assigns more weights to important features while suppressing features that are not useful for the current task. The outputs of these three classification models are averaged to obtain the final classification results.

Before the fully connected classification layer, the generalized-mean (GEM) pooling (Radenović et al., 2018) is applied to the learned features, which is defined as f ( g ) = [ f 1 ( g ) … f k ( g ) … f K ( g ) ] T ,   w h e r e ,   f k ( g ) = ( 1 | X k | ∑ x ∈ X k x p k ) 1 p k (1) where X is the input features and f ( g ) represents the output vector, respectively. p k ∈ [1, +∞), when p k = 1, the equation denotes the average pooling, when p k tends to positive infinity, the equation represents the max pooling, respectively. K denotes the channel of the output feature vector f ( g ) . Following (Radenović et al., 2018), this work set p k to 3.

After the GEM operation, age information of patient should also be taken into account (Reni et al., 2017; Ostrom et al., 2018). The reason is that glioma subtypes have different age distributions, where astrocytoma is more often seen in young men while glioblastoma is more frequently distributed in the older groups (Tamimi and Juweid, 2017). In our work, age is used as additional feature, which is concatenated with the extracted image features for the final classification task. Finally, these combined features are adopted as input to a classification branch and a regression branch. The two branches adopt a cross-entropy (L_BCE) and a smooth L1 loss (L_loc), respectively, to achieve a more robust brain tumor classification. L B C E = − ∑ l [ ( y l ⁡ log y l ^ ) + ( 1 − y l ) log ( 1 − y l ^ ) ] (2) s m o o t h L 1 = { 0.5 ( y l ^ − y l ) 2                     i f   | y l ^ − y l | < 1 , | y l ^ − y l | − 0.5                                   o t h e r w i s e .   (3)   L l o c = ∑ l [ s m o o t h L 1 ( y l ^ − y l ) ] (4) where y l and y l ^ denote the ground-truth and predicted labels for the l ^th sample, respectively.

For the 3D MRI model, the input images are four types of MRI sequence images, including T1, T1-Gd, T2, and FLAIR. Each 3D MRI image is cropped to 192 × 192 × 128 pixels. 3D ResNet is adopted as a backbone to learn the residual representation between the input and output, which has become the basic feature extraction network in the computer vision community. Then, global average pooling and fully connected layers are used to generate the final classification results. We also employ the cross-entropy and smooth L1 loss functions. Finally, the output probabilities of the 2D WSI and 3D MRI models are averaged to drive the final classification results.

We utilize a public dataset released from the MICCAI 2020 CPM-RadPath Challenge ¹ , which is proposed for three subtypes of glioma classification, including Glioblastoma, Oligodendroglioma, and Astrocytoma. It includes paired radiology scans and digitized histopathology images with global image-level labels. The CPM-RadPath Challenge splits these data into three parts: training (270 cases), validation (35 cases), and testing (73 cases).

It is noted that only the annotations of training data are released. Thus, our model is developed depending on the training set. The validation and test sets have not released their labels, which can be regarded as two unseen test sets. All these data are collected using multi-parametric MRI (mp-MRI) and digital pathology scanners at 16 international institutions.

Specifically, all mp-MRI scans are acquired using 1–3T scanners. The MRI scans are provided as the NIFTI files (. nii.gz). Each case includes four sequences: T1, T1-Gd, T2, and FLAIR. All MRI images have been preprocessed, co-registered to the same anatomical template, and interpolated to the same resolution (1 cubic mm) in all three directions. Annotations are generated by board-certified neuroradiologists, neurosurgeons, and neuropathologists with at least 4 years of experience.

Meanwhile, tissue specimens are made from tissues removed from the patient during surgery and then stained with hematoxylin and eosin (H&E), which are scanned at 20× or 40× magnification to generate digital histopathology images called WSIs. The color and intensity of WSIs varied across images due to different acquisition times, image fading, or image acquisition artifacts. All WSIs are stored in tiled tiff format. Sixteen professional neuropathologists with at least 4 years of experience annotate these cases by referring to the 2016 WHO classification scheme.

The algorithmic performance is evaluated using three metrics, including F1-score, balanced accuracy, and Cohen’s kappa. The F1 score is a weighted average of accuracy and recall, which takes into account both false positives and false negatives. The balanced accuracy score is a more appropriate metric to evaluate data with imbalanced categories, which is defined as the arithmetic mean of the proportion of correct predictions in each category. Cohen’s kappa is used for consistency testing and can also be used to measure classification accuracy. A higher kappa coefficient means that the classifier is more effective.

It is assumed that true positive (TP), false positive (FP), false negative (FN) is the number of true positives, false positives, and false negatives respectively. The three metrics can be computed as follows. S e n s i t i v i t y = T P T P + F N (5) P r e c i s i o n = T P T P + F P (6) F 1 _ S c o r e = 2 ∗ ( P r e c i s i o n ∗ S e n s i t i v i t y ) P r e c i s i o n + S e n s i t i v i t y (7) B a l a n c e d   A c c u r a c y = ∑ C l a s s = 1 c S e n s i t i v i t y / c (8) K a p p a = p 0 − p e 1 − p e (9)

In the balanced accuracy score, c is the number of classes. In Cohen’s kappa, p 0 is obtained by dividing the sum of the number of correctly classified samples by the total number of samples. It is assumed that the true number of samples for each category is a 1 , a 2 ,..., a c , and the predicted number of samples for each category is b 1 , b 2 ,..., b c . The total number of samples is denoted by n. In Kappa, the p e is calculated by p e = a 1 × b 1 + a 2 × b 2 + … + a c × b c n × n (10)

Our training data are augmented by horizontal flipping, vertical flipping, random scaling and rotation, and random jitter. We use the ImageNet pre-trained weights to initialize the 2D WSI and 3D MRI models, and the weights of the decoder part are initialized randomly. We used the Adam optimizer (Kingma and Ba, 2015) with an initial learning rate of 3 × 10⁻⁴ for all experiments. The learning rate decreases by 10 times at the 50^th and 80^th epochs. The training batch size is set to 24. All networks are implemented based on the PyTorch framework and trained using four NVIDIA Tesla P40 GPU cards.

We used 5-fold cross-validation based only on the training set to find the optimal network parameters for each deep learning model. The best-performing fold is taken as the final training model for the given architecture. It is noted that multiple models with different backbones are trained at each stage of the 2D WSI model, and their ensemble is used to obtain the final results.

Table 1 shows the results of our ablation experiments conducted on the validation data in the CPM-RadPath challenge. Table 2 shows the contribution of different components in our 2D classification framework. These networks are trained using the same parameter settings as described in the previous section.

As the performances shown in Table 1, our two-stage classification strategy contributes to gaining higher accuracy on three evaluation metrics. Moreover, the classification model of MRI and WSI data can complement each other to obtain more robust and accurate results. We conducted a corresponding ablation study to verify the benefits of our proposed multimodal framework and two-stage strategy. We use the same training schedule and parameter settings as the full method described in Section 4.3 for comparison with the ablation method (i.e., with the corresponding components removed).

The experiments demonstrate that our proposed multimodal framework can help improve the accuracy of glioma subtype classification. To validate the importance of multimodal image information, we use MRI or WSI alone as the training set under the most ideal conditions (two-stage), and the balanced accuracy is reduced by 6.7% or even 15.6%. It is well known that in such tasks with unbalanced data, the balanced accuracy provides a better measure of the performance of the algorithm (Brodersen et al., 2010). In addition, there is also a greater than 5% drop in F1 score in this case and a greater than 10% drop in Kappa. In contrast, the multimodal complementary model can reach the best results with 88.9% of the balanced accuracy on the validation set.

Also, we demonstrate the benefit of our two-stage scheme. Another set of experiments is conducted in this study by comparing our two-stage approach with the single-stage alternative. The first stage network is used directly to classify the three glioma subtypes. As shown in the corresponding rows of Table 1, compared to the two-stage approach, the single-stage strategy results in a 3.3 and 5.5% performance decrease in the balanced accuracy in the MRI validation set and WSI validation set, respectively. This result further demonstrates the advantages of our proposed two-stage solution. From the results, it can be seen that our network can integrate feature information from MRI and WSI data.

In addition, we also validate the advantages of the multi-task learning strategy in the 2D WSI model. We use three highly correlated classification networks to train the 2D WSI model: EfficientNet-B2, EfficientNet-B3, and SEResNeXt101. As shown in Table 2, with both classification branch (cls), regression branch (reg), and a fully connected layer (gem), we obtain better results than other methods.

Table 3 lists the results of the top six teams in the MICCAI 2020 CPM-RadPath Challenge. As shown in Table 3, our method has obtained the best classification performance in the testing phase.

In particular, we outperform the second-place team, Tabulo, by 8.8% in terms of balanced accuracy. According to a report submitted to the challenge organizer, Tabulo utilizes two additional publicly available datasets: the Multimodal Brain Tumor Segmentation Challenge 2019 (BraTS-2019) ² and MoNuSAC ³ . While our method does not use any external data, all models on the classification pipeline are not pre-trained on medical data. Plmoer, a team from the University of Pittsburgh Medical Center, which processes noisy labels in many patches from WSI of each category. It is similar to us while our two-stage strategy brings an obvious improvement. The Marvinler team used multi-instance learning for the WSI processing, which tends to have high memory requirements. In addition, the Hanchu team requires an additional image segmentation task, yet this does not improve the model performance. The Azh2 team trains a densely connected network for a specific configuration of the DenseNet, which also performs much worse than our model.

It should be noted that we are not able to upload the 3D MRI model in time, so this result is only available for the 2D WSI model with the two-stage strategy. However, it can be seen from Table 1 that the 2D WSI model combined with the 3D MRI model will improve the result substantially. So, we believe that better results will be obtained if the ensemble model is uploaded. Since the test dataset is not yet open, we would like to include the results of the ensemble model with the two-stage strategy if the test dataset is open.

Table 4 summarizes relevant work on the MICCAI 2019 and 2020 CPM-RadPath Challenge. These results are all obtained on the validation set and are all in the one-stage classification framework.

The early work (Chan et al., 2019; Pei et al., 2019) first segmented the tumors before classifying them in 2019, however, the final classification results are not satisfactory. In the MICCAI 2020 CPM-RadPath Challenge, (Pei et al., 2020; Yin et al., 2020; Zhao et al., 2020), still use the segmentation before the classification framework, with a significant improvement over last year’s results. The performance of these methods is affected by the segmentation results. Methods in (Hamidinekoo et al., 2020; Lerousseau et al., 2020) and our solutions do not require the segmentation process.

Without segmenting the images, a direct comparison between our method and (Lerousseau et al., 2020) is feasible because of the similar performance on the validation set. However, in the MICCAI 2020 CPM-RadPath Challenge, method in (Lerousseau et al., 2020) performs not well on the test set and the model does not have strong generalization ability. Similarly, although methods in (Yin et al., 2020; Zhao et al., 2020) perform equal or better than our method on the validation set, our method obtains first place on the final competition test set.

We visualize the output of our classification model. We show here the visualization of pathological images of three different glioma subtypes on a classification model. The input patch size is 1,024 × 1,024 pixels, and then each patch is put into our model to calculate the predicted probability for each glioma subtype. Each patch is labeled with the color that represents its probability. Finally, every patch is integrated to obtain the overall probability of the WSI belonging to glioblastoma or oligodendroglioma, or astrocytoma. As shown in Figure 5, probabilities are converted to color maps in logarithmic form. Therefore, those lower probabilities (<50%) are shown as very light colors, and higher probabilities (≥50%) are shown as dark colors.

Two pathologists were invited to view these pathological images and their corresponding patches, and they confirmed the presence of the typical glioma subtype patches and normal patches, which are shown in the right two columns of Figure 5.