Method and system for detecting introgression based on local principal component configuration

Abstract

The application belongs to the technical field of computer-aided genetic analysis, and discloses a gene introgression detection method and system based on local principal component configuration. A genome is divided into multiple windows, and a standardized genotype matrix of each window is constructed and principal component decomposition is performed to obtain a local coordinate matrix of individual x K-dimensional principal components. Stable anchor samples are used to estimate orthogonal transformation. For each candidate donor, the directionality geometric feature and neighborhood topological feature of the recipient-donor relationship relative to the recipient-background relationship are calculated, and are robustly summarized to obtain a window geometric score. The classical population genetic statistical evidence of the window geometric score is fused to obtain a window detection score. Isolated states that do not meet the minimum number of continuous windows are filtered, and candidate introgression segments are output. The application directly analyzes an individual x principal component numerical matrix instead of image pixels, and can provide candidate introgression segments with clear donor direction, repeatable calculation and convenience for classical statistical review.

Description

Technical Field

[0001] This invention relates to the fields of population genomics, bioinformatics, and computer-aided genetic analysis, specifically to a gene introgression detection method and system based on local principal component configurations. The method utilizes geometric configuration changes, neighborhood topological changes, and chromosome continuity of the local principal component numerical matrix to detect candidate gene introgression segments between populations. Background Technology

[0002] Gene introgression is the process by which genetic material from different species, subspecies, varieties, or populations hybridizes and then backcrosses to enter another population. Common detection methods include the D statistic and its derived statistics, the f statistic, relative node depth, nucleotide differences, population differentiation, local ancestry inference, hidden Markov models, and local gene trees or topological weighting analysis. These methods rely on pre-defined four-population topologies, allele frequencies, haplotype reference panels, population history models, or window-by-window tree construction, respectively.

[0003] Local principal component analysis (PCA) can represent population relationships across different genomic windows as low-dimensional coordinate matrices. Existing local PCA methods are mostly used to detect differences in recombination rates, inversions, or abnormal population structure, or for visualization through window distances. When eigenvalues ​​are distinct, the sign of a single eigenvector remains non-unique; when two or more eigenvalues ​​are equal or very close, the basis in the corresponding eigenspace can undergo orthogonal rotation, reflection, axis exchange, or reordering. The label sets and covariance matrices of different windows are also different, therefore, PC axes of the same name do not have a naturally consistent biological orientation. Directly comparing scatter plots or comparing coordinates column-by-column can easily misinterpret the uncertainty of numerical decomposition as biological changes. Converting scatter plots to images and then using image recognition frameworks introduces dependencies on resolution, coordinate range, point coverage, and sample size, which is detrimental to accurate reproduction.

[0004] Furthermore, a single background population may contain infiltration, sample imbalance, or misjudgment of its relationship with the recipient; independent thresholding window by window can produce discontinuous, isolated, high-scoring false positives; and candidate boundaries obtained from fixed-length windows should not be mistaken for precise breakpoints. Therefore, a detection scheme is needed that is centered on a numerical matrix, eliminates uncertainty in the direction of local principal components, can represent the donor orientation, supports multiple background populations, utilizes chromosome continuity, and can be cross-validated with classical statistical methods. Summary of the Invention

[0005] This invention provides a gene introgression detection method and system based on local principal component configurations, aiming to solve the problems of incomparable local PCA configurations, loss of numerical information in image analysis, unstable results under single background, discontinuous high scores in the window, and lack of classical statistical verification in new methods. It also provides the boundary uncertainty range of candidate segments rather than spurious single-base precision.

[0006] Specific technical solutions:

[0007] A gene introgression detection method based on local principal component conformation includes the following steps:

[0008] S1, acquire data including recipient population, at least one candidate donor population, at least one background population, and individual genotypes and marker locations;

[0009] The minimum number of background groups is 1. When the number of background groups is greater than 1, the receptor-background control is calculated for each background group, and a consensus value is formed by the median or truncated mean. At the same time, the dispersion between background groups is output to identify unstable windows caused by background group selection.

[0010] S2, the genome is divided into multiple windows according to the marker positions, and a standardized genotype matrix is ​​constructed for each window;

[0011] The window can be a fixed physical length window, a fixed number of markers window, a window based on the recombined graph, or a window based on linked unbalanced blocks, and adjacent windows may or may not overlap.

[0012] S3, perform principal component decomposition on each standardized genotype matrix to obtain the local coordinate matrix of individual × K-dimensional principal components;

[0013] K is greater than or equal to 2, and the local principal component numerical matrix is ​​retained instead of the scatter image as the decision input.

[0014] S4. Use pre-determined stable anchor samples to estimate the orthogonal transformation and align each local coordinate matrix to the reference coordinate matrix;

[0015] The stable anchoring samples are derived from one or more background populations or outgroups that do not act as recipient-donor pairs in the candidate infiltration relationship; the orthogonal transformation is obtained by minimizing the Frobenius norm between the local coordinates of the anchoring samples and the reference coordinates, without performing isotropic scaling that would eliminate variations in population dispersion.

