Anti-noise bi-clustering method and system for gene expression data

Abstract

According to the anti-noise bi-clustering method and system for the gene expression data, the thought that a seed list is obtained through a longest common subsequence searching method is utilized, the longest approximate common subsequence is searched through an element exchange permitting strategy, interference of noise in a data set on sequence matching can be effectively reduced, and the anti-noise bi-clustering method and system for the gene expression data are provided. The algorithm can still accurately identify the similarity in the gene expression mode in a complex noise environment, so that the accuracy of a bi-clustering result is improved; the binary search algorithm is adopted to determine the insertion interval, the method can fully consider the characteristics of noise data, and compared with a traditional expansion mode, the method can more reasonably process the influence of noise and avoid error expansion caused by the noise, so that the bi-clustering structure is more accurately constructed.

Description

Technical Field

[0001] The present invention belongs to the technical field related to gene expression data analysis, and in particular relates to a noise-resistant biclustering method and system for gene expression data. Background Art

[0002] The statements in this section merely provide background information related to the present invention and do not necessarily constitute prior art.

[0003] With the rapid development of sequencing technology in recent years, a vast amount of gene expression data has been generated, posing new opportunities and challenges for its analysis. The introduction of whole-genome transcriptome sequencing has made it possible to detect gene modules from gene expression data. In the context of bioinformatics, gene modules are defined as groups of genes with similar expression profiles. These genes are often functionally related and co-regulated. Because transcriptional regulation is highly context-specific, detecting gene modules under specific conditions is crucial. The advent of high-throughput sequencing technologies has made it possible to simultaneously measure whole-genome expression data from many more samples, enabling more efficient inference of gene modules under specific conditions. For example, scRNA-seq data reveals cellular heterogeneity, enabling the inference of cell-type-specific gene modules. However, high-throughput sequencing often generates a significant amount of systematic background noise. This noise can obscure important biological information and interfere with the inference of gene modules, making the identification of gene modules challenging.

[0004] Inferring conditional gene modules can be accomplished by identifying a set of genes that exhibit significant similarity under certain conditions within complex gene expression data. From a computational biology perspective, this process can be formulated as a biclustering problem. Biclustering algorithms can simultaneously cluster rows representing genes and columns representing experimental conditions in a gene expression matrix, effectively capturing co-expression relationships between genes under specific conditions and providing a viable approach for discovering biologically meaningful gene modules. However, the high noise levels in high-throughput gene expression data can interfere with biclustering analysis results, so biclustering algorithms applied to such data must exhibit strong noise immunity.

[0005] Since the algorithm proposed by Cheng and Church was applied to gene expression data in 2000, numerous biclustering algorithms have been developed. Some of these algorithms focus on addressing the problem of high noise in the data to improve the accuracy of biclustering analysis. The QUBIC2 algorithm uses the LTMG algorithm to discretize expression data to eliminate the effects of noise. However, this algorithm has certain limitations. On the one hand, the discretization process is computationally complex, resulting in long processing times. On the other hand, its performance is sensitive to the discretization level parameter set by the user, and different parameter settings can lead to significantly different analysis results. The ARBic algorithm claims to be robust to noise and can theoretically maintain good performance in noisy environments. However, in practice, when processing large-scale expression data, especially high-throughput data, the algorithm's operational efficiency is low and the runtime is significantly increased, which to some extent limits its application in large-scale data scenarios.

[0006] In view of the shortcomings of existing algorithms in dealing with noise, operating efficiency and parameter sensitivity, the development of a new biclustering algorithm has important practical significance. Summary of the Invention

[0007] To overcome the deficiencies of the above-mentioned prior art, the present invention provides a noise-resistant biclustering method and system for gene expression data. The scheme of the present invention can complete calculations within a reasonable time range, effectively deal with the interference of noise on the biclustering results, and at the same time has strong robustness to model parameters, ensuring stable performance under different parameter settings, providing a more reliable and efficient method for the analysis of gene expression data.

[0008] In order to achieve the above object, the present invention adopts the following technical solutions:

[0009] In a first aspect, the present invention provides a noise-resistant biclustering method for gene expression data, comprising:

[0010] Normalizing the gene expression matrix to obtain an index matrix sorted in ascending order, and calculating the coefficient of variation of each row in the index matrix;

[0011] Based on the longest common subsequence algorithm, in the index matrix, if the elements at two index positions are within the noise perturbation range and the error range determined by the corresponding coefficient of variation, elements are allowed to be exchanged to find the longest approximate common subsequence, and the positive and negative correlations between genes are considered to determine the initial seed;

[0012] The initial seeds are biclustered and the rows and columns of the biclusters are expanded to obtain the noise-resistant biclustering results of the gene expression matrix.

