Image processing method and system based on in-situ hybridization technology and medium

Abstract

The invention discloses an image processing method and system based on an in-situ hybridization technology and a medium, and relates to the technical field of biological information, and the method comprises the following steps: converting a cell DAPI dyeing result into two-dimensional data from three-dimensional data by utilizing Z-axis maximum intensity projection, identifying and separating a single cell from the two-dimensional data by utilizing cell segmentation, and when the cell segmentation is used for processing an overlapping region, identifying and separating the single cell from the two-dimensional data. Using a registration algorithm Ashlar to calculate an error between adjacent visual fields of the same round to obtain cell position information; an imaging result after hybridization of the fluorescent probe and the gene is subjected to multiple rounds of registration through a registration algorithm Ashlar to generate a spliced image, local maximum values of the spliced image under different rounds are marked by using a fluorescent dot recognition algorithm, so that position information of each fluorescent dot is obtained, the position information of each fluorescent dot corresponds to the gene in a transcript, and the position information of each fluorescent dot corresponds to the gene in the transcript. Obtaining gene position information; distributing genes into cells by utilizing the cell position information and the gene position information to obtain a cell gene matrix and displaying the cell gene matrix; according to the image processing method and system and the medium, accurate matching of sequencing data and spatial information is achieved, and then the distribution rule of gene expression in space is better revealed.

Description

Technical Field

[0001] The present invention relates to the field of bioinformation technology, and in particular to an image processing method, system and medium based on in situ hybridization technology. Background Art

[0002] Spatial transcriptomics technology is an important biological research method that combines traditional single-cell sequencing technology, in situ hybridization technology and other omics technologies to obtain gene expression data while preserving the spatial location information of the sample. Currently, there are two major technical paths for the implementation of spatial transcriptomics technology: sequencing-based and imaging-based. Compared with sequencing-based technology, imaging-based methods have higher requirements for equipment and operation, but have high resolution and can easily achieve single-cell accuracy. At the same time, their cost is relatively low, and the cost per sample is diluted as the number of samples increases, thus having the advantages of high precision and low cost. All imaging-based spatial transcriptomics technology processes include three parts: sample preprocessing, microscopic imaging and image processing. Among them, the image processing system is responsible for accurately extracting, parsing and visualizing the spatial information of gene expression from complex and large image data.

[0003] Currently, data processing pipelines for sequencing-based spatial transcriptomics technologies are relatively well-developed. However, the development of data processing pipelines for imaging-based spatial transcriptomics technologies currently faces more challenges. This is because imaging data requires complex image processing, including image enhancement, denoising, segmentation, and feature extraction. Furthermore, the data generated by imaging-based spatial transcriptomics technologies typically contain multiple dimensions, such as spatial location, gene type, and gene expression levels. These data require precise alignment, registration, and quantification to ensure the accuracy of subsequent analysis.

[0004] The current data processing pipeline for imaging-based spatial transcriptomics includes steps such as Z-axis maximum intensity projection, inter-round registration, cell segmentation, point identification, decoding, splicing results, gene assignment, and quality control. When imaging a sample using a microscope, the microscope scans each field of view in a Z-shaped or serpentine pattern. However, during this movement, due to the precision of the microscope's motion module, the actual imaging result may be horizontally offset from the defined field of view. This offset affects the accuracy of subsequent splicing of single-field cell segmentation results and point identification results at the edges of the field of view. Specifically, problems such as cell overlap, cell loss, and abnormal cell morphology may occur at the edges of the field of view, leading to inaccurate final results. Furthermore, when the field of view is cut, cells located along the cut line are split into multiple pieces, resulting in subsequent cell segmentation errors. This further affects the accuracy of imaging cells at the edges of the field of view. Summary of the Invention

[0005] Based on the technical problems existing in the background technology, the present invention proposes an image processing method, system and medium based on in situ hybridization technology to achieve accurate matching of sequencing data and spatial information, thereby better revealing the spatial distribution pattern of gene expression.

[0006] The image processing method based on in situ hybridization technology proposed by the present invention includes:

[0007] The DAPI staining results of cells were converted from three-dimensional to two-dimensional data using Z-axis maximum intensity projection. Individual cells were identified and separated from the two-dimensional data using cell segmentation. When processing overlapping areas, the Ashlar registration algorithm was used to calculate the error between adjacent fields of view in the same round. The cell segmentation results of each field of view were spliced ​​based on the calculated offset to obtain cell position information.

[0008] The imaging results after hybridization of the fluorescent probe and the gene are registered multiple times using the Ashlar registration algorithm to generate a spliced ​​image. The local maximum of the spliced ​​image at different rounds is marked using the fluorescent spot recognition algorithm to obtain the position information of each fluorescent spot. The position information of each fluorescent spot is then matched with the gene in the transcript to obtain the gene position information.