[0016] S5. For each candidate donor, calculate the directional geometric features and neighborhood topological features of the acceptor-donor relationship relative to the acceptor-background relationship, and robustly summarize the features obtained from multiple background groups to obtain the window geometric score.

[0017] The directional geometric features and neighborhood topological features include at least two of the following: Euclidean distance difference and Mahalanobis distance difference between the population centroids; energy distance, maximum mean difference and slice Wasserstein distance; k-nearest neighbor donor connection ratio and minimum spanning tree cross-population edge ratio; Procrustes residual, centered kernel alignment difference, Grassmann subspace distance, covariance trace and anisotropy.

[0018] S6, The window geometric score is used alone or fused with classical population genetic statistical evidence of the same window to obtain the window detection score;

[0019] Furthermore, the Patterson D statistic is used as the classic baseline; the adopted statistic is robustly standardized in a consistent direction to form a classic evidence score, which is then fused with the window geometric score with preset weights, while the scores before and after fusion are saved separately.

[0020] S7 identifies candidate infiltration states based on chromosome continuity of adjacent window detection scores, filters out isolated states that do not meet the minimum number of consecutive windows, and outputs candidate infiltration segments.

[0021] Furthermore, a hidden Markov model containing background state and candidate infiltration state is used to perform state decoding on the window detection score; the influence of isolated high-scoring windows and isolated low-scoring windows is reduced by state transition probability; only candidate segments that reach the preset minimum number of consecutive windows are retained, and the midpoint of the first and last positive windows is used as the boundary point for estimation, and the corresponding window range is used as the boundary uncertainty interval, or a second round of boundary refinement is performed using a smaller window within this interval.

[0022] A gene infiltration detection system based on local principal component configuration is used to implement the gene infiltration detection method based on local principal component configuration. The system includes: a data acquisition module, a local dimensionality reduction module, an anchoring alignment module, a feature calculation module, a score fusion module, a continuous segment identification module, and an output and verification module.

[0023] The data acquisition module is used to acquire individual genotype data and marker location information of the recipient population, at least one candidate donor population, and at least one background population, and to divide the genome into windows and construct a normalized genotype matrix for each window based on the marker location information; the local dimensionality reduction module is used to perform principal component decomposition on the normalized genotype matrix of each window to generate a local coordinate matrix of individual × K-dimensional principal components; the anchoring and alignment module is used to estimate orthogonal transformation based on stable anchored samples and align each local coordinate matrix to a reference coordinate matrix; the feature calculation module is used to calculate the directional geometric features and neighborhood topological features of the recipient-donor relationship relative to the recipient-background relationship for each candidate donor, and to calculate the features obtained from multiple background populations. Features are robustly aggregated to obtain window geometric scores; the score fusion module is used to use the window geometric scores as window detection scores alone, or to perform robust standardization and weighted fusion with classical population genetic statistical evidence of the same window in a consistent direction to obtain window detection scores; the continuous segment identification module is used to decode candidate infiltration states based on the chromosome continuity of adjacent window detection scores, filter out isolated states that do not meet the preset minimum number of continuous windows, and identify candidate infiltration segments; the output and verification module is used to output candidate infiltration segments, boundary point estimates, and boundary uncertainty intervals, and save window geometric scores, classical evidence scores, and fusion scores for verification by classical population genetic statistical evidence or known ground truth data.

[0024] The technical effects of this invention are as follows:

[0025] (1) Directly process the numerical matrix of individual × principal components to avoid information loss and plotting parameter dependence caused by image rasterization;

[0026] (2) Eliminate the uncertainty of sign, rotation and reflection of local PCs by orthogonally aligning stable anchor populations, while preserving biological variations in population scale and discreteness;

[0027] (3) Distinguish between general local structural abnormalities and candidate infiltrations directed by the donor by comparing the directionality of the receptor relative to the donor and the background;

[0028] (4) Supports separate calculation and robust aggregation for multiple background groups, and can explicitly report background selection instability;

[0029] (5) Use HMM and minimum number of consecutive windows to suppress isolated false positives and provide a boundary uncertainty range that matches the window resolution;

[0030] (6) Retain pure geometric fractions, classical consensus fractions, and fusion fractions to facilitate the interpretation of evidence sources;

[0031] (7) No need to build gene trees window by window, matrix operations can be divided into blocks, parallelized and ported to Linux high-performance computing environment.

[0032] The key technology of this invention does not lie in individually proposing PCA, Procrustes, distance statistics, or HMM, but in establishing a mutually constrained data processing chain for local gene infiltration detection: after the local PCA numerical matrix is ​​aligned with a stable population anchor, the directional geometry and neighborhood topological changes of the recipient-donor relationship relative to the recipient-background are calculated, followed by robust summarization of multiple backgrounds, parallel or fusion of classical evidence, and identification of continuous chromosome segments. This combination allows coordinate uncertainty, donor orientation, background sensitivity, window continuity, and boundary resolution to be explicitly handled in the same computational flow. See Table 1 for details.

