BINGO: a blind unmixing algorithm for ultra-multiplexing fluorescence images

Abstract Motivation Spectral imaging is often used to observe different objects with multiple fluorescent labels to reveal the development of the biological event. As the number of observed objects increases, the spectral overlap between fluorophores becomes more serious, and obtaining a “pure” picture of each fluorophore becomes a major challenge. Here, we propose a blind spectral unmixing algorithm called BINGO (Blind unmixing via SVD-based Initialization Nmf with project Gradient descent and spare cOnstrain), which can extract all kinds of fluorophores more accurately from highly overlapping multichannel data, even if the spectra of the fluorophores are extremely similar or their fluorescence intensity varies greatly. Results BINGO can isolate up to 10 fluorophores from spectral imaging data for a single excitation. nine-color living HeLa cells were visualized distinctly with BINGO. It provides an important algorithmic tool for multiplex imaging studies, especially in intravital imaging. BINGO shows great potential in multicolor imaging for biomedical sciences. Availability and implementation The source code used for this paper is available with the test data at https://github.com/Xinyuan555/BINGO_unmixing


trial step calculation：Update W and H
According to linear mixing model of spectra recognized by fluorescence imaging over the years (Zimmermann 2005, Zimmermann et al. 2014, Rakhymzhan et al. 2017), the total detected fluorescence signal for every channel  can be written as I(λ), I(λ) =  1  1 () +  2  2 () + ⋯ +     (). (1) () represents the spectral contribution of the fluorophores to every channel, and Ax represents the abundances of the fluorophores in the measured spot.
On the other hand, BINGO was designed with NMF, the decomposition structure of NMF is where, X = [I(λ 1 ), I(λ 2 ), ⋯ , I(λ  )], X is a flatten multi-channel image matrix with a scale of  × , m denotes pixel number of each image and n denotes number of channels.[I(λ 1 ), I(λ 2 ), ⋯ , I(λ

Supplementary Note 2 NNDSVD Initialization
Algorithm 2 shows the steps of NNDSVD Initialization.Input the image matrix X and the number of fluorophores to be solved, r, to output the initialization value W0 and H0.At first, it performs a singular value decomposition of the image matrix X with r as the rank, and the goal of the second SVD is to calibrate all elements of the matrix to non-negative, which mainly based on equation (4).
The principle of the second SVD is as follows.The elements in   and   include both positive and negative elements.u+ is written as all the elements of ≥ 0 in vector u, and the rest are replaced by 0. u-is written as all the elements of < 0 in vector u, and the rest are replaced by 0. Similarly for vector v, at this point, After ignoring  − , all elements less than 0 can be corrected. � ± ∶=  ± / � ± �， � ± ∶=  ± / � ± �，  ± = � ± �� ± �，here [ � + ,  � − ] and [ � + ,  � − ] became unitary matrix. ± were singular value.writing  + with equation (4)， Since rank(C)=1, rank( + )is at most 2, we can approximate  + by the larger term on the righthand side of Eq. ( 6).The initial non-negative matrices W0 and H0 are obtained after a low-rank approximation for each of the latter terms.Although NNDSVD gives up part of the similarity in the process of calibrating negative elements, the features extracted by SVD are retained to the maximum extent, which can be compensated in the subsequent iterative optimization.
where p is the number of pixels in each image, m is the number of image channels, and this reference value is proposed in the hyperspectral decomposition to be less than 0.1.

SAD
For blind source separation (BSS), the extracted spectrum can be compared with the actual measured spectrum and the SAD (spectral angular distance) is calculated to assess the similarity of the two spectra.
The larger the difference between the two spectra, the larger the corresponding spectral angular distance value, and the smaller the value indicates that the two spectra are more similar.For two vectors a and b, For each fluorophore, a SAD value is calculated, and the average value of all fluorophores SAD ������ is often used to measure the overall capability of the spectral unmixing algorithm during a single experiment.

Correlation coefficient
The 2D correlation coefficient is usually used to examine the similarity between two images, which can be calculated by the corrcofe function in MATLAB.The two-dimensional correlation numbers between different channel images are plotted as a heat map (Fig1b).The horizontal and vertical coordinates of the heat map all correspond to fluorophore types, and are represented by 1, 2, 3, ..., which can visualize the crosstalk between images.the 2D correlation between each channel image and itself is 1, and the value of 2D correlations with other channels should be close to 0. Calculating the 2D correlation coefficients between different channel images can quantify the degree of crosstalk between all channel images.

SSIM
For data that can provide ground truth, the SSIM (Structural Similarity) between the unmixing result and the standard image can be calculated, which can be calculated by the ssim function in MATLAB.
SSIM is a widely used image quality evaluation indicator that evaluates the similarity between two images based on the assumption that the human eye extracts structured information from them when viewing an image.As with SAD, an SSIM value is calculated for each fluorophore, and the average value of all fluorophores  ������� is generally used to measure the overall capability of the spectral unmixing algorithm.
the SSIM takes a value between [0,1], and the closer to 1 means the more similar the image is, the more accurate the unmixing result is.

Dice
After binarizing the fluorophore grayscale image, structural information such as area and volume of the target can be obtained by counting the number of pixels.The Dice coefficient is an ensemble similarity measure function that is usually used to calculate the similarity of two samples and takes the value of [0,1].
Unlike SSIM, which can measure similarity for grayscale images, the Dice coefficient can only be used for binarized images.If the Dice coefficient between the unmixing result after binarization and the standard image is large, it means that the result of the later image analysis is closer to the true value and the result is more reliable.
the paper, take Fluo 1(blue solid line) and Fluo 2(green solid line) for example, both spectra have been normalized itself respectively.Peak gap means the distance between the wavelengths corresponding to the their emission peak.Spectral overlap means the ratio between the area of the overlapping region(diagonal stripe region) and the area of the spectrum itself(blue region).bEmission spectrum of Fluo 2 in Fig. 2 and Fig. 3 before(solid line) and after(dotted line) spectral shift.Simulation included peak shift and boarden, the corresponding wavelength of the peak shifted to shorter wavelength about 5 nm, and the FWHM of it boarden to about 1.2 times.c Single-channel image of ground truth, sythetic images, and the results of different algorithms.Scale bar: 100 μm.
fluorophore were all set about 0~ 8 nm to shorter or longer wavelength randomly, and broaden about 1~1.5 times randomly.c Detailed parameters used in Fig. 3. d Synthetic image of 3 color image.e SAD ���� value while increasing the number of fluorophores from 3 to 11. BINGO achieved the smaller SAD ���� value compared to the other algorithms, indicating that its estimated spectra are closer to those after the artificial spectral shift.f Single-channel images and 2D-corr heat map of three fluorophores in Fig. 3ai after unmixing with different algorithms.

Supplementary
W and H are called the abundance and feature matrices,  ∈  ×   ∈  × , respectively.What r represents is the number of fluorophores in the sample,  ≤ min (, ).The abundance matrix W, as the name suggests, W = [A 1 , A 2 , ⋯ , A  ] , represents the spatial concentration distribution of each fluorophore, and the corresponding image of each fluorophore is obtained after rearranging the matrix to an image; H = [ 1 ,  2 ,   ]  , the feature matrix H reflects the intensity distribution of each fluorophore in each channel, corresponding to its estimated spectrum.