Automated analysis of scattering-based light sheet microscopy images of anal squamous intraepithelial lesions

We developed an algorithm for automatically analyzing scattering-based light sheet microscopy (sLSM) images of anal squamous intraepithelial lesions. We developed a method for automatically segmenting sLSM images for nuclei and calculating seven features: nuclear intensity, intensity slope as a function of depth, nuclear-to-nuclear distance, nuclear-to-cytoplasm ratio, cell density, nuclear area, and proportion of pixels corresponding to nuclei. 187 images from 80 anal biopsies were used for feature analysis and classifier development. The automated nuclear segmentation method provided reliable performance with the precision of 0.97 and recall of 0.91 when compared with the manual segmentation. Among the seven features, six showed statistically significant differences between high-grade squamous intraepithelial lesion (HSIL) and non-HSIL (non-dysplastic or low-grade squamous intraepithelial lesion, LSIL). A classifier using linear support vector machine (SVM) achieved promising performance in diagnosing HSIL versus non-HSIL: sensitivity of 90%, specificity of 70%, and area under the curve (AUC) of 0.89 for per-image diagnosis, and sensitivity of 90%, specificity of 80%, and AUC of 0.92 for per-biopsy diagnosis.


Introduction
Anal cancer is an under-served disease with rapidly increasing incidence and mortality.The incidence of anal squamous carcinoma increased by 2.7% per year between 2001 and 2015.Within this period, both incidence of higher stage disease and anal cancer mortality increased [1].The trend continued, with the anal cancer mortality increasing by 5.7% per year between 2014 and 2020 [2].The recent ANCHOR (Anal Cancer HSIL Outcomes Research) trial showed that treatment of anal precancer can significantly reduce the risk of anal cancer [3].However, the standard methods for screening of anal cancer have limitations.Anal cytology via Pap smear is the first-line screening method, however, it provides moderate diagnostic performance, 81% sensitivity and 63% specificity for detecting anal precancer (high-grade squamous intraepithelial lesion, HSIL) [4].Therefore, patients with abnormal anal cytology results need to undergo high-resolution anoscopy-guided biopsy to make the final diagnosis.High-resolution anoscopy provides magnified views of the anal tissue, however, it cannot visualize microscopic details of the tissue that are critical in accurately diagnosing anal lesions.
In vivo optical diagnostic methods can non-invasively examine microscopic details of the tissue.Previously, light scattering spectroscopy methods were evaluated for diagnosing HSIL of cervix tissues [5], which are biologically and histomorphologically similar to anal tissues [6].In light scattering spectroscopy, differences in nuclear morphologic features (e.g.nuclear diameter) between HSIL and non-HSIL (non-dysplastic or low-grade intraepithelial lesion, LSIL) generate differences in light scattering spectra.A previous study evaluating light scattering spectroscopy demonstrated 77% sensitivity and 62% specificity for diagnosing cervical HSIL versus non-HSIL [5].Angle-resolved low coherence interferometry (a/LCI) advanced light scattering spectroscopy further to investigate depth-resolved nuclear morphologic features.Recent in vivo imaging studies showed that a/LCI could provide great diagnostic performance, 100% sensitivity and 71-82% specificity for diagnosing HSIL versus non-HSIL in cervix [7].Confocal microscopy is another promising in vivo optical diagnostic approach because it can directly visualize nuclei.In previous studies, quantitative analysis of confocal microscopy images of cervix provided a high sensitivity (100%) and specificity (92-100%) for diagnosing HSIL versus non-HSIL [8,9].High-resolution micro-endoscopy (HRME) is a low-cost in vivo microscopy approach that can also directly visualize nuclei.HRME showed a promising sensitivity of 84-97% and specificity of 54-74% for diagnosing high-grade cervical lesions in low-resource settings [10,11].More recently, HRME was evaluated for examining anal lesions in vivo.A deep learning-based model was used to analyze HRME images, which demonstrated great diagnostic performance, sensitivity of 91% and specificity of 87% [12,13].
Scattering-based light sheet microscopy (sLSM) is another in vivo optical diagnostic approach that can visualize nuclei with a high resolution (lateral = 1-2 µm; axial = 5-6 µm) using intrinsic scattering contrast [14,15].In sLSM, the tissue is illuminated with light sheet at 45°relative to the tissue surface, and scattered light is detected at 90°from the illumination beam with low-NA microscope optics.In a recent study, a bench sLSM device was evaluated for imaging anal biopsies ex vivo [16].This study showed that sLSM could visualize key nuclear morphologic features used during standard of care histopathologic diagnosis of anal squamous intraepithelial lesions.Additionally, 11 pathologists made diagnosis based on sLSM images, which provided 91% sensitivity and 85% specificity for diagnosing HSIL versus non-HSIL [16].The feasibility of integrating sLSM into a compact probe was also demonstrated recently [17].
One potential use case of sLSM is to guide anal biopsy.In this approach, the anoscopist first examines the anal mucosa with high-resolution anoscopy to identify areas suspicious of HSIL and images the suspicious areas with the sLSM probe.sLSM images are then analyzed to determine the probability of HSIL.This information can be used to ensure biopsy of likely HSIL tissues and avoid biopsy of obviously non-HSIL (non-dysplastic or LSIL) tissues.One critical task in this sLSM-guided biopsy approach is to interpret sLSM images at the time of anoscopy.Our previous study showed that pathologists were able to make accurate diagnosis based on sLSM images after a short training session of reviewing only 8 pairs of sLSM and hematoxylin and eosin (H&E) images.Therefore, training pathologists to read sLSM images might not be a significant hurdle.Once trained, the pathologist can analyze sLSM images in real time either at the anoscopy clinic or remotely.
However, relying on pathologist's interpretation could hamper the adoption of sLSM due to their limited availability.Automated analysis of sLSM images can mitigate this issue by providing real-time diagnostic information that the anoscopist can act upon.In this paper, we present the development of an algorithm for automatically analyzing sLSM images of anal squamous intraepithelial lesions.The processes for automatically segmenting sLSM images for nuclei and generating image features are presented.Results from feature analysis and classifier development are also presented.

Scattering-based light sheet microscopy
The working principle of sLSM was described previously [15].Briefly, light from an LED (center wavelength = 640 nm; bandwidth = 25 nm) and a slit (width = 5 µm; length = 3 mm) was focused as light sheet on the tissue by relay optics with the demagnification of 2 (Fig. 1).As the light sheet intersects the tissue, nuclei scatter a portion of the illumination light.The scattered light is captured by the detection objective lens (focal length = 20 mm; water immersion; NA = 0.3) and focused to generate a magnified image of the tissue.The bench sLSM device used for imaging anal biopsies ex vivo had the lateral resolution of 1.4 µm, and axial resolution of 5.2 µm.The anal biopsies were treated with a low-concentration acetic acid before sLSM imaging to enhance the nuclear contrast.The biopsy sample was translated laterally by a motorized stage with the step size of 100 µm.An sLSM image was acquired at each lateral location.This made the vertical distance between two neighboring sLSM images along the tissue depth direction also 100 µm.Neighboring sLSM images did not have any overlaps.We used sLSM images of anal biopsies collected from our previous ex vivo study [16].187 sLSM images obtained from 80 biopsies were used in image analysis.Remaining 30 biopsies out of 110 biopsies from our previous study were not included in the automated image analysis because these 30 biopsies showed columnar epithelium, denuded tissue, or unclear nuclear contrast.Among the 80 biopsies, 14 were non-dysplastic, 36 were LSIL, and 30 were HSIL per final pathology report.Representative sLSM images are shown in Fig. 2. Nuclei are visualized as bright dots in these images.

Image pre-processing
We used MATLAB 2024a (Mathworks) installed on a PC with the following specifications: Intel Core i9-12900F CPU, 32GB DDR5 RAM, NVIDIA GeForce RTX 3080 GPU with 10GB VRAM.sLSM images were pre-processed before nuclear segmentation and feature extraction.Each image was cropped to a size of 1064 × 1064 pixels, which corresponds to 500 µm × 500 µm on the tissue space (Fig. 3(a)).Bright signal at the top generated by specular reflection of the tissue surface and/or imaging window was masked out using Gaussian filtering followed by binary thresholding (Fig. 3(b)).Next, the tissue surface was identified using adaptive binarization based on Otsu's method [18] and edge detection (Fig. 3(c)).The detected tissue surface was used as the reference to perform column shifting, flattening the tissue surface (Fig. 3(d)).The column-shifted image underwent another round of surface detection to identify the row index of the tissue surface (Fig. 3(e)).Most sLSM images were acquired with an exposure time of 0.1 seconds, but different exposure times were used for a small number of biopsies to adjust the brightness displayed during the image acquisition.The actual exposure time was used to standardize the image intensity to the reference exposure time of 0.1 seconds (Fig. 3(f)).

Automated segmentation of nuclei
Our overall approach is to segment sLSM images for nuclei and analyze morphometric features of the segmented nuclei.Each sLSM image visualizes nuclei at multiple depth levels, and the nuclear contrast and intensity change as a function of depth within the image.Therefore, we developed a segmentation method that binarizes each row (or each depth) of the image at a time.At each row (yellow line in Fig. 4(a)), an intensity profile was generated (Fig. 4(b)).The background signal of the intensity profile was calculated by using a one-dimensional adaptation of the rolling ball background subtraction method [19] with the ball radius of 25 pixels.The background signal curve was then subtracted from the intensity profile (Fig. 4(c)).The binarization threshold (red dotted line in Fig. 4(c)) was calculated as the average of the 91 st percentile intensity and 15 th percentile intensity.The reference percentiles for threshold calculation were optimized via grid search to have a high precision and recall for segmenting nuclei as discussed later.Finally, the intensity profile was binarized by the threshold as shown in Fig. 4(d).This process was repeated for every row to generate a two-dimensional binary mask showing nuclei.The segmentation performance was evaluated by comparing automatically-segmented nuclei to manually-segmented nuclei.850 manually-segmented nuclei from 9 sLSM images (3 nondysplastic, 3 LSIL, and 3 HSIL images) were used for evaluation (Figs.5(c, d)).The row-by-row binarization method (Figs.5(e, f)) provided good segmentation performance, precision of 0.97 (95% confidence interval, CI = 0.95-0.98)and recall of 0.91 (95% CI = 0.88-0.92).The corresponding F1 score was 0.94.
We compared our row-by-row binarization approach to a deep learning-based cell segmentation method, StarDist [20].The sLSM image (Figs.5(a),(b)) was background-subtracted using the twodimensional rolling ball method [19] with a radius of 25 pixels.The background-subtracted image was segmented by the Fiji plugin implementation of StarDist.Various probability thresholds were evaluated for segmentation performance via grid search.When the threshold was set to produce the highest F1 score, 0.90, the precision was 0.90 (95% CI = 0.88-0.92)and recall was 0.91 (95% CI = 0.88-0.92).Figures 5(g),(h) show StarDist-segmented images at this probability threshold.It is notable that the StarDist-segmented HSIL image (Fig. 5(h)) shows some of the nuclei at deeper regions are not segmented, and segmented nuclei are smaller than the manually-segmented nuclei.The precision was 0.94 and recall was 0.86 for this StarDist-segmented image, while the precision was 0.97 and recall was 0.90 using the row-by-row binarization method (Fig. 5(f)).While the performance of StarDist might be further improved by fine-tuning, we decided to use the row-by-row binarization approach because of its better performance especially for HSIL images.

Feature calculations
Seven features were calculated from the automatically-segmented nuclei within the depth range of 20-100 µm.The lower bound of 20 µm was used to exclude any residual contribution of the surface specular reflection.The upper bound of 100 µm was used because nuclear features were reliably visible up to the depth of 100 µm in most sLSM images and also because previous confocal  (c, d) manually-segmented images, (e, f) automatically-segmented images using the row-by-row binarization, (g, h) automatically-segmented images using StarDist, (i, j) nuclear-to-nuclear distance estimation using Delaunay neighbors, and (k, l) nuclear-to-cytoplasmic ratio estimation using Voronoi diagrams.
microscopy studies showed that the imaging depth of 50 µm revealed nuclear morphologic features relevant for diagnosis of HSIL versus non-HSIL [8,9].
The first two features were related to nuclear intensity.The sLSM image before background subtraction (Figs.5(a),(b)) was multiplied with the binary nuclear mask (Figs.5(e),(f)).Intensity values of non-zero pixels, also the pixels corresponding to nuclei, were averaged to calculate the representative nuclear intensity for each image.Nuclear intensity versus depth curve was generated using every row, and the slope of this curve was calculated by linear fitting.The linear fitting performance was acceptable with the average R 2 value of 0.59 for all 187 images.
The next two features, nuclear-to-nuclear distance and nuclear-to-cytoplasmic ratio (N/C ratio), were estimated based on the approach previously developed for analyzing confocal microscopy images of cervical biopsies [9].Nuclear-to-nuclear distance was calculated as the average of distances between a nucleus to its three closest Delaunay neighbors (Figs.5(i),(j)).N/C ratio was estimated by finding centroids of nuclei, generating a Voronoi diagram using the centroids (Figs.5(k),(l)), calculating the area of each Voronoi polygon (estimated cellular area) and the area of its enclosed nucleus (estimated nuclear area), and dividing the nuclear area by the difference between the cellular area and nuclear area (estimated cytoplasmic area).
The last three features were cell density, nuclear area, and proportion of nuclear pixels.Cell density was calculated by dividing the number of nuclei by the image area in the unit of µm 2 .Nuclear area was calculated by averaging areas of nuclei.Proportion of nuclear pixels was calculated by dividing the number of non-zero pixels in the binary nuclear mask (Figs.5(e),(f)) by the number of all pixels of the image.
Statistical significance was assessed using one-way analysis of variance (ANOVA) and Bonferroni post-hoc multiple comparisons.The statistical analysis was conducted using the MATLAB 2024a Statistics and Machine Learning Toolbox.

Classifier development and evaluation
We evaluated a linear support vector machine (SVM) classifier for diagnosing HSIL versus non-HSIL (non-dysplastic or LSIL).Classifications were made into two categories (non-HSIL, HSIL) instead of three (non-dysplastic, LSIL, HSIL).This is because HSIL finding would indicate immediate biopsy or treatment, while non-HSIL finding would indicate monitoring over time.
All seven features were used during the classifier development.Before training and testing, each feature was normalized by applying median interquartile range normalization.The dataset was then divided into training and testing sets.We evaluated classifier performance through a five-fold cross validation.The data was randomly divided into five folds.The data stratification was based on biopsies rather than images.This ensured that all feature data from each biopsy belonged to the same fold, making data from each biopsy used for either training or testing but not for both.Each fold included feature data from approximately 37 images obtained from 3 non-dysplastic, 7 LSIL, and 6 HSIL biopsies.Four folds (approximately 150 images from 12 non-dysplastic, 28 LSIL, and 24 HSIL) were used to train the classifier, and the remaining fold was used for testing.This process was repeated five times to test each image exactly once.
Classification scores for all 187 images were aggregated from the five folds.Classification scores and their corresponding histopathologic diagnoses were used to generate a receiver operating characteristic (ROC) curve for diagnosing HSIL versus non-HSIL.This ROC curve showed the performance of per-image diagnosis.We also evaluate the performance of per-biopsy diagnosis.For each biopsy, the minimum of classification scores was chosen as the representative score for the biopsy.We hypothesized that choosing the minimum score would reduce false positives and increase the specificity of diagnosing HSIL.The ROC curve for per-biopsy diagnosis was generated by comparing 80 classification scores from 80 biopsies with the corresponding histopathologic diagnoses.The area under the curve (AUC) was calculated for each ROC curve.In order to compensate for the imbalance in data (more non-HSIL data than HSIL data), the weight of false negatives in the cost function was optimized to maximize the AUC of the ROC curve via grid search.The operating threshold was set to achieve the sensitivity of 90%.

Feature analysis
The processing time for pre-processing, nuclear segmentation, and feature calculation of each image was 0.55 (± 0.19) seconds.Figure 6 shows a set of two-dimensional scatter plots of two related features.As shown in Fig. 6(a), nuclear intensity was higher, and intensity slope was more negative in HSIL than in non-dysplastic or LSIL.The more negative slope, or more signal attenuation over depth, in HSIL can be explained by the stronger nuclear scattering and higher cell density in HSIL.Nuclear-to-nuclear distance was smaller, and N/C ratio was larger in HSIL (Fig. 6(b)), which is in good agreement with the pathologists' assessment of sLSM images of HSIL: enlarged and crowded nuclei.While cell density was higher in HSIL, nuclear area did not reveal noticeable difference between three categories (y-axis of Fig. 6(c)).Proportion of nuclear pixels, which is essentially the product of cell density and nuclear area, was larger in HSIL than in non-dysplastic or LSIL (Fig. 6(d)).
Figure 7 shows a set of scatter plots for each feature along with the average of the feature for each diagnostic category.The trends observed in two-dimensional scatter plots (Fig. 6) are also noted in Fig. 7: higher intensity, more negative intensity slope, smaller nuclear-to-nuclear distance, larger N/C ratio, higher cell density, and larger proportion of nuclear pixels in HSIL images.Nuclear area did not show statistically significant differences between the three diagnostic categories.Four of the morphologic features were analyzed in the previous confocal microscopy study [9] (nuclear-nuclear distance, N/C ratio, cell density, nuclear area) .Three of these features showed similar trends to those found in confocal microscopy images of HSIL: smaller nuclear-to-nuclear distance, larger N/C ratio, and higher cell density than those in non-dysplastic or LSIL images.While the previous confocal microscopy study found that nuclear area of non-dysplastic was significantly different from that of LSIL or HSIL, the same trend was not observed in sLSM images.The same confocal microscopy study found no significant differences in nuclear area between LSIL and HSIL, which agreed with our finding.Results from one-way ANOVAs and multiple comparisons are summarized in Table 1.

Classifier performance
The processing time for generating the classification score from the seven features of each image was short, 4.0 × 10 −5 (± 2.3 × 10 −5 ) seconds.The ROC curve for per-image diagnosis is shown in Fig. 8(a).At the operating point providing the sensitivity of 90% (red circle in Fig. 8(a)), the specificity was 70%, and accuracy 77%.The AUC was 0.89.The ROC curve for per-biopsy diagnosis is shown in Fig. 8(b).At the operating point with the sensitivity of 90%, the specificity was 80%, and accuracy 84%.The AUC for per-biopsy diagnosis was 0.92.The performance of the automated classification was comparable to manual diagnosis of sLSM (blue squares) and H&E images (green squares) as summarized in Table 2 [16].Of the 80 biopsies included in sLSM image analysis, 35 were submitted for p16 immunohistochemistry staining per standard of care.p16 staining is used when H&E-stained images are morphologically ambiguous in diagnosing anal intraepithelial neoplasia 2 (AIN2), a type of HSIL, versus LSIL or reactive metaplasia [6].Among the 35 biopsies submitted for p16 staining, 17 were non-HSIL, and 18 were HSIL per final pathology report.We evaluated the differences of classification scores between biopsies that underwent p16 staining and those that did not.As shown in Fig. 9, classification scores for non-HSIL biopsies with p16 staining (black circles) were higher than the scores for those without (p-value = 9.6 × 10 −7 ).Conversely, HSIL biopsies analyzed with p16 staining (black circles) had lower classification scores compared to those without (p-value = 1.28 × 10 −2 ).The overlap in classification scores between non-HSIL and HSIL images was greater for biopsies with p16 staining.The d' value, calculated by dividing the difference between the means of non-HSIL and HSIL scores by the square root of their average variance, was 1.05 for the p16 group and 2.14 for non-p16 group.Classification score HSIL Non-HSIL Fig. 9. Classification scores for non-HSIL and HSIL images.black circles -scores for images of biopsies that went through p16, and white circles -scores for images without p16.

Discussion and conclusion
In this paper, we presented the development of an automated algorithm for analyzing sLSM images of anal squamous intraepithelial lesions.We developed a rule-based, row-by-row binarization method for automatically segmenting sLSM images for nuclei.The segmentation method provided a high precision and recall when compared with manual segmentation.Among the seven features extracted from segmented nuclei, six showed statistically significant difference between HSIL and non-dysplastic or LSIL.The diagnostic performance of the trained classifier was promising.The linear SVM classifier demonstrated good per-image diagnostic performance (sensitivity = 90%, specificity = 70%) and also good per-biopsy diagnostic performance (sensitivity = 90%, specificity = 80%).The per-biopsy sensitivity and specificity were similar to those of manual analysis of sLSM images by 11 pathologists, 91% and 85%, respectively.Additionally, the per-biopsy diagnostic performance was comparable to that of other optical imaging methods successfully evaluated for in vivo diagnosis of HSIL, including a/LCI (sensitivity = 100%, specificity =71-82%) and HRME (sensitivity = 91%, specificity = 87%) [7,12,13].
In this work, use of the rule-based method of segmenting nuclei and extracting features showed promising diagnostic performance.However, we expect that the rule-based method will face challenges when there is debris above the tissue surface, the surface is fragmented or folded, or the keratinized layer is thick on the surface.Additionally, the rule-based method was slow, 0.55 seconds for calculating seven features per image, and might not have captured some of the key morphologic features pathologists use.Deep learning-based methods could achieve high diagnostic performance for these potentially challenging cases, decrease the processing speed, and improve overall diagnostic performance beyond what is achievable with rule-based methods.In the future, we will explore various deep learning methods of using a convolutional neural network (CNN) followed by a fully-connected network [12,21] for analyzing sLSM images.
The image analysis algorithm and its diagnostic performance were optimized for the use case of guiding biopsy during high-resolution anoscopy.The biopsy tissues imaged by sLSM were those deemed suspicious per expert anoscopist's assessment of high-resolution anoscopy images.These tissues, especially those determined non-dysplastic per final pathology report, might show different sLSM features than the tissues appearing normal under high-resolution anoscopy.Therefore, for a different use case such as imaging the entire anal canal with a miniaturized sLSM probe as a first-line screening method, a new set of sLSM images will need to be obtained from normal-looking tissue areas and used to refine the image analysis algorithm.
One challenge in using sLSM for guiding anal biopsy can be its small FOV.In this study, 500 µm was used as the FOV, which is similar to the high-power field diameter used during histopathologic analysis of H&E slides.However, the sLSM FOV is significantly smaller than the diameter of the anal tissue accessible through the anoscope, 15 mm.Our current plan is to have the sLSM probe manually maneuvered by the anoscopist.Therefore, any tissue locations deemed benign or low-grade by the anoscopist will not be imaged by sLSM.This can be problematic because there is a shortage of trained anoscopists and diagnostic performance of high-resolution anoscopy needs improvement [22][23][24].A potential solution to address this issue is systematic sLSM imaging of pre-determined tissue locations (e.g. at every quadrant of the circular tissue area accessible for imaging at the distal end of the anoscope) in addition to the tissue locations deemed high-grade per anoscopist's assessment.Another direction is using deep learning-based methods to analyze high-resolution anoscopy images and identify high-risk tissue areas [25] and image the deep learning-identified tissue areas with sLSM .This approach can be helpful in guiding anoscopists, especially those with less experience in high-resolution anoscopy.
Funding.National Cancer Institute (P30CA023074); National Institute of Biomedical Imaging and Bioengineering

Fig. 4 .
Fig. 4. Row-by-row binarization of sLSM image.(a) sLSM image with the yellow line indicating the row for intensity profile, (b) intensity profile before background subtraction, (c) intensity profile after background subtraction, and (d) binarized intensity profile.

Fig. 5 .
Fig. 5. Automated nuclear segmentation and analysis of nuclear morphometric features for representative LSIL and HSIL images.(a, b) sLSM images,(c, d) manually-segmented images, (e, f) automatically-segmented images using the row-by-row binarization, (g, h) automatically-segmented images using StarDist, (i, j) nuclear-to-nuclear distance estimation using Delaunay neighbors, and (k, l) nuclear-to-cytoplasmic ratio estimation using Voronoi diagrams.

Fig. 6 .
Fig. 6.Two-dimensional scatter plots of feature distributions for three diagnostic categories: non-dysplastic (green), LSIL (yellow), and HSIL (red).(a) Nuclear intensity versus intensity slope, (b) nuclear-to-nuclear distance versus N/C ratio, (c) cell density versus nuclear area, and (d) cell density versus proportion of nuclear pixels.

Fig. 8 .
Fig. 8. ROC curves for diagnostic performance evaluation.(a) ROC curve for per-image diagnosis, and (b) ROC curve for per-biopsy diagnosis.red circles -operating points for automated classification, green square -operating point for manual diagnosis of H&E images, and blue square -operating point for manual diagnosis of sLSM images.