Method and system for measuring transcriptome similarity

Abstract

The invention relates to the technical field of bioinformatics, in particular to a transcriptome similarity measurement method and system. The invention discloses a method for calculating an elastic network weighted transcriptome feature distance (EN-wTSD). The invention further discloses a transcriptome similarity analysis method based on the calculation method. A system for implementing the analysis method is also disclosed.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics technology, specifically to a method and system for measuring transcriptome similarity. Background Technology

[0002] In the fields of precision medicine and new drug development, accurately assessing the similarity between models (such as experimental animal models, organoids, and organ-on-a-chip) and real human physiological and pathological states at the physiological / pathological transcriptome level is crucial for improving the accuracy of disease research and drug prediction. Traditional animal models have significant limitations in predicting human drug responses, primarily due to the vast differences in gene expression patterns and regulatory mechanisms among species. Despite the continuous development of novel microphysiological systems (MPS) such as organ-on-a-chips, organoids, and assemblies, there is currently a lack of standardized, cross-model applicable, and highly sensitive computational methods to quantify the distance or similarity between models and human tissues at the transcriptome level, whether for traditional animal models, in vitro models, or emerging MPS.

[0003] Existing cross-species or cross-model comparison methods typically rely on simple differential gene overlap or simplistic measures such as Pearson correlation coefficients and Euclidean distances. These methods struggle to capture the complex regulatory network differences between species or models. While Manatakis et al. (2020) proposed Transcriptomic Signature Distance (TSD), offering a novel perspective, their original algorithm relies on the Glasso (Graphical Least Absolute Shrinkage and Selection Operator) algorithm for weight calculation. This results in excessively high baseline noise (background noise) for similar samples (e.g., human-to-human samples) and insufficient resolution for subtle transcriptomic drift, leading to low algorithm efficiency. Furthermore, in large-scale cross-species or cross-model data integration analysis, systematic solutions for removing batch effects, extracting conserved key feature genes, and constructing robust evaluation systems remain lacking.

[0004] References

[0005] 1. Manatakis DV, VanDevender A, Manolakos E S. An information-theoretic approach for measuring the distance of organ tissue samples using their transcriptomic signatures [J]. Bioinformatics, 2020, 36(21): 5194-204.

[0006] 2. Bailey P, Emad A. wCorr: Weighted correlations [J]. R packageversion, 2023, 1(8).

[0007] 3. Love M, Anders S, Huber W. Differential analysis of count data–theDESeq2 package [J]. Genome Biol, 2014, 15(550): 10-1186. Summary of the Invention

[0008] To address the technical problems existing in the prior art, this invention innovatively introduces the ElasticNet algorithm to improve the weighted Transcriptomic Signature Distance (TSD), constructing the Elastic Net weighted Transcriptomic Signature Distance (EN-wTSD), which significantly reduces baseline noise and improves the ability to identify subtle feature differences.

[0009] This invention also discloses a system for assessing cross-species transcriptome similarity. It can be widely applied to cross-species transcriptome similarity analysis between various biological models such as organoids, organ-on-a-chip, and experimental animal models and humans. It can be used for quantitative assessment of human physiological / pathological similarity and to assist in drug screening decisions, providing a reliable basis for the evaluation of organ-on-a-chip and organoid models. The system includes multi-source data standardization, feature gene optimization, and automated analysis processes, solving the problems of high baseline noise, low sensitivity to weak signals, and poor computational efficiency in traditional measurement methods for cross-species comparisons. The system can be implemented as a computer program (e.g., an R package).

[0010] This invention provides the following technical solutions:

[0011] ①. A method for analyzing the transcriptome similarity between two transcriptome data samples from a first biological sample and a second biological sample, respectively, wherein the first biological sample and the second biological sample each include at least two transcriptome data samples, the method comprising: calculating a weighted transcriptome distance between the two transcriptome data samples from the first biological sample and the second biological sample, respectively, and determining the transcriptome similarity between the two transcriptome data samples;

[0012] The method for calculating the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample, respectively, includes the following steps:

[0013] (1) Obtain the expression matrix of characteristic genes in the transcriptome of each sample of the first biological sample. Where n1 is the number of transcriptome data samples of the first biological sample, p1 is the number of characteristic genes of the first biological sample, and the elements in matrix X1 are the expression levels of each characteristic gene in each transcriptome data sample of the first biological sample; using the same method, a matrix of expression levels of characteristic genes in each transcriptome data sample of the second biological sample is obtained. Where n2 is the number of transcriptome data samples of the second biological sample, p2 is the number of characteristic genes of the second biological sample, and the elements in the matrix X2 are the expression levels of each characteristic gene in each transcriptome data sample of the second biological sample.

[0014] (2) Take any one of the p1 genes in the matrix X1. The vector of expression levels as response variables The expression matrices of the remaining p1-1 genes were used as predictor variables. A resilient network regression model is constructed, where the nodes of the resilient network are each gene. The regression coefficient vector of the nodes is minimized through the following objective function. The objective function is:

[0015]

[0016] in For mixed parameters, As the penalty coefficient, the same calculation is performed on the remaining p1-1 genes except for gene k to obtain their regression coefficient vectors. A symmetric matrix A1 is constructed using all the regression coefficient vectors. The symmetric matrix A1 is a covariance matrix or a precision matrix. Using the same method, a symmetric matrix A2 is obtained from the matrix X2. The symmetric matrix A2 is a covariance matrix or a precision matrix.

[0017] (3) Calculate the square of the Mahalanobis distance of the characteristic genes of any transcriptome data sample i among the n1 transcriptome data samples of the first biological sample. The calculation formula is:

[0018]

[0019] in Transcriptome data samples The vector of expression levels of characteristic genes; the same calculation is performed on any transcriptome data sample j from the n2 transcriptome data samples of the second biological sample to obtain the square of its Mahalanobis distance; further, the specific weight of transcriptome data sample i is calculated. and the specific weights of the transcriptome data sample j The calculation formula is:

[0020] ;

[0021] The normalized specific weights of the transcriptome data sample i are further calculated. and the normalized specific weights of the transcriptome data sample j The calculation formula is:

[0022]

[0023] ;

[0024] (4) Calculate the weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j. i-j The weighted transcriptome distance EN-wTSD is calculated using the following formula:

[0025] EN-wTSD

[0026] in, For weighted Jensen-Shannon divergence The weighted correlation distance is calculated based on the weighted Pearson correlation coefficient, where γ is the coefficient.

[0027] in The calculation formula is:

[0028]

[0029] in and These are the probabilistic vector distributions of gene expression levels for transcriptome data sample i and transcriptome data sample j, respectively. The function is , For the sample in the gene The expression probability distribution on M is and Constructed hybrid distribution:

[0030]

[0031] in The calculation formula is:

[0032]

[0033] in The weighted Pearson correlation coefficient;

[0034] The method for determining the transcriptomic similarity between the two transcriptomic data samples is as follows:

[0035] The weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j i-j The closer to 0, the higher the transcriptome similarity between transcriptome data sample i and transcriptome data sample j.

[0036] In the implementation scheme, the mixing parameter α is in the range of 0-1, preferably 0.5.

[0037] In the implementation scheme, the range of γ is 0-1, preferably 0.5.

[0038] ②. A method for analyzing transcriptomic similarity between two transcriptomic data samples from a biological sample, the biological sample comprising at least three transcriptomic data samples, the method comprising: calculating a weighted transcriptomic distance between two transcriptomic data samples from the biological sample, and determining transcriptomic similarity between the two transcriptomic data samples;

[0039] The method for calculating the weighted transcriptomic distance between two transcriptomic data samples from the biological sample includes the following steps:

[0040] (1) Obtain the expression matrix of characteristic genes in the transcriptome of each of the biological samples. , where n is the number of transcriptome data samples of the biological sample, p is the number of characteristic genes of the biological sample, and the elements in the matrix X are the expression levels of each characteristic gene in each transcriptome data sample of the biological sample.

[0041] (2) Take any one of the p genes in the matrix X. The vector of expression levels as response variables The expression matrices of the remaining p-1 genes are used as predictor variables. A resilient network regression model is constructed, where the nodes of the resilient network are each gene. The regression coefficient vector of the nodes is minimized through the following objective function. The objective function is:

[0042]

[0043] in For mixed parameters, The same calculation is performed on the remaining p-1 genes except for gene k, using the penalty coefficient, to obtain their regression coefficient vectors respectively. A symmetric matrix A is constructed using all the regression coefficient vectors, where the symmetric matrix A is the covariance matrix or the precision matrix.