[0013] In a second aspect, the present invention provides a noise-resistant biclustering system for gene expression data, comprising:

[0014] a processing module configured to: perform normalization processing on the gene expression matrix to obtain an index matrix sorted in ascending order of rows, and calculate the coefficient of variation of each row in the index matrix;

[0015] The calculation module is configured to: based on the longest common subsequence algorithm, if the elements at two index positions in the index matrix are within the noise disturbance range and the error range determined by the corresponding coefficient of variation, allow the elements to be exchanged, find the longest approximate common subsequence, and consider the positive and negative correlations between genes to determine the initial seed;

[0016] The biclustering and expansion module is configured to: perform biclustering on the initial seeds, perform row and column expansion on the biclusters, and obtain noise-resistant biclustering results of the gene expression matrix.

[0017] In a third aspect, the present invention provides an electronic device comprising a memory and a processor, and computer instructions stored in the memory and executed on the processor, wherein the computer instructions, when executed by the processor, perform the method described in the first aspect.

[0018] In a fourth aspect, the present invention provides a computer-readable storage medium for storing computer instructions, wherein when the computer instructions are executed by a processor, the method described in the first aspect is performed.

[0019] One or more of the above technical solutions have the following beneficial effects:

[0020] In the present invention, the longest common subsequence search method is used to obtain the seed list, and the element exchange strategy is allowed to find the longest approximate common subsequence, which can effectively reduce the interference of noise in the dataset on sequence matching, so that the algorithm can still accurately identify the similarity in gene expression patterns in a complex noisy environment, thereby improving the accuracy of the biclustering results.

[0021] In the present invention, when performing biclustering, not only the positive regulatory relationships between genes are taken into account, but also the negative regulatory relationships between genes, so that each bicluster contains positively correlated and negatively correlated genes; and the scheme of the present invention can perform efficient biclustering inference on high-throughput gene expression data with a scale of more than 10,000 rows within a reasonable running time; at the same time, the scheme has strong robustness to model parameters, reducing the result fluctuations caused by differences in parameter settings, and ensuring the reliability and consistency of the analysis results.

[0022] Advantages of additional aspects of the present invention will be given in part in the following description and in part will be obvious from the following description, or will be learned through practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] The accompanying drawings, which constitute a part of the present invention, are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute improper limitations on the present invention.

[0024] Figure 1 This is a block diagram of the noise-resistant biclustering method for gene expression data in Example 1 of the present invention;

[0025] Figure 2 This is a performance comparison of different algorithms in Example 1 of the present invention on a simulated data set;

[0026] Figure 3 The performance comparison of different algorithms in Example 1 of the present invention on a batch gene expression dataset is shown below:

[0027] Figure 4 This is a performance comparison of different algorithms in Example 1 of the present invention on a single-cell gene expression dataset;

[0028] Figure 5 Schematic diagram of the running time and memory status of the first embodiment of the present invention on simulated data. DETAILED DESCRIPTION

[0029] It should be noted that the following detailed descriptions are exemplary and intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present invention belongs.

[0030] It should be noted that the terms used herein are for describing particular embodiments only and are not intended to limit the exemplary embodiments according to the present invention.

[0031] In the absence of conflict, the embodiments of the present invention and the features thereof may be combined with each other.

[0032] Example 1

[0033] This embodiment discloses a noise-resistant biclustering method for gene expression data, comprising:

[0034] Step 1: Normalize the gene expression matrix to obtain an index matrix sorted in ascending order, and calculate the coefficient of variation of each row in the index matrix.

[0035] This embodiment can be used for gene expression data of various scales, such as microarray data, RNA-seq data, scRNA-seq data, etc. This embodiment uses the expression matrix M from a microarray, RNA sequencing experiment, or single-cell RNA sequencing experiment as an example of input data.

[0036] When the number of rows in a microarray data matrix exceeds 20,000, it is filtered using the coefficient of variation (CV) score, and the top 20,000 rows with the highest CV scores are retained. Input matrices for RNA sequencing data are normalized using the FPKM method. When the input matrix comes from a single-cell RNA sequencing experiment, highly variable genes are selected, and missing values ​​in the single-cell RNA sequencing data are corrected using the MAGIC method.

[0037] The processed gene expression matrix M is then normalized using the minimum-maximum normalization method:

[0038]

[0039] Among them, M′ i. is the normalized expression vector of the gene, M i. Represents the i-th row vector.

[0040] After that, we get the index matrix A that sorts the input matrix in ascending order of value. That is, the index matrix A is sorted in ascending order of value within each row. Each element of the index matrix is ​​stored as a structure (idx, val, val O ) in the form of (idx,val,val O ) represent the column index of the element in the matrix M, the current value and initial value of the element, val and val respectively. O Both are initialized to M i The value corresponding to the idx column in ′. For example:

