A method, apparatus, and medium for transcriptome abnormality detection

Abstract

The present application relates to a kind of methods, equipment and medium for transcriptome anomaly detection, the method obtains gene expression, splicing event and isomer ratio data;Negative binomial, beta and normal distribution are fitted to each feature, the best distribution is selected based on MLE / AIC, and statistical anomaly score is calculated accordingly;Pseudo-normal samples are screened by sliding window, used to train WGAN-GP and obtain adversarial reconstruction anomaly score;Again, cross fitting reconstruction is carried out using attention autoencoder, and the robust Z score of gene and sample dimension is calculated based on pseudo-normal sample error, to obtain autoencoder reconstruction anomaly score;Three types of scores are weighted and fused into consensus score, converted into p value and corrected by FDR, and finally determine abnormal event.The method realizes the unified detection of three types of transcriptome anomalies through distribution guidance, pseudo-normal screening, adversarial reconstruction, attention cross fitting and weighted consensus.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics and transcriptome data analysis, specifically to a method, device, and medium for detecting transcriptome anomalies. Background Technology

[0002] Transcriptome anomaly detection is a key tool for analyzing the functional impact of rare genetic variations. Increasing evidence suggests that rare genetic variations often exert functional effects through aberrant transcriptome phenotypes (manifested as statistically significant outliers in RNA-seq), thus directly linking genetic variations to disease etiology. Recognized transcriptome anomalies are mainly classified into three categories: aberrant gene expression, aberrant splicing, and isoform ratio anomalies (i.e., differentially expressed transcripts).

[0003] Aberrant gene expression refers to extreme overexpression or underexpression after covariate adjustment, typically reflecting rare, detrimental variations that interfere with transcriptional regulation, RNA processing, or stability. In existing techniques, aberrant gene expression is often modeled using a negative binomial distribution, which captures the inherent mean-variance dependence and excessive dispersion of RNA-seq counting data. This assumption forms the basis of widely used RNA-seq analysis tools, including DESeq2 and frameworks specifically designed for outlier detection, such as OUTRIDER, which combines likelihood-based modeling with latent factor correction to robustly detect outliers in the presence of biological and technical confounding factors.

[0004] Aberrant splicing, including exon skipping, intron retention, and occult splice site activation, stems from splice site disruption and is a common mechanism in human disease. For aberrant splicing, the focus of analysis has shifted from total counts at the gene level to proportions at the junction site level. The beta-binomial distribution provides a natural framework that accommodates both binomial sampling of splicing events and hyperbinomial variations arising from biological heterogeneity. In existing techniques, tools such as FRASER apply this distribution through latent factor correction based on autoencoders, enabling sensitive detection of small but functionally significant splicing alterations.

[0005] Isomer ratio outliers capture variations in the relative abundance of different transcripts from the same gene, often leading to increases or losses of protein domains and functional alterations. Isomer-level analysis significantly improves sensitivity compared to gene-level pooling analysis. However, isomer-level quantification relies on probabilistically assigning reads to individual isomers, which often results in distributions that deviate from the canonical parameter family, exhibit multimodality, or uncertainty. Existing methods such as DRIMSeq employ Dirichlet-multinomial distributions to capture variability in isomer use, while integrated frameworks like IsoformSwitchAnalyzeR further link differential transcript use to functional consequences. However, the uncertainty in quantification and the heterogeneity of distributions make reliable isomer-level anomaly detection extremely challenging.

[0006] Beyond the statistical distribution-based methods mentioned above, machine learning approaches are being increasingly explored. Unsupervised architectures such as autoencoders and variational autoencoders learn latent transcriptomic representations and utilize reconstruction errors to identify outliers; generative adversarial networks simulate the distribution of normal variations and detect biases exceeding this learning baseline. Deep learning models trained on sequence context, including convolutional neural networks and Transformer-based architectures, extend detection to the areas of splicing and isoform regulation by directly linking genomic sequences to transcriptomic results.

[0007] However, existing technologies still have the following shortcomings: First, detection methods for the three types of anomalies—gene expression, splicing, and isoform ratio—remain fragmented. Each method is based on specific distribution assumptions, lacking a unified framework, resulting in poor cross-category comparability and hindering joint analysis. Second, while statistical models provide interpretability and a theoretical basis, they lack adaptability to heterogeneous distributions (especially the multimodal distribution of isoform ratios); while machine learning methods offer flexibility and scalability, they often come at the cost of sacrificing transparency and a fundamental distributional basis. Third, existing methods fail to effectively eliminate potential outlier interference in samples when constructing a reference normal distribution, leading to contamination of the reference baseline and affecting detection specificity. Fourth, for isoform ratio outliers, there is no mature and widely accepted statistical calibration method, resulting in low detection reliability of existing schemes.

[0008] Therefore, developing a distribution-aware detection method that can uniformly handle gene expression, splicing, and isoform ratio outliers, possess both statistical interpretability and machine learning adaptability, and establish a robust reference benchmark is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0009] This invention provides a method, device, and medium for transcriptome abnormality detection. Its purpose is to solve the technical problems in the prior art that it is impossible to uniformly detect three types of abnormalities: gene expression, splicing, and isoform ratio, and it is difficult to balance statistical interpretability and machine learning adaptability.

[0010] To achieve the above objectives, the first aspect of the present invention provides a method for detecting transcriptome abnormalities, comprising the following steps: Acquire transcriptome data, which includes gene expression data, splicing event data, and isoform ratio data; For each feature in the transcriptome data, a negative binomial distribution, a beta distribution, and a normal distribution are fitted respectively. Based on maximum likelihood estimation and the Akaike information criterion, the distribution with the highest goodness of fit is selected as the modeling distribution for that feature. If the feature is isoform ratio data and the goodness of fit is lower than a preset threshold, then no distribution-based statistical anomaly score is calculated. Calculate the statistical anomaly score for each feature-sample pair based on the modeled distribution; Pseudo-normal samples were filtered out from the expression values ​​of each gene using a sliding window; The pseudo-normal samples are used to train a Wasserstein generative adversarial network with gradient penalty, and the transcriptome data is input into the trained Wasserstein generative adversarial network with gradient penalty to obtain first reconstructed data. Based on the first reconstructed data and the first reconstructed error between the first reconstructed data and the transcriptome data, an adversarial reconstruction anomaly score is generated. The transcriptome data is input into the attention autoencoder for cross-fit reconstruction, and the second reconstruction error corresponding to the pseudo-normal sample is used as a standardization benchmark to obtain the second reconstruction data. Based on the second reconstruction error between the second reconstruction data and the transcriptome data, the robust Z-score for the gene and the robust Z-score for the sample are calculated respectively, and the abnormality score matrix is ​​combined to obtain the abnormality score reconstructed by the autoencoder. The consensus score is obtained by weighting and fusing the statistical anomaly score, the adversarial reconstruction anomaly score, and the autoencoder reconstruction anomaly score. The consensus score is converted into a p-value, and the error detection rate is controlled using a Benjamini-Hochberg process. Abnormal events are determined based on the corrected p-value.

[0011] Furthermore, after acquiring transcriptome data, methods for selecting the distribution with the highest goodness of fit as the modeling distribution for this feature based on maximum likelihood estimation and the Akaike Information Criterion include: For each feature, its observed data vector is denoted as... Fit the negative binomial distribution, beta distribution, and normal distribution respectively; The parameters of each distribution are estimated by maximizing the log-likelihood function, where maximizing the log-likelihood function is:

[0012] in, This represents the observed data vector of a certain feature across all samples. The total number of samples; For this feature in the first Measurements in a sample This is the probability density function or mass function of the corresponding distribution; The set of parameters representing the selected candidate distribution; This represents maximizing the log-likelihood function; Calculate the Akaike information content criterion values ​​for each distribution based on the estimated parameters:

[0013] in, This represents the Akaike Information Content Criterion value. This indicates the number of free parameters in the candidate distribution; These are the parameter estimates obtained by maximizing the log-likelihood function; Compare the Akaike information content criterion values ​​of each distribution, and select the distribution with the smallest Akaike information content criterion value as the modeling distribution for this feature.

[0014] Furthermore, methods for filtering pseudo-normal samples from the expression values ​​of each gene using a sliding window include: For each gene, its expression value in all samples is subjected to variance stabilization transformation, and then standardized to obtain the standardized expression value; The standardized expression values ​​are sorted in ascending order to obtain the sorted sequence and the corresponding sample order. For each position in the sorted sequence, construct a symmetrical window containing a fixed number of samples centered on that position; Within each symmetrical window, the local mean and local standard deviation are calculated, and the local bias score is calculated based on the standardized expression value of the sample at the center of the symmetrical window and the local mean and local standard deviation.

[0015] in, The local deviation score; It is a local mean; Local standard deviation; Standardized expression values; Samples with local deviation scores less than a preset threshold are labeled as pseudo-normal samples for that gene.

[0016] Furthermore, methods for obtaining the first reconstruction data include: The transcriptome data were subjected to logarithmic transformation and robust scaling normalization to obtain the preprocessed input matrix; The input matrix is ​​divided into multiple blocks according to features, and a generative adversarial network is trained independently for each block; For each block, a binary mask matrix is ​​constructed based on the pseudo-normal samples, where the position with a mask value of 1 indicates that the corresponding sample is a pseudo-normal sample. A generator is constructed, which adopts an encoder-decoder-encoder structure: the first encoder maps the input to the latent space, the decoder reconstructs the data from the latent vectors, and the second encoder encodes the reconstructed data again to obtain the latent vectors; Construct a discriminant whose output scalar score serves as the Wasserstein discriminant; The generator and the discriminator are trained using Wasserstein loss with gradient penalty, wherein the loss of the discriminator includes the score of the real sample, the score of the generated sample, and the gradient penalty term, and the loss of the generator includes adversarial loss, mask reconstruction loss, and latent space consistency loss. The trained generator is used to reconstruct the input data for each block, and the reconstruction results of each block are merged to obtain the first reconstructed data.