[0044] (3) Calculate the square of the Mahalanobis distance of the characteristic genes of any transcriptome data sample i among the n transcriptome data samples of the biological sample. The calculation formula is:

[0045]

[0046] in It is a vector of the expression levels of characteristic genes in transcriptome data sample i; the same calculation is performed on any transcriptome data sample j among the n transcriptome data samples of the biological sample to obtain the square of its Mahalanobis distance;

[0047] Further calculate the specific weights of the transcriptome data sample i. and the specific weights of the transcriptome data sample j The calculation formula is:

[0048] ;

[0049] The normalized specific weights of the transcriptome data sample i are further calculated. and the normalized specific weights of the transcriptome data sample j The calculation formula is:

[0050]

[0051] ;

[0052] (4) Calculate the weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j. i-j The weighted transcriptome distance EN-wTSD is calculated using the following formula:

[0053] EN-wTSD

[0054] in, For weighted Jensen-Shannon divergence The weighted correlation distance is calculated based on the weighted Pearson correlation coefficient, where γ is the coefficient.

[0055] in The calculation formula is:

[0056]

[0057] in and These are the probabilistic vector distributions of gene expression levels for transcriptome data sample i and transcriptome data sample j, respectively. The function is , For the sample in the gene The expression probability distribution on M is and Constructed hybrid distribution:

[0058]

[0059] in The calculation formula is:

[0060]

[0061] in The weighted Pearson correlation coefficient;

[0062] The method for determining the transcriptomic similarity between the two transcriptomic data samples is as follows:

[0063] The weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j i-j The closer to 0, the higher the transcriptome similarity between transcriptome data sample i and transcriptome data sample j.

[0064] In the implementation scheme, the mixing parameter α is in the range of 0-1, preferably 0.5.

[0065] In the implementation scheme, the range of γ is 0-1, preferably 0.5.

[0066] ③. A method for analyzing transcriptomic similarity between a first biological sample and a second biological sample, wherein the first biological sample and the second biological sample each comprise at least two transcriptomic data samples, the method comprising:

[0067] (1) Based on the method described in ① for calculating the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample respectively, calculate the EN-wTSD between all pairwise transcriptome data samples, and further calculate their average EN-wTSD. 平均 Or the median value EN-wTSD 中位 The two transcriptome data samples mentioned therein are either transcriptome data samples of the first biological sample or transcriptome data samples of the second biological sample.

[0068] The EN-wTSD 平均 Or EN-wTSD 中位 The closer to 0, the higher the transcriptomic similarity between the first biological sample and the second biological sample;

[0069] (2) Based on the method described in ① for calculating the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample respectively, calculate the EN-wTSD between all pairwise transcriptome data samples of the first biological sample, and construct its set D. HH Calculate the quantile threshold T sim , ,in Represents the quantile function;

[0070] Based on the method described in ① for calculating the weighted transcriptome distance between two transcriptome data samples from the first and second biological samples respectively, the EN-wTSD between pairwise transcriptome data samples from all different biological samples is calculated, and the average EN-wTSD is further calculated. 平均 Or the median value EN-wTSD 中位 The pairwise transcriptome data samples of the different biological samples are either transcriptome data samples of the first biological sample or transcriptome data samples of the second biological sample.

[0071] If EN-wTSD 平均 ≤T sim or EN-wTSD 中位 ≤T sim If the transcriptome of the second biological sample is similar to that of the first biological sample, then the transcriptome of the second biological sample is determined to be similar; if EN-wTSD 平均 >T sim or EN-wTSD 中位 >T sim If so, it is determined that the transcriptome of the second biological sample is not similar to the transcriptome of the first biological sample; or

[0072] (3) Based on the method described in ① for calculating the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample respectively, obtain the EN-wTSD between all pairwise transcriptome data samples of the first biological sample, and construct its set D. HH Calculate the quantile threshold T sim , ,in Represents the quantile function;

[0073] Based on the method described in ① for calculating the weighted transcriptome distance between two transcriptome data samples from the first and second biological samples respectively, calculate the weighted transcriptome distance EN-wTSD between all pairwise transcriptome data samples from different biological samples, and construct its set D. HM The set is calculated where each pair of transcriptome data samples from different biological samples is either a transcriptome data sample from the first biological sample or a transcriptome data sample from the second biological sample. middle The number of pairwise transcriptome data samples Z was used to further calculate the transcriptome similarity ratio between the first biological sample and the second biological sample. The calculation formula is: ,in For set The total number of pairwise transcriptome data samples;

[0074] Similarity ratio The higher the value, the greater the transcriptome similarity between the first biological sample and the second biological sample.

[0075] ④. The method of any one of ①-③, wherein the transcriptome data samples are derived from organ-on-a-chip, organoids, assemblies, animal tissues or ex vivo human tissues.

[0076] ⑤. A system for analyzing transcriptomic similarity between two transcriptomic data samples from a first biological sample and a second biological sample, respectively, wherein the first biological sample and the second biological sample each comprise at least two transcriptomic data samples, the system comprising:

[0077] 1) Input terminal, used to receive transcriptome data from each transcriptome data sample, including gene name and gene expression level;

[0078] 2) Matrix construction module: Constructs matrices of the expression levels of characteristic genes in the transcriptomes of each transcriptome data sample from the first biological sample. Where n1 is the number of transcriptome data samples of the first biological sample, p1 is the number of characteristic genes of the first biological sample, and the elements in matrix X1 are the expression levels of each characteristic gene in each transcriptome data sample of the first biological sample; using the same method, a matrix of expression levels of characteristic genes in each transcriptome data sample of the second biological sample is obtained. Where n2 is the number of transcriptome data samples of the second biological sample, p2 is the number of characteristic genes of the second biological sample, and the elements in the matrix X2 are the expression levels of each characteristic gene in each transcriptome data sample of the second biological sample.

[0079] Take any one of the p1 genes in the matrix X1. The vector of expression levels as response variables The expression matrices of the remaining p1-1 genes were used as predictor variables. A resilient network regression model is constructed, where the nodes of the resilient network are each gene. The regression coefficient vector is minimized through the following objective function. The objective function is:

[0080]

[0081] in For mixed parameters, As the penalty coefficient, the same calculation is performed on the remaining p-1 genes except gene k to obtain their regression coefficient vectors. A symmetric matrix A1 is constructed using all the regression coefficient vectors. The symmetric matrix A1 is a covariance matrix or a precision matrix. Using the same method, a symmetric matrix A2 is obtained from the matrix X2. The symmetric matrix A2 is a covariance matrix or a precision matrix.

[0082] 3) Weight calculation module: Calculates the squared Mahalanobis distance of the characteristic genes of any transcriptome data sample i among the n1 transcriptome data samples of the first biological sample. The calculation formula is:

[0083]

[0084] in It is a sample The standardized vector of expression levels of characteristic genes; the same calculation is performed on any transcriptome data sample j in the n2 samples of the second biological sample to obtain the square of its Mahalanobis distance;

[0085] Further calculate the specific weights of the transcriptome data sample i. and the specific weights of the transcriptome data sample j The calculation formula is:

[0086] ;

[0087] The normalized specific weights of the transcriptome data sample i are further calculated. and the normalized specific weights of the transcriptome data sample j The calculation formula is:

[0088]

[0089]

[0090] 4) Transcriptome distance calculation module: Calculates the weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j. i-j The weighted transcriptome distance EN-wTSD is calculated using the following formula:

[0091] EN-wTSD

[0092] in, For weighted Jensen-Shannon divergence; The weighted correlation distance is calculated based on the weighted Pearson correlation coefficient, where γ is the coefficient.

[0093] in The calculation formula is:

[0094]

[0095] in and These are vectors representing the probabilistic gene expression levels of transcriptome data sample i and transcriptome data sample j, respectively. The function is , For the sample in the gene The expression probability distribution on M is and Constructed hybrid distribution:

[0096]

[0097] in The calculation formula is:

[0098]

[0099] The weighted Pearson correlation coefficient;

[0100] 5) Output terminal, output the EN-wTSD.

[0101] ⑥. A system for analyzing transcriptomic similarity between two transcriptomic data samples from a biological sample, said biological sample comprising at least three transcriptomic data samples, said system comprising:

[0102] 1) Input terminal, used to receive transcriptome data from each transcriptome data sample, including gene name and gene expression level;

[0103] 2) Matrix construction module: Constructs a matrix of the expression levels of characteristic genes in the transcriptome of each transcriptome data sample of the biological sample. Where n is the number of transcriptome data samples of the biological sample, p is the number of characteristic genes of the biological sample, and the matrix The elements in the table represent the expression levels of each characteristic gene in each transcriptome data sample of the biological sample.