[0041] M′ i· =(0.8,0,0.25,0.2,1,0.6)

[0042] A i· =((2,0,0),(4,0.2,0.2),(3,0.25,0.25),(6,0.6,0.6),(1,0.8,0.8),(5,1,1))

[0043] M′ i· is the normalized expression vector of the gene, A i· is the corresponding index vector.

[0044] Calculate the CV or coefficient of variation of each row of the index matrix A:

[0045]

[0046] Among them, μ i and σ i They represent the mean and variance of the values ​​in the i-th row respectively.

[0047] In the subsequent steps, the index matrix A will be clustered.

[0048] Step 2: Based on the longest common subsequence algorithm, in the index matrix, if the elements at two index positions are within the noise perturbation range and the error range determined by the corresponding coefficient of variation, elements are allowed to be exchanged to find the longest approximate common subsequence. Taking into account the positive and negative correlations between genes, the initial seeds are determined, the initial seeds are biclustered, and the rows and columns of the biclusters are expanded to obtain the noise-resistant biclustering results of the gene expression matrix.

[0049] After data preprocessing, the obtained index matrix A is clustered according to the following process:

[0050] Step 21: Find seeds.

[0051] The index matrix A is divided into k groups. In this embodiment, k is 4. The grouping principle is the same as the UniBic algorithm. In each group, the longest approximate common subsequence is found for each pair of rows. If the length of the longest approximate common subsequence exceeds a certain threshold and the repetition degree meets the conditions, the corresponding two rows are used as seeds and stored in the priority queue Q in the form of (x1, x2, len). seed In the example, the priority is len; x1, x2 are the corresponding row indices, and len is the length of the longest approximate common subsequence to be sought.

[0052] Specifically, this embodiment proposes the UniBic algorithm, which draws on the idea of ​​using the longest common subsequence search method to obtain a seed list in the existing UniBic algorithm and makes significant improvements and optimizations on this basis. The UniBic algorithm designs a longest approximate common subsequence algorithm. During the solution process of the classic longest common subsequence algorithm, the two elements being compared are allowed to be exchanged with each other under certain conditions, thereby finding the longest approximate common subsequence within the allowable error range.

[0053] Taking the two rows A1 and A2 of the index matrix A as an example, the longest approximate common subsequence of rows A1 and A2 is calculated according to the column index idx, and the numerical value is used to assist in determining whether to exchange elements.

[0054] The dynamic programming state transition equation for solving LACS (Longest Approximate Common Subsequence) is as follows:

[0055]

[0056] Among them, C[i,j] represents A 1[0..i] and A 2[0..j] The longest approximate common subsequence length, A 1i .idx represents the list corresponding to the i-th element of row A1, A 2j .idx represents the list corresponding to the j-th element in row A2.

[0057] When comparing the i-th element of A1 and the j-th element of A2, for A 1i .idx and A 2j The comparison of .idx has the following four results:

[0058] Case 1: A 1i .idx and A 2j .idx is the same, no need to exchange. Otherwise, calculate A 1i .idx and A 2j .idx are the penalties for swapping on A1 and A2 respectively.

[0059] Use id 11 ,id 12 Respectively represent A 1i .idx and A 2j .idx is the subscript position in A1, i.e. id 11 =i,id 12 satisfy id 21 ,id 22 Respectively represent A 1i .idx and A 2j .idx is the position in A2.

[0060] A 1i and A 2j In A1, the exchange is Exchange separately and and in, refer to Corresponding val O value, refer to Corresponding val O Value, A 1i and A 2j In A1, commutativity must satisfy the following conditions:

[0061] 1)id xy There exists, x, y∈{1,2} and id 12 >id 11 ,id 21 >id 22 ,

[0062] 2) exist and within the noise disturbance range, and Also there and within the noise disturbance range.

[0063] For any given value v aand v b ,

[0064] v a In v b Within the noise disturbance range Among them, α is provided by the user, the default is 0.2, and CV refers to v a and v b The CV score of the current row, and fabs is the fabs function.

[0065] exist and Within the noise disturbance range, that is and CV1 refers to and The CV score of the current row.

[0066] Also there and Within the noise disturbance range, that is

[0067]

[0068] and CV1 refers to and The CV score of the current row.

[0069] If both conditions 1) and 2) above are satisfied, then swap A in A1. 1i and A 2j The penalty points are:

[0070]

[0071] Among them, CV1 refers to and The CV score of the current row.

[0072] If it is not exchangeable, the penalty is infinite, and A is exchanged in A2. 1i and A 2j The penalty score 2 is calculated similarly.

[0073]

[0074] Among them, CV2 refers to and The CV score of the current row.