[0033] Table 1. Technical solutions of the present invention to address potential technical problems.

[0034] PC symbol flipping, axis swapping, or near-degenerate subspace rotation Misjudging the same configuration as a window change Orthogonal Procrustes alignment with stable background / outgroup anchors, while preserving Procrustes residuals. When the anchored population itself is abnormal, it cannot be forcibly corrected. Convert PCA scatter plot to image Affected by resolution, coordinate range, point coverage and drawing parameters Directly analyze the individual × K-dimensional PC numerical matrix; the image is for illustrative purposes only. The input sample order and labeling quality control must remain fixed. Local structural anomalies are often mistaken for infiltration. Inversion, selection, or batch effects can also produce anomalous PCA. Calculate the directional contrast between receptor-donor and receptor-background and reference configurations. Verification with classical statistical and biological evidence is still needed. Incorrect background selection or background bleeding. Candidate signals depend on specific background Calculate the values ​​for multiple backgrounds separately, sum them up using the median, and output the standard deviation between backgrounds. Error topology backgrounds should not be simply merged. Imbalanced sample size and a small number of outliers Distance, nearest neighbor, and graph indices shifted. Within the population, normalized mean, robust standardization, and sensitivity analysis that allows for equal resampling. Even very small populations can still cause instability in covariance and kNN. Singular or ill-conditioned covariance matrix Mahalanobis distance cannot be stably inverted. Add a λI regularization term to the merged covariance matrix λ should be pre-defined and sensitivity analysis performed. Per-window threshold generates isolated high scores or internally isolated low scores This can lead to discontinuous false positives or interruption of the true range. HMM state transition constraints, minimum number of consecutive windows, and physical large-interval chain breakage. Conservative parameters will reduce the sensitivity to weak penetration. Fixed window misinterpreted as precise breakpoint Creating false single-base boundaries Output the first and last positive points in the window and the uncertain interval of the complete window boundary. Fine-grained boundaries require nested small window or singlet methods. The source of the combined evidence is not transparent. Unable to determine whether classical statistics or geometric features are driving the process. Store pure geometry, classical consensus, and fusion scores separately, and mark evidence conflicts. The fusion weights must be predetermined using independent data. Building trees window by window involves large computational costs and results in unstable topology. Difficult to scale to high-density whole-genome data Computation is performed using block-based linear algebra, distance graphs, and neighborhood graphs, eliminating the need for window-by-window tree construction. This method is complementary to, but not equivalent to, local tree evidence.

[0035] Orthogonal alignment only eliminates coordinate representation differences that can be determined by stable anchored samples, without performing isotropic scaling or eliminating real local configuration changes that cannot be explained by orthogonal transformation; the latter is further characterized by Procrustes residuals, CKA differences, Grassmann distance, and donor orientation features. Detailed Implementation

[0036] The following embodiments are used to explain the present invention, not to limit the scope of protection.

[0037] The system used in this embodiment is a gene infiltration detection system based on local principal component configuration. The system includes a data acquisition module for acquiring individual genotype data and marker location information from the recipient population, at least one candidate donor population, and at least one background population, and for dividing genomic windows and constructing standardized genotype matrices for each window based on the marker location information; a local dimensionality reduction module for performing principal component decomposition on the standardized genotype matrices of each window to generate individual × K-dimensional principal component local coordinate matrices; an anchoring and alignment module for estimating orthogonal transformations based on stable anchored samples and aligning each local coordinate matrix to a reference coordinate matrix; and a feature calculation module for calculating the directional geometric features and neighborhood topology of the recipient-donor relationship relative to the recipient-background relationship for each candidate donor. The system identifies and robustly aggregates features from multiple background populations to obtain window geometric scores. A score fusion module is used to either use the window geometric score as a separate window detection score or to robustly normalize and weight it with classical population genetic statistical evidence from the same window to obtain a window detection score. A continuous segment identification module decodes candidate infiltration states based on the chromosomal continuity of adjacent window detection scores, filters out isolated states that do not meet the preset minimum number of continuous windows, and identifies candidate infiltration segments. An output and verification module outputs candidate infiltration segments, boundary point estimates, and boundary uncertainty intervals, and saves the window geometric score, classical evidence score, and fusion score for verification by classical population genetic statistical evidence or known ground truth data.

[0038] Unless otherwise specified, the genotype dose is 0, 1, or 2, and missing values ​​are filled in within the window using the mean of the markers; all population labels are fixed before analysis, and candidate segments are not included in weighting.

[0039] Group roles and anchored group definitions:

[0040] The population names in this specification indicate the analytical role of the sample in the current test, as shown in Table 2, rather than newly added biological classifications. The same natural population may play different roles in different pre-designed tests, but the role within a single test should be fixed before analysis.

[0041] Table 2. Group names indicate the analytical role of the sample in the current test.