[0104] Take any one of the p genes in the matrix X. The vector of expression levels as response variables ,the remaining A matrix of gene expression levels was used as a predictor variable. A resilient network regression model is constructed, where the nodes of the resilient network are each gene. The regression coefficient vector of the nodes is minimized through the following objective function. The objective function is:

[0105]

[0106] in For mixed parameters, The same calculation is performed on the remaining p-1 genes except for gene k, using the penalty coefficient, to obtain their regression coefficient vectors respectively. A symmetric matrix A is constructed using all the regression coefficient vectors, where the symmetric matrix A is the covariance matrix or the precision matrix.

[0107] 3) Weight calculation module: Calculates the squared Mahalanobis distance of the characteristic genes of any transcriptome data sample i among the n transcriptome data samples of the biological sample. The calculation formula is:

[0108]

[0109] in It is a vector of the expression levels of characteristic genes in transcriptome data sample i; the same calculation is performed on any transcriptome data sample j among the n transcriptome data samples of the biological sample to obtain the square of its Mahalanobis distance;

[0110] Further calculate the specific weights of the transcriptome data sample i. and the specific weights of the transcriptome data sample j The calculation formula is:

[0111] ;

[0112] The normalized specific weights of the transcriptome data sample i are further calculated. and the normalized specific weights of the transcriptome data sample j The calculation formula is:

[0113]

[0114] ;

[0115] 4) Transcriptome distance calculation module: Calculates the weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j. i-j The weighted transcriptome distance EN-wTSD is calculated using the following formula:

[0116] EN-wTSD

[0117] in, For weighted Jensen-Shannon divergence The weighted correlation distance is calculated based on the weighted Pearson correlation coefficient, where γ is the coefficient.

[0118] in The calculation formula is:

[0119]

[0120] in and These are the probabilistic vector distributions of gene expression levels for transcriptome data sample i and transcriptome data sample j, respectively. The function is , For the sample in the gene The expression probability distribution on M is and Constructed hybrid distribution:

[0121]

[0122] ,in The calculation formula is:

[0123]

[0124] in The weighted Pearson correlation coefficient;

[0125] 5) Output terminal, output the EN-wTSD.

[0126] In the implementation scheme, the mixing parameter α is in the range of 0-1, preferably 0.5.

[0127] In the implementation scheme, the range of γ is 0-1, preferably 0.5.

[0128] ⑦. A system for analyzing transcriptomic similarity between a first biological sample and a second biological sample, wherein the first biological sample and the second biological sample each comprise at least two transcriptomic data samples, the system comprising:

[0129] 1) Input terminal, used to receive transcriptome data from each transcriptome data sample, including gene name and gene expression level;

[0130] 2) Matrix construction module, as defined in ⑤;

[0131] 3) Weight calculation module, as defined in ⑤;

[0132] 4) Transcriptome distance calculation module, as defined in ⑤;

[0133] 5) Similarity index calculation module: Calculates quantile thresholds based on EN-wTSD values. and similarity ratio The similarity ratio As defined in ③, the quantile threshold As defined in ③;

[0134] 6) Output terminal: Output the EN-wTSD and quantile threshold values. and similarity ratio .

[0135] ⑧. The system described in any one of ⑤-⑦, wherein the transcriptome data samples are derived from organ-on-a-chip, organoids, assemblies, animal tissues, or ex vivo human tissues.

[0136] 9. A computer device comprising at least one processor and a memory storing computer instructions executable on the processor, the instructions, when executed by the processor, implementing the steps of any one of the methods described in ①-④.

[0137] ⑩. A computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of any one of ①-④.

[0138] The elastic network-weighted transcriptome distance (EN-wTSD) used in this invention has significant technical advantages:

[0139] 1. Significantly reduced baseline noise, establishing a more accurate biological baseline. By introducing a regularization mechanism using an elastic network, this invention can more effectively distinguish between "biological signals" and "statistical noise." Validation data shows that in self-comparisons of similar samples (such as human-to-human controls), the distance values ​​calculated by this invention show a significant decreasing trend compared to traditional methods (baseline reduction of more than 3 times). This means that a purer "biological zero point" can be established, effectively reducing false positive results caused by algorithmic errors, making the similarity assessment between the model and the human body more objective.

[0140] 2. Improved sensitivity for detecting subtle transcriptomic drift. Thanks to more robust weight estimation, this invention significantly improves the resolution of difference measures (statistical resolution improved by more than 10 orders of magnitude). When dealing with small batch effects or drug-induced early transcriptomic changes, this invention can capture subtle feature drifts that are difficult to identify using traditional methods. Experiments show that, on the same dataset, this invention can significantly improve the statistical significance of subtle differences, thus providing a more sensitive quantitative indicator for drug screening and efficacy evaluation.

[0141] 3. Enhanced robustness to collinear gene data. Gene expression data naturally exhibits high collinearity. Traditional Glasso algorithms are prone to instability in weight estimation when processing such high-dimensional features. This invention utilizes the properties of elastic networks combined with L1 and L2 regularization, making it suitable for processing gene network data with group effects. It can handle strongly correlated gene clusters while preserving key features, thereby obtaining weight allocations that are more consistent with biological realities and improving the robustness of the evaluation system.

[0142] 4. The algorithm converges quickly, significantly improving computational efficiency. The node regression strategy employed in this invention boasts superior computational performance in its algorithm implementation. Compared to the computationally complex Glasso iteration process in traditional methods, this invention significantly reduces the computation time for large-scale transcriptome data while maintaining accuracy, making it more suitable for real-time analysis of large-scale cross-species databases.

[0143] 5. The algorithm of this invention has strong versatility and breaks through the limitation of single application. This computational framework does not depend on a specific experimental model and can be applied to evaluate the transcriptome similarity between any biological experimental model, such as organ-on-a-chip, organoids, and 3D cultured tissues, and the target species, thus having broad applicability. Attached Figure Description

[0144] Figure 1 This is a flowchart of transcriptome similarity analysis using the EN-wTSD method.

[0145] Figure 2 Figure A shows the clustering diagram of principal component analysis performed in Example 1. Figure B shows the clustering diagram before preprocessing and the clustering diagram after preprocessing. PC1 represents the first principal component, PC2 represents the second principal component, and the disease group and healthy group represent the diseased samples and healthy samples in each dataset, respectively.

[0146] Figure 3 This is a heatmap of pairwise distances between EN-wTSD samples obtained from the cross-species analysis of group 1 in Example 1.

[0147] Figure 4 The box plots are for cross-species analysis results of group 1 using different algorithms in Example 1, where human 1, human 2, chip, mouse 1 and mouse 2 indicate the datasets in Table 1. Figure A shows the results using the EN-wTSD algorithm, Figure B shows the results using the GL-wTSD algorithm, Figure C shows the results using the unweighted TSD algorithm, and Figure D shows the results using the ED algorithm.

[0148] Figure 5 The box plots are for cross-species analysis results of group 2 in Example 1 using different algorithms, where human 1, human 2, chip, mouse 1 and mouse 2 indicate the datasets in Table 1. Figure A shows the results using the EN-wTSD algorithm, Figure B shows the results using the GL-wTSD algorithm, Figure C shows the results using the unweighted TSD algorithm, and Figure D shows the results using the ED algorithm.

[0149] Figure 6 Figure A shows the relationship between the relative fluctuation ratio of node degree and network density under different edge strength thresholds. Figure B shows the relationship when the edge strength threshold is 0.001, Figure C shows the relationship when the edge strength threshold is 0.005, and Figure C shows the relationship when the edge strength threshold is 0.01. Detailed Implementation

[0150] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention. In this application, characteristic genes refer to a number of genes selected for a specific purpose. For example, when evaluating the applicability of an animal model for a specific disease, genes related to that specific disease can be selected as specific genes.

[0151] In this application, the target species is a species selected based on a specific purpose. For example, humans are the target species when evaluating the applicability of animal models for specific human diseases.

[0152] Biological samples are biological substances obtained directly from or derived from organisms, their tissues, cells, or body fluids.

[0153] A transcriptome data sample is a sample of transcriptome data that includes samples of a specific organism.

[0154] Transcriptome data refers to any RNA expression information in a specific biological sample.

[0155] In this application, biological samples include organ-on-a-chip, organoids, assemblies, animal models, or human models, wherein the human models are derived from ex vivo human tissue.