[0075] Case 2: score1 and score2 are both infinite, then A 1i and A 2j is not commutative in both A1 and A2.

[0076] Case 3: If score2 ≥ score1, then swap A in A1 1i and A 2j .

[0077] Otherwise, there is Case 4: If score2 < score1, swap A in A2 1i and A 2j .

[0078] The longest approximate common subsequence between rows can be solved through the dynamic programming state transition equation. The longest approximate common subsequence calculated for A1 and A2 represents the positive correlation relationship between gene 1 and gene 2. Considering that there are two cases of positive and negative co-expression between genes, take the negative value of each numerical value in A2 and re-sort it in ascending order of numerical values, and the obtained sequence is denoted as -A′2. At the same time, calculate the LACS between -A′2 and row A1, which represents the negative correlation relationship between gene 1 and gene 2. The LACS calculated between A1 and -A′2 and row A2 respectively, and take the longest approximate common subsequence with the longest length among them to update the LACS of the two sequences, that is:

[0079] LACS(A1,A2)←length_max{LACS(A1,A2),LACS(A1,-A′2)}

[0080] Select two rows with the longest approximate common subsequence as the initial seeds. During the clustering process, the longest approximate common subsequence algorithm can be used to expand the rows of the submatrix, so that the possible errors caused by data being perturbed by noise can be fully considered from the perspective of the matrix rows.

[0081] In addition, define the repetition rate of a sequence as the frequency of the most frequently occurring numerical value divided by the total length of the sequence. When calculating the longest approximate common subsequence of two sequences for seed selection, it is necessary to calculate the repetition rate of the obtained longest approximate common subsequence. If the repetition rate exceeds a certain threshold, the corresponding two sequences cannot be used as initial seeds.

[0082] Step 22: Strict clustering.

[0083] Select seeds for clustering according to the storage order of the priority queue. Assume S = Q seed .top, that is, the initial seed with the longest approximate common subsequence.

[0084] Initialize the bicluster B: (x1, x2, row in , row de , column), B.x1 = S.x1, B.x2 = S.x2, S.x1 represents the row index of one of the rows in the initial seed, represents the corresponding row in A, and for each row A of the index matrix Ar , to consider the case of negative correlation, take the negative of each value of A r , and denote the sequence obtained by re - sorting them in ascending order of values as -A′ r , and at the same time take A r , -A′ r and the row Find the longest approximate common subsequence, and store the solution result with the larger value of the obtained longest approximate common subsequence as the r - th row and the S.x1 - th row in Ans[r]. Its structure is (row, len, sig, P), where row represents the row index to be found, len represents the length of the longest approximate common subsequence to be found, sig is a 01 variable, taking 1 means the result comes from the original sequence indicating positive correlation between the two rows, 0 means the result comes from the sequence after taking the negative and reversing, meaning negative correlation between the two rows to be found, and P records the elements of the longest approximate common subsequence to be found.

[0085] Sort Ans in descending order according to len. If Ans[0].sig = 1, then add Ans[0].row to B.row in , otherwise add it to B.row de . Update B.column with Ans[0].P, delete the elements in Ans[i].P that do not belong to Ans[0].P, update Ans[i].len according to the size of Ans[i].P after deletion, finally assign 0 to Ans[0].len, re - sort Ans in descending order according to len, and then iteratively update Ans and B until the termination condition is reached.

[0086] The termination condition of this embodiment is set as:

[0087] 1) min(size(B.row<00,00071>)) + B.row in ), size(B.column)) is less than before iteration, or

[0088] 2) size(B.column) < S.len * δ, where δ is provided by the user, with a default value of 0.15.

[0089] Step 23: Expand the rows and columns of the bicluster obtained in Step 22 within the error range.

[0090] First, expand the columns of the bicluster B. Take each column in the index matrix A that does not belong to column as a candidate column, use the binary search algorithm for column expansion to calculate the insertion position of the candidate column. If the candidate column with the highest obtained frequency satisfies f max ≥ (size(B.row de +B.row in)*η), the candidate column is considered eligible for insertion. Next, row expansion is performed on the expanded biclusters. The longest approximate common subsequence is calculated for each row in the index matrix A that does not belong to B and the row B.x1 of the bicluster. If the length of the obtained longest approximate common subsequence is greater than size(B.column)*η, the row is considered eligible for insertion, where η is provided by the user and defaults to 0.85.

[0091] Specifically, when performing column expansion on the desired bicluster, the possible value range of the candidate column is determined, and a binary search is performed in the row direction for the upper and lower bounds of the range to determine the corresponding insertion interval, thereby effectively considering the impact of noise disturbance at the numerical level.