[0009] The cell position information and the gene position information are used to distribute genes to cells, and a cell gene matrix is ​​obtained and displayed.

[0010] Furthermore, the registration algorithm Ashlar is used to calculate the error between adjacent fields of view in the same round, and the cell segmentation results of each field of view are spliced ​​according to the calculated offset to obtain the cell position information, specifically:

[0011] The two adjacent fields of view are joined after the cells in the overlapping area of ​​the two adjacent fields of view are unioned;

[0012] The spliced ​​visual field is gap-processed to obtain the spliced ​​cells;

[0013] The cell ID offsets before and after splicing were exported for subsequent gene assignment.

[0014] Furthermore, the stitched visual field is subjected to gap processing to obtain the stitched cells, specifically:

[0015] Merge at least partially aligned cells into the same cell;

[0016] For cells that are completely misaligned, the average area of ​​cells within a certain range around the cells is calculated. If the cell is greater than or equal to the average area of ​​the surrounding cells, it is merged into one cell; if the cell is less than or equal to the average area of ​​the surrounding cells, it is discarded.

[0017] Furthermore, the imaging results after hybridization of the fluorescent probe and the gene are subjected to multiple rounds of registration using the registration algorithm Ashlar to generate a spliced ​​image, specifically:

[0018] The different fields of view in the first round are stitched together, and the Ashlar registration algorithm is used to stitch the fields of view together for the overlapping areas;

[0019] The fields of view of subsequent rounds are sequentially registered with the fields of view completed by stitching in the first round of imaging;

[0020] Merge all corrected image positions to generate a stitched image.

[0021] Furthermore, the fluorescent spot recognition algorithm is used to mark the local maximum of the stitched image at different rounds, thereby obtaining the position information of each fluorescent spot, specifically:

[0022] Clip filtering: The clip operation limits the brightness of all fluorescent spot signal images in a round to a set threshold [Clip Min ,Clip Max ] to achieve the purpose of denoising, where Clip Min is the lower limit of brightness, Clip Max is the upper limit of brightness;

[0023] Finding local maxima: Define an a×b window. For each pixel, compare the brightness of each pixel with all other pixels in the corresponding a×b window neighborhood. If the brightness of the pixel is higher than that of all neighboring pixels, it is marked as a local maximum. Repeat all rounds to obtain the position information of each fluorescent point. a and b are the length and width of the window.

[0024] Furthermore, the Clip Min Set to the average brightness of all pixels in the full field of view within a round, Clip Max It is set to twice the 99% quantile of the brightness of all fluorescent spots in M ​​random fields of view within a round.

[0025] Furthermore, the position information of each fluorescent spot is matched with the gene in the transcript to obtain the gene position information, specifically:

[0026] In repeating all rounds to obtain the position information of the fluorescent spot, in each round, the fluorescent spot generates an intensity vector;

[0027] The intensity range of all fluorescent spots in different rounds was calculated and normalized;

[0028] After comprehensive standardization, the multidimensional vector of each point is calculated, the L2 norm of the multidimensional vector is normalized, and the Euclidean distance of the codebook is calculated;

[0029] The code with the closest Euclidean distance is assigned to the fluorescent point. If the signal of the fluorescent point in other rounds is higher than the set intensity threshold, the decoding is abandoned.

[0030] The fluorescent spots are filtered according to the Euclidean distance of each point, and then the duplicate spots are removed. Finally, the duplicate spots in the overlapping areas of adjacent fields of view are removed to obtain the gene location information.

[0031] Furthermore, the intensity threshold is 1 / 4 of the sum of the two brightest intensities of each fluorescent spot as the threshold.

[0032] An image processing system based on in situ hybridization technology includes a cell location information generation module, a gene location information generation module, and a cell gene matrix generation module;

[0033] The cell position information generation module is used to convert the cell DAPI staining results from three dimensions into two dimensions using Z-axis maximum intensity projection, identify and separate individual cells from the two-dimensional data using cell segmentation, and use the Ashlar registration algorithm to calculate the error between adjacent fields of view in the same round when processing overlapping areas of Z-axis maximum intensity projection and cell segmentation. The cell segmentation results of each field of view are spliced ​​according to the calculated offset to obtain cell position information;

[0034] The gene position information generation module is used to perform multi-round registration of the imaging results after hybridization of the fluorescent probe and the gene using the registration algorithm Ashlar to generate a spliced ​​image, and use the fluorescent point recognition algorithm to mark the local maximum of the spliced ​​image in different rounds to obtain the position information of each fluorescent point, and correspond the position information of each fluorescent point to the gene in the transcript to obtain the gene position information;

[0035] The cell gene matrix generation module is used to use the cell position information and the gene position information to distribute genes to cells, obtain a cell gene matrix and display it.