[0156] In this application, the model includes biological models (e.g., mouse, cynomolgus monkey, etc.) or experimental models (e.g., organoids, organoid assemblies); for reference models and comparative models, for example, when comparing the physiological / pathological similarity between a mouse model and a human (to evaluate the feasibility of using the mouse model for screening drugs for humans), the human is the target species, and in the comparison, the human constitutes the reference model and the mouse model is the comparative model.

[0157] A preferred embodiment of the transcriptome similarity analysis method of the present invention comprises the following three parts:

[0158] (1) Cross-species data collection and standardization

[0159] To construct high-quality computational input, the multi-source transcriptomics data underwent standardization and cleaning, specifically including:

[0160] Multi-source data fusion: Collect transcriptome data from biological samples and construct a unified expression matrix across species using homologous gene mapping tools (such as biomaRt).

[0161] Preprocessing (including normalization and debaterization): The expression matrix is ​​normalized and debaterized (e.g., using tools such as ComBat, TMM, DESeq2, logCPM, etc.). The debaterized result is then transformed to be suitable for distance calculation (e.g., variance stabilization transformation) to obtain a unified expression matrix.

[0162] Preferably, gene screening can also be performed based on a reference gene set to obtain characteristic genes for constructing the matrix. The reference gene set can be used to construct a gene set for the sample's gene signature (i.e., a set of genes with co-expression patterns), used for subsequent calculation of the probability distribution of gene expression and the gene ranking vector. The reference gene set can be derived from tissue / organ-specific signature gene sets in authoritative databases, or a custom reference gene set can be constructed according to the research direction. Relevant classic biomedical databases include The Human Protein Atlas (HPA) and CellMarker.

[0163] (2) Calculation of transcriptome feature distance based on elastic network weighting (EN-wTSD)

[0164] This section is the core of the invention, aiming to calculate the transcriptomic feature distance between the model and the target species. Traditional TSD algorithms rely on graphical lasso estimation of the accuracy matrix, which is sensitive to noise in high-dimensional data (where the number of genes is much greater than the number of samples) and prone to overfitting. This invention introduces an elastic network algorithm to reconstruct the weight learning mechanism, with the specific steps as follows:

[0165] Step 1: Construct a weight learning model based on elastic networks. This invention employs a node-wise regression strategy combined with elastic network regularization to estimate the accuracy matrix between genes. This precision matrix reflects the partial correlation structure between genes and is used to subsequently evaluate the specific weights of samples, thus reflecting the confidence level of the samples.

[0166] Input: Based on the preprocessed genes in part (1), obtain the gene expression matrix. ,in For the sample size, The number of characteristic genes.

[0167] Each gene is considered a node in the network. Conditional dependencies between nodes are estimated using resilient network node regression, thereby determining the network's structural relationships. The expression level of any gene k among the p genes in the matrix is ​​used as the response variable. ,the remaining The expression level of each gene was used as a predictor variable. The following elastic network regression function is used to construct an elastic network regression model, and the regression coefficient vector is solved by minimizing the following objective function. :

[0168]

[0169] in: This is a vector of regression coefficients; For genes The vector of expression levels across all samples. In addition to genes The expression levels of other genes in all samples.

[0170] For the mixed parameter, α ranges from 0 to 1, and is preferably set as follows: This setup combines the sparsity feature selection capability of LASSO regression (L1 regularization) with the group effect capability of ridge regression (L2 regularization) in handling collinear features, making it more robust to handling strongly correlated gene clusters in gene regulatory networks compared to Glaasso.

[0171] The penalty coefficient is adopted. ,pass Cross-validation automatically selects the corresponding value with the smallest mean square error. .

[0172] Output: Based on the above elastic network regression model, solve for the regression coefficient vector. They are arranged into a symmetric matrix A, which is the covariance matrix. or precision matrix (When the sample size is small (p is 10 times n or more), the precision matrix is ​​better; when the sample size is large (p is less than 10 times n), the covariance matrix is ​​better.)

[0173] When the precision matrix is ​​used, it can be a sparse precision matrix. The obtained regression coefficient vector is filled into the precision matrix to obtain the initial sparse precision matrix. Then, perform symmetry processing on it to obtain The formula for symmetry processing is:

[0174]

[0175] When it is a covariance matrix, the obtained sub-regression coefficient vectors are filled into the covariance matrix to obtain the covariance matrix.

[0176] Step 2: Calculate the specificity weights of the samples. Using matrix A obtained in Step 1, calculate the squared Mahalanobis distances of the characteristic genes in the transcriptome data samples. The calculation results were then used to construct sample-specific weights, thereby reducing the impact of outlier samples on cross-species transcriptome distance metrics. The specific calculation formula is as follows:

[0177] ,

[0178] in It is a vector of the normalized expression levels of characteristic genes in transcriptome data sample i (i.e., 1×p). i The matrix, p i (where A is the number of characteristic genes in sample i). or .

[0179] Then, the square of the Mahalanobis distance between the transcriptome data samples is converted into sample weights. This weight reflects the statistical "typicality" of the sample. A larger weight value indicates that the transcriptome data sample is closer to the distribution center, thus reducing the interference of outliers or low-quality data on similarity assessment. The weight of transcriptome data sample i is calculated using an exponential mapping function, and its calculation formula is as follows:

[0180] ,

[0181] The exponential mapping function described herein adopts a form consistent with the exponential kernel of the normal distribution, and .

[0182] In pairwise sampling (transcriptome data samples) With transcriptome data samples When comparing, the final weights are calculated through normalization:

[0183]

[0184] .

[0185] Step 3: Combine sample weights and gene ranking features to calculate the final EN-wTSD distance. EN-wTSD consists of two parts: weighted Jensen-Shannon divergence (wJSD) and weighted ranking correlation distance (wRCD).

[0186] Weighted Jensen-Shannon divergence (wJSD) is used to measure the difference in gene expression probability distributions. The generalized wJSD known in the art is expressed as:

[0187]

[0188] in The expression probability distributions of two transcriptome data samples to be compared (transcriptome data samples i and j) are shown. For sample weights, It is the first to participate in the mixing There are several probability distributions.

[0189] Traditional JSD assumes that the weights of the two distributions are equal, while this invention uses the weights calculated above. and Constructing a hybrid distribution :

[0190]

[0191] in and These are the vector distributions of probabilistic gene expression levels for two transcriptome data samples, respectively. The improved weighted Jensen-Shannon divergence of this invention is defined as the two distributions relative to a mixed distribution. Sum of weighted Kullback-Leibler (KL) divergences:

[0192] ,

[0193] in, This refers to the KL divergence (Kullback-leibler divergene). , Represents a specific characteristic gene k. Representative transcriptome data samples in genes The expression probability distribution is shown. After introducing the weights calculated by the elastic network, wJSD can effectively suppress baseline noise.

[0194] Weighted Relevance Distance (wRCD): Used to measure the consistency of the ranking of relative gene expression abundance across different samples.

[0195]

[0196] in The weighted Pearson correlation coefficient is calculated using the method described in (Bailey P, Emad A, 2023), specifically as follows: For each transcriptome data sample... With transcriptome data samples The expression levels of each gene are sorted to obtain the corresponding ranking vector. and Subsequently, based on the combined gene weights As observation weights, for and Calculate the weighted Pearson correlation coefficient The weighted correlation is to... By directly incorporating the estimation process of weighted mean, weighted covariance, and weighted variance, we obtain... The specific formula for wRCC is: ,in