[0092] Then, statistics are collected on the candidate column insertion positions of all rows in the bicluster, and the column label with the highest frequency is used as the insertion position of the candidate column. in ,row de , column) to store the biclustering submatrix, where x1 and x2 store the initial seed row labels for generating the bicluster, and row in ,row de The expanded rows that are positively correlated with x1 and the rows that are negatively correlated with x1 are stored separately, and column stores all the columns of the biclustering submatrix.

[0093] Take any column to be inserted For row de ,row in For each row in M′[r][c].val, r∈row de ,row in The value range of the element to be inserted into column c is the value range of the element in the corresponding row of the bicluster, and the insertion interval is binary searched in the rth row of the bicluster. For the rth row, extract M′[r][l].val, l∈column, that is, the elements of the rth row in the bicluster, to form a new vector C, and calculate the coefficient of variation score of C, recorded as cv r , using cv r =cv r *γ to calculate the value range of M′[r][c].val [val l ,val r ], where γ is provided by the user and defaults to 0.2:

[0094]

[0095] Then the insertion interval [low, up] of M′[r][c] in the rth row of the biclustering submatrix, that is, M′[r][l], l∈column) can be respectively calculated by val l,val r Call the lower_bound and upper_bound functions of C++ to get the index position of the left and right elements of the interval. Count the [low,up] intervals for each row and calculate the index position col of the column with the highest frequency and the corresponding frequency f. max Write it down.

[0096] Step 24: Output the number of biclusters required by the user.