[0017] Furthermore, the method for generating an adversarial reconstruction anomaly score based on a first reconstruction error between the first reconstructed data and the transcriptome data includes: After inverse transformation back to the original scale, the absolute error between the transcriptome data and the first reconstructed data is calculated as the first reconstruction error; For each feature, the median and median absolute deviation are calculated using the first reconstruction error of the pseudo-normal sample corresponding to that feature, and then the robust Z-score of the feature dimension is calculated for each feature-sample pair. For each sample, the median and median absolute deviation are calculated using the first reconstruction error of all features within that sample, and then a robust Z-score for the sample dimension is calculated for each feature-sample pair. The robust Z-scores of the above feature dimensions and sample dimensions are geometrically averaged to obtain the adversarial reconstruction anomaly score.

[0018] Furthermore, the method for obtaining the second reconstructed data by inputting the transcriptome data into an attention autoencoder for cross-fitting reconstruction and using the second reconstruction error corresponding to the pseudo-normal sample as a standardization benchmark includes: The transcriptome data were subjected to logarithmic transformation and robust scaling normalization to obtain a preprocessed input matrix, which was then transposed into sample-feature form. An attention autoencoder is constructed, comprising: a learnable feature attention layer for weighting input features; an encoder for mapping the weighted input features to latent representations; and a decoder for reconstructing the original input from the latent representations. The attention autoencoder is trained in a denoising manner, that is, Gaussian noise is added to the input features. The training loss includes smooth L1 reconstruction loss, cosine similarity loss, latent representation regularization term and attention weight regularization term. K-fold cross-fitting is used for reconstruction, so that each sample is reconstructed by a model that has never seen the sample before. The reconstruction results are inversely transformed back to the logarithmically transformed scale to obtain the second reconstructed data.

[0019] Furthermore, based on the second reconstruction error between the second reconstructed data and the transcriptome data, robust Z-scores for genes and robust Z-scores for samples are calculated respectively, and combined to obtain an anomaly score matrix. This method for reconstructing anomaly scores using an autoencoder includes: The absolute error between the transcriptome data and the second reconstructed data is calculated as the second reconstruction error; For each feature, the median and median absolute deviation are calculated using the second reconstruction error of the pseudo-normal sample corresponding to that feature. Then, for each feature-sample pair, the robust Z-score of the feature dimension is calculated, retaining only the positive deviation. For each sample, the median and median absolute deviation are calculated using the second reconstruction error of all features within that sample. Then, for each feature-sample pair, a robust Z-score for the sample dimension is calculated, retaining only the positive bias. The robust Z-score of the feature dimension and the robust Z-score of the sample dimension are geometrically averaged to obtain the basic anomaly score. Calculate the average attention weight of each feature in the attention autoencoder, multiply the base anomaly score by a scaling factor composed of the average attention weight, and obtain the autoencoder reconstructed anomaly score matrix, which is used as the autoencoder reconstructed anomaly score.

[0020] Furthermore, the method for weightedly fusing the statistical anomaly score, the adversarial reconstruction anomaly score, and the autoencoder reconstruction anomaly score to obtain the consensus score includes: For each feature and each sample, the following three standardized anomaly scores are obtained: statistical anomaly score from distribution modeling, adversarial reconstruction anomaly score from generative adversarial networks, and autoencoder reconstruction anomaly score from attention autoencoders. The consensus score is calculated using the following weighted summation formula:

[0021] in, For consensus score, As a feature, As a sample, To count outlier scores, To combat the reconstruction of abnormal scores, To reconstruct anomalous scores for the autoencoder, , , The weight coefficients are non-negative and satisfy the following conditions: .

[0022] To achieve the above objectives, a second aspect of the present invention provides an electronic device including a memory and a processor, the memory being used to store a program supporting the processor in executing the method for transcriptome abnormality detection, the processor being configured to execute the program stored in the memory.

[0023] To achieve the above objectives, a third aspect of the present invention provides a computer-readable storage medium storing a computer program, which, when executed by a processor, performs the steps of the method for detecting transcriptome anomalies.

[0024] The beneficial effects of this invention are: Compared with existing technologies, this invention provides a method, device, and medium for transcriptome anomaly detection. It automatically selects the most suitable statistical distribution for each feature using maximum likelihood estimation and the Akaike information criterion, and marks isoform ratio features with low fit as "uncertain" to avoid forced distribution modeling, thus preserving statistical interpretability. Then, a sliding window density filter screens pseudo-normal samples for each gene, and these pseudo-normal samples are used to train a Wasserstein generative adversarial network with gradient penalty to obtain adversarial reconstruction anomaly scores. Simultaneously, an attention autoencoder is used for cross-fit reconstruction, and robust Z-scores are calculated using the error of pseudo-normal samples as a standardized benchmark, thus establishing a less contaminated and more reliable reference benchmark. Finally, the statistical anomaly score, adversarial reconstruction anomaly score, and autoencoder reconstruction anomaly score are weighted and fused to obtain a consensus score, which is then converted into a p-value and controlled by the Benjamini-Hochberg error detection rate to determine the anomaly event. This achieves unified detection of three types of anomalies across gene expression, splicing, and isoform ratio, balancing the principles of statistical models with the adaptability of machine learning models, significantly improving the robustness and accuracy of detection. Attached Figure Description

[0025] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below.

[0026] Figure 1 This is a flowchart of a method for detecting transcriptome abnormalities disclosed in an embodiment of the present invention.

[0027] Figure 2 This is a diagram showing the effect of sliding window filtering on the fidelity of GANs, as disclosed in an embodiment of the present invention. (a) shows that for gene expression features, the sliding window method achieves negative values ​​in all three distribution alignment indices compared to full sample training, indicating that the reference distribution quality is better; (b) shows that under the splice inclusion feature, the sliding window method outperforms the full sample method in all three distribution alignment indices; (c) shows that under the isoform ratio feature, the sliding window method outperforms the full sample method in all three distribution alignment indices.

[0028] Figure 3 This is an adversarial fitting plot under three modalities disclosed in an embodiment of the present invention, wherein (a) shows an adversarial fitting plot under gene expression; (b) shows an adversarial fitting plot under splicing inclusion; and (c) shows an adversarial fitting plot under isoform usage.

[0029] Figure 4 This is a Spearman correlation plot of sample-level outlier load in a GTEx whole blood cohort at the expression, transcript, and splice layer levels, as disclosed in an embodiment of the present invention. Detailed Implementation

[0030] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0031] According to embodiments of the present invention, it should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the following methods, in some cases the steps shown or described may be executed in a different order than that shown here.

[0032] Before detailing the implementation steps of this invention, it is necessary to explain the technical defects of the prior art and why this invention can overcome these defects. As mentioned earlier, gene expression, splicing events, and isoform ratios are three types of transcriptomic abnormalities that are closely related in biology. However, existing detection methods are isolated from each other, based on different distribution assumptions such as negative binomial, beta-binomial, or Dirichlet-multinomial, which makes it impossible to compare and integrate cross-category abnormal signals within the same framework. At the same time, although traditional statistical models are interpretable, they are difficult to flexibly adapt to non-normal, multimodal distributions such as isoform ratios. Pure machine learning methods (such as autoencoders and generative adversarial networks) can learn complex data manifolds, but they often lack distribution priors, leading to increased false positives or poor anomaly interpretability. In addition, when constructing a "normal" reference baseline, existing methods often use all samples directly, where potential outliers contaminate the reference distribution and reduce detection specificity. Especially for isoform ratio abnormalities, there is currently no mature and statistically calibrated detection scheme.

[0033] To avoid the aforementioned problems, this invention provides a method for transcriptome anomaly detection. This method automatically selects the optimal distribution for each feature or categorizes it as "uncertain" using maximum likelihood estimation and the Akaike information criterion, thus preserving statistical interpretability and adapting to heterogeneous distributions. A sliding window density filter screens pseudo-normal samples for each gene, and these are used to train a Wasserstein generative adversarial network and an attention autoencoder, significantly improving the purity of the reference benchmark and the calibration of reconstruction errors. Finally, by weighted fusion of three complementary anomaly scores and controlling the error detection rate, a unified, robust, and interpretable detection of three types of anomalies—gene expression, splicing, and isoform ratio—is achieved. The following detailed description of the method, with specific steps, is as follows. Figure 1 As shown, the method includes: Step S100: Obtain transcriptome data, which includes gene expression data, splicing event data, and isoform ratio data; Transcriptome data comprises three main categories: gene expression data, splicing event data, and isoform ratio data. These data typically originate from RNA sequencing (RNA-seq) experiments on biological samples (such as fibroblasts, skin tissue, and whole blood samples), and are obtained after data alignment, quantification, and standardization.

[0034] Specifically, gene expression data is provided in the form of a gene expression count matrix, where each row represents a gene, each column represents a sample, and the values ​​in the matrix represent the read count of that gene detected in that sample. To avoid noise introduced by low-expressed genes and to ensure the reliability of statistical modeling, the gene expression data needs to be filtered after acquisition.

[0035] In a preferred embodiment, the raw read counts are converted to counts per million (CPM) by dividing the raw read count of each gene for a given sample by the total number of valid reads in that sample and then multiplying by one million. CPM can eliminate the influence of sequencing depth differences between different samples. Subsequently, only genes with expression levels greater than 1 CPM in at least 5% of the samples are retained to exclude sparse or near-zero expression characteristics, allowing subsequent analysis to focus on stable expression genes suitable for negative binomial modeling.