[0197] Gene weight vector It is calculated based on the variance stability of the ranking of characteristic genes in each transcriptome data sample. The specific calculation method is as follows: for any gene Based on its expression level, the ranking of gene k among all genes is obtained, and the mean ranking of gene k in each transcriptome data sample is calculated. ) and standard deviation ( Assuming the distribution of gene rankings approximately follows a normal distribution. For a specific transcriptome data sample genes (The specific ranking is as follows) The likelihood of ranking is calculated using the probability density function of the normal distribution. To smooth the numerical values ​​and construct an additive model, a logarithmic transformation was performed on the likelihood to obtain the gene weights of transcriptome data sample i. :

[0198]

[0199] Likelihood of ranking gene k in transcriptome data sample j Perform the same calculations to obtain the gene weights of transcriptome data sample j. To further obtain the combined gene weights It is composed of the weighted sum of transcriptome data sample i and transcriptome data sample j:

[0200]

[0201] The EN-wTSD is calculated using the following formula:

[0202] EN-wTSD

[0203] γ is a coefficient, 0 < γ < 1, preferably 0.5.

[0204] The closer the EN-wTSD value is to 0, the more similar the transcriptomic features of different transcriptomic data samples are. The method for obtaining EN-wTSD described above is named the EN-wTSD algorithm.

[0205] (3) Evaluation of transcriptome similarity

[0206] After obtaining the EN-wTSD distance values ​​among all samples, a similarity assessment is further performed. The purpose is usually to evaluate the substitutability or equivalence of the model in simulating the physiological processes of the target species (e.g., humans), but other purposes are not excluded, such as providing data for subsequent research. There are three methods for similarity assessment: Method 1, assessing similarity based on the mean or median of the EN-wTSD; Method 2, assessing similarity based on the mean or median of the EN-wTSD and the quantile threshold T. sim Method 3: Evaluate similarity based on similarity ratio. Evaluate similarity.

[0207] Method 1 includes the following steps:

[0208] Calculate the inter-model transcriptome distance between any human transcriptome data sample and any mouse transcriptome data sample, obtain the transcriptome distance for all the inter-model samples, and further calculate their average EN-wTSD. 平均Or the median value EN-wTSD 中位 The EN-wTSD 平均 Or EN-wTSD 中位 The closer to 0, the higher the transcriptome similarity between the first biological sample and the second biological sample;

[0209] Method 2 includes the following steps:

[0210] ① Select the target species (usually human), obtain the dataset of the target species, calculate the EN-wTSD distance between each sample in the reference model group H (e.g., human-to-human EN-wTSD), and construct a distance set. ; Calculate the EN-wTSD distance between the comparison model group M and the reference model group H (e.g., EN-wTSD of human-chip).

[0211] ②Based on the aforementioned distance set Calculate the quantile threshold It is used to characterize the baseline fluctuation range of the reference model group. ,in This represents the quantile function, indicating that the proportion of a sample not exceeding a certain threshold is at least 1 / 3. The minimum threshold, Possible values ​​include 0.75, 0.90, 0.95, 0.99, or 0.995, etc.

[0212] ③ Similarity evaluation based on the threshold: When the mean or median EN-wTSD distance between the comparison model group M and the reference model group H is... When the transcriptomes of the comparison model M and the reference model group H are determined to be similar; when At that time, it was determined that the transcriptomes of the comparison model M and the reference model group H were not similar.

[0213] Method 3 includes the following steps:

[0214] ① Calculate the comparison model group (Selected from organ-on-a-chip models, organoid models, assembly-like models, or animal models, etc.) and reference model set (e.g., human models) cross-group EN-wTSD, and construct their set. ;

[0215] ② Calculate the similarity index (similarity ratio) for model selection and ranking. Obtain the set. middle Given the number of sample pairs Z, calculate the similarity ratio between different models or different datasets. The calculation formula is: , Represents a set The total number of sample pairs. The value ranges from 0 to 1, with values ​​closer to 1 indicating higher similarity and values ​​closer to 0 indicating lower similarity.

[0216] Figure 1 The architecture and process of a system for evaluating the similarity of transcriptomes from different datasets are presented.

[0217] Example 1:

[0218] Cross-species physiological suitability assessment of organ-on-a-chip for alcoholic liver disease based on elastic network weighting.

[0219] This embodiment utilizes the computational framework proposed in this invention to perform cross-species transcriptome similarity quantification assessment on ALD organ-on-a-chip models, traditional mouse models, and real human pathological tissues, taking alcoholic liver disease (ALD) as an example. This embodiment aims to verify the technical advantages of the improved transcriptome feature distance (EN-wTSD) algorithm described in this invention in reducing baseline noise, improving sensitivity to subtle differences, and increasing computational efficiency.

[0220] Alcoholic liver disease (ALD) is a leading cause of chronic liver disease worldwide, with a 5-year mortality rate exceeding 50%, and even higher in later stages of the disease. Late-stage ALD often goes undetected clinically before symptoms appear, highlighting the importance of prevention, treatment, and prognosis. This embodiment focuses on an organ-on-a-chip study of ALD, conducting a cross-species analysis of its physiological compatibility. The specific steps include:

[0221] S1: Cross-species data collection. Organ microarray data, human clinical data on alcoholic liver disease, and transcriptomic data on mouse alcoholic liver disease were collected from the GEO database, as shown in Table 1 below:

[0222] Table 1. Data Information

[0223]

[0224] Construct a matrix,

[0225] ,in, The number of transcriptome data samples for biological samples. The number of characteristic genes in a biological sample, elements Indicates biological sample Middle characteristic genes The standardized expression value.

[0226] S2: Data preprocessing, mainly including normalization. The DESeq2 software is used with a size factor normalization method. The original gene expression levels are normalized, and further stabilization and batch effect removal can be performed, such as logarithmic transformation. Figure 2 The cluster diagrams of principal component analysis before and after preprocessing are shown. It can be seen that the preprocessing better reflects the similarity and differences between groups, indicating that the preprocessing achieved its intended effect.

[0227] S3: Gene screening. Obtain the specific gene set corresponding to the organ tissue from the Human Protein Atlas database as the feature gene set; then take the intersection of the feature gene set with the genes in the above dataset to determine the final feature genes used for distance calculation.

[0228] S4: WTSD distance calculation based on elastic networks. Distance metrics are measured for "human-human," "human-chip," and "human-mouse" pairs. The regression coefficient vector of the minimized nodes is obtained using the above method (mixture parameter α=0.5, optimal value automatically selected through cross-validation). Construct a sparse precision matrix :

[0229]

[0230] Then construct its symmetric matrix: .

[0231] Weight learning: Input the feature gene matrix and calculate the specific weights between each sample. Compared with the GL-wTSD algorithm, this step effectively solves the problem of weight estimation distortion caused by gene collinearity.

[0232] Distance synthesis: Combining weighted Jensen-Shannon divergence (wJSD) and weighted rank-related distance (wRCD), the EN-wTSD value between each sample is output. Furthermore, a transcriptome distance (TSD) heatmap and a transcriptome distance box plot can be generated; see [link to relevant documentation]. Figures 3-5 .

[0233] S4: Transcriptome similarity assessment, with specific implementation details including:

[0234] Based on the EN-wTSD values ​​obtained from the above calculations, a heatmap is created (see [link]). Figure 3The diagram visually illustrates the distance between transcriptome samples. Humans and the chip are relatively closer; mice 1 and 2 are closer, while the chip is further away from both types of mice, showing greater separation. For specific transcriptome data pairs, such as the EN-wTSD value of human 9 and chip 2 (0.22) and the EN-wTSD value of human 9 and mice 1-5 (0.33), the transcriptome similarity between chip 2 and human 9 is higher than that between mice 1-5 and human 9.

[0235] Based on TSD box plot ( Figure 4 and 5 As can be seen, in the transcriptome feature space, the distance between ALD organ microarray and human pathological samples is significantly smaller than that of mouse models (EN-wTSD). 人-芯片 <EN-wTSD 人-小鼠 This indicates that the organ-on-a-chip has higher physiological adaptability than traditional mouse models in simulating the pathological features of human ALD. Furthermore, compared with other algorithms, the TSD algorithm with elastic network weighting calculates the smallest transcriptomic distance between human samples and can distinguish differences between mice in different datasets (the transcriptomic distances between human 1-mouse 1 and human 1-mouse 2 differ significantly in the EN-wTSD algorithm), which verifies the stability and sensitivity of the EN-wTSD algorithm.

[0236] To verify the technological advancement of the proposed EN-wTSD algorithm, it was compared with Euclidean distance (ED), unweighted TSD algorithm, and GL-wTSD algorithm relying on Glasso. The data was divided into two groups (Group 1 includes data from chips, human 1, mouse 1, and mouse 2 from Table 1; Group 2 includes data from chips, human 2, mouse 1, and mouse 2 from Table 1) for analysis. The Euclidean distance algorithm is described in (Bailey P and Emad A, 2023), and the unweighted TSD algorithm and GL-wTSD algorithm are described in (Manatakis et al., 2020).

[0237] Ideally, the transcriptome distance between similar samples (human-to-human) should be close to 0, and a transcriptome distance greater than 0 represents baseline noise. As shown in Tables 2 and 3, the mean human-to-human group distance calculated by the unweighted TSD algorithm is 0.261, and the mean human-to-human group transcriptome distance calculated by GL-wTSD is 0.270, both exhibiting high systematic baseline noise. In contrast, the EN-wTSD algorithm of this invention reduces the mean human-to-human group transcriptome distance to 0.097 (Group 1) and 0.084 (Group 2), with significantly lower baseline noise than other algorithms, reducing baseline noise by approximately 64%-69%. This establishes a more accurate evaluation benchmark and effectively avoids false positives caused by algorithm errors.

[0238] according to Figure 2 , 3 The results in 4 suggest that, compared with the mouse model, organ-on-a-chip can better simulate the overall condition of diseased tissues in the human body; and compared with other algorithms, the transcriptome distance between human samples calculated by the EN-wTSD method is the smallest (i.e., the human-to-human transcriptome distance value is smaller), which is more in line with the actual situation.

[0239] Compared to other algorithms, the EN-wTSD algorithm of this invention significantly improves the sensitivity of difference resolution, demonstrating higher statistical power when distinguishing subtle differences across species or batches. The Wilcoxon rank-sum test was used to examine the significance of differences. As shown in Table 5, when comparing subtle differences between the "human 2-mouse 1" and "human 2-mouse 2" groups, the significance P-value of GL-wTSD was 0.019 (P<0.05), while the P-value of the EN-wTSD algorithm of this invention reached 5.65 × 10⁻⁶. -15 That is, the statistical significance of this invention is improved by more than 10 orders of magnitude. This indicates that in scenarios requiring high sensitivity, such as drug screening, this invention can capture subtle transcriptome drift that is difficult to identify using traditional methods. The fact that the ED algorithm and the unweighted TSD algorithm could not detect the subtle differences between the "human-mouse 1" and "human-mouse 2" mouse data (statistical significance was not significant) also demonstrates the advantage of EN-wTSD. Although GL-wTSD has a smaller standard deviation in some results, this often stems from "oversmoothing" caused by its forced positive definite process. This invention's EN-wTSD, while reducing the mean, retains a moderate degree of dispersion (e.g., in Table 2, for human 1-chipset, the GL standard deviation is 0.024, and the EN standard deviation is 0.057), which better reflects the objective laws of individual differences in biological samples and avoids "overfitting" masking true biological heterogeneity.

[0240] The EN-wTSD algorithm is more computationally efficient. As shown in Table 6, when processing datasets of the same size, the GL-wTSD algorithm takes approximately 970 seconds, while the EN-wTSD algorithm of this invention only takes about 200 seconds, representing a speedup of approximately 4.5 to 4.9 times (an efficiency improvement of approximately 80%). This is because the elastic network uses the coordinate descent method for solving, avoiding the complex maintenance of matrix positive definiteness during iteration in Glasho, making it more suitable for rapid analysis of large-scale omics data.

[0241] S5: Similarity assessment, used to comprehensively evaluate the equivalence or substitutability of the model under evaluation compared to the human model. The similarity assessment is performed based on the results of the EN-wTSD analysis, using quantile thresholds. Conduct a similarity evaluation.

[0242] Median of EN-wTSD and Similarity evaluation: quantiles are taken Based on the set of distances within the reference group (person1-person1) 0.95 quantile threshold For each group to be evaluated, the distance set of "person 1-model" is... Take the median as the group-level summary statistic and compare it with... A qualitative judgment was made by comparison, and the results are shown in Table 7. Figure 4 and 5 As shown. The results are shown in the summary table: Person 1-Person 2 The median is 0.20, and the number of chips is 1. The median is 0.24, which satisfies all conditions. Therefore, it was judged as "similar"; conversely, human 1-mouse 1 was similar. The median was 0.33, and the median for human 1-mouse 2 was... A median of 0.30 does not meet the threshold condition. Therefore, it was determined to be "dissimilar". This conclusion indicates that the homologous human datasets (human 1-human 2) and the organ-on-a-chip model (human 1-chip) can all reach the similarity standard defined by the 0.95 quantile threshold of the reference group at the transcriptome level. However, the transcriptome differences between the two mouse models and human tissues exceed this similarity threshold, and they are not equivalent to the reference human tissues as a whole.

[0243] Similarity is evaluated based on similarity ratio: quantiles are used. Based on the obtained quantile threshold Evaluations were conducted on any pair of samples to be evaluated (e.g., "human-chip", "human-mouse 1", "human-mouse 2", etc.), and the results are shown in Table 8. Figure 4 and 5 As shown, the similarity ratio between human 1 and human 2 is 0.88, indicating that the similarity between different datasets from the same source is high, which is consistent with biological norms. The similarity ratio between the chip and human 1 is about 0.41, which is significantly higher than the similarity ratios of other groups (human 1-mouse 1, human 1-mouse 2), indicating that the similarity between the organ-on-a-chip model and human tissue is significantly higher than that between the mouse model and human tissue.

[0244] In summary, this embodiment demonstrates through data that the proposed resilient network-based wTSD algorithm significantly outperforms existing technologies in baseline accuracy, detection sensitivity, and computational efficiency. The EN-wTSD algorithm provides a powerful tool for the standardized evaluation of cross-species experimental models. Furthermore, it verifies the effectiveness of ALD organ-on-a-chip as a human disease model.

[0245] Table 2. Mean and standard deviation of transcriptome distance for each algorithm in Group 1

[0246]

[0247] Table 3. Mean and standard deviation of transcriptome distance for each algorithm in Group 2

[0248]

[0249] Table 4. Statistical significance test results (p-value and significance) of transcriptome distance results for each algorithm in Group 1.

[0250]

[0251] ns indicates no significant difference (p>0.05), *** indicates p<0.001, ** indicates p<0.01, and * indicates p<0.05.

[0252] Table 5. Statistical significance test results (p-value and significance) of transcriptome distance results for each algorithm in Group 2.

[0253]

[0254] ns indicates no significant difference (p>0.05), *** indicates p<0.001, ** indicates p<0.01, and * indicates p<0.05.

[0255] Table 6. Comparison of operating efficiency between EN-wTSD and GL-wTSD

[0256]

[0257] Table 7. Transcriptome similarity determination between each model and human tissues ( )

[0258]

[0259] Table 8. Transcriptome similarity ratio between each model and human tissue

[0260]

[0261] Thanks to the technological advantages of elastic networks in reducing baseline noise and improving sensitivity, the transcriptomic similarity distance between each model and real human pathological tissues was successfully quantified. The calculation results intuitively and significantly demonstrate that this organ-on-a-chip has a high degree of similarity to human pathological tissues in gene expression profiles and regulatory networks, and in some aspects, its physiological simulation efficacy is superior to traditional experimental mouse models.

[0262] Example 2

[0263] In practical applications, gene expression data may be subject to rejection or replenishment due to quality control failures of individual samples (i.e., changes in a small number of samples). To evaluate the robustness of the model, the fluctuation of key features of the output network is assessed under the condition that only a small number of samples change. Specifically, when evaluating or testing the robustness of the EN-wTSD method in handling gene collinearity data, a resampling method is used for comparison. The principle is as follows: from the same dataset, a certain proportion of samples are randomly drawn or sampled multiple times, and the proportion drawn each time is consistent. Gene relationship networks are constructed using both the GL-wTSD method and the EN-wTSD method of this application. To ensure fair comparison, the density of the networks output by the two methods is forced to be the same in each drawing or sampling, that is, the number of connections between nodes and the overall connectivity are completely consistent, avoiding bias in conclusions due to differences in network sparsity. Subsequently, the network features obtained by the two methods in multiple samplings are compared to see if there are fluctuations; the smaller the fluctuation, the stronger the robustness.

[0264] The data used is from Table 1 in Example 1. Unordered repeated sampling was employed to simulate collinearity data processing, evaluating the robustness of the EN-wTSD method for handling gene collinearity data. Specifically, gene relationship networks were constructed using both the GL-wTSD method and the EN-wTSD method of this application. From a fixed sample set, sampling was repeated 50 times, with 80% of the samples randomly selected in each sampling to form a subset. Compared to the overall sample, each sampled subset differs by only a small number of samples (e.g., 20%). Therefore, different sampling iterations effectively simulate input perturbations caused by small increases or decreases in the number of samples.

[0265] Both the GL-wTSD method and the EN-wTSD method of this application output a gene network on each subset of the dataset, and further output a matrix of edge strength. For the GL-wTSD method, the edge strength is given by the absolute value of the off-diagonal elements of the accuracy matrix; for the EN-wTSD regression method, the edge strength is obtained by symmetrically combining the absolute values ​​of the bidirectional regression coefficients. The relative fluctuation ratio of node degree is obtained from the edge strength matrix to evaluate the robustness of the two methods.

[0266] To ensure a fair comparison, both methods retain the same number of edges in each subset sampling to avoid bias caused by inconsistent network densities. Each node (gene) in the network is connected by several edges. The number of edges (i.e., node degree) represents the number of edges associated with a gene in the network. If the algorithm is robust, the node degree of the same gene should remain relatively stable after 50 samplings, without significant fluctuations. The variation in node degree for each gene across the 50 samplings is statistically analyzed, and then divided by the average node degree of that gene across the 50 samplings to obtain the relative fluctuation ratio of node degree. The smaller this ratio, the more stable the node degree of the gene, and the more stable the algorithm.

[0267] First, a threshold for edge strength (k) is set (e.g., 0.01, 0.005, or 0.001). This threshold is used to remove very weak connections that are close to noise, retaining only edges with edge strength greater than the threshold. The smaller value between the edge strength values ​​retained by the two methods is used to set the edge strength for both methods, avoiding bias caused by different network densities. Then, the relative fluctuation ratio of node degree is calculated using the method described above. The results are shown in Table 9 and... Figure 6 As shown, under different thresholds, the mean fluctuation ratio of gene node degrees obtained by the EN-wTSD method is significantly lower than that of the GL-wTSD method. The relative fluctuation ratio of gene node degrees obtained by the EN-wTSD method is also more concentrated (compared to the GL-wTSD method). This indicates that the elastic network provides a more stable estimate of the connectivity of gene network nodes when the samples change slightly, and the EN-wTSD method has stronger robustness. Collinearity data often exists in gene expression data, manifested as highly synchronized changes in multiple genes (possibly due to interactions between genes). For processing this type of collinearity data, the GL-wTSD model estimation may lead to frequent changes in certain connectivity relationships between different sampling times, resulting in significant fluctuations. However, the EN-wTSD method, by simultaneously introducing L1 and L2 regularization penalty terms, results in smaller fluctuations in regression coefficients when the input variables are highly similar (i.e., collinear), thus better maintaining the stability of node connectivity and exhibiting stronger robustness.

[0268] Table 9

[0269]

[0270] The system constructed using the above method simplifies the complex cross-species analysis process and has significant promotional value and technical universality in organoid / organ-on-a-chip quality control, drug screening, and mechanism research.

Claims

1. A method for analyzing transcriptomic similarity between two transcriptomic data samples, respectively from a first biological sample and a second biological sample, wherein the first biological sample and the second biological sample each comprise at least two transcriptomic data samples, the method comprising: Calculate the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample, respectively, and determine the transcriptome similarity between the two transcriptome data samples; The method for calculating the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample, respectively, includes the following steps: (1) Obtain the expression matrix of characteristic genes in the transcriptome of each sample of the first biological sample. Where n1 is the number of transcriptome data samples of the first biological sample, p1 is the number of characteristic genes of the first biological sample, and the elements in matrix X1 are the expression levels of each characteristic gene in each transcriptome data sample of the first biological sample; using the same method, a matrix of expression levels of characteristic genes in each transcriptome data sample of the second biological sample is obtained. Where n2 is the number of transcriptome data samples of the second biological sample, p2 is the number of characteristic genes of the second biological sample, and the elements in the matrix X2 are the expression levels of each characteristic gene in each transcriptome data sample of the second biological sample. (2) Take any one of the p1 genes in the matrix X1. The vector of expression levels as response variables The expression matrices of the remaining p1-1 genes were used as predictor variables. A resilient network regression model is constructed, where the nodes of the resilient network are each gene. The regression coefficient vector of the nodes is minimized through the following objective function. The objective function is: in For mixed parameters, As the penalty coefficient, the same calculation is performed on the remaining p1-1 genes except for gene k to obtain their regression coefficient vectors. A symmetric matrix A1 is constructed using all the regression coefficient vectors. The symmetric matrix A1 is a covariance matrix or a precision matrix. Using the same method, a symmetric matrix A2 is obtained from the matrix X2. The symmetric matrix A2 is a covariance matrix or a precision matrix. (3) Calculate the square of the Mahalanobis distance of the characteristic genes of any transcriptome data sample i among the n1 transcriptome data samples of the first biological sample. The calculation formula is: in Transcriptome data samples The vector of expression levels of characteristic genes; the same calculation is performed on any transcriptome data sample j from the n2 transcriptome data samples of the second biological sample to obtain the square of its Mahalanobis distance; further, the specific weight of transcriptome data sample i is calculated. and the specific weights of the transcriptome data sample j The calculation formula is: ； The normalized specific weights of the transcriptome data sample i are further calculated. and the normalized specific weights of the transcriptome data sample j The calculation formula is: ； (4) Calculate the weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j. i-j The weighted transcriptome distance EN-wTSD is calculated using the following formula: EN-wTSD in, For weighted Jensen-Shannon divergence The weighted correlation distance is calculated based on the weighted Pearson correlation coefficient, where γ is the coefficient. in The calculation formula is: in and These are the probabilistic vector distributions of gene expression levels for transcriptome data sample i and transcriptome data sample j, respectively. The function is , For the sample in the gene The expression probability distribution on M is and Constructed hybrid distribution: in The calculation formula is: in The weighted Pearson correlation coefficient; The method for determining the transcriptomic similarity between the two transcriptomic data samples is as follows: The weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j i-j The closer to 0, the higher the transcriptome similarity between transcriptome data sample i and transcriptome data sample j.

2. A method for analyzing transcriptomic similarity between two transcriptomic data samples from a biological sample, said biological sample comprising at least three transcriptomic data samples, the method comprising: Calculate the weighted transcriptome distance between two transcriptome data samples from the biological sample, and determine the transcriptome similarity between the two transcriptome data samples; The method for calculating the weighted transcriptomic distance between two transcriptomic data samples from the biological sample includes the following steps: (1) Obtain the expression matrix of characteristic genes in the transcriptome of each of the biological samples. , where n is the number of transcriptome data samples of the biological sample, p is the number of characteristic genes of the biological sample, and the elements in the matrix X are the expression levels of each characteristic gene in each transcriptome data sample of the biological sample. (2) Take any one of the p genes in the matrix X. The vector of expression levels as response variables The expression matrices of the remaining p-1 genes are used as predictor variables. A resilient network regression model is constructed, where the nodes of the resilient network are each gene. The regression coefficient vector of the nodes is minimized through the following objective function. The objective function is: in For mixed parameters, The same calculation is performed on the remaining p-1 genes except for gene k, using the penalty coefficient, to obtain their regression coefficient vectors respectively. A symmetric matrix A is constructed using all the regression coefficient vectors, where the symmetric matrix A is the covariance matrix or the precision matrix. (3) Calculate the square of the Mahalanobis distance of the characteristic genes of any transcriptome data sample i among the n transcriptome data samples of the biological sample. The calculation formula is: in It is a vector of the expression levels of characteristic genes in transcriptome data sample i; the same calculation is performed on any transcriptome data sample j among the n transcriptome data samples of the biological sample to obtain the square of its Mahalanobis distance; Further calculate the specific weights of the transcriptome data sample i. and the specific weights of the transcriptome data sample j The calculation formula is: ； The normalized specific weights of the transcriptome data sample i are further calculated. and the normalized specific weights of the transcriptome data sample j The calculation formula is: ； (4) Calculate the weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j. i-j The weighted transcriptome distance EN-wTSD is calculated using the following formula: EN-wTSD in, For weighted Jensen-Shannon divergence The weighted correlation distance is calculated based on the weighted Pearson correlation coefficient, where γ is the coefficient. in The calculation formula is: in and These are the probabilistic vector distributions of gene expression levels for transcriptome data sample i and transcriptome data sample j, respectively. The function is , For the sample in the gene The expression probability distribution on M is and Constructed hybrid distribution: in The calculation formula is: in The weighted Pearson correlation coefficient; The method for determining the transcriptomic similarity between the two transcriptomic data samples is as follows: The weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j i-j The closer to 0, the higher the transcriptome similarity between transcriptome data sample i and transcriptome data sample j.

3. A method for analyzing transcriptomic similarity between a first biological sample and a second biological sample, wherein the first biological sample and the second biological sample each comprise at least two transcriptomic data samples, the method comprising: (1) According to the method for calculating the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample, as described in claim 1, calculate the EN-wTSD between all pairwise transcriptome data samples, and further calculate their average EN-wTSD. 平均 Or the median value EN-wTSD 中位 The two transcriptome data samples mentioned therein are either transcriptome data samples of the first biological sample or transcriptome data samples of the second biological sample. The EN-wTSD 平均 Or EN-wTSD 中位 The closer to 0, the higher the transcriptomic similarity between the first biological sample and the second biological sample; (2) According to the method for calculating the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample, as described in claim 1, calculate the EN-wTSD between all pairwise transcriptome data samples of the first biological sample and construct its set D. HH Calculate the quantile threshold T sim , ,in Represents the quantile function; According to the method for calculating the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample, as described in claim 1, the EN-wTSD between pairwise transcriptome data samples from all different biological samples is calculated, and the average EN-wTSD is further calculated. 平均 Or the median value EN-wTSD 中位 The pairwise transcriptome data samples of the different biological samples are either transcriptome data samples of the first biological sample or transcriptome data samples of the second biological sample. If EN-wTSD 平均 ≤T sim or EN-wTSD 中位 ≤T sim If the transcriptome of the second biological sample is similar to that of the first biological sample, then the transcriptome of the second biological sample is determined to be similar; if EN-wTSD 平均 >T sim or EN-wTSD 中位 >T sim If so, it is determined that the transcriptome of the second biological sample is not similar to the transcriptome of the first biological sample; or (3) According to the method for calculating the weighted transcriptome distance between two transcriptome data samples from the first biological sample and the second biological sample, as described in claim 1, the EN-wTSD between all pairwise transcriptome data samples of the first biological sample is obtained, and its set D is constructed. HH Calculate the quantile threshold T sim , ,in Represents the quantile function; According to the method for calculating the weighted transcriptome distance between two transcriptome data samples from a first biological sample and a second biological sample, as described in claim 1, the weighted transcriptome distance EN-wTSD between pairwise transcriptome data samples from all different biological samples is calculated, and its set D is constructed. HM The set is calculated where each pair of transcriptome data samples from different biological samples is either a transcriptome data sample from the first biological sample or a transcriptome data sample from the second biological sample. middle The number of pairwise transcriptome data samples Z was used to further calculate the transcriptome similarity ratio between the first biological sample and the second biological sample. The calculation formula is: ,in For set The total number of pairwise transcriptome data samples; Similarity ratio The higher the value, the greater the transcriptome similarity between the first biological sample and the second biological sample.

4. The method of any one of claims 1-3, wherein the transcriptome data samples are derived from organ-on-a-chip, organoids, assemblies, animal tissues, or ex vivo human tissues.

5. A system for analyzing transcriptomic similarity between two transcriptomic data samples, respectively from a first biological sample and a second biological sample, wherein the first biological sample and the second biological sample each comprise at least two transcriptomic data samples, the system comprising: 1) Input terminal, used to receive transcriptome data from each transcriptome data sample, including gene name and gene expression level; 2) Matrix construction module: Constructs matrices of the expression levels of characteristic genes in the transcriptomes of each transcriptome data sample from the first biological sample. Where n1 is the number of transcriptome data samples of the first biological sample, p1 is the number of characteristic genes of the first biological sample, and the elements in matrix X1 are the expression levels of each characteristic gene in each transcriptome data sample of the first biological sample; using the same method, a matrix of expression levels of characteristic genes in each transcriptome data sample of the second biological sample is obtained. Where n2 is the number of transcriptome data samples of the second biological sample, p2 is the number of characteristic genes of the second biological sample, and the elements in the matrix X2 are the expression levels of each characteristic gene in each transcriptome data sample of the second biological sample. Take any one of the p1 genes in the matrix X1. The vector of expression levels as response variables The expression matrices of the remaining p1-1 genes were used as predictor variables. A resilient network regression model is constructed, where the nodes of the resilient network are each gene. The regression coefficient vector is minimized through the following objective function. The objective function is: in For mixed parameters, As the penalty coefficient, the same calculation is performed on the remaining p-1 genes except gene k to obtain their regression coefficient vectors. A symmetric matrix A1 is constructed using all the regression coefficient vectors. The symmetric matrix A1 is a covariance matrix or a precision matrix. Using the same method, a symmetric matrix A2 is obtained from the matrix X2. The symmetric matrix A2 is a covariance matrix or a precision matrix. 3) Weight calculation module: Calculates the squared Mahalanobis distance of the characteristic genes of any transcriptome data sample i among the n1 transcriptome data samples of the first biological sample. The calculation formula is: in It is a sample The standardized vector of expression levels of characteristic genes; the same calculation is performed on any transcriptome data sample j in the n2 samples of the second biological sample to obtain the square of its Mahalanobis distance; Further calculate the specific weights of the transcriptome data sample i. and the specific weights of the transcriptome data sample j The calculation formula is: ； The normalized specific weights of the transcriptome data sample i are further calculated. and the normalized specific weights of the transcriptome data sample j The calculation formula is: 4) Transcriptome distance calculation module: Calculates the weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j. i-j The weighted transcriptome distance EN-wTSD is calculated using the following formula: EN-wTSD in, For weighted Jensen-Shannon divergence; The weighted correlation distance is calculated based on the weighted Pearson correlation coefficient, where γ is the coefficient. in The calculation formula is: in and These are vectors representing the probabilistic gene expression levels of transcriptome data sample i and transcriptome data sample j, respectively. The function is , For the sample in the gene The expression probability distribution on M is and Constructed hybrid distribution: in The calculation formula is: The weighted Pearson correlation coefficient; 5) Output terminal, output the EN-wTSD.