[0097] For each bicluster, its area is calculated as size(B.column)*size(B.row de +B.rowin, sort all the obtained biclusters in descending order of area. The overlapping area of ​​two biclusters is the area of ​​the overlapping submatrix. If biclusters b1 and b2 overlap and satisfy area overlap ≥f*min(area b1 ,area b2 ), the bicluster with the smallest area is deleted. The parameter f is provided by the user and is 0.8 by default.

[0098] This embodiment proposes the NoiBic algorithm for biclustering. This algorithm draws on the existing UniBic algorithm's approach of using the longest common subsequence search method to obtain a seed list, and significantly improves and optimizes this approach, resulting in its own unique algorithmic features: First, the NoiBic algorithm proposes a strategy that allows element swapping to find the longest approximate common subsequence. This strategy effectively reduces the interference of noise in the dataset on sequence matching, enabling the algorithm to accurately identify similarities in gene expression patterns even in complex noisy environments, thereby improving the accuracy of biclustering results. Second, during the column expansion phase, the NoiBic algorithm uses a method that calculates insertable intervals. This method fully accounts for the characteristics of noisy data. Compared to traditional expansion methods, it can more rationally handle the impact of noise, avoid erroneous expansion caused by noise, and thus more accurately construct bicluster structures. Third, when performing biclustering, NoiBic considers not only positive regulatory relationships between genes, but also negative regulatory relationships. As a result, each bicluster obtained by the algorithm contains both positively and negatively correlated genes. Fourth, the NoiBic algorithm exhibits excellent computational performance and parameter stability. It can perform efficient biclustering inference on high-throughput gene expression data exceeding 10,000 rows within a reasonable runtime. Furthermore, the algorithm is highly robust to model parameters, reducing fluctuations in results caused by differences in parameter settings and ensuring the reliability and consistency of analysis results.

[0099] This example evaluates the proposed solution from the following two aspects: (1) Comparison of this solution with other algorithms on simulated data, and (2) Comparison of this solution with other algorithms on real gene expression data.

[0100] (1) Generation and evaluation criteria of simulation data.

[0101] Because there are no real-world gold-standard datasets for benchmarking biclustering algorithms, algorithms are typically validated on simulated datasets. Five 10,000×100 matrices are generated, each with randomly generated elements from a standard normal distribution. Submatrices are generated by randomly selecting a predetermined number of rows and columns, and their elements are rearranged to preserve the trend. Overlap between submatrices from the same initial matrix is ​​allowed. To test the model's ability to eliminate the effects of noise, varying levels of noise are added to matrices drawn from a standard normal distribution N(0,1). The noise level represents the proportion of randomly selected rows and columns to which noise is added.

[0102] Using Amela The matching score proposed, M1, M2 are two groups of biclusters, b1, b2 are two biclusters from M1 and M2 respectively. Then, the matching score between the two groups of biclusters is as follows:

[0103]

[0104] J(b1,b2) in the above formula represents the Jaccard similarity coefficient between b1 and b2, which is defined as follows:

[0105]

[0106] where T and P are the set of true biclusters embedded in the simulated data matrix and the set of predicted biclusters obtained by the biclustering method, respectively. The matching scores S(T,P) and S(P,T) are the recovery score and correlation score, respectively, which are the evaluation criteria for the simulated data.

[0107] (2) Evaluation criteria for true expression data.

[0108] For microarray and RNA sequencing data, the gene modules found by the algorithm were validated. The gene modules in each bicluster were evaluated using the F1 score mentioned in ARBic. A bicluster is considered enriched if and only if the gene module is enriched in at least one pathway in the KEGG database (P value ≤ 5 after the Benjamini-Hochberg test provided by the R toolkit clusterProfiler). If an algorithm finds k biclusters, of which l biclusters are enriched, the F score is calculated as follows:

[0109]

[0110] For a bicluster B = (I, J), where I represents the gene set and J represents the sample set, the gene enrichment score g is calculated as follows:

[0111]

[0112] Among them, K i ,i=1,2,…,p is the pathway enriched in bicluster B.

[0113] Then, the G score is calculated as follows:

[0114]

[0115] The F1 score used to evaluate the algorithm's ability to detect gene modules is defined as follows:

[0116]

[0117] For scRNA-Seq data, the cell type was used as the true label, and the cell clustering ability of the algorithm was evaluated by calculating the adjusted Rand index (ARI) and normalized mutual information (NMI) scores.

[0118] NoiBic is robust to noise on simulated datasets:

[0119] Because the biclustering problem lacks a gold-standard real-world dataset to validate the results, we benchmarked our algorithms on generated simulated data. To verify their robustness to data noise, we ran them on simulated data with varying noise levels. We generated four 10,000-by-100 data matrices with noise levels of 0.1, 0.2, 0.3, and 0.4, embedding three 1,000-by-20 trend-preserving submatrices within each matrix. We then averaged the recovery and correlation scores calculated for the four matrices at varying noise levels to measure the accuracy of the model in identifying biclusters. The detailed calculation methods for the recovery and correlation scores are described in the Methods section.

[0120] like Figure 2 As shown in Figure A, NoiBic outperforms other algorithms at different noise levels. As the noise level increases, the performance of other algorithms shows a clear downward trend. For example, the scores of RecBic and ARBic drop from close to 0.9 to 0.5. However, as the noise level increases, our method is able to maintain a level above 0.75, which proves that our method can effectively eliminate the noise effect in the data. Figure 2As shown in B, the index position of the embedding submatrix is ​​disrupted and the benchmark is re-performed. NoiBic still outperforms other algorithms.

[0121] like Figure 4 As shown in the figure, on 10,000 rows of simulated data, NoiBic maintains an accuracy above 0.75 while keeping the running time and memory usage at an average level of 19 minutes and 4.34GB. Each point in the scatter plot represents a simulated data set.

[0122] NoiBic effectively identifies functional gene modules in bulk gene expression data: To evaluate its ability to identify functional gene modules, we compared NoiBic with other methods on six bulk transcriptome datasets. These datasets include two from Escherichia coli, three from humans, and one from mice. The HerBC dataset was obtained from breast cancer patients using RNA sequencing, while the other five datasets were obtained using microarray technology. Details of these datasets can be found in Table 1.

[0123] The F1 score mentioned in ARBic was used to evaluate the functional gene modules in the biclusters found by the algorithm. The calculation method of the F1 score has been described in the evaluation scheme section. In order to minimize the bias caused by the selection of tool parameters, a grid search was performed on the parameters of eight algorithms including NoiBic, QUBIC, QUBIC2, UniBic, FABIA, ISA, OPSM and CPB. These eight algorithms were run on 40 to 60 parameter combinations on six data sets. The detailed information on the parameter combination selection is provided in Table 3. Figure 3 In each boxplot, each possible parameter combination corresponds to a data point. Since ARBic and RecBic take a long time to run on the scale of real datasets, it is not possible to perform grid search on their parameters, so they are not included in the comparison. Figure 3 As shown in the figure, NoiBic outperformed other software on all six data sets. Notably, on the E. coli data, NoiBic achieved an F1 score of 1 under various parameter combinations, significantly higher than the next closest competitor, UniBic.

[0124] NoiBic is able to effectively classify cell clusters on scRNA-seq data: NoiBic is compared with other algorithms on eleven datasets of different sizes, as detailed in Table 2. These datasets include seven datasets from humans and four datasets from mice. Since other algorithms are not applicable to single-cell data, this example is only compared with QUBIC2. In order to make a fair comparison, this example performed a grid search on the parameters of QUBIC2 and NoiBic. Fifteen parameter combinations were used for each method, as detailed in Table 3. Taking into account that there are now many single-cell data with cell type labels, this example calculates the adjusted Rand index (ARI) and normalized mutual information (NMI) by comparing the cluster labels of the cells with the true cell type labels of the cells. This requires that there is no overlap between the two clusters, so this example sets the f parameters of QUBIC and NoiBic to zero.

[0125] like Figure 4 As shown in the figure, NoiBic performed better than QUBIC2 in both ARI and NMI scores in both human and mouse data. In three human datasets and four mouse datasets, the Wilcoxon rank sum test showed that NoiBic's ARI scores were significantly higher than QUBIC2 at a significance level below 0.05. Furthermore, in four human datasets and four mouse datasets, the Wilcoxon rank sum test showed that NoiBic's NMI scores were significantly higher than QUBIC2 at a significance level below 0.05.

[0126] Table 1: Batch gene expression dataset details

[0127]

[0128] Table 2: Single-cell gene expression dataset details

[0129]

[0130]

[0131] Table 3: Parameter selection of each algorithm

[0132]

[0133]

[0134] Example 2

[0135] The purpose of this embodiment is to provide a noise-resistant biclustering system for gene expression data, including:

[0136] a processing module configured to: perform normalization processing on the gene expression matrix to obtain an index matrix sorted in ascending order of rows, and calculate the coefficient of variation of each row in the index matrix;

[0137] A calculation module is configured to: based on the longest common subsequence algorithm, in the index matrix, if the elements at two index positions are within the noise perturbation range and the error range determined by the corresponding coefficient of variation, allow the elements to be exchanged, find the longest approximate common subsequence, and consider the positive and negative correlations between genes to determine the initial seed;

[0138] The biclustering and expansion module is configured to: perform biclustering on the initial seeds, perform row and column expansion on the biclusters, and obtain noise-resistant biclustering results of the gene expression matrix.

[0139] In further embodiments, there is also provided:

[0140] An electronic device includes a memory and a processor, and computer instructions stored in the memory and executed by the processor. When the computer instructions are executed by the processor, the method described in Example 1 is performed. For the sake of brevity, no further details are given here.

[0141] It should be understood that in this embodiment, the processor may be a central processing unit (CPU), or may be other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), off-the-shelf field-programmable gate arrays (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any conventional processor, etc.

[0142] The memory may include a read-only memory and a random access memory, and provides instructions and data to the processor. A portion of the memory may also include a non-volatile random access memory. For example, the memory may also store information about the device type.

[0143] A computer-readable storage medium is used to store computer instructions, and when the computer instructions are executed by a processor, the method described in embodiment 1 is performed.

[0144] The method in Example 1 can be directly implemented as being executed by a hardware processor, or by a combination of hardware and software modules within the processor. The software module can be located in a storage medium well-established in the art, such as random access memory, flash memory, read-only memory, programmable read-only memory, electrically erasable programmable memory, or registers. The storage medium is located in the memory, and the processor reads the information in the memory and, in conjunction with its hardware, completes the steps of the above method. To avoid repetition, a detailed description is not given here.

[0145] A computer program product includes a computer program, and when the computer program is executed by a processor, the method described in embodiment 1 is implemented.

[0146] The present invention also provides at least one computer program product tangibly stored on a non-transitory computer-readable storage medium. The computer program product includes computer-executable instructions, such as instructions contained in program modules, which are executed in a device on a real or virtual processor of a target to perform the process / method described above. Generally, program modules include routines, programs, libraries, objects, classes, components, data structures, etc. that perform specific tasks or implement specific abstract data types. In various embodiments, the functionality of program modules can be combined or divided between program modules as needed. The machine-executable instructions for the program modules can be executed in local or distributed devices. In distributed devices, program modules can be located in local and remote storage media.

[0147] The computer program code for implementing the method of the present invention can be written in one or more programming languages. These computer program codes can be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing device so that when the program code is executed by the computer or other programmable data processing device, the functions / operations specified in the flow chart and / or block diagram are implemented. The program code can be executed entirely on a computer, partially on a computer, as an independent software package, partially on a computer and partially on a remote computer, or entirely on a remote computer or server.

[0148] In the context of the present invention, computer program code or related data can be carried by any appropriate carrier to enable a device, apparatus, or processor to perform the various processes and operations described above. Examples of carriers include signals, computer-readable media, and the like. Examples of signals include electrical, optical, radio, acoustic, or other forms of propagated signals, such as carrier waves, infrared signals, and the like.

[0149] Those skilled in the art will appreciate that the units and algorithm steps of the various examples described in conjunction with this embodiment can be implemented in electronic hardware or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professionals and technicians can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0150] Although the above describes the specific embodiments of the present invention in conjunction with the accompanying drawings, it is not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art on the basis of the technical solution of the present invention without any creative work are still within the scope of protection of the present invention.

