Mouse mammary gland tumor tissue samples were harvested from PyMT (WT) and PyMT/TGFβRII^LysM (KO) mice as described in our previous work where we found increased ECM deposition in PyMT/TGFβRII^LysM samples (Stringer et al., 2021). Human lung adenocarcinoma tissue samples were collected from 14 patients after surgery under IRB approved protocol 000616 at Vanderbilt University Medical Center. Using Computer-Aided Nodule Assessment and Risk Yield (CANARY) software we stratified the patients into two risk groups (Indolent and Aggressive) regarding behavior of adenocarcinoma as we showed in (Senosain et al., 2021). Two TMA cores from each tissue sample were processed.

Tissues were fixed in 10% neutral buffered formalin for 24 h at room temperature. Fixed samples were embedded in paraffine blocks, and 5 μm sections were cut. Tissue microarrays (TMA) were generated from patient’s material. All cores were examined by a pathologist to exclude tissue with massive necrosis areas, large blood vessels and artifacts. H&E staining was used to evaluate core quality. For visualization of collagen in the tissue samples we used picrosirius red staining and CHP–Collagen Hybridizing Peptide (3Helix, United States) staining labeled with Cy3 dye. Multiplex immunofluorescent (IF) staining was used to combine CHP staining with a set of cellular markers (Senosain et al., 2021). Epithelial/tumor cells were visualized with anti-PanCK-Alexa Fluor 488 antibodies (eBioscience, United States), immune cells with - anti-CD45 antibodies (Biolegend, United States) + Fab anti-mouse-Cy3 antibodies (Jackson ImmunoResearch, United States), and T-cells with anti-CD3 antibodies (Agilent Dako, United States) + anti-rabbit-Cy5 antibodies (Thermo Fisher Scientific, United States). Nuclei were visualized by Gold Antifade Mountant with DAPI (Thermo Fisher Scientific, United States). Images from Picrosirius red stained samples were obtained using the Keyence BZ-X710 microscope (KEYENCE, United States) using ×20 objective and Texas Red filter (excitation wavelength 560/40 nm, emission wavelength 630/75 nm). Whole image scanning was performed on the Apiro Versa 200 (Leica, United States) with Texas Red filter (excitation wavelength 560/55 nm, emission wavelength 645/75 nm), Spectrum Green filter (excitation wavelength 480/30 nm, emission wavelength 535/40 nm), Cy5/AlexaFluor 647 filter (excitation wavelength 620/50 nm, emission wavelength 690/50 nm), Cy7 filter (excitation wavelength 720/75 nm, emission wavelength 810/90 nm) and DAPI filter (excitation wavelength 350/50 nm, emission wavelength 460/50 nm).

The KNIME (KNIME 4.1.2) analytical platform was used as an environment for integration of our Python package for fiber segmentation and spatial analysis with workflow for single-cell analysis of multiplexed fluorescent stained images. Segmentation, labeling and classification of cells was performed as described in our previous work (Senosain et al., 2021). Quantitative information from single-cell data (such as X, Y coordinates etc.) was used for spatial analysis using developed Python package (Python 3.7). Python nodes were incorporated in this pipeline for collagen fiber segmentation, labeling and analysis of geometrical features and spatial analysis of cells and fibers.

Set of mouse and human cancer tissue samples was used to implement performance tests and validation of the algorithm. Breast cancer image size: 1920 x 1440 pixels, human lung adenocarcinoma (whole slide scans) image size: 5120 x 5120 pixels. Collagen fibers were manually annotated by pathologist and their length was measured using ImageJ. These results were used as a ground truth for further tests. Number of highlighted fibers and their length were chosen as only features for ground truth due to restrictions of manual measurements of other parameters (width, angle etc.). Mean value of the parameters was calculated for each image. The results were presented as average for the set of images. Accuracy test was performed as a comparison to ground truth and to results obtained using CT-FIRE algorithm (Bredfeldt et al., 2014a). Performance time test was realized using a set of 50 images. The results were compared to the CT-FIRE algorithm (Bredfeldt et al., 2014a). Neighborhood outline for breast cancer samples was set as 100 pixels and for lung cancer tissue as 250 pixels and was equal to ∼50 μm and ∼100 μm in diameter respectively.