[0036] A computer-readable storage medium stores a plurality of classification programs, wherein the plurality of classification programs are used to be called by a processor and execute the image processing method described above.

[0037] Those skilled in the art will understand that all or part of the steps of implementing the above-mentioned method embodiment can be completed by hardware related to program instructions, and the aforementioned program can be stored in a computer-readable storage medium. When the program is executed, it executes the steps of the above-mentioned method embodiment; and the aforementioned storage medium includes: ROM, RAM, disk or optical disk, etc. Various media that can store program codes.

[0038] The advantages of the image processing method, system and medium based on in situ hybridization technology provided by the present invention are: they can quickly and accurately process the cell data and gene point data obtained by imaging, eliminate the imaging problem of the edge of the field of view after field of view splicing, and finally obtain spatial transcriptome information, realize the accurate matching of sequencing data and spatial information, and thus better reveal the spatial distribution pattern of gene expression; additionally introduce the local maximum point recognition algorithm (i.e., fluorescent point recognition algorithm) to overcome the traditional LOG (Laplacian of The algorithm solves the problem of long running time of the Gaussian algorithm by greatly improving experimental efficiency and making the processing of large-scale data sets more feasible. The brightness threshold range is limited by the clip filtering operation, and a set window is defined to find the local maximum. This algorithm can significantly reduce calculation time while maintaining detection accuracy. The registration algorithm Ashlar is used to calculate the relative offset between adjacent fields of view in the same round, and then the cell segmentation results of each field of view are adjusted and stitched based on these offsets. The registration algorithm Ashlar can eliminate the horizontal offset between adjacent fields of view due to the accuracy limitation of the microscope movement module, which causes cell overlap, loss or morphological abnormalities at the edge of the field of view in the stitched image. In particular, when processing overlapping areas, partially aligned cells are merged, and the retention of small cells that cannot be aligned is determined based on the average area of ​​the surrounding cells, thereby effectively solving the edge problem of repeated fields of view. This not only improves the accuracy of the final imaging results, but also ensures that subsequent steps such as gene assignment can be performed based on more accurate cell position information, further improving the reliability of the overall experiment. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 It is a schematic diagram of the process of the present invention;

[0040] Figure 2 It is a flowchart corresponding to the image processing method;

[0041] Figure 3 Schematic diagram of the cell segmentation processing strategy using the registration algorithm Ashlar;

[0042] Figure 4 Schematic diagram of 5 types of splicing;

[0043] Figure 5 Schematic diagram of gene quality quantification results;

[0044] Figure 6 Schematic diagram comparing high-quality gene processing results with the Allen brain atlas. DETAILED DESCRIPTION

[0045] The technical solutions of the present invention are described in detail below through specific embodiments. Numerous specific details are set forth in the following description to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways than those described herein, and those skilled in the art may make similar modifications without departing from the scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0046] like Figures 1 to 5 As shown, the image processing method based on in situ hybridization technology proposed by the present invention includes:

[0047] Step 1: Convert the cell DAPI staining results from 3D to 2D data using Z-axis maximum intensity projection. Use cell segmentation to identify and separate individual cells from the 2D data. When processing overlapping areas between Z-axis maximum intensity projection and cell segmentation, use the Ashlar registration algorithm to calculate the error between adjacent fields of view in the same round. Use the calculated offset to stitch the cell segmentation results of each field of view to obtain cell position information.

[0048] Step 2: The imaging results after hybridization of the fluorescent probe and the gene are registered multiple times using the Ashlar registration algorithm to generate a spliced ​​image. The local maximum of the spliced ​​image at different rounds is marked using the fluorescent point recognition algorithm to obtain the position information of each fluorescent point. The position information of each fluorescent point is then matched with the gene in the transcript to obtain the gene position information;

[0049] Step 3: Using the cell location information and the gene location information, genes are assigned to cells to obtain and display a cell gene matrix.

[0050] This embodiment can quickly and accurately process the cell data and gene point data obtained by imaging, and eliminate the imaging problem of the edge of the field of view after field of view stitching, and finally obtain spatial transcriptome information, realize the precise matching of sequencing data and spatial information, and thus better reveal the spatial distribution pattern of gene expression.