Claims

1. A noise-resistant biclustering method for gene expression data, characterized in that: include: Normalizing the gene expression matrix to obtain an index matrix sorted in ascending order, and calculating the coefficient of variation of each row in the index matrix; Based on the longest common subsequence algorithm, in the index matrix, if the elements at two index positions are within the noise perturbation range and the error range determined by the corresponding coefficient of variation, elements are allowed to be exchanged to find the longest approximate common subsequence, and the positive and negative correlations between genes are considered to determine the initial seed; The initial seeds are biclustered and the rows and columns of the biclusters are expanded to obtain the noise-resistant biclustering results of the gene expression matrix.

2. The noise-resistant biclustering method for gene expression data according to claim 1, wherein: The gene expression matrix is ​​normalized using a minimum-maximum normalization method to obtain an index matrix sorted in ascending order of rows; wherein each element in the index matrix includes the column index of the element in the gene expression matrix, the current value of the element, and the initial value of the element.

3. The noise-resistant biclustering method for gene expression data according to claim 1, wherein: The element A in the index matrix 1i and element A 2j The following conditions must be met to allow the swap in row A1: Condition 1: id xy Exists, x,y∈{1,2} and id 12 >id 11 ,id 21 >id 22 ; Among them, id 11 ,id 12 Respectively represent A 1i .idx and A 2j .idx is the subscript position in A1, id 21 ,id 22 Respectively represent A 1i .idx and A 2j .idx is the position of element A in A2 1i is the i-th element in row A1 of the index matrix, element A 2j is the jth element in row A2 of the index matrix; A 11 .idx and A 2j .idx represents the column index corresponding to the i-th element in row A1 and the column index corresponding to the j-th element in row A2; Condition two: exist and within the noise disturbance range, and Also there and within the error range; If conditions one and two are satisfied, calculate the penalty score1 for swapping A in A1 1i and A 2j and calculate the penalty score2 for swapping A in A2 1i and A 2j If score2 ≥ score1, then swap A in A1 1i and A 2j ; if score2 < score1, then swap A in A2 1i and A 2j ; A 1i and A 2j Swapping in A1 means separately swap with and where (idx, val, val O ) respectively represent the column index, current value, and initial value of the element in the gene expression matrix.