[0036] Splicing event data are presented at the splice junction level, quantified using the percentage of splice inclusion (PSI) or an equivalent splice proportion (e.g., ψ5 or ψ3). PSI represents the proportion of a given exon included in a mature transcript, with values ​​ranging from 0 to 1. For each splicing event, the PSI value is calculated as the ratio of the number of split-reads spanning that junction to the total number of reads covering that region. Strict filtering is required when acquiring splicing event data to ensure quality. Specifically, only expressed junctions with sufficient splice read coverage in the target tissue (e.g., fibroblasts) are retained. Furthermore, the analysis is restricted to sites exhibiting variable splice site usage, characterized by multiple donors or recipients sharing a common splice site, thus focusing on events with biological variation. Junctions with PSI values ​​that remain nearly constant across samples (typically concentrated near the 0 or 1 boundary value) should be excluded to mitigate boundary effects and enhance the numerical stability of subsequent distribution-aware models.

[0037] Isoform ratio data reflects the relative abundance of different transcripts within the same gene. Transcript-level quantification of RNA-seq data was performed using tools such as RSEM to obtain the TPM value for each transcript. TPM, while eliminating the influence of sequencing depth and transcript length, provides a more intuitive representation of the relative expression ratio of transcripts. For each gene, the TPM value of each transcript is divided by the sum of the TPM values ​​of all transcripts of that gene to obtain the isoform utilization ratio, which ranges from 0 to 1. Furthermore, the sum of the isoform utilization ratios of all transcripts within the same gene is 1, constituting a constitutive constraint. To ensure the reliability of quantification in large-scale cohorts, after data acquisition, the isoform utilization ratio ψ value is limited to the range of 5% to 95%. This excludes isoforms with extremely low utilization (<5%) or near-complete dominance (>95%), as these isoforms often involve high quantitative uncertainty or a lack of information. This filtering strategy eliminates unstable or uninformative isoform configurations, ensuring that downstream analysis targets reliably quantified, biologically significant isoform-level variation.

[0038] It should be noted that the gene expression data, splicing event data, and isoform ratio data mentioned above can come from different analysis levels within the same cohort, or from different cohorts; this invention does not impose any restrictions on this. After completing data acquisition and necessary pre-filtering, the subsequent distribution modeling and anomaly detection steps can be performed.

[0039] Step S200: For each feature in the transcriptome data, fit a negative binomial distribution, a beta distribution, and a normal distribution respectively. Based on maximum likelihood estimation and the Akaike information criterion, select the distribution with the highest goodness of fit as the modeling distribution for that feature. If the feature is isoform ratio data and the goodness of fit is lower than a preset threshold, then no distribution-based statistical anomaly score calculation is performed. In step S200, distribution modeling is performed on each feature in the acquired transcriptome data (i.e., each gene expression feature, each splicing event feature, or each isoform ratio feature) to determine the most suitable statistical distribution type for that feature. Since the three types of outlier data have different inherent statistical characteristics—gene expression counts typically exhibit excessive dispersion with a mean-variance correlation, making them suitable for a negative binomial distribution; splicing event PSI values ​​are limited to the 0-1 range and exhibit super-binomial variation, making them suitable for a beta distribution; and while isoform ratios are also limited to the 0-1 range, they often exhibit multimodal, heterogeneous, or even uncertain distribution patterns, and may not conform to any parameter family—this invention employs a data-driven model selection strategy, rather than manual specification.

[0040] Specifically, for a certain feature, its observations across all samples are denoted as a vector. ,in Let be the total number of samples. The vector is fitted using three candidate distributions: the negative binomial distribution (suitable for counting non-negative integers, but can also be used approximated after continuous processing), the beta distribution (suitable for continuous proportions in the 0-1 interval), and the normal distribution (as a benchmark for general continuous features). For each candidate distribution, its parameters are solved using maximum likelihood estimation (MLE). The core idea of ​​maximum likelihood estimation is to find a set of parameters... This allows the current data to be observed under this parameter. The probability (i.e., the likelihood function) is maximized when the probability of a given sample point is maximized. The likelihood function is defined as the product of the probability densities (or masses) of each sample point, calculated by taking the logarithm:

[0041] in, For this feature in the first Measurements in a sample For the corresponding distribution, it is the probability density function (for continuous distributions) or the probability mass function (for discrete distributions). The set of parameters representing the selected candidate distribution; This represents maximizing the log-likelihood function; finding the log-likelihood function through optimization algorithms (such as gradient ascent or analytical solution) is crucial. Maximum parameter estimate .

[0042] After obtaining parameter estimates for each distribution, the Akaike Information Criterion (Akaike Information Criterion) is used. To compare the goodness of fit of the three candidate distributions. A penalty for model complexity is introduced based on the likelihood function, and its calculation formula is as follows:

[0043] in, This represents the Akaike Information Content Criterion value. This indicates the number of free parameters in the candidate distribution (for example, the negative binomial distribution usually has 2 parameters, the beta distribution has 2 parameters, and the normal distribution has 2 parameters). These are the parameter estimates obtained by maximizing the log-likelihood function; The smaller the value, the better the distribution achieves a balance between goodness of fit and model simplicity. Therefore, for each feature, the values ​​of the negative binomial distribution, beta distribution, and normal distribution are calculated separately. Value, select The distribution with the smallest value is used as the modeling distribution for this feature.

[0044] To further quantify the relative strength of evidence for model selection, Δ can also be calculated. And Akaike weights. Assume three candidate distributions... The values ​​are respectively Let the minimum value be... Then each distribution Δ for:

[0045] in, For the first The difference between the AIC value of each distribution and the minimum AIC value; For the first Akaike information criterion values ​​for each candidate distribution; It is the minimum value of AIC among all candidate distributions.

[0046] Akaike Weight Defined as:

[0047] in, For Akaike weights; For the first Distribution Value (defined as above); is the total number of candidate distributions.

[0048] The Akaike weights can be approximated as the posterior probability that the distribution represents the true distribution given the data. A smaller Δ... The large Akaike weights indicate that this distribution is strongly supported by the data.

[0049] Specifically, for isomer ratio data, since its distribution often deviates from the canonical parameter family, the goodness of fit of all three candidate distributions may be poor. For example, when the Akaike weights of all candidate distributions are below a certain preset threshold (e.g., 0.2) or the maximum log-likelihood value is extremely low, it indicates that the feature cannot be effectively described by any simple parametric distribution. In this case, the present invention marks the feature as "uncertain" and does not perform distribution-based statistical anomaly score calculation. This means that in the subsequent step S300, statistical anomaly scores are only calculated for features successfully assigned to negative binomial, beta, or normal distributions, while for features in the "uncertain" category (usually complex isomer ratio features), their statistical anomaly scores will be set to a default value or skipped directly, and anomaly detection will instead rely on subsequent machine learning components (adversarial reconstruction and autoencoder reconstruction). This design avoids false positives or false negatives caused by forcibly using incorrect distributions.

[0050] Through this step, each gene expression feature and splicing event feature is assigned the best-matching statistical distribution, while the isoform ratio feature is intelligently diverted to either statistical modeling or machine learning modeling paths based on its fit quality.

[0051] Step S300: Calculate the statistical anomaly score for each feature-sample pair based on the modeled distribution; After determining the optimal modeling distribution (negative binomial, beta, or normal distribution) for each feature and estimating the corresponding distribution parameters in step S200, this step further calculates the statistical anomaly score corresponding to each feature-sample pair. This score is used to quantify the degree of deviation between the observed value of that feature in a specific sample and its expected value under a given distribution.

[0052] Specifically, for features determined to follow a negative binomial distribution (such as gene expression data), whose observed values ​​are non-negative counts or normalized expression levels, the statistical anomaly score reflects the position of the observed value relative to a reference distribution defined by the mean and overdispersion parameter of the negative binomial distribution, for example, the degree to which it is located at the tail of the distribution. For features determined to follow a beta distribution (such as the PSI value of splicing events), whose observed values ​​are between 0 and 1, the statistical anomaly score, based on the density function or cumulative distribution function of the beta distribution, measures the degree to which the PSI value deviates from the typical proportion; the greater the deviation, the higher the score. For features determined to follow a normal distribution, the statistical anomaly score is the deviation value of the observed value after standardization by the mean and standard deviation. The calculation of all statistical anomaly scores relies on the distribution parameters obtained by maximum likelihood estimation in step S200, without using any additional cross-sample information or posterior adjustment, thus ensuring the statistical interpretability of the scores.

[0053] It is important to note that for features marked as "uncertain" in step S200 (common in complex isomer ratio data), since they cannot be effectively fitted by any candidate distribution, no statistical anomaly score is calculated for them in this step. In the subsequent weighted fusion in step S700, these features rely solely on the adversarial reconstruction anomaly score and the autoencoder reconstruction anomaly score. Through step S300, each feature-sample pair (limited to features that were successfully modeled) obtains a standardized statistical anomaly score. .

[0054] Step S400: Screen for pseudo-normal samples from the expression values ​​of each gene using a sliding window; In this step, for each gene (i.e., each gene expression characteristic), a set of "pseudo-normal" samples is selected from the expression values ​​of all its samples. These pseudo-normal samples are not pre-known normal controls, but rather samples automatically identified using an unsupervised, local density-based sliding window algorithm that are within the typical range of expression levels. They will be used in subsequent steps S500 and S600 to train the generative adversarial network and as a benchmark for normalizing reconstruction errors, thereby avoiding the contamination of the reference distribution by outliers that may be mixed in from the full sample.