[0042] Receptor population R The target population to be determined whether a certain segment accepts foreign genetic components Evaluate whether it is relatively close to the candidate donor in terms of background. Candidate donor group D Potential genetic source populations pre-selected based on breeding history, geographical origin, or phylogenetic relationships. Provides the direction of penetration; the detection results support the proximity relationship with the candidate donor. Background Group B A control group with a well-defined relationship to the recipient, used to represent the comparison relationship when no target donor is introduced. Construct a receptor-background control; calculate multiple backgrounds separately and then perform a robust summary. Outer group O Outside of the predefined population topology, a more distantly related population used to determine ancestry and derived allele status. Primarily used for four-group statistics such as Patterson D; not all geometric analyses require its use. Anchored sample set A A pre-selected set of reference individuals from one or more stable background groups and / or outgroups. It is only used to estimate the rotation, reflection, and translation from the local PCA to the reference PCA, and does not directly represent the donor orientation.

[0043] "Anchored group" is a shorthand for the group from which the sample is anchored; what is actually used in the calculation is the sample set. Anchor samples should have sufficient effective genotypes in each window, the number of samples should be greater than the number of principal components retained, and there should be no known strong infiltration, large fragment inversions, or systematic deletions in the region to be examined. Anchor samples and their source populations should be fixed before whole-genome scanning and should not be changed window by window based on the results. The recipient and candidate donor used in this test are generally not used as primary anchor samples to avoid their true proximity changes being partially absorbed by the alignment process.

[0044] Step 1: Input, window, and local principal component matrix;

[0045] The input must include at least the SampleID, Group label, genotype dose for each marker, and SNP number, chromosome, and physical location. The SampleID must be unique and consistent across all files. (Click window) Obtain One sample and 1 valid token, of which Display window The number of markers that passed quality control; subscript This indicates that the quantity varies depending on the window size. (Mark) The dosage is standardized as follows: In the formula: For individual i, label the dose of j as 0, 1, or 2 within window w; This serves as the reference allele frequency for the marker; Standardized dosage. Markers with a frequency of 0 or 1, excessively high missing rates, or insufficient variance were removed.

[0046] For the standardized matrix Perform singular value decomposition and take the first-to-last value. Each component is used as a local coordinate: ;

[0047] In the formula: For n× Standardized genotype matrix; and It is a left and right singular vector matrix; It is a singular value diagonal matrix; Let K be the coordinate matrix of the individual principal components, which is n×K; K is usually taken as 2 to 10 and is determined by the sample size, the number of labels, and the cumulative explanatory power.

[0048] This invention preserves Instead of using PC1-PC2 scatter plot pixels as model input, the images are used for manual verification only. Windows can overlap; smaller windows have higher boundary resolution, but lower-marked windows have increased variance.

[0049] Step 2: Stabilize and orthogonally align;

[0050] The sign of a PCA eigenvector is not unique; when corresponding eigenvalues ​​are equal or close, the same eigenvector... The principal component subspace can also be represented by different orthogonal bases and may undergo rotation, reflection, or axis exchange under numerical perturbation. Different windows are formed by different label sets; therefore, unprocessed local coordinate columns cannot be directly compared column-by-column across windows. This embodiment constructs reference coordinates. Reference coordinates can be calculated using sparse whole-genome labeling, pre-determined neutral labeling, or a "leave one chromosome" label. "Leave one chromosome" means that when analyzing a chromosome, only the labels of the remaining chromosomes are used to establish reference PCA coordinates, thereby avoiding the entry of actual infiltration signals from the chromosome under test into the reference configuration and causing information leakage; each chromosome under test establishes its own corresponding reference. The aforementioned pre-determined stable anchoring sample set is used. Solve for each window: In the formula: Q is a K×K orthogonal matrix; I is the identity matrix; and These are the submatrices of the anchored samples in the local and reference coordinates, respectively; the double vertical lines in the formula represent the Frobenius norm; This is the optimal rotation or reflection transformation.

[0051] Center the local coordinates of all individuals and map them to the reference coordinate system: In the formula: The coordinates are the aligned local coordinates; 1 represents an n-dimensional vector of all 1s. and These are the centroids of the anchored sample in the local and reference coordinates, respectively. This embodiment does not perform isotropic scale normalization to avoid eliminating population diffusion or contraction information.

[0052] The Procrustes transformation is fitted only to the anchored sample; therefore, the residual measures the difference between the local configuration of the anchored sample and the reference configuration that cannot be explained by rotation or reflection. Infiltration of the recipient population as a whole or some individuals towards the donor does not directly increase the residual calculated based on the anchored sample; such changes are captured by recipient-donor orientation features and individual neighborhood features. The anchoring residual may only increase when the anchored population itself undergoes strong infiltration, inversion, intrapopulation heterogeneity, or a large number of missing anchored samples within the window. This residual primarily serves as a quality control measure and an indicator of non-specific configurational anomalies, and cannot independently prove infiltration.

[0053] Therefore, anchored Procrustes residuals are not a direct measure of infiltration completeness. If we consider the local morphology of the recipient population itself, incomplete infiltration can increase intrapopulation dispersion and anisotropy due to only some individuals or haplotypes moving closer to the donor; conversely, when donor genetic material undergoes near-complete and relatively uniform substitution in the recipient, internal recipient deformation may actually decrease. Regardless of whether internal recipient deformation increases or decreases, the directional distance and neighborhood indices between the recipient and donor relative to the recipient and background are still used to capture signals of their approach to the donor.