4. The noise-resistant biclustering method for gene expression data according to claim 3, wherein: exist and The error range is as follows: and Also there and The noise disturbance range is as follows: and CV1 refers to and The CV score of the current row, fabs is the fabs function, and α is the coefficient.

5. The noise-resistant biclustering method for gene expression data according to claim 1, wherein: Combined with the longest approximate common subsequence and considering the positive and negative correlations between genes, the initial seed is determined as follows: Calculate positive and negative correlations between genes; The longest approximate common subsequence among the positive and negative correlations between genes is taken as the initial seed.

6. The noise-resistant biclustering method for gene expression data according to claim 1, wherein: Perform row and column expansion on the double clustering, specifically: Perform column expansion on the biclusters, taking each column in the index matrix that does not belong to the bicluster columns as a candidate column, and using the binary search algorithm of column expansion to calculate the insertion position of the candidate column. If the candidate column with the highest frequency meets the conditions, then the candidate column can be inserted; After column expansion, row expansion is performed on the biclusters. The longest approximate common subsequence is calculated for each row in the index matrix that does not belong to the bicluster and one of the rows of the bicluster. If the length of the obtained longest approximate common subsequence is greater than the length of the bicluster column, the row of the bicluster can be inserted.

7. The noise-resistant biclustering method for gene expression data according to claim 1, wherein: After column expansion of the row-expanded biclusters using a binary search algorithm, the method includes: performing overlapping filtering on the two biclusters after column expansion, and deleting the bicluster with the smallest area among the two biclusters if the overlapping area is greater than a set multiple of the minimum area of ​​the two biclusters, thereby obtaining a noise-resistant biclustering result for the expression data.

8. A noise-resistant biclustering system for gene expression data, characterized in that: include: a processing module configured to: perform normalization processing on the gene expression matrix to obtain an index matrix sorted in ascending order of rows, and calculate the coefficient of variation of each row in the index matrix; The calculation module is configured to: based on the longest common subsequence algorithm, if the elements at two index positions in the index matrix are within the noise disturbance range and the error range determined by the corresponding coefficient of variation, allow the elements to be exchanged, find the longest approximate common subsequence, and consider the positive and negative correlations between genes to determine the initial seed; The biclustering and expansion module is configured to: perform biclustering on the initial seeds, perform row and column expansion on the biclusters, and obtain noise-resistant biclustering results of the gene expression matrix.

9. An electronic device, characterized in that: The method comprises a memory and a processor, and computer instructions stored in the memory and executed on the processor, wherein when the computer instructions are executed by the processor, the method according to any one of claims 1 to 7 is completed.

10. A computer-readable storage medium, characterized in that Used to store computer instructions, which, when executed by a processor, complete the method according to any one of claims 1 to 7.

Images