[0055] Specifically, for each gene Suppose that it is in The original expression value vector in each sample is ,in, For genes In the sample The original expression value in Using the transpose symbol, the row vector is transposed to a column vector. First, a variance-stabilizing transformation is performed on these expression values. Because RNA-seq counting data exhibits a mean-variance relationship (smaller variance at low expression, larger variance at high expression), directly using the original values ​​would lead to the variance of highly expressed genes dominating subsequent calculations. One transformation is... The purpose of multiplying by 100 and adding 1 is to adjust the counting scale to an appropriate range while avoiding taking the logarithm of zero. The index of the sample indicates the first... A sample. This transformation makes the variance approximately constant throughout the dynamic range. The transformed value is denoted as... .

[0056] Next, the transformed expression values ​​of each gene are standardized to eliminate differences in baseline expression levels between genes. The standardization formula is:

[0057] in, For genes In the sample The standardized expression value of the Chinese character. For genes In the sample The original expression value in and These are the mean and standard deviation of the transformed expression value of the gene across all samples. A very small positive number (e.g.) (), used to prevent division by zero. Standardized With zero mean and unit variance, the expression levels of different genes can be compared relatively fairly.

[0058] Obtain the standardized expression value vector Then, sort them in ascending order, and record the original sample index corresponding to each value. Let the sorted sequence be... ,in It is the sorted sample index mapping function. For genes After sorting, the first The standardized representation value of each position. .

[0059] For each position in the sorted sequence (Corresponding rank is) Construct a symmetrical window centered on this location. The window contains a fixed number of samples, denoted as the window size. In a preferred embodiment, take A symmetrical window means that the center of the window is the location. The window extends to the left and right sides by approximately We have a limited number of samples, but we need to handle boundary cases: the window may be asymmetrical near the ends of the sequence; in this case, we should use only the actually usable samples. Within each window... Within the calculation window, the local mean and local standard deviation of all standardized expression values ​​are calculated.

[0060] Subsequently, the window center sample (i.e., the rank of) is calculated. The local bias score of a sample is defined as the absolute value of the difference between the standardized expression value of that sample and the local mean of the window, divided by the local standard deviation.

[0061] in, The local deviation score; It is a local mean; Local standard deviation; This is the standardized expression value; the essence of this local bias score is the degree of standardized deviation of the central sample relative to its neighboring samples. If the central sample expresses similarly to its surrounding samples, then... Smaller; if the central sample is a local outlier (e.g., excessively high or low expression), then Relatively large.

[0062] Set a threshold The local deviation score is less than The sample is labeled as a pseudo-normal sample of this gene. In a preferred embodiment, take... This is because in a standard normal distribution, approximately 95% of the data points fall within the range of the mean ± 2 standard deviations. This indicates that the expression value of this sample is not significantly different from that of its neighboring samples and can be considered "normal". For each gene... The set of indices of its pseudo-normal samples is denoted as :

[0063] It is worth noting that the screening of pseudo-normal samples is performed independently on a gene-by-gene basis. That is, each gene may select a different set of samples as its own pseudo-normal samples. This is because different genes have different expression patterns, and the abnormally high expression of certain genes in a specific sample may not affect the reference distribution of other genes. This gene-specific screening of pseudo-normal samples is more refined and reasonable than global sample filtering.

[0064] From a computational complexity perspective, sorting each gene requires... The time complexity is O(n), while the calculation of window statistics can be completed in a single linear scan, with a complexity of O(n log n). Because the number of genes is usually very large (tens of thousands), but the number of samples is small. Typically small (tens to hundreds), the algorithm is efficient and feasible in practice.

[0065] In step S400, a set of highly purified pseudo-normal samples is obtained for each gene. These samples will be used as training references for the Wasserstein generative adversarial network in step S500 (for constructing the mask matrix) and as a benchmark for normalizing the reconstruction error of the attention autoencoder in step S600 (for calculating the median and median absolute bias). Compared to directly using all samples, this sliding window-based local density filtering strategy can effectively eliminate potential outliers from contaminating the reference distribution, thereby improving the specificity and robustness of subsequent anomaly detection.

[0066] Step S500: Train a Wasserstein generative adversarial network with gradient penalty using the pseudo-normal samples, and input the transcriptome data into the trained Wasserstein generative adversarial network with gradient penalty to obtain first reconstructed data, and generate an adversarial reconstruction anomaly score based on the first reconstruction error between the first reconstructed data and the transcriptome data. A Wasserstein generative adversarial network (WGAN-GP) with gradient penalty was trained using the obtained pseudo-normal samples. This network was then used to reconstruct transcriptome data, and an adversarial reconstruction anomaly score was calculated based on the reconstruction error. This step aims to learn a generative model capable of capturing the dominant structure of the transcriptome data (i.e., normal variation patterns), aligning the reconstructed samples with the real samples in terms of distribution, thereby identifying sample-feature pairs that cannot be well reconstructed as anomalies. Unlike traditional generative adversarial networks, WGAN-GP significantly improves training stability and generation quality through Wasserstein distance and gradient penalty mechanisms.

[0067] Specifically, the input transcriptome data matrix is ​​first preprocessed. Let the original data matrix be... ,in For the sample size, These are the feature numbers (including gene expression, splicing events, and isoform ratio features). To reduce data skewness and stabilize variance, a logarithmic transformation is applied to each element:

[0068] in, For the first The sample, the first The values ​​of each feature after logarithmic transformation Adding 1 to the original expression level (such as CPM, PSI, or isoform ratio) is to avoid taking the logarithm of zero. Next, RobustScaler is used to standardize each feature (i.e., each column). RobustScaler scales based on quantiles (typically using the median and interquartile range), which is more robust to outliers than standard normalization (based on the mean and standard deviation).

[0069] Since directly inputting high-dimensional transcriptome data (with tens of thousands of features) into a single generative adversarial network would lead to excessive computational burden and the risk of overfitting, this invention divides the features into several chunks of similar size. Each chunk contains a set of features, and each chunk trains its own generative adversarial network independently. This chunking strategy reduces memory and computational requirements while allowing each chunk to learn the joint distribution of its internal features.

[0070] For each block, the pseudo-normal sample index is selected according to the criteria in step S400 for each gene (feature). Construct a binary mask matrix ,in The number of features within this block. The number of samples. Elements of the mask matrix. If and only if the sample In features The above were marked as pseudo-normal samples (i.e. , Features Corresponding genes The set of pseudo-normal sample indices (from step S400)); otherwise The purpose of this mask matrix is ​​to tell the generator that during training, only the locations with a mask value of 1 need to be reconstructed, while the locations with a mask value of 0 (which may contain outliers) are partially ignored, thereby preventing the model from learning anomalous patterns.

[0071] The generator adopts a three-stage structure of "encoder-decoder-encoder". Let the input vector be... For the sample Feature values ​​within the current block (after preprocessing). First encoder. Will Mapping to a low-dimensional latent space yields latent vectors. Decoder Will Reconstructed This refers to the reconstructed sample output by the generator. The second encoder... Reconstruct the sample Encode again to obtain The entire generator can be formally represented as:

[0072]

[0073]

[0074] in, For latent vectors, For the first encoder, For the input vector, This represents all trainable parameters of the generator; To reconstruct the sample; For decoders; This is the latent vector obtained by the second encoder re-encoding the reconstructed sample. This is the second encoder. The second encoder is introduced to impose latent space consistency constraints, i.e., requiring... and The closer the input is to the source, the more the generator is encouraged to retain the core semantic information of the input.

[0075] Discriminator As a Wasserstein discriminator, it outputs a scalar score instead of the probability value found in traditional GANs. This score can be interpreted as a continuous measure of the "realism" of the input sample: the higher the score for a real sample, the lower the score for a generated sample. The discriminator also uses a multilayer perceptron structure, but it does not require the use of a sigmoid activation function in the output layer.

[0076] The training process alternately optimizes the discriminator and generator. The objective of the discriminator is to maximize the Wasserstein distance estimate between real and generated samples. Its loss function is:

[0077] in, The loss function; The number of samples; For discriminator networks; This is the gradient penalty coefficient; The gradient penalty term consists of three parts: the first term encourages the discriminator to give high scores to real samples, the second term encourages it to give low scores to generated samples, and the third term is the gradient penalty term, used to enforce the Lipschitz constraint (the theoretical requirement of Wasserstein GAN). The formula for calculating the gradient penalty is:

[0078] here These are the linear interpolation points between real samples and generated samples. From uniform distribution Random sampling in the middle; For interpolation points The gradient operator. The gradient penalty requires the gradient norm of the discriminator at the interpolation point to be as close to 1 as possible, thus implicitly satisfying the 1-Lipschitz constraint. In a preferred embodiment, the gradient penalty coefficient is... .

[0079] The generator's loss function consists of a weighted sum of three parts:

[0080] in, The loss function for the generator; To counteract the loss of weight; Loss weights are used to reconstruct the mask; Weights for latent space consistency loss; To combat the losses; Loss due to mask reconstruction; This represents the consistency loss in the latent space.

[0081] Combating losses This encourages the generation of samples that can deceive the discriminator, i.e., make the discriminator give the generated samples a high score. Defined as:

[0082] in, For the discriminator to generate samples The output score.

[0083] Mask Reconstruction Loss The generator is forced to accurately reconstruct the original input at the locations of pseudo-normal samples. This loss only calculates the mask matrix. The positions marked as 1 are ignored as they may be abnormal. Defined as:

[0084] in, It is a small constant that prevents division by zero. The nth in the mask matrix The sample, the first The value of each feature; For the generator to the first The sample, the first Reconstructed values ​​of each feature; The original input value; This represents the number of features in the current block. The denominator is the sum of all mask elements, which serves as an average.