[0051] By introducing an additional local maximum point recognition algorithm (i.e., fluorescent point recognition algorithm), it aims to overcome the traditional LOG (Laplacian of The problem of long running time of the Gaussian algorithm is greatly improved, making the processing of large-scale data sets more feasible. The brightness threshold range is limited by the clip filtering operation, and a set window is defined to find the local maximum. The algorithm can significantly reduce the calculation time while maintaining detection accuracy; the registration algorithm Ashlar is used to calculate the relative offset between adjacent fields of view in the same round, and then the cell segmentation results of each field of view are adjusted and spliced ​​according to these offsets (see step (a9) for details). The registration algorithm Ashlar can eliminate the horizontal offset between adjacent fields of view due to the accuracy limitation of the microscope movement module, which causes the problem of cell overlap, loss or morphological abnormalities at the edge of the field of view in the spliced ​​image; in particular, when processing overlapping areas, partially aligned cells will be merged, and small cells that cannot be aligned will be retained based on the average area of ​​the surrounding cells, thereby effectively solving the edge problem of repeated fields of view, not only improving the accuracy of the final imaging results, but also ensuring that subsequent steps such as gene assignment can be performed based on more accurate cell position information, further improving the reliability of the overall experiment.

[0052] In this example, three functional modules are mainly set up: 1> Processing the cell DAPI (4',6-diamidino-2-phenylindole) staining results to obtain cell location information. 2> Processing the imaging results after the fluorescent probe and gene hybridization to obtain gene location information. 3> Using the singleton mode combined with the Snakemake framework to synchronize and control the execution of 1> and 2>, and combining the cell location information with the gene location information, assigning genes to cells, obtaining a cell gene matrix and displaying the data quality.

[0053] In one embodiment, 1> the DAPI (4',6-diamidino-2-phenylindole) staining results of the cells are processed to obtain cell location information, that is, step 1 is specifically:

[0054] Various spatial transcriptomics technologies process DAPI staining results in two steps: Z-axis maximum intensity projection (Z-axis maximum intensity projection) and cell segmentation. Z-axis maximum intensity projection converts the original three-dimensional imaging data (i.e., the three-dimensional image of the DAPI staining results) into two-dimensional data. The cell segmentation algorithm accurately identifies and isolates individual cells in the DAPI staining image, facilitating subsequent gene assignment and other operations. This example builds on these two steps by adding offset correction and splicing for the cell segmentation results.

[0055] When processing overlapping fields of view, the cell segmentation results in the overlapping areas will have overlapping, abnormal shapes or missing problems. These problems are collectively referred to as edge problems with overlapping fields of view. The edge problem is caused by two reasons. One is that when the field of view is cut, the cells located on the cutting line will be cut into multiple parts. Then, incorrect segmentation will occur during cell segmentation. The second reason is that the original process did not take into account the error during the movement of the imaging stage, which will cause offsets between adjacent fields of view, which will aggravate the edge problem. In order to solve the edge problem, this embodiment proposes an offset correction and stitching algorithm.

[0056] Specifically, we first need to calculate the relative offset of the field of view to eliminate the error when the imaging stage moves. Figure 3 As shown, this embodiment uses the registration algorithm Ashlar to calculate the error between adjacent fields of view in the same round. Then, the cell segmentation results of each field of view are spliced ​​according to the calculated offset, specifically:

[0057] Since the cell segmentation results of the overlapping areas of two adjacent fields of view may not be exactly the same, edge problems may occur at the junction of the union and the two fields of view. The last step is to perform gap processing, and finally export the changes in cell IDs before and after splicing for subsequent gene assignment.

[0058] Specifically, the registration algorithm Ashlar is as follows: the cells in the overlapping area of ​​two adjacent fields of view are unioned and then the two adjacent fields of view are spliced ​​together; the spliced ​​fields of view are gap processed to obtain spliced ​​cells; and the cell ID offsets before and after splicing are derived for subsequent gene allocation.

[0059] In this embodiment, the stitched field of view is gap-processed to obtain stitched cells, specifically: cells that are at least partially aligned are merged into the same cell; for cells that are completely unaligned, the average area of ​​cells within a certain range around the cell is calculated, and if the cell is greater than or equal to the average area of ​​the surrounding cells, it is merged into one cell; if the cell is less than or equal to the average area of ​​the surrounding cells, it is discarded.

[0060] In one embodiment, 2> processing the imaging result after hybridization of the fluorescent probe and the gene to obtain gene location information, that is, step 2 specifically includes:

[0061] The process of processing the imaging results after hybridization of fluorescent probes and genes is mainly divided into four parts: multi-round registration, fluorescent spot identification, decoding, and offset correction and field of view stitching for fluorescent staining results. Specifically:

[0062] 2.1 Multi-round Registration: In practice, factors such as sample handling, slide movement, and optical system variations can cause positional deviations between data from different rounds. In this example, the Ashlar algorithm was used for multi-round registration. The process is divided into the following three parts.

[0063] 2.1.1> Stitch the different fields of view in the first round, and use the registration algorithm Ashlar to stitch the fields of view in the overlapping areas.

[0064] 2.1.2> Sequentially align the fields of view of subsequent rounds with the field of view completed by the first round of imaging stitching.

[0065] 2.1.3>Merge all corrected image positions to generate a stitched image.

[0066] Through 2.1.1> to 2.1.3>, multi-run registration was used to ensure that data acquired from different experimental runs could be accurately aligned and compared.

[0067] 2.2>The fluorescent spot recognition algorithm can calculate the position information of each fluorescent spot (points*position matrix) based on the fluorescent spot signal image.

[0068] The current mainstream point recognition algorithm is the LOG (Laplacian of Gaussian) algorithm. The basic idea of ​​the LOG algorithm is to first apply a Gaussian filter to the image to remove noise, and then apply the Laplace transform to highlight edges and points in the image. The general steps are as follows: First, a series of different scales (standard deviation σ) are selected. For each scale, a Gaussian-Laplacian filter is calculated at that scale. The input image (i.e., the stitched image) is then convolved with these filters at different scales, resulting in a series of scale-space response maps. Local maxima are then found in the response map at each scale, and weak response points are filtered out using thresholding.

[0069] The LOG algorithm combines the advantages of Gaussian smoothing and Laplace transform, and can effectively detect point structures at different scales. However, in actual use, the LOG algorithm takes too long to run, and the data of thousands of fields of view and thousands of genes need to be processed for 7 days, resulting in a poor actual use experience. Therefore, this embodiment develops a point recognition algorithm based on local maximum (i.e., fluorescent point recognition algorithm), which removes noise by clip screening and searches for maximum points in local areas. For a batch of 1120 fields of view, 46 rounds, and 3-channel data, running on a single core on an Intel Xeon CPU E7-8860 v4 2.20GHZ, the fluorescent point recognition algorithm can efficiently process large-scale data within 24 hours, and the efficiency is 7 times higher than that of the LOG algorithm.

[0070] The fluorescent spot recognition algorithm includes the following steps 2.2.1> to 2.2.2>:

[0071] 2.2.1> Clip screening: The clip operation limits the brightness of all fluorescent spot signal images in a round to a set threshold [Clip Min ,Clip Max ] to achieve the purpose of denoising. Among them, Clip Min is the lower limit of brightness, Clip Max Is the brightness limit, Clip Min Set to the average brightness of all pixels in the full field of view within a round, Clip Max It is set to twice the 99% quantile of the brightness of all fluorescent spots in M ​​random fields (for example, M=30) within a round.

[0072] 2.2.2> Finding local maxima: The present invention defines an a×b (e.g., 5×5) window and compares the brightness of each pixel with all other pixels in the corresponding a×b window neighborhood. If the brightness of the pixel is higher than that of all neighboring pixels, it is marked as a local maximum. All rounds are repeated to obtain the position information of each fluorescent point, where a and b are the length and width of the window.

[0073] Using the same computing resources, this embodiment only takes 24 hours to process the 1,000-view 1,000-gene data using the fluorescent spot recognition algorithm. Compared to the LOG algorithm, this embodiment's fluorescent spot recognition algorithm shortens the spot recognition process time by 7 times.

[0074] 2.3> The decoding algorithm is responsible for matching the identified points with the genes in the transcript and calculating the gene position information (gene*position matrix). The decoding algorithm has the following steps:

[0075] 2.3.1> In the process of repeating all rounds to obtain the position information of the fluorescent spot, in each round, the fluorescent spot generates an intensity vector, and the intensity range of the spot in different rounds is calculated and normalized;

[0076] 2.3.2> After comprehensive standardization, calculate the L2 norm of each point and normalize the multidimensional vector and the Euclidean distance of the codebook;

[0077] 2.3.3> Take the code with the closest Euclidean distance and assign it to the fluorescent point. If the signal of other rounds of the fluorescent point is higher than the set intensity threshold, the decoding is abandoned.

[0078] 2.3.4> Filtering: Filter the fluorescent spots according to the Euclidean distance of each spot, then remove duplicate spots, and finally remove duplicate spots in the overlapping areas of adjacent fields of view to obtain gene location information.

[0079] When processing the imaging results after hybridization of fluorescent probes and genes, we also face the problem of offset between adjacent fields of view mentioned in the process of generating cell position information. Therefore, the offset between adjacent fields of view calculated in the same round is also used in the process of generating gene position information, and the decoded gene position information is shifted according to this offset.

[0080] In one embodiment, 3> creates a global variable corresponding to the global resource in the instance of the singleton pattern, and then uses the workflow of the Snakemake framework to combine the global variable to control the execution order of 1> and 2>, and combines the cell position information with the gene position information, assigns the gene to the cell, obtains the cell gene matrix and displays the data quality. That is, step three is specifically as follows:

[0081] After obtaining gene location information and cell location information, this embodiment distributes genes to cells to obtain a cell gene matrix. This embodiment adds a function for calculating indicators of basic spatial transcriptome technology. The following indicators can be displayed:

[0082] 3.1>Number of genes;

[0083] 3.2>Number of genes after decoding filtering / Number of genes before filtering;

[0084] 3.3>Gene quality calculated based on the above indicators;

[0085] In view of the use of the above three main functions, this embodiment has developed a command line system to provide the above functions in order to more efficiently utilize computing resources and improve operating speed. The development details are as follows: Due to the need to control resource access and implement global state management, the singleton mode is chosen, and the development is based on the Snakemake framework, and the project initialization is completed in the form of incoming JSON. A core data structure Bassdata and two main threads are defined to process DAPI staining results and fluorescent point signals respectively. Bassdata records the main data path, the algorithm parameters entered by the user, and can also synchronize the execution progress between the two main threads, allocate available computing resources, and regulate the execution logic of the entire process.

[0086] Example 1

[0087] (a1) Initialize Bassdata and Z-axis maximum projection: Input the DAPI staining results of cells and the imaging results after hybridization of fluorescent probes with genes, the designed codebook, and the configuration JSON file. The configuration JSON file contains the input and output paths and key parameters required for running the entire pipeline: set the size to 2150*2150 pixels; set the overlap of adjacent fields of view to 5%; set the DAPI channel to 4001, the fluorescent spot signal channels to 5501 and 7401; set the number of fields of view to 1120; and set the number of rounds to 46. Then, create the specified folder according to the output path in the configuration JSON. Then, execute the Z-axis maximum projection. The JSON file sets the valid Z-axis imaging range to stack2 to stack14. Then, within stack2-stack14, the maximum value of each pixel at position (i, j) along the Z axis is calculated and converted to uint16 format, where (i, j) is the horizontal and vertical coordinate of the pixel point. During execution, multi-threaded computing is performed using the parallelism specified in the JSON file, with a default setting of 20 threads. 1120 views takes 20 hours.

[0088] (a2) Registration Algorithm: Ashlar multi-pass registration: First, rename the data from different channels and different passes of the same field of view according to the Ashlar format and enter a specific path. Then, write an Ashlar script for each field of view path and execute it. During execution, multi-threaded operations are performed using the Ashlar parallelism specified in the JSON file, with a default setting of 20 threads. Registering 1120 fields of view takes approximately 4 hours.

[0089] (a3) Cell segmentation: The cell segmentation algorithm (deepcell) configuration JSON is used to configure the parameters required for deepcell operation: set the segmentation type to nuclear; set the input data resolution image_mpp value to between 0.2 and 0.4, the cell smoothness interior_model_smooth value to between 0.2 and 0.4, the model estimate of the conservative value of the cell and background interior_threshold value to between 0.5 and 0.7, and the cell radius value radius value to between 0.2 and 0.4. Then, the built-in trained deep convolutional neural network (DCNN) is used to automatically learn the mapping relationship between the input image and the output segmentation map, thereby completing the cell segmentation and outputting the cell segmentation results in the specified path. Each cell ID uses the cell's grayscale value as a unique identifier. The default setting is 20 threads of parallelism. It takes approximately 4 hours to perform cell segmentation on 1120 fields of view.

[0090] (a4) Clip screening: Calculate the average brightness of all pixels in all fields of view within a round and set this value as Clip Min ; Randomly select n fields of view to perform point recognition, and then calculate twice the 90% to 99% percentile of the brightness of all points in these n fields of view as Clip Max . Then the clip interval of each round [Clip Min ,Clip Max ] is saved to clip.csv. Then, for each image in each round, the corresponding clip interval is read and its pixel value is limited to the clip interval range. Finally, the clip-filtered image is saved in the specified path.

[0091] (a5) SpaceTX format conversion: The SpaceTX format is a data format designed by Starfish technology. Since Starfish will be used for subsequent fluorescence spot signal decoding operations, the data format needs to be converted to SpaceTX format. First, rename all images to the specified format of "field number-round number-channel number.tif", then extract the file coordinate information based on the image file name and save it to the file coordinates.csv. Then, based on the information in coordinates.csv, call Starfish's format_structured_dataset formatting function to generate files such as experiment.json, nuclei.json, and primary.json. These files will be used in the subsequent decoding process.

[0092] (a6) Local maximum point identification: For each pixel of all images after clip filtering, compare its brightness value with that of all other pixels in its 8x8 neighborhood. If the brightness value of the pixel is higher than that of all its neighbors, it is marked as a local maximum. Set the deduplication radius sz value. If multiple local maxima are within sz / 2 jump, only the brightest one is retained. Then create an image containing only the maximum points, and around each point The maximum value within the region is found, and the remaining values ​​are set to zero, leaving only the strongest peak within each region. The data is then converted to (x,y) format and saved in the blob_result.csv file. The default setting is 20 threads in parallel. Point recognition for 46 rounds, 2 channels, and 1120 fields of view takes approximately 20 hours.

[0093] (a7) Decoding: First, load the fluorescent spot identification results from the blob_result.csv file. Then traverse all the point signal images, and for each identified fluorescent spot in each image, calculate the average brightness of its center point and the surrounding neighborhood. Since a total of 46 rounds of imaging are performed, a 46-dimensional intensity vector is generated for each fluorescent spot. Then decode the fluorescent spot according to the code book. After comprehensive standardization, calculate the L2 norm of the multidimensional vector and normalize it as well as the Euclidean distance of the code book, and assign the code with the closest Euclidean distance to the fluorescent spot. Take 1 / 4 of the sum of the intensities of the brightest rounds of each point as the intensity threshold. If the signal of other rounds of the point is higher than the intensity threshold, the decoding is abandoned. Finally, save the decoding results to the gene name.parquet file.

[0094] (a8) Filtering: For each gene, load the corresponding parquet file. First, use the Gaussian Mixture Model (GMM) to model the Euclidean distance column in the parquet file, then find the number of clusters that satisfies the minimum Bayesian Information Criterion (BIC), and then fit the model. Then, for each cluster, determine whether the average distance of the cluster is greater than 0.5. If it is greater than 0.5, the points of the cluster are considered to be background points; if it is less than 0.5, the points of the cluster are considered to be foreground points and saved to the filtered .parquet file. Since there are overlapping areas between adjacent fields of view, it is necessary to deduplicate the points in the overlapping areas. First, traverse all fields of view and calculate the global coordinates of each fluorescent point in the entire image. Then, according to the provided field of view position table fov_position, mark the fluorescent points that need to be deleted in the overlapping area of ​​adjacent fields of view to avoid duplication. Then, it filters out these marked fluorescent points and saves the processed data to the deduplicated .parquet file.

[0095] (a9) Offset correction and stitching: First, all fields of view are divided into five stitching types (e.g. Figure 4), which are the five cases where the field of view to be spliced ​​is located to the right, left, bottom, lower right, and lower left of the spliced ​​field of view. Then, different Ashlar scripts are used according to different splicing types to calculate the relative offset between all fields of view. Then, the absolute position of the corrected field of view is calculated based on the relative offset between fields of view and the size of a single field of view. For the fluorescent point signal results, you can just move according to the absolute position; for the cell segmentation results, you need to take the union of the cells in the overlapping areas between adjacent fields of view, and then merge the partially aligned cells into the same cell; for cells that are completely unaligned, the average area of ​​the cells within a certain range around them will be calculated. If the cell is greater than or equal to the average area of ​​the surrounding cells, it will be retained; if the cell is less than or equal to the average area of ​​the surrounding cells, it will be discarded. The cell id offset before and after splicing is then derived for subsequent gene allocation. (such as Figure 3 The default setting is 20 threads in parallel. Recognizing points in 46 rounds and 1120 fields of view takes about 5 hours.

[0096] (a10) Gene assignment: Use the unique function to obtain the area and position of each cell from the spliced ​​cell segmentation results, and then assign the gene information in the .parquet to the cells at the corresponding positions.

[0097] (a11) Effect evaluation: Gene quality is calculated by counting the ratio of the number of all genes to the number of genes after filtering / the number of genes before filtering. Output a scatter plot, such as Figure 5 As shown, the x-axis is the percentage of genes retained after filtering in the original data, and the y-axis is the number of genes. After calculating the top 10 horizontal / vertical axes and taking the intersection, we can get the highest quality genes and display them in different colors. Figure 5 At the same time, this embodiment can also display some genes with obvious spatial specificity on this map, using different colored dots as markers, such as Figure 5 shown. Figure 6 The spatial expression results of the Allen Brain Atlas and the high-quality gene Gria1, both from the Allen Brain Atlas and this example, are shown. The left graph shows the expression results of the Allen Brain Atlas, while the right graph shows the expression results of this example. It can be seen that the results of this example are almost identical to those of the Allen Brain Atlas.

[0098] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with the technical field, within the technical scope disclosed by the present invention, who makes equivalent replacements or changes based on the technical solution and inventive concept of the present invention, should be covered by the scope of protection of the present invention.

Claims

1. An image processing method based on in situ hybridization technology, characterized in that: include: The DAPI staining results of cells were converted from three-dimensional to two-dimensional data using Z-axis maximum intensity projection. Individual cells were identified and separated from the two-dimensional data using cell segmentation. When processing overlapping areas, the Ashlar registration algorithm was used to calculate the error between adjacent fields of view in the same round. The cell segmentation results of each field of view were spliced ​​based on the calculated offset to obtain cell position information. The imaging results after hybridization of the fluorescent probe and the gene are registered multiple times using the Ashlar registration algorithm to generate a spliced ​​image. The local maximum of the spliced ​​image at different rounds is marked using the fluorescent spot recognition algorithm to obtain the position information of each fluorescent spot. The position information of each fluorescent spot is then matched with the gene in the transcript to obtain the gene position information. The cell position information and the gene position information are used to distribute genes to cells, and a cell gene matrix is ​​obtained and displayed.

2. The image processing method according to claim 1, wherein: The registration algorithm Ashlar is used to calculate the error between adjacent fields of view in the same round, and the cell segmentation results of each field of view are spliced ​​according to the calculated offset to obtain the cell position information, specifically: The two adjacent fields of view are joined after the cells in the overlapping area of ​​the two adjacent fields of view are unioned; The spliced ​​visual field is gap-processed to obtain the spliced ​​cells; The cell ID offsets before and after splicing were exported for subsequent gene assignment.

3. The image processing method according to claim 2, wherein: The stitched visual field is subjected to gap processing to obtain stitched cells, specifically: Merge at least partially aligned cells into the same cell; For cells that are completely misaligned, the average area of ​​cells within a certain range around the cells is calculated. If the cell is greater than or equal to the average area of ​​the surrounding cells, it is merged into one cell; if the cell is less than or equal to the average area of ​​the surrounding cells, it is discarded.

4. The image processing method according to claim 1, wherein: The imaging results after hybridization of the fluorescent probe and the gene are subjected to multiple rounds of inter-registration using the registration algorithm Ashlar to generate a spliced ​​image, specifically: The different fields of view in the first round are stitched together, and the Ashlar registration algorithm is used to stitch the fields of view together for the overlapping areas; The fields of view of subsequent rounds are sequentially registered with the fields of view completed by stitching in the first round of imaging; Merge all corrected image positions to generate a stitched image.

5. The image processing method according to claim 4, characterized in that The fluorescent spot recognition algorithm is used to mark the local maximum of the stitched image at different rounds, thereby obtaining the position information of each fluorescent spot, specifically: Clip filtering: The clip operation limits the brightness of all fluorescent spot signal images in a round to a set threshold [Clip Min ,Clip Max ] to achieve the purpose of denoising, where Clip Min is the lower limit of brightness, Clip Max is the upper limit of brightness; Finding local maxima: Define an a×b window. For each pixel, compare the brightness of each pixel with all other pixels in the corresponding a×b window neighborhood. If the brightness of the pixel is higher than that of all neighboring pixels, it is marked as a local maximum. Repeat all rounds to obtain the position information of each fluorescent point. a and b are the length and width of the window.

6. The image processing method according to claim 5, characterized in that The Clip Min Set to the average brightness of all pixels in the full field of view within a round, Clip Max It is set to twice the 99% quantile of the brightness of all fluorescent spots in M ​​random fields of view within a round.

7. The image processing method according to claim 5, characterized in that: The position information of each fluorescent spot is matched with the gene in the transcript to obtain the gene position information, specifically: In repeating all rounds to obtain the position information of the fluorescent spot, in each round, the fluorescent spot generates an intensity vector; The intensity range of all fluorescent spots in different rounds was calculated and normalized; After comprehensive standardization, the multidimensional vector of each point is calculated, the L2 norm of the multidimensional vector is normalized, and the Euclidean distance of the codebook is calculated; The code with the closest Euclidean distance is assigned to the fluorescent point. If the signal of the fluorescent point in other rounds is higher than the set intensity threshold, the decoding is abandoned. The fluorescent spots are filtered according to the Euclidean distance of each point, and then the duplicate spots are removed. Finally, the duplicate spots in the overlapping areas of adjacent fields of view are removed to obtain the gene location information.

8. The image processing method according to claim 7, wherein: The intensity threshold is 1 / 4 of the sum of the two brightest intensities of each fluorescent spot as the threshold.

9. An image processing system based on in situ hybridization technology, characterized in that: It includes a cell position information generation module, a gene position information generation module and a cell gene matrix generation module; The cell position information generation module is used to convert the cell DAPI staining results from three dimensions into two dimensions using Z-axis maximum intensity projection, identify and separate individual cells from the two-dimensional data using cell segmentation, and use the Ashlar registration algorithm to calculate the error between adjacent fields of view in the same round when processing overlapping areas of Z-axis maximum intensity projection and cell segmentation. The cell segmentation results of each field of view are spliced ​​according to the calculated offset to obtain cell position information; The gene position information generation module is used to perform multi-round registration of the imaging results after hybridization of the fluorescent probe and the gene using the registration algorithm Ashlar to generate a spliced ​​image, and use the fluorescent point recognition algorithm to mark the local maximum of the spliced ​​image in different rounds to obtain the position information of each fluorescent point, and correspond the position information of each fluorescent point to the gene in the transcript to obtain the gene position information; The cell gene matrix generation module is used to use the cell position information and the gene position information to distribute genes to cells, obtain a cell gene matrix and display it.

10. A computer-readable storage medium, characterized in that The computer-readable storage medium stores a plurality of classification programs, which are used to be called by a processor and execute the image processing method according to any one of claims 1 to 8.

Images