[0054] Step 3: Summary of background group size and robustness;

[0055] The minimum available configuration is one receptor R, one candidate donor D, and a background with a well-defined relationship to the receptor. If we calculate D and Furthermore, outgroup O is also required. The minimum configuration can calculate directional contrast, but cannot determine whether the results depend on a specific background. The recommended configuration is two or more independent background populations that are known to be related to the receptor and have no obvious candidate infiltrations. The program does not directly merge multiple backgrounds, but calculates them one by one and then summarizes them, as shown in Table 3.

[0056] Table 3. Background Group Size and Robustness Summary

[0057] Minimum configuration 1 Exploratory positioning Background selection sensitivity cannot be assessed Recommended configuration ≥2 Stable positioning and consistency of direction Calculate for each background, summarize the medians, and report the standard deviation. Multiple backgrounds but different topologies ≥2 Sensitivity analysis Simple merging should be avoided; stratification should be based on pre-defined phylogenetic relationships. Background appears to be self-infiltration. any Not suitable as the main background Excluded or used solely as evidence of conflict

[0058] For feature k and H background populations, principality and heterogeneity are defined as follows: In the formula: The k-th feature is calculated for window w relative to the h-th background group; A consensus based on the median across multiple contexts; This represents the standard deviation between background values. When H=1, only a single background value is output, without interpreting background heterogeneity.

[0059] Step 4: Distance, similarity, and topological metrics;

[0060] The indicators used are divided into four layers according to their function. Referring to the implementation, all available indicators in Table 4 are calculated simultaneously. After directional unification, robust standardization, and preset weighting, a geometric score is formed, and the original value of each indicator is retained for interpretation; therefore, these indicators are not for users to choose arbitrarily based on preference afterward. As an optional implementation of this invention, when the sample size is insufficient, computational resources are limited, or parameters have been predetermined by independent data, a pre-specified subset of indicators can also be used, but it should at least include one donor direction indicator and one non-Euclidean distribution or neighborhood indicator, and should not be temporarily selected based on the results of the section to be inspected. Each indicator directly acts on the aligned local coordinate matrix.

[0061] Table 4 Available Indicators

[0062] Geometric direction Centroid Euclidean distance difference The group center is relatively close Affected by outliers and scale Geometric direction Mahalanobis distance difference Center proximity after considering covariance Small samples require covariance regularization. Distribution differences Energy distance Overall difference in arbitrary dimension distribution Resampling is advisable when sample sizes are imbalanced. Distribution differences RBF-MMD Kernel mean embedding differences Bandwidth needs to be fixed or use the median rule. Distribution differences Sliced ​​Wasserstein Multi-directional one-dimensional transportation distance The random projection direction should have a fixed seed. Neighborhood topology kNN donor connection ratio Donor enrichment in local neighbors of the receptor k must be less than the number of valid samples. Neighborhood topology MST cross-group ratio Mixture of two groups on the minimum spanning tree Parallel distances require stable handling Overall configuration Procrustes residuals Changes in local and reference dimensions that cannot be explained by orthogonal transformations Non-donor specific Overall configuration Linear CKA Difference Inter-sample inner product structural similarity Centralization is required Overall configuration Grassmann distance Principal component space angle K needs to be consistent across windows. Group morphology Covariance trace, anisotropy Receptor population diffusion and stretching Direction cannot be stated alone.

[0063] Taking the Euclidean distance between the centroids as an example, the donor orientation comparison is as follows: In the formula: , , These represent the centroids of the receptor, background, and donor in window w, respectively; the subscript 0 indicates the reference configuration; and d represents the Euclidean distance. An increase indicates that the receptor is closer to the donor and farther from the background than the reference.

[0064] The distance between the centroids of Mahalanobis is: In the formula: and Let the centroids of groups a and b be denoted as 'a' and 'b', respectively. Let be the combined covariance matrix of the two groups; λ is the regularization coefficient; I is the identity matrix. The construction of the donor orientation comparison is the same as that in equation (6).

[0065] Energy distance is used to compare two sets of multidimensional coordinate distributions: In the formula: X and X′ are independent coordinate vectors from the same population distribution; Y and Y′ are from another population; ‖·‖ is the Euclidean norm. In actual calculations, the mean of pairwise distances between samples is used instead of the expected value.

[0066] The maximum mean difference of the RBF kernel is: In the formula: and Let k(x,y) = exp(-‖xy‖^2 / (2 ^2)) is an RBF core; The bandwidth is defined as the kernel bandwidth. In this embodiment, the bandwidth is determined by the median of the non-zero distance between the merged samples.

[0067] The Wasserstein distance of the slice is: In the formula: K is the unit direction; L is the fixed number of random directions; The distance is a one-dimensional, first-order Wasserstein distance. A fixed random seed ensures reproducibility across different runs.