6. A system for analyzing transcriptomic similarity between two transcriptomic data samples from a biological sample, said biological sample comprising at least three transcriptomic data samples, said system comprising: 1) Input terminal, used to receive transcriptome data from each transcriptome data sample, including gene name and gene expression level; 2) Matrix construction module: Constructs a matrix of the expression levels of characteristic genes in the transcriptome of each transcriptome data sample of the biological sample. Where n is the number of transcriptome data samples of the biological sample, p is the number of characteristic genes of the biological sample, and the matrix The elements in the table represent the expression levels of each characteristic gene in each transcriptome data sample of the biological sample. Take any one of the p genes in the matrix X. The vector of expression levels as response variables ,the remaining A matrix of gene expression levels was used as a predictor variable. A resilient network regression model is constructed, where the nodes of the resilient network are each gene. The regression coefficient vector of the nodes is minimized through the following objective function. The objective function is: in For mixed parameters, The same calculation is performed on the remaining p-1 genes except for gene k, using the penalty coefficient, to obtain their regression coefficient vectors respectively. A symmetric matrix A is constructed using all the regression coefficient vectors, where the symmetric matrix A is the covariance matrix or the precision matrix. 3) Weight calculation module: Calculates the squared Mahalanobis distance of the characteristic genes of any transcriptome data sample i among the n transcriptome data samples of the biological sample. The calculation formula is: in It is a vector of the expression levels of characteristic genes in transcriptome data sample i; the same calculation is performed on any transcriptome data sample j among the n transcriptome data samples of the biological sample to obtain the square of its Mahalanobis distance; Further calculate the specific weights of the transcriptome data sample i. and the specific weights of the transcriptome data sample j The calculation formula is: ； The normalized specific weights of the transcriptome data sample i are further calculated. and the normalized specific weights of the transcriptome data sample j The calculation formula is: ； 4) Transcriptome distance calculation module: Calculates the weighted transcriptome distance EN-wTSD between transcriptome data sample i and transcriptome data sample j. i-j The weighted transcriptome distance EN-wTSD is calculated using the following formula: EN-wTSD in, For weighted Jensen-Shannon divergence The weighted correlation distance is calculated based on the weighted Pearson correlation coefficient, where γ is the coefficient. in The calculation formula is: in and These are the probabilistic vector distributions of gene expression levels for transcriptome data sample i and transcriptome data sample j, respectively. The function is , For the sample in the gene The expression probability distribution on M is and Constructed hybrid distribution: in The calculation formula is: in The weighted Pearson correlation coefficient; 5) Output terminal, output the EN-wTSD.

7. A system for analyzing transcriptomic similarity between a first biological sample and a second biological sample, wherein the first biological sample and the second biological sample each comprise at least two transcriptomic data samples, the system comprising: 1) Input terminal, used to receive transcriptome data from each transcriptome data sample, including gene name and gene expression level; 2) A matrix construction module, as defined in claim 5; 3) Weight calculation module, as defined in claim 5; 4) Transcriptome distance calculation module, as defined in claim 5; 5) Similarity index calculation module: Calculates quantile thresholds based on EN-wTSD values. and similarity ratio The similarity ratio The quantile threshold as defined in claim 3 As defined in claim 3; 6) Output terminal: Output the EN-wTSD and quantile threshold values. and similarity ratio .

8. The system of any one of claims 5-7, wherein the transcriptome data samples are derived from organ-on-a-chip, organoids, assemblies, animal tissues, or ex vivo human tissues.

9. A computer device comprising at least one processor and a memory storing computer instructions executable on the processor, the instructions, when executed by the processor, implementing the steps of the method according to any one of claims 1-4.

10. A computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of the method according to any one of claims 1-4.

Images