[0085] Latent space consistency loss : Requires the latent coding of the original sample The re-encoding of the reconstructed samples They are close in Euclidean distance, thus maintaining the stability of the latent representation. Defined as:

[0086] In a preferred embodiment, the loss weight is set as follows: , , A higher reconstruction loss weight emphasizes accurate reconstruction of pseudo-normal samples.

[0087] During training, the Adam optimizer is used, and the learning rate is set to... In each iteration step, the discriminator is updated multiple times (e.g., 5 times), and then the generator is updated once to maintain the discriminator's discriminative ability. The total number of training epochs is set to 300, and the training can be stopped early based on the smoothness of the loss curve.

[0088] After training, the generator is used to reconstruct all samples (including those not used for training pseudo-normal samples) to obtain the first reconstructed data. The reconstruction results are inversely transformed back to the scale after logarithmic transformation (i.e., (space), then calculate the absolute reconstruction error:

[0089] in, In order to combat the reconstruction module, features ,sample The absolute reconstruction error; These are the original values ​​after preprocessing (after logarithmic transformation); The value is to counteract the reconstruction of the generative network.

[0090] In order to transform reconstruction errors into statistically interpretable outlier scores, this invention employs a two-dimensional Z-score calibration method based on robust statistics.

[0091] First, for each feature The pseudo-normal sample index set obtained in step S400 is used. Extract the corresponding values ​​of these samples in the first reconstruction error, and denote them as the set. Calculate the median of this set. and median absolute deviation Then, for each feature-sample pair Calculate the robust Z-score for the feature dimension:

[0092] Among them, factor It is a proportionality constant that calibrates MAD to be consistent with the standard deviation (assuming the error follows a normal distribution). It is a small constant that prevents division by zero; This is the robust Z-score for the feature dimension. The Z-score represents the degree of deviation of the reconstruction error of the current sample from the error of normal samples under this feature; the larger the positive value, the more significant the anomaly.

[0093] Secondly, for each sample Considering the reconstruction error of all features within the sample, i.e., the set .calculate the median of and median absolute deviation Then calculate the robust Z-score for the sample dimension:

[0094] in, The robust Z-score for the sample dimension.

[0095] This score measures the performance in the sample. Internal features The degree of deviation of the reconstruction error from the overall error level of the sample.

[0096] To combine the information from the two dimensions mentioned above, this invention performs a geometric average of the two and adds a very small constant. To ensure numerical stability, the basic adversarial reconstruction anomaly score is obtained:

[0097] in, Based on the adversarial reconstruction of abnormal scores, take This is to retain only positive bias (excessive reconstruction error indicates anomaly), while negative bias (reconstruction error below normal level) does not have biological significance as an indication of anomaly, so it is set to zero.

[0098] Through step S500, the present invention obtains an adversarial reconstruction anomaly score matrix based on WGAN-GP, which captures the degree of deviation between sample-feature pairs and the learned normal data manifold. Due to the guidance of pseudo-normal samples and the constraints of mask reconstruction, the score has high specificity and robustness to potential outliers.

[0099] Step S600: Input the transcriptome data into the attention autoencoder for cross-fit reconstruction, and use the second reconstruction error corresponding to the pseudo-normal sample as a standardization benchmark to obtain the second reconstruction data. Based on the second reconstruction error between the second reconstruction data and the transcriptome data, calculate the robust Z-score for the gene and the robust Z-score for the sample respectively, and combine them to obtain the anomaly score matrix as the autoencoder reconstruction anomaly score. In this step, a lightweight attention-based autoencoder is used to reconstruct the transcriptome data, and a cross-fitting strategy is employed to reduce the model's optimistic bias towards training samples. Unlike the generative adversarial network in step S500, this autoencoder focuses on capturing the statistical properties of reconstruction errors from both the feature and sample dimensions, thereby generating a complementary anomaly score matrix. An autoencoder is an unsupervised neural network architecture where the encoder compresses the input data into a low-dimensional latent representation, and the decoder recovers the original input from the latent representation. The training objective is to make the reconstructed output as close as possible to the original input. The attention mechanism allows the model to assign different weights to different features, enabling the autoencoder to focus on features more important to the reconstruction task.

[0100] Specifically, the transcriptome data matrix is ​​first preprocessed. Let the original data matrix be... ,in For the sample size, The matrix contains the number of features (including gene expression, splicing events, and isoform ratio features). Missing values ​​are filled with zeros, and the matrix is ​​then transposed into sample-feature form, where each row corresponds to a sample and each column corresponds to a feature. A logarithmic transformation is then applied to each element. To stabilize variance and reduce data skewness, the transformed data is standardized by feature using RobustScaler, which scales based on the median and interquartile range, providing good robustness to outliers.

[0101] The structure of an attention autoencoder is as follows: For the input sample vector (i.e., the first) (The values ​​of each sample across all features). The model first applies a learnable feature attention mechanism. Let the original attention parameters be... Through the Sigmoid function Map it to The interval is used to obtain the attention weight vector. Then, the input vector is multiplied element-wise by the attention weights to obtain the weighted input. : ,in This represents element-wise multiplication. For the first The input feature vectors of each sample. This attention mechanism enables the model to automatically learn which features are more critical for the reconstruction task and assign them higher weights.

[0102] The weighted input is fed into the encoder. Mapping to a low-dimensional latent space yields a latent representation. Subsequently, the decoder Reconstructing the latent representation back to the original feature space: Encoders and decoders typically employ a multilayer perceptron (MLP) structure with a small number of layers to ensure the model's lightweight nature.

[0103] To improve the model's robustness to noise and outliers, this step employs a denoising training strategy. Specifically, Gaussian noise is added to the input data during training, forcing the autoencoder to learn feature representations insensitive to input perturbations, thereby enhancing the stability of the reconstruction. The training loss function consists of four weighted components: SmoothL1 Reconstruction Loss: Measures the reconstruction output With the original input The difference between them. SmoothL1 loss (also known as Huber loss) behaves like L2 loss when the error is small and like L1 loss when the error is large. It is more robust to outliers than mean squared error (MSE).

[0104] Cosine similarity loss: defined as ,in Calculate the cosine similarity between two vectors. This loss encourages the reconstructed vector to maintain the same orientation as the original vector, which helps preserve the relative proportions of features.

[0105] Latent representation regularization term: Calculate latent representations The mean of the absolute values ​​(i.e., L1 regularization), multiplied by a coefficient of 0.01, is used to constrain the complexity of the latent space and prevent overfitting.

[0106] Attention weight regularization term: Calculates the attention weight vector mean and The differences between features are multiplied by a coefficient of 0.01 to encourage a more dispersed distribution of attention weights and avoid over-concentration on a few features.

[0107] To reduce the optimism bias generated during the self-reconstruction process (i.e., the model's reconstruction error for training samples is systematically low because the model has already memorized these samples), this step employs a K-fold cross-fit reconstruction strategy. The specific steps are as follows: All samples are randomly divided into... A subset of samples of similar size. For each subset, use all samples from all other subsets (i.e., the remaining samples). Train an attention autoencoder on a subset of samples, then use that model to reconstruct samples in the current subset. Repeat this process. This process ensures that each sample is reconstructed by a model that has never seen that sample during training. Cross-fit reconstruction effectively eliminates the model's overfitting bias to its own training samples, making the reconstruction error more reflective of the sample's "unexpectedness" rather than a memory effect. In a preferred embodiment, taking... or 10.

[0108] After completing the cross-fit reconstruction, the reconstruction results are inversely transformed back to the scale after logarithmic transformation (i.e., (space), to obtain the second reconstruction data For each sample and characteristics Calculate the absolute reconstruction error:

[0109] in, For absolute reconstruction error, These are the original values ​​after preprocessing. These reconstruction errors constitute a... The error matrix.

[0110] Next, the pseudo-normal sample index set obtained by screening for each feature (gene) in step S400 is used. These error values ​​are then subjected to robust statistical calibration, and robust Z-scores for the feature dimension and the sample dimension are calculated respectively.

[0111] Robust Z-score for each feature dimension: The reconstruction errors of all pseudo-normal samples under this feature are collected and denoted as set. Calculate the median of this set. And the median absolute deviation (MAD). MAD is a very robust measure of dispersion and is insensitive to outliers. Then, for each feature-sample pair... Calculate the robust Z-score for the feature dimension:

[0112] in, Robust Z-score for the feature dimension; The set of reconstruction errors for pseudo-normal samples; coefficients It is a proportionality constant that calibrates MAD to be consistent with the standard deviation (assuming the error follows a normal distribution). It is a very small positive number used to prevent division by zero. The Z-score reflects the position of the reconstruction error of the current sample relative to the error distribution of pseudo-normal samples under this feature: the larger the positive value, the more significantly the reconstruction error of the sample on this feature is higher than the normal level, that is, the higher the degree of abnormality; a negative value indicates that the reconstruction error is lower than the normal level, which usually does not have a biologically significant abnormality indication.

[0113] Robust Z-score for each sample dimension: The reconstruction errors of all features within the sample are collected and denoted as the set. .calculate the median of and median absolute deviation Then calculate the robust Z-score for the sample dimension:

[0114] in, For the sample dimension, a robust Z-score is used; this Z-score measures the robustness of the sample dimension. Internal features The Z-score is the degree of deviation of the reconstruction error from the overall error level of the sample. If the reconstruction error of a certain feature in the sample is much larger than the reconstruction errors of most features in the sample, the Z-score of the sample dimension will be larger.

[0115] In order to integrate the information from the two dimensions mentioned above into a unified anomaly score, this step uses a geometric mean to combine them and only retains positive biases (because only excessively high reconstruction errors indicate anomalies).

[0116] Specifically, first, take the positive parts of the two Z-scores respectively: and Then, the geometric mean of the two is calculated, and a very small constant is added to ensure numerical stability:

[0117] in, Reconstructing anomaly scores based on the autoencoder (without attention scaling); To retain only positive deviations, negative values ​​are set to zero.

[0118] The geometric mean can penalize errors that are insignificant in both dimensions simultaneously, and the baseline outlier score will be high only when the Z-scores in both dimensions are positive and large. This is more selective than simple summation or averaging.

[0119] Finally, attention weights are introduced as a gentle scaling factor. For each feature... Calculate its average attention weight across all samples. The average attention weight is scaled to the range of 0.5 to 1.0, and then multiplied by the base anomaly score to obtain the final autoencoder reconstruction anomaly score:

[0120] in, Reconstruct the anomalous score for the final autoencoder; This represents the average attention weight.

[0121] The scaling factor serves to: optimize the scaling for high-attention features automatically learned by the model (i.e., ... For features with larger anomalies, their outlier scores will be moderately amplified because these features are more important for the reconstruction task, and their reconstruction errors are more noteworthy; while for features with low attention, their outlier scores will be compressed to reduce their interference with the final result. All sample-feature pairs... Together they form the autoencoder reconstruction anomaly score matrix, which has the same dimension as the statistical anomaly score in step S300 and the adversarial reconstruction anomaly score in step S500.

[0122] Step S700: The statistical anomaly score, the adversarial reconstruction anomaly score, and the autoencoder reconstruction anomaly score are weighted and fused to obtain a consensus score; In this step, anomaly scores from three complementary detection modules are integrated into a unified consensus score. These three modules characterize transcriptomic anomalies from different perspectives: the statistical anomaly score provided in step S300 (denoted as...). Based on probability distribution assumptions (negative binomial, beta, or normal distribution), it reflects the degree of deviation of gene expression or splicing characteristics from their theoretical distribution; the adversarial reconstruction anomaly score provided in step S500 (denoted as...) Based on a Wasserstein generative adversarial network, the deviation between sample-feature pairs and the learned normal data manifold is measured; step S600 provides an autoencoder to reconstruct anomaly scores (denoted as...). Based on attention autoencoders and cross-fit reconstruction, reconstruction errors are calibrated from both the feature and sample dimensions. These three score matrices have the same dimension and are all non-negative (or positively truncated); a higher score indicates a higher degree of anomalousness. The consensus score is a comprehensive score obtained by combining multiple heterogeneous anomaly scores through a weighted summation. Its core idea is that only events jointly identified as anomalous by multiple independent detectors are more likely to be genuine biological anomalies, thereby reducing false positives or false negatives caused by a single model.

[0123] Specifically, for each feature and each sample Standardized anomaly scores are obtained from the three modules respectively. , and Among them, outlier scores were statistically analyzed. Only features successfully assigned to negative binomial, beta, or normal distributions are defined; for features marked as "uncertain" in step S200 (common in complex heterogeneous proportion data), the score is set to a default value (e.g., 0), and subsequent fusion relies only on the other two scores. Adversarial reconstruction anomaly scores and autoencoder reconstruction anomaly scores are defined for all features.

[0124] The consensus score is calculated using a weighted linear combination:

[0125] in, For consensus score, As a feature, As a sample, To count outlier scores, To combat the reconstruction of abnormal scores, To reconstruct anomalous scores for the autoencoder, , , The weight coefficients are non-negative and satisfy the following conditions: The weights are chosen to balance statistical calibration, resistance to manifold sensitivity, and robust reconstruction evidence, ensuring that no single component dominates the final decision.

[0126] In a preferred implementation, adaptive weighting can be applied based on the independent performance of each module on the validation set (e.g., AUPRC), or an equal-weight strategy can be used as the default setting. For the heterogeneity ratio feature, due to the lack of statistical anomaly scores, only the adversarial reconstruction score and the autoencoder reconstruction score are used in actual fusion, at which point the weights can be renormalized. The weighted fusion operation makes the three detectors complementary: the statistical anomaly score provides the theoretical basis for distribution awareness, the adversarial reconstruction score is sensitive to the shift in the overall data manifold, and the autoencoder reconstruction score provides robust local calibration in both sample and feature dimensions. By requiring consensus among independent and complementary detectors, the consensus mechanism effectively reduces the bias of individual models, improves robustness against noise, distribution misjudgments, and reconstruction artifacts, and provides a stable and interpretable comprehensive anomaly index for subsequent statistical inference. After calculation, each feature-sample pair Each obtained a consensus score This score will be used as input for significance evaluation in step S800.

[0127] Step S800: Convert the consensus score into a p-value, and use the Benjamini-Hochberg process to control the error detection rate, and determine the abnormal event based on the corrected p-value.

[0128] In this step, the consensus score obtained in step S700 is... The results are converted into statistically significant p-values, and the final anomalous events are identified through false discovery rate (FDR) control. The p-value refers to the probability of observing the current or more extreme results under the assumption that the null hypothesis is true; in this invention, the null hypothesis is set as "the feature-sample pair does not contain anomalies," meaning the consensus score only reflects normal random fluctuations. The false discovery rate (FDR) is the expected proportion of false rejections (false positives) among rejected null hypotheses (i.e., events judged as anomalous). The Benjamini-Hochberg process (BH process for short) is an FDR control method widely used in high-throughput biological data. It detects as many true anomalies as possible while controlling the false positive rate by sorting p-values ​​and setting adaptive thresholds.

[0129] Specifically, firstly, the consensus scores of all feature-sample pairs calculated in step S700 are... Convert to p-value. Since the consensus score is a standardized comprehensive anomaly indicator, a larger value indicates a higher degree of anomaly.

[0130] Through the above steps, this invention achieves statistical calibration and reliable identification of anomalous events at three levels: expression, splicing, and isomer ratio. This avoids a large number of false positives caused by multiple hypothesis testing and ensures the statistical reliability of the final anomaly detection results.

[0131] After completing the anomaly detection process in steps S100 to S800 above, this embodiment further conducts a systematic performance evaluation and experimental verification of the proposed method on multiple independent RNA-seq cohorts to demonstrate its effectiveness and superiority in detecting gene expression, splicing, and isoform ratio anomalies. The datasets used for evaluation include: (1) 119 fibroblast samples from the Kremer rare disease dataset as a rare disease benchmark, which contains previously reported anomalous expression and splicing events; (2) 606 sun-exposed skin samples from GTEx v8 as a large-scale health reference for benchmark variant modeling of gene expression; and (3) whole blood samples from GTEx v8 for assessing compositional variation at the isoform level in a unique tissue context. By applying the multi-stage process of this invention, including distribution-guided branching, adversarial reconstruction, autoencoder-based anomaly scoring, and consensus integration, anomalous events at the gene, splicing, and isoform levels are systematically quantified.

[0132] To systematically evaluate the impact of sliding window filtering (step S400), controlled anomaly injection experiments were constructed on the GTEx v8 skin dataset. Experiments were conducted at four effect sizes (…). Artificial outliers were introduced, and the sliding window size was varied for each setting (range: 100, 200, 300, 400, 500 samples). Two complementary evaluation criteria were used: (i) the gene-level average Jaccard index, which quantifies the average overlap between the predicted and actual outlier sets across genes; and (ii) the gene-level median Jaccard index, which is more robust to gene-specific variations. The Jaccard index is defined as the ratio of the intersection size to the union size of two sets, and is used to measure the similarity between the predicted and actual outlier sets. The evaluation results are shown in Table 1. Table 1. Sliding window filtering performance under different anomaly intensities (Z) and window sizes.

[0133]

[0134] At all Z thresholds, the Jaccard exponent remained consistently high (both the mean and median were greater than 97%), demonstrating the robustness of the proposed sliding window method to different window sizes. Under stronger perturbations (Z=4, 6), the Jaccard value exceeded 99%, indicating near-perfect recovery of manually inserted outliers. This demonstrates that the sliding window procedure effectively preserves genuine outlier signals while reducing noise introduced by extreme samples.

[0135] In addition to Jaccard-based evaluation, three distributional dissimilarity metrics were used to compare the generative adversarial networks trained with a "sliding window with reserved samples" approach with those trained with a "full sample" approach: one-dimensional Fréchet distance, center (10% truncation) Wasserstein distance, and two-sample Kolmogorov-Smirnov p-value. The cross-gene mean was consistently negative in terms of expression, splicing ratio, and isoform usage, indicating that the mean-variance structure and center alignment were more consistent when using a sliding window. These results confirm that sliding window pre-filtering significantly improves the quality of the reference distribution for adversarial generation, thus achieving more stable and reliable anomaly detection.

[0136] like Figure 2 As shown, for each feature type—gene expression, splicing inclusion (…) The ratio of isoforms—a one-dimensional Wasserstein generative adversarial network (WGAN) was trained on a sliding window pseudo-normal subset. Numerical values ​​in... Variance stabilization is performed across scales, and standardization is applied based on features. Both the generator and discriminator employ multilayer perceptrons (MLPs) with LeakyReLU activation functions. Batch normalization and dropout are used in the generator to stabilize training and mitigate overfitting. The optimization process follows the RMSProp algorithm, with each generator step corresponding to five discriminator steps. Lipschitz constraints are implemented through weight clipping, a fixed random seed is set, and early stopping is implemented when the loss function reaches a plateau.

[0137] Complementary distribution diagnostic metrics are used to quantify the fidelity of the generated data: (i) the two-sample Kolmogorov-Smirnov (KS) statistic (sensitive to global distribution shifts); (ii) the 10% truncated Wasserstein (Earth-Mover) distance (robust to extreme tails); and (iii) the integrated absolute error (IAE) between kernel density estimates (KDEs) (sensitive to differences in probability distribution shape). Furthermore, the stability of moment alignment and generator / discriminator loss is monitored to prevent mode collapse.

[0138] Sliding window training consistently improved GAN fidelity across all three modalities (gene expression, splicing ratio, and isoform usage). Negative values ​​indicate that sliding window training makes the "pseudo-normal" distribution closer to the real data in terms of mean-variance matching and overall alignment; while small or near-zero values ​​reflect stable discriminability. (The image shows a superimposed plot of the true and generated densities of representative features from the three modalities.) Figure 3 As can be seen, the reference distribution generated adversarially closely replicates the empirical distribution in terms of location, scale, and shape, including the shoulder and tail behaviors used by the isomers. Among the representative features, the KS test generally fails to reject the equality hypothesis at the 5% significance level, while the truncated Wasserstein and IAE values ​​remain at low levels, indicating minimal distributional differences. These observations support using generative adversarial network references as high-fidelity priors in the masked autoencoder stage, thereby reducing spurious residuals and improving the specificity of anomaly detection.

[0139] In the Kremer fibroblast cohort ( The detection performance of gene expression outliers was evaluated using controlled artificially synthesized perturbations. Based on a filtered counting matrix, outliers were injected into 1% of all gene-sample entries. For each gene, the injected value was... The data is generated spatially based on the gene-specific mean and standard deviation, and then transformed back to the original counting scale. Four perturbation amplitudes were considered ( This covers biases from weak to strong. OutSingle, OUTRIDER-PCA, OUTRIDER-Base, and the method of this invention (MAE-GUIDE) are compared using precision-recall (PR) curves and receiver operating characteristic (ROC) curves, combined with AUPRC (area under the precision-recall curve), AUROC (area under the ROC curve), the best F1 score, and the precision at the best F1 threshold to summarize performance.

[0140] The results show that performance steadily improves with increasing perturbation strength. At Z=2, the methods are largely equivalent in AUPRC, with the method of this invention achieving the best AUROC (0.9599). Starting from Z=3, the method of this invention becomes the best performing method among the main summary metrics, with AUPRC / AUROC of 0.8658 / 0.9956 at Z=3, 0.9477 / 0.9984 at Z=4, and 0.9741 / 0.9991 at Z=6, achieving the highest overall best F1 score under moderate and strong perturbation levels. These results indicate that once the bias becomes moderately separable, the method of this invention provides the most reliable overall detection performance while remaining competitive under weak perturbations.

[0141] Further comparison of the outlier sets identified by the method of this invention and the OUTRIDER reveals that outliers detected solely by the method of this invention tend to appear more frequently at relatively high average expression levels, while outliers detected solely by the OUTRIDER are more concentrated in the low to moderate abundance range, with common outliers occupying the middle distribution. This indicates that the two methods do not simply repeat the same expression anomalies, but rather exhibit different detection preferences across the entire expression profile. The method of this invention complements the OUTRIDER by capturing an additional subset of anomalous events in more robustly expressed genes.

[0142] Anomaly detection at the joinotome level was evaluated in a Kremer fibroblast cohort using a controlled PSI (splice inclusion percentage) perturbation scheme. Anomalies at the joinotome level were introduced by modifying the PSI value while preserving local splice competition around the same donor or recipient site to maintain empirical joinotome usage structure. Three perturbation mechanisms were considered: weak (0.2–0.3), moderate (0.3–0.7), and strong (0.7–1.0). The method of this invention was compared with FRASER-AE and FRASER-PCA.

[0143] Under all perturbation mechanisms, the method of this invention achieved the strongest overall performance. Under weak perturbations, the method of this invention achieved an AUPRC of 0.7208 and an AUROC of 0.9701, outperforming FRASER-AE (0.7007 / 0.9676) and FRASER-PCA (0.6554 / 0.9602). Under medium perturbations, the method of this invention achieved an AUPRC of 0.9231, an AUROC of 0.9943, and a best F1 score of 0.8877, significantly exceeding FRASER-PCA (0.8766 AUPRC, 0.9908 AUROC) and FRASER-AE (0.7987 AUPRC, 0.9938 AUROC). Under strong perturbations, the method of this invention remained the best performing method, with an AUPRC of 0.9234, an AUROC of 0.9964, and a best F1 score of 0.9060. The PR and ROC curves further demonstrate that the method of the present invention consistently maintains excellent ranking behavior under weak, medium, and strong perturbation settings, with the most significant advantage in the medium-effect mechanism.

[0144] Isoform usage (iso-pct) outlier detection was evaluated on a GTEx whole blood cohort. The method involved perturbing 1% of the isoform-sample pairs for each gene and renormalizing the matching isoforms within the same gene to maintain compositional constraints. Three perturbation mechanisms were considered: weak (0.2–0.5), medium (0.5–0.7), and strong (0.7–1.0). The method of this invention (MAE-GUIDE) was compared with that of autoencoders-only (AE-Only), generative adversarial networks-only (GAN-Only), and other methods.

[0145] Among all three mechanisms, the method of this invention ranks first overall. At the weak perturbation level, the method of this invention achieves the best AUPRC (0.3560), AUROC (0.9189), and best F1 (0.4517), outperforming AE-Only and GAN-Only. In the medium perturbation mechanism, the method of this invention achieves 0.8604 AUPRC, 0.9859 AUROC, and a best F1 of 0.8119. Under strong perturbation, the method of this invention remains the best performing model, with 0.9818 AUPRC, 0.9989 AUROC, and a best F1 of 0.9418. These results demonstrate that integrating the autoencoder and generative adversarial network components produces more robust and accurate heterogeneous-level anomaly detection than using either component alone.

[0146] This embodiment further examines whether the proposed framework can recover abnormalities independently identified in clinical and functional genetics. Using the same analytical workflow as the benchmark experiments, this method successfully re-detected several disease-related expression aberrations in the Kremer rare disease fibroblast cohort, including TIMMDC1, CLPP, TALDO1, ALDH18A1, MCORN1, COASY, and PANK2. These genes are involved in mitochondrial complex assembly, protein homeostasis, pentose phosphate pathway metabolism, proline biosynthesis, lysosomal ion transport, and coenzyme A biosynthesis, which is highly consistent with previous reports and relevant clinical context of this cohort.

[0147] In addition to total expression levels, this framework also recovered identified splicing-level aberration signals. In the Kremer cohort, several known pathogenic splicing defects were re-identified at the joinpoint level, including events involving CLPP, TIMMDC1, MCOLLN1, and aberrant donor use events in TAZ. These results demonstrate that the method is capable of recovering clinically relevant splicing aberrations in real-world rare disease samples, and is effective not only in simulated perturbation settings.

[0148] At the isoform usage level, evaluation of the GTEx whole blood cohort further demonstrated that the ensemble model consistently outperformed its individual components under weak, moderate, and strong perturbation mechanisms, indicating that combined reconstruction and adversarial strategies can robustly detect constitutive abnormalities in transcript usage. In summary, these results demonstrate that the framework has been successfully transferred from controlled benchmark experiments to biologically and clinically meaningful real-world scenarios, retrieving established pathogenic signals across expression, splicing, and isoform usage levels.

[0149] To investigate the relationship between gene expression outliers, splicing outliers, and transcriptional outliers, a matched multilevel dataset was constructed from the GTEx whole blood cohort. Starting with the cleaned gene count matrix, transcript TPM matrix, and splicing junction matrix, the three analysis layers were linked within a shared sample space, and then the corresponding workflows described above were used to identify outliers. For each sample, the number of detected outliers at the expression, transcript, and splicing levels was counted, and their consistency was assessed using Spearman correlation analysis.

[0150] Analysis results as follows Figure 4 The results showed that the three types of outliers were positively correlated, but the consistency was only moderate. The strongest association was observed between transcript levels and splicing outliers. This indicates that transcriptomic abnormalities are more tightly coupled with splicing variations than with total gene expression shifts. In contrast, expression abnormalities are more closely coupled with transcriptomic abnormalities (…). ) and splicing outliers ( The correlations between the three analytical layers are weak. These results indicate that the three analytical layers capture relevant but non-redundant aspects of transcriptomic dysregulation: splicing abnormalities are more directly reflected at transcript usage levels, while many such events may only be faintly visible in total gene expression. This analysis supports the value of multilayer integration, as signals at expression, splicing, and transcript levels provide complementary perspectives on aberrant molecular states, rather than being interchangeable measures.

[0151] This invention implements a unified distribution-aware framework for detecting outliers in gene expression, splicing, and isoform ratios. The framework utilizes maximum likelihood estimation / Akaike information criterion for distribution-guided branch selection, employs a sliding window strategy to filter pseudo-normal samples, establishes robust reference standards through a Wasserstein generative adversarial network with gradient penalties, uses a lightweight attention autoencoder for denoising and K-fold cross-fit reconstruction, and ensures the reliability of the results through a weighted consensus layer combined with Benjamini-Hochberg FDR control. Experimental results demonstrate that, compared to existing methods, this invention exhibits superior AUPRC and higher consistency with empirical distributions, successfully identifying known disease-related anomalies and providing a scalable, cross-modal, and statistically calibrated new approach for the discovery of transcriptomic anomalies.

[0152] According to another aspect of the embodiments of this application, an electronic device is also provided, including a processor and a memory, wherein the processor is configured to implement the steps of the method when executing a computer program stored in the memory.

[0153] In the above embodiments of the present invention, the descriptions of each embodiment have different focuses. For parts not described in detail in a certain embodiment, please refer to the relevant descriptions of other embodiments.

[0154] In the several embodiments provided in this application, it should be understood that the disclosed technical content can be implemented in other ways. The device embodiments described above are merely illustrative; for example, the division of units can be a logical functional division, and in actual implementation, there may be other division methods. For instance, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the displayed or discussed mutual coupling, direct coupling, or communication connection may be through some interfaces; the indirect coupling or communication connection between units or modules may be electrical or other forms.

[0155] Furthermore, the functional units in the various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit.

[0156] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, read-only memory (ROM), random access memory (RAM), portable hard drives, magnetic disks, or optical disks.

[0157] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for detecting transcriptome abnormalities, characterized in that, Includes the following steps: Acquire transcriptome data, which includes gene expression data, splicing event data, and isoform ratio data; For each feature in the transcriptome data, a negative binomial distribution, a beta distribution, and a normal distribution are fitted respectively. Based on maximum likelihood estimation and the Akaike information criterion, the distribution with the highest goodness of fit is selected as the modeling distribution for that feature. If the feature is isoform ratio data and the goodness of fit is lower than a preset threshold, then no distribution-based statistical anomaly score is calculated. Calculate the statistical anomaly score for each feature-sample pair based on the modeled distribution; Pseudo-normal samples were filtered out from the expression values ​​of each gene using a sliding window; The pseudo-normal samples are used to train a Wasserstein generative adversarial network with gradient penalty, and the transcriptome data is input into the trained Wasserstein generative adversarial network with gradient penalty to obtain first reconstructed data. Based on the first reconstructed data and the first reconstructed error between the first reconstructed data and the transcriptome data, an adversarial reconstruction anomaly score is generated. The transcriptome data is input into the attention autoencoder for cross-fit reconstruction, and the second reconstruction error corresponding to the pseudo-normal sample is used as a standardization benchmark to obtain the second reconstruction data. Based on the second reconstruction error between the second reconstruction data and the transcriptome data, the robust Z-score for the gene and the robust Z-score for the sample are calculated respectively, and the abnormality score matrix is ​​combined to obtain the abnormality score reconstructed by the autoencoder. The consensus score is obtained by weighting and fusing the statistical anomaly score, the adversarial reconstruction anomaly score, and the autoencoder reconstruction anomaly score. The consensus score is converted into a p-value, and the error detection rate is controlled using a Benjamini-Hochberg process. Abnormal events are determined based on the corrected p-value.

2. The method for detecting transcriptome abnormalities as described in claim 1, characterized in that, After acquiring transcriptome data, methods for selecting the distribution with the highest goodness of fit as the modeling distribution for this feature based on maximum likelihood estimation and the Akaike information criterion include: For each feature, its observed data vector is denoted as... Fit the negative binomial distribution, beta distribution, and normal distribution respectively; The parameters of each distribution are estimated by maximizing the log-likelihood function, where maximizing the log-likelihood function is: in, This represents the observed data vector of a certain feature across all samples. The total number of samples; For this feature in the first Measurements in a sample This is the probability density function or mass function of the corresponding distribution; The set of parameters representing the selected candidate distribution; This represents maximizing the log-likelihood function; Calculate the Akaike information content criterion values ​​for each distribution based on the estimated parameters: in, This represents the Akaike information content criterion value. This indicates the number of free parameters in the candidate distribution; These are the parameter estimates obtained by maximizing the log-likelihood function; Compare the Akaike information content criterion values ​​of each distribution, and select the distribution with the smallest Akaike information content criterion value as the modeling distribution for this feature.

3. The method for detecting transcriptome abnormalities as described in claim 1, characterized in that, Methods for filtering pseudo-normal samples from the expression values ​​of each gene using a sliding window include: For each gene, its expression value in all samples is subjected to variance stabilization transformation, and then standardized to obtain the standardized expression value; The standardized expression values ​​are sorted in ascending order to obtain the sorted sequence and the corresponding sample order. For each position in the sorted sequence, construct a symmetrical window containing a fixed number of samples centered on that position; Within each symmetrical window, the local mean and local standard deviation are calculated, and the local bias score is calculated based on the standardized expression value of the sample at the center of the symmetrical window and the local mean and local standard deviation. in, The local deviation score; It is a local mean; Local standard deviation; Standardized expression values; Samples with local deviation scores less than a preset threshold are labeled as pseudo-normal samples for that gene.

4. The method for detecting transcriptome abnormalities as described in claim 1, characterized in that, Methods for obtaining the first reconstruction data include: The transcriptome data were subjected to logarithmic transformation and robust scaling normalization to obtain the preprocessed input matrix; The input matrix is ​​divided into multiple blocks according to features, and a generative adversarial network is trained independently for each block; For each block, a binary mask matrix is ​​constructed based on the pseudo-normal samples, where the position with a mask value of 1 indicates that the corresponding sample is a pseudo-normal sample. A generator is constructed, which adopts an encoder-decoder-encoder structure: the first encoder maps the input to the latent space, the decoder reconstructs the data from the latent vectors, and the second encoder encodes the reconstructed data again to obtain the latent vectors; Construct a discriminant whose output scalar score serves as the Wasserstein discriminant; The generator and the discriminator are trained using Wasserstein loss with gradient penalty, wherein the loss of the discriminator includes the score of the real sample, the score of the generated sample, and the gradient penalty term, and the loss of the generator includes adversarial loss, mask reconstruction loss, and latent space consistency loss. The trained generator is used to reconstruct the input data for each block, and the reconstruction results of each block are merged to obtain the first reconstructed data.

5. The method for detecting transcriptome abnormalities as described in claim 4, characterized in that, A method for generating an adversarial reconstruction anomaly score based on a first reconstruction error between the first reconstructed data and the transcriptome data includes: After inverse transformation back to the original scale, the absolute error between the transcriptome data and the first reconstructed data is calculated as the first reconstruction error; For each feature, the median and median absolute deviation are calculated using the first reconstruction error of the pseudo-normal sample corresponding to that feature, and then the robust Z-score of the feature dimension is calculated for each feature-sample pair. For each sample, the median and median absolute deviation are calculated using the first reconstruction error of all features within that sample, and then a robust Z-score for the sample dimension is calculated for each feature-sample pair. The robust Z-scores of the above feature dimensions and sample dimensions are geometrically averaged to obtain the adversarial reconstruction anomaly score.

6. The method for detecting transcriptome abnormalities as described in claim 1, characterized in that, The method for obtaining the second reconstructed data by inputting the transcriptome data into an attention autoencoder for cross-fit reconstruction and using the second reconstruction error corresponding to the pseudo-normal sample as a standardization benchmark includes: The transcriptome data were subjected to logarithmic transformation and robust scaling normalization to obtain a preprocessed input matrix, which was then transposed into sample-feature form. An attention autoencoder is constructed, comprising: a learnable feature attention layer for weighting input features; an encoder for mapping the weighted input features to latent representations; and a decoder for reconstructing the original input from the latent representations. The attention autoencoder is trained in a denoising manner, that is, Gaussian noise is added to the input features. The training loss includes smooth L1 reconstruction loss, cosine similarity loss, latent representation regularization term and attention weight regularization term. K-fold cross-fitting is used for reconstruction, so that each sample is reconstructed by a model that has never seen the sample before. The reconstruction results are inversely transformed back to the logarithmically transformed scale to obtain the second reconstructed data.

7. The method for detecting transcriptome abnormalities as described in claim 1, characterized in that, Based on the second reconstruction error between the second reconstructed data and the transcriptome data, robust Z-scores for genes and robust Z-scores for samples are calculated separately, and combined to obtain an anomaly score matrix. This method for reconstructing anomaly scores using an autoencoder includes: The absolute error between the transcriptome data and the second reconstructed data is calculated as the second reconstruction error; For each feature, the median and median absolute deviation are calculated using the second reconstruction error of the pseudo-normal sample corresponding to that feature. Then, for each feature-sample pair, the robust Z-score of the feature dimension is calculated, retaining only the positive deviation. For each sample, the median and median absolute deviation are calculated using the second reconstruction error of all features within that sample. Then, for each feature-sample pair, a robust Z-score for the sample dimension is calculated, retaining only the positive bias. The robust Z-score of the feature dimension and the robust Z-score of the sample dimension are geometrically averaged to obtain the basic anomaly score. Calculate the average attention weight of each feature in the attention autoencoder, multiply the base anomaly score by a scaling factor composed of the average attention weight, and obtain the autoencoder reconstructed anomaly score matrix, which is used as the autoencoder reconstructed anomaly score.

8. The method for detecting transcriptome abnormalities as described in claim 1, characterized in that, The method for weighted fusing the statistical anomaly score, the adversarial reconstruction anomaly score, and the autoencoder reconstruction anomaly score to obtain the consensus score includes: For each feature and each sample, the following three standardized anomaly scores are obtained: statistical anomaly score from distribution modeling, adversarial reconstruction anomaly score from generative adversarial networks, and autoencoder reconstruction anomaly score from attention autoencoders. The consensus score is calculated using the following weighted summation formula: in, For consensus score, As a feature, As a sample, To count outlier scores, To combat the reconstruction of abnormal scores, To reconstruct anomalous scores for the autoencoder, , , The weight coefficients are non-negative and satisfy the following conditions: .

9. An electronic device, comprising a memory and a processor, characterized in that, The memory is used to store programs that support a processor in executing the method for transcriptome abnormality detection as described in any of claims 1-8, wherein the processor is configured to execute the programs stored in the memory.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, The computer program, when run by a processor, performs the steps of the method for transcriptome abnormality detection as described in any one of claims 1-8.

Images