[0068] The k-nearest neighbor donor affinity of the recipient sample is: In the formula: Number of receptor samples; (i) represents the k nearest neighbors of individual i in the merged sample of recipient, donor, and background; I is the indicator function. The actual change is formed by subtracting the reference value from the window value.

[0069] The minimum spanning tree has the following cross-population edge ratio: In the formula: The edge set of the Euclidean spanning tree formed by the minimum spanning tree of the recipient and comparison group samples; | represents the number of sides; The group to which individual i belongs. The increase in the receptor-donor crossover ratio relative to the receptor-background crossover ratio supports donor affinity.

[0070] The linearly centered kernel alignment similarity between the local configuration and the reference configuration is: In the formula: and They are respectively and The linear Gram matrix after sample centering; tr is the matrix trace. The program uses 1-CKA as the configurational variance measure.

[0071] Step 5: Geometric fractions, classical evidence fusion, and conflict marking;

[0072] Robust standardization was performed on each directional feature using the median and median absolute deviation, and after limiting extreme values, a pure geometric fraction was formed: In the formula: This is a set of donor orientation features; It is a set of non-specific configurational variation characteristics; and For the corresponding characteristic number; is the robust standard score; clip is the truncation function. The coefficient 0.2 is a preset value for the example and can be determined in advance through independent simulation.

[0073] Within the same window and for the same population definition, the Patterson D statistic is used as the primary classic baseline. For sensitivity analyses requiring multiple evidence, a separate statistic can be calculated. The dXY difference and Hudson FST difference were analyzed, and under a pre-defined scheme, all directions were unified to increase the numerical value to support the construction of classical consensus after the recipient-donor proximity. This is an optional supplement and not a necessary prerequisite for choosing this method. In the formula: z represents robust standardization; and Calculated according to the four-group order (((B,R),D),O); and These represent absolute differences between the receptor and background, and between the receptor and donor, respectively; the FST subscript has the same meaning.

[0074] The final fusion score is: In the formula: The fusion score is α, and the classical evidence weight is α=0.25. In this embodiment, α=0 represents pure geometric detection. α must be determined through independent simulation, pre-registered parameters, or reserved chromosomes, and is not optimized on the region to be detected.

[0075] Here, "fusion" refers to combining the PCA local configurational geometric scores of this invention with classic population genetic statistical evidence such as Patterson D according to preset weights, rather than simply equating the aforementioned multiple geometric distances with multiple independent methods. The program simultaneously outputs the pure geometric score, the classic evidence score, and the fusion score. When the pure geometric score is consistent with the direction of classic evidence, it is marked as jointly supported; when the geometric score is high and the classic score is low, it is marked as a configuration-specific candidate, which may reflect complex topology, inversion, or donor mismatch; when the classic score is high and the geometric score is low, it is marked as a frequency-sharing candidate, which requires examination of minority locus drivers, sample size, and reference population. This parallel output avoids the inability to trace the source of evidence after fusion.

[0076] Step Six: HMM Continuity, Discontinuity False Positive Correction and Boundaries;