DELL Precision T7920, two Intel Xeon Gold 6136 Processors (24 Core, 3.0 Ghz, 24.75 M Cache each), 128 Gb RAM.

The Microsa source code, including trial dataset and user manual, is available at https://github.com/VGeorgii/Microsa.

Results were presented as mean ± SEM. Two-group comparison was performed using two-sample t tests or Wilcoxon Rank-Sum test as appropriate. P values less than 0.05 were considered statistically significant. Hierarchical cluster analysis was performed using Ward’s method for Euclidian distances.

Implementation of computational image analysis for scientific and clinical purposes is a robust, effective, and precise approach. The complexity of tissue structure requires novel methods and technologies for image acquisition and processing that will extend the dimensionality of extracted information. Tissues are characterized by a combination of objects with certain morphology, number, and spatial distribution. Cells in tissue can interact by close (cell junctions) and distant (cytokines, chemokines, metabolites etc.) types of communication. The distance between objects plays a crucial role in the determining the effects of such interactions. Our goal was to develop and test computational methods which would be able to combine analysis of cellular and fibrous constituents of tissue structure in a spatially dependent manner (Figure 1A). Spatial analysis can be implemented to identify of structural patterns in tissue that help to decipher disease pathogenesis or that can be used for outcome prediction. There are several approaches that can be used for spatial analysis: distance calculation (Parra, 2021), spatial homogeneity assessment (Holmes et al., 2009), and clustering or graph algorithms (Roeder et al., 2012; Heindl et al., 2015; Parra, 2021). Our tool implements distance calculation. This approach can be realized in two ways. The first approach (“number”) is based on the analysis of neighborhood outlines around processed cells or other objects in the tissue. The radius of outline is set by the user and determined by the type of interactions the researcher is interested in (Figure 1B). A small radius allows the researchers to calculate objects in close contact to the processing object, while a bigger radius can be used to calculate the number of objects in distant contact. The second approach (“distance”) involves calculating the distance between the processed object and other objects of interest. For instance, the distance between a tumor cell and the closest T-cell, or between a T cell and a blood vessel etc. After that, extracted information can be passed to K-nearest neighbors, N-distance, or self-organizing map algorithms for identification of neighborhoods with similar features (Roeder et al., 2012; Heindl et al., 2015; Parra, 2021).

To perform segmentation, labeling and geometrical feature analysis of fibers we used the Python programming language. This decision was determined by the flexibility of Python and diversity of packages designed for image processing that demonstrate good performance and efficacy. The developed Python package was designed for assessment of fiber features and combination of extracted information with cell-based datasets for spatial analysis. Our package uses several popular python libraries: Pandas (https://pandas.pydata.org/), SciPy (Virtanen et al., 2020), NumPy (Harris et al., 2020) and Scikit-image (van der Walt et al., 2014), Matplotlib (Hunter, 2007). The package consists of four modules: i) fiber segmentation and labeling, ii) fiber geometrical feature calculation and spatial analysis, iii) cell-based spatial analysis, and iv) visualization (Figure 1C).

Before implementing of these modules of our Python package, the user needs to perform foreground-background separation. By our experience, the most efficient methods are simple thresholding and Frangi filter-based thresholding (Frangi et al., 1998; Hemler et al., 2004; Shi and Yang, 2009; Park et al., 2013; Comin et al., 2014). Simple thresholding efficiently processes images with small numbers of fibers and low levels of background. The Frangi filter method was designed specifically for objects that form net-like structures (blood vessels, fibers etc), which can better identify densely packed fibers (Frangi et al., 1998).

ECM fibers in the tissue often have complex curve-like structure with different lengths, widths, and linearities, and often have additional branches. To perform single-fiber analysis, we decided to process each fiber as a curve with no branches, which reflects general geometrical features and does not have intersections with another fibers. We found that the most efficient way to generate a simplified model of the fibers is to use thinning or skeletonization algorithms (Comin et al., 2014). However, because of the complexity of collagen fiber shapes the algorithm we used for skeletonization generated multiple artifacts that appeared as spikes or additional small branches. We therefore implemented an extra process, which helps to remove these artifacts from the fiber skeleton. The visualization module can be used to overlay segmented fibers with the original image to control quality of segmentation and labeling (Figure 2A). After fiber segmentation and labeling is complete, the geometrical feature calculation module can be used to extract fiber properties.

This parameter represents the length of the longest branch of the fiber. The value is calculated as length of the curve between two endpoints.

This parameter represents average distance from pixels of the fiber’s skeleton to the outline of fiber’s mask. T h i c k n e s s = ( ∑ i = 0 N d i ) / N where d–distance from pixel to fiber’s mask outline, N number of the pixel in the fiber’s skeleton.

This parameter represents the angle between the X axis and a straight line between the fiber’s curve endpoints. Angle can take values between 0 and 180°.

This parameter represents average intensity of the fiber’s skeleton pixels I n t e n s i t y = ( ∑ i = 0 N I i ) / N where I–pixel intensity (from 0 to 255), N–total number of pixels in the fiber’s skeleton.

This coefficient represents how close the shape of the fiber is to straight line. It takes the values from 0 to 1. S t r a i g h t n e s s = l / L where l–length of fiber’s curve, L–Euclidian distance between endpoints.

This coefficient reflects the alignment of processing fiber and its neighboring fibers. It can take values from 0 (not aligned) to 1 (absolutely aligned). The number of neighboring fibers depends on radius of the neighborhood outline which can be regulated. A l i g n m e n t = ( ∑ i = 0 n cos ( α i ) ) / n where α – angle between processing fiber and neighboring fiber i , n–total number of neighboring fibers.

The output of these modules is a dataset which contains information about each fiber’s location (XY coordinates) and geometrical features. Implementation of cell-based spatial analysis gives an opportunity to investigate features of fibers co-localized with processing cells.

To test the performance of fiber segmentation and fiber features analysis algorithms we used a set of images obtained from mouse mammary gland tumor samples stained with Picrosirius red (Figures 2Ba). The results of our algorithm have been compared to other software designed for fiber analysis–CT-FIRE (Stein et al., 2008; Bredfeldt et al., 2014a). This approach was implemented in multiple works (Bredfeldt et al., 2014a; Bredfeldt et al., 2014b). Our team was inspired by idea that is realized in CT-FIRE–to track and analyze geometrical features of each fiber in the tissue image. CT-FIRE uses different method for segmentation and labeling of fibers and we were interested to compare performance of our method and CT-FIRE. Both algorithms were executed using parallel computation mode (13 cores). We found that our approach demonstrated better results (Figure 2Bb). To assess the accuracy of our developed algorithm we used a set of samples manually annotated by pathologist and compared with results generated by our approach and CT-FIRE. We found that number of fibers segmented and labeled by our method is closer to the ground truth (Figure 2Bc)). Both approaches demonstrated good accuracy estimating length of segmented fibers (Figure 2Be). Additional analysis demonstrated that pruning is the most time-consuming step in the fiber segmentation algorithm and needs further amendments or revision (Figure 2Bd).

This package was designed to combine single-cell and single-fiber datasets in spatially dependent manner. The generated data can be used for statistics and machine learning methods. Cell segmentation and labeling must be performed using appropriate software or algorithms (QuPath (Bankhead et al., 2017), KNIME (Dietz and Berthold, 2016)) and generated single-cell datasets can be utilized as an input for spatial analysis algorithm. To generate single-fiber datasets the user can use images with fibers visualized by brightfield staining (Picrosirius red, Trichrome blue, IHC) or IF staining (antibodies against collagen CHP/3Helix etc.) techniques and apply modules for fiber segmentation and fiber geometrical feature analysis from our package, or utilize other methods or software designed for that purpose. Using modules for spatial analysis, the user can generate additional information for each labeled cell or fiber based on objects (cells or fibers) localized in a neighborhood outline of a given radius (Figure 2C).

To validate the accuracy and performance of our developed package and test its integration with other algorithms and software, we used several sets of tissue section slides stained by different techniques. Geometrical features, deposition of collagen fibers and their visualization methods can significantly affect algorithm performance. In our current work, we used Picrosirius red and collagen hybridizing peptide (CHP/3Helix) staining to visualize fibers.

In our previous work we have used the PyMT mouse model (spontaneous mammary gland tumor formation) (Vasiukov et al., 2020). This work demonstrated that lack of TGFβ signaling in myeloid cells (PyMT/TGFβRII^LysM) is related to reduced tumor growth and increased collagen deposition. Here, for visualization of collagen fibers we have used Picrosirius red staining with subsequent calculation of collagen deposition. One of the advantages of Picrosirius red staining is the ability to use polarized light to identify immature and mature collagen fibers. In our work we demonstrated that in the KO group tumor tissue had higher amounts of immature collagen that indicate more intensive ECM remodeling. Using the fiber segmentation module in our developed package, we detected increased number of collagen fibers in tumor tissue in the PyMT /TGFβRII^LysM group. This data confirms results we obtained in our previous work using different calculation methods. In addition, analysis of geometrical features of collagen fibers demonstrated tendency that collagen fibers in the PyMT /TGFβRII^LysM group have increased straightness and alignment, but decreased thickness in comparison to PyMT/TGFβRII^WT group (Figure 3A).

For collagen detection in human breast cancer tissue samples, we used CHP staining. CHP visualizes collagen well with low background noise. Moreover, this method can be used for multiplex IF staining that allows simultaneous detection of different biological markers. We visually chose three different areas that were characterized by different collagen deposition and fiber spatial distribution. Analysis demonstrated that Area 1 had the highest number of collagen fibers, whereas Area three had the lowest number. Investigation of fibers co-localization demonstrated similar results. Area 1 showed the highest rate of fiber co-localization (Figure 3B).

As was mentioned above, CHP staining is suitable for multiplex IF that allows simultaneous visualization of cellular markers and collagen fibers and can be used for integrated analysis to reveal relationships between them. In our previous work (Senosain et al., 2021), we used a panel of antibodies for visualizing of epithelial/tumor cells (PanCK), immune cells (CD45) and T-cells (CD3) (Figure 4A) to reveal relationships between anticancer immune response and behavior (indolent or aggressive) of human lung adenocarcinomas. To perform cell’s segmentation, labeling and classification we utilized workflow developed in KNIME. In current work, we additionally performed CHP (3Helix) staining to identify collagen fibers. To analyze tumor tissue from patients with indolent and aggressive lung adenocarcinoma, we used dataset generated in our previous work (Senosain et al., 2021). First, we performed analysis of basic parameters of tumor tissue (total number of different type of cells, total number of collagen fibers and their morphology). We did not find a significant difference in cellular composition between indolent and aggressive groups of patients. Calculation of collagen fibers showed that there is no difference in number of fibers and their geometrical features (Figure 4B). Next, we performed spatial analysis and used tumor cell as a processing object. Neighborhood outline was equal to 100 μm. Cells and collagen fibers within the neighborhood outline was counted and geometrical features of fibers were estimated (Figure 4C). The analysis revealed that tumor cells from the aggressive group were co-localized with lower number of CD3⁻ immune cells. In addition, tumor cells from aggressive lung adenocarcinomas were co-localized with lower number of collagen fibers and these fibers generally had smaller length (Figures 4D,E).