[0077] Sort each chromosome by window position, and set background state 0 and candidate infiltration state 1. Let... ( () represents the emission density in state s. The state transition probabilities are calculated using Viterbi recursion: In the formula: (s) is the maximum log probability that the window ends in w and the state ends in s; r is the previous state; For the transition probability; This refers to emission density. A larger one... and This allows isolated high or low scores to be corrected when the context does not support them.

[0078] Viterbi decodes only retains at least Status 1 window. This embodiment... If a single low-scoring window appears during actual penetration, the HMM can fill it due to the persistence probability of state 1; if there is only one or a few isolated high-scoring windows, the state transition cost and They should be filtered out. For large physical gaps caused by marker gaps, the transfer probability should be adjusted according to the actual distance or the chromosome strand should be broken; blind connection should be avoided.

[0079] For the first positive window and the last positive window The boundary point estimates and uncertainty intervals are as follows: In the formula: and Estimate the left and right boundary points respectively; and For the boundary uncertainty interval; and This represents the start and end positions of the window. This range clearly reflects the resolution of the first round of windows.

[0080] When more precise boundaries are required, and Within the extended intervals on both sides, the features and states are recalculated with smaller windows and step sizes until the label density or computational cost is limited; the final boundaries are still reported as intervals. If phase haplotype is present, the local ancestor or haplotype breakpoint method can be used again within the candidate segment.

[0081] Step 7: Calculation of classical statistics;

[0082] To ensure that the new method and the traditional method are compared on the exact same window, this embodiment calculates the D statistic according to (((P1 background, P2 recipient), P3 donor), O outgroup): In the formula: and The expected schema weights of j, calculated based on allele frequencies across the four populations, are used to label j; a positive D indicates that P2 shares more alleles with P3 than P1 shares alleles with P3. When outpopulations are not fixed, the complementarity schema is also considered.

[0083] window The statistic is written as: In the formula: , , , The allele frequencies of P1, P2, P3 and outgroups at marker j are respectively. =max( , The larger donor orientation frequency is defined for each site. Windows with denominators close to 0 are considered missing.

[0084] The absolute nucleotide differences between the two groups a and b are: In the formula: The number of valid markers for the window; and Allele frequencies in the two populations. The program uses... - As a comparison of donor directions.

[0085] Hudson FST calculates the ratio by summing the numerator and denominator of the site, and then forms... - contrast: In the formula: and This indicates the effective number of alleles of j in populations a and b; other symbols are the same as before. An increase in this contrast indicates that the recipient is relatively closer to the donor.

[0086] Step 8: Compare the simulation with the known truth value with the classical method;

[0087] To verify the program's operability and the usability of the detection scores, 120 diploid individuals were constructed, including 30 background individuals (B), 30 recipient individuals (R), 30 donor individuals (D), and 30 outgroups (O). The chromosome length was 12 Mb, containing 2400 markers. Known introgression regions from D to R were set at 4–6 Mb. The window length was 500 kb, the step size was 250 kb, K=3, and the minimum number of consecutive positive windows was 4. Three settings were used: weak, medium, and strong, with recipient carrier ratios of 0.25, 0.50, and 0.75, and donor allele substitution probabilities in the introgression regions of 0.35, 0.50, and 0.70, respectively. Each setting was replicated 10 times using independent random seeds, with an additional 10 replicates as a negative control without introgression.

[0088] All methods used the same window and truth labels. ROC-AUC measures window ranking, average precision measures the concentration of positive windows, and Top-K base IoU measures the cross-union ratio between the union of selected windows and the true segment. The results are shown in Table 5 below.

[0089] Table 5 Comparison Results

[0090] weak LPGA Geometry 0.640±0.133 0.409±0.137 0.222 weak Patterson D 0.657±0.099 0.399±0.139 0.211 middle LPGA Geometry 0.911±0.063 0.869±0.046 0.630 middle Patterson D 0.781±0.106 0.562±0.182 0.336 powerful LPGA Geometry 0.932±0.040 0.909±0.041 0.635 powerful Patterson D 0.912±0.059 0.847±0.094 0.584

[0091] The results show that under weak penetration, the ROC-AUCs of LPGA geometry and Patterson D are 0.640 and 0.657, respectively, with average precisions of 0.409 and 0.399, indicating that both methods have limited ability to distinguish weak signals. Under moderate penetration, the ROC-AUCs are 0.911 and 0.781, respectively, with average precisions of 0.869 and 0.562. Under strong penetration, the ROC-AUCs are 0.932 and 0.912, respectively, with average precisions of 0.909 and 0.847. This embodiment uses Patterson D as the primary classical baseline, and the comparison aims to evaluate whether local PCA geometric evidence provides supplementary window ranking capabilities, rather than claiming a complete replacement for classical methods.

[0092] HMM decoding and setting using pure geometric fractions In moderate infiltration, continuous segments were detected in 9 out of 10 replicates, with an average precision of 0.809, recall of 0.874, and IoU of 0.783. In strong infiltration, all 10 replicates detected continuous segments, with an average precision of 0.849, recall of 0.999, and IoU of 0.848. Weak infiltration failed to form segments meeting the continuity requirement. The 10 negative controls without infiltration did not output any candidate segments. These results indicate that the continuity rule reduces isolated high-scoring false positives, but it also sacrifices sensitivity for weak signals.

[0093] In the default single-shot example, the true range is 4.00–6.00 Mb, and the program outputs the original candidate window union, which is approximately 3.51–6.01 Mb. The uncertain range for the left boundary is approximately 3.51–4.01 Mb, and the uncertain range for the right boundary is approximately 5.51–6.01 Mb, both covering the true boundary. This example verifies the meaning of the boundary ranges and shows that points within the window should not be interpreted as precisely reconstructed breakpoints.

[0094] The above simulations are for program functionality and relative comparison verification and do not represent all population histories, LD, ILS, or mismatch donor scenarios. For formal applications and research use, common ancestor or forward simulations, real known hybridization materials, sample size imbalance, deletion rate, recombination rate, background mismatch, and outpopulation mismatch tests should be added, and independent software should be used to verify D or .

[0095] Step Nine: Evaluation Indicators; In the formula: Let i be the detection score for window i; For truth labels; and These represent the number of positive and negative windows, respectively; rank is the merged score rank. AUC is the probability that the score of a random positive window is higher than that of a random negative window. In the formula: windows are arranged in descending order of score; y_((k)) is the label of the k-th window; Precision(k) is the positive rate of the first k windows; N is the number of valid windows. In the formula: S is the set of genomic regions merged by the selected window; T is the set of actual infiltrated regions; |·| represents the base length. This indicator penalizes both overly wide boundaries and omissions.

Claims

1. A gene introgression detection method based on local principal component conformation, characterized in that, Includes the following steps: S1, acquire data including recipient population, at least one candidate donor population, at least one background population, and individual genotypes and marker locations; S2, the genome is divided into multiple windows according to the marker positions, and a standardized genotype matrix is ​​constructed for each window; S3, perform principal component decomposition on each standardized genotype matrix to obtain the local coordinate matrix of individual × K-dimensional principal components; S4. Use pre-determined stable anchor samples to estimate the orthogonal transformation and align each local coordinate matrix to the reference coordinate matrix; S5. For each candidate donor, calculate the directional geometric features and neighborhood topological features of the acceptor-donor relationship relative to the acceptor-background relationship, and robustly summarize the features obtained from multiple background groups to obtain the window geometric score. S6, The window geometric score is used alone or fused with classical population genetic statistical evidence of the same window to obtain the window detection score; S7 identifies candidate infiltration states based on chromosome continuity of adjacent window detection scores, filters out isolated states that do not meet the minimum number of consecutive windows, and outputs candidate infiltration segments.

2. The gene introgression detection method based on local principal component conformation according to claim 1, characterized in that, In S1, the minimum number of the background population is 1; when the number of the background population is greater than 1, the receptor-background control is calculated for each background population, and a consensus value is formed by the median or truncated mean. At the same time, the dispersion between backgrounds is output to identify unstable windows caused by background population selection.

3. The gene introgression detection method based on local principal component conformation according to claim 1, characterized in that, In S2, the window is a fixed physical length window, a fixed number of markers window, a window based on the recombined graph, or a window based on the chain of unbalanced blocks, and adjacent windows may overlap or not overlap.

4. The gene introgression detection method based on local principal component conformation according to claim 1, characterized in that, In S3, K is greater than or equal to 2, and the local principal component numerical matrix is ​​retained instead of using the scatter image as the decision input.

5. The gene introgression detection method based on local principal component conformation according to claim 1, characterized in that, In S4, the stable anchoring samples are derived from one or more background populations or outgroups that do not act as recipient-donor pairs in the candidate infiltration relationship; the orthogonal transformation is obtained by minimizing the Frobenius norm between the local coordinates and the reference coordinates of the anchoring samples, without performing isotropic scaling that would eliminate variations in population dispersion.

6. The gene introgression detection method based on local principal component conformation according to claim 1, characterized in that, In S5, the directional geometric features and neighborhood topological features include at least two of the following categories: Euclidean distance difference and Mahalanobis distance difference between the population centroids; energy distance, maximum mean difference and slice Wasserstein distance; k-nearest neighbor donor connection ratio and minimum spanning tree cross-population edge ratio; Procrustes residual, centered kernel alignment difference, Grassmann subspace distance, covariance trace and anisotropy.

7. The gene introgression detection method based on local principal component conformation according to claim 1, characterized in that, In S6, the Patterson D statistic is used as the classic baseline. The classic evidence score is formed by performing a robust standardization of the adopted statistic in a consistent direction, and then fused with the window geometric score with preset weights. At the same time, the scores before and after fusion are saved separately.

8. The gene introgression detection method based on local principal component conformation according to claim 1, characterized in that, In S7, a hidden Markov model containing background state and candidate infiltration state is used to decode the window detection score; the influence of isolated high-scoring windows and isolated low-scoring windows is reduced by the state transition probability. Only candidate segments that reach the preset minimum number of consecutive windows are retained, and the midpoint of the first and last positive windows is used as the boundary point for estimation, and the corresponding window range is used as the boundary uncertainty interval, or a second round of boundary refinement is performed using a smaller window within this interval.

9. A gene introgression detection system based on local principal component configuration, characterized in that, The system is used to implement the gene infiltration detection method based on local principal component configuration as described in any one of claims 1 to 8; the system includes: a data acquisition module, a local dimensionality reduction module, an anchoring alignment module, a feature calculation module, a score fusion module, a continuous segment identification module, and an output and verification module; The data acquisition module is used to acquire individual genotype data and marker location information of the recipient population, at least one candidate donor population, and at least one background population, and to divide the genome into windows and construct a normalized genotype matrix for each window based on the marker location information; the local dimensionality reduction module is used to perform principal component decomposition on the normalized genotype matrix of each window to generate a local coordinate matrix of individual × K-dimensional principal components; the anchoring and alignment module is used to estimate orthogonal transformation based on stable anchored samples and align each local coordinate matrix to a reference coordinate matrix; the feature calculation module is used to calculate the directional geometric features and neighborhood topological features of the recipient-donor relationship relative to the recipient-background relationship for each candidate donor, and to calculate the features obtained from multiple background populations. Features are robustly aggregated to obtain window geometric scores; the score fusion module is used to use the window geometric scores as window detection scores alone, or to perform robust standardization and weighted fusion with classical population genetic statistical evidence of the same window in a consistent direction to obtain window detection scores; the continuous segment identification module is used to decode candidate infiltration states based on the chromosome continuity of adjacent window detection scores, filter out isolated states that do not meet the preset minimum number of continuous windows, and identify candidate infiltration segments; the output and verification module is used to output candidate infiltration segments, boundary point estimates, and boundary uncertainty intervals, and save window geometric scores, classical evidence scores, and fusion scores for verification by classical population genetic statistical evidence or known ground truth data.