Method and system for type 1 diabetes risk prediction based on eQTL and HLA typing

Abstract

The application discloses a type 1 diabetes risk prediction method and system based on eQTL and HLA typing, obtains high-resolution HLA typing and genome SNP genotype data of a specific population queue, combines single-cell transcriptome data to screen an immune cell-specific cis-eQTL site set, integrates eQTL sites, specific population risk HLA alleles and public clinical phenotype SNP features to obtain an initial genetic feature set, refines features through a backward elimination iterative algorithm, adopts an L-BFGS-B quasi-Newton algorithm with L2 regularization and coefficient boundary constraints to optimize weights, and constructs a GRS genetic risk score prediction model. The genotype of a to-be-detected subject is detected, and a genetic risk score is calculated. The application is matched with a corresponding prediction system, improves the risk prediction accuracy of early-onset type 1 diabetes, and is suitable for genetic screening of high-risk populations.

Description

Technical Field

[0001] This invention relates to the fields of molecular genetics, bioinformatics, and disease risk assessment, and in particular to a method and system for predicting the risk of type 1 diabetes based on eQTL and HLA typing. Background Technology

[0002] Type 1 diabetes is a T-cell-mediated autoimmune disease in which the body's immune system mistakenly attacks pancreatic beta cells, resulting in an absolute deficiency of insulin secretion. Patients require lifelong insulin therapy. This disease commonly affects children and adolescents, with some patients exhibiting early-onset characteristics. Early-onset type 1 diabetes patients typically develop the disease before the age of 18, experience faster disease progression, and have a higher risk of complications. Type 1 diabetes has a high genetic susceptibility, and genetic factors play a crucial regulatory role in the development and progression of the disease. Conducting genetic risk prediction can enable early screening of high-risk groups, providing a reference for early clinical intervention and disease monitoring.

[0003] Meanwhile, most existing GRS only utilize the statistical association between SNPs obtained from GWAS and diseases, ignoring the expression regulation effect of variants in specific immune cells and lacking cell-specific mechanisms; conventional SNP chips have insufficient HLA typing accuracy, making it difficult to resolve the fine alleles of key genes such as HLA-DRB1 and HLA-DQA1, resulting in biased risk assessment.

[0004] In summary, developing a type 1 diabetes risk prediction scheme that is adapted to the Chinese population and integrates immune cell-specific eQTL regulatory mechanisms and high-resolution HLA typing is of great practical significance for improving the accuracy of disease risk prediction. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of existing traditional genetic risk prediction models, which are mostly based on European and American populations and have insufficient predictive efficacy in Chinese / Asian populations; ignore immune cell-specific regulatory information; have insufficient HLA typing resolution; have insufficient integration of multi-source genetic features; and are prone to overfitting of model weights. This invention provides a method and system for predicting the risk of type 1 diabetes based on eQTL and HLA typing, which can achieve accurate assessment of the risk of developing type 1 diabetes and / or early-onset type 1 diabetes in the test subjects. It is especially suitable for the genetic background of the Chinese population and helps to screen high-risk groups in China.

[0006] To achieve the above-mentioned objectives, the first aspect of this invention provides a method for predicting the risk of type 1 diabetes based on eQTL and HLA typing, comprising the following steps: S1. Obtain biological samples from the cohort population to obtain high-resolution HLA typing data and genomic SNP genotype data of the cohort population; S2. Based on the genomic SNP genotype data and immune cell single-cell transcriptome sequencing data of the cohort population, a set of immune cell-specific eQTL sites is screened. S3. Integrate the set of immune cell-specific eQTL sites, HLA allele characteristics associated with the risk of type 1 diabetes, and SNP characteristics of publicly disclosed clinical phenotypes of type 1 diabetes to obtain an initial set of genetic characteristics. S4. Perform feature screening and weight coefficient optimization on the initial genetic feature set to obtain the final feature set and the risk weight coefficients corresponding to each feature, and construct a genetic risk score (GRS) prediction model for assessing the risk of developing type 1 diabetes. S5. Obtain biological samples from the subject to be tested, and detect the genotype data of the subject to be tested corresponding to the final feature set; S6. Input the genotype data of the subject to be tested into the GRS prediction model to calculate the genetic risk score of the subject to be tested; S7. Based on the genetic risk score, assess the risk of the subject developing type 1 diabetes and / or early-onset type 1 diabetes.

[0007] Furthermore, this invention also specifies the sample source, HLA typing resolution, eQTL analysis type and analysis window, complete eQTL screening sub-process, immune cell type, correction covariate, HLA allele screening method, feature set composition, feature screening algorithm, weight optimization algorithm, and GRS calculation formula.

[0008] A second aspect of the present invention provides a type 1 diabetes risk prediction system based on eQTL and HLA typing, wherein the method comprises: The sample genotype acquisition module is used to acquire biological samples from the cohort population and the subjects to be tested and obtain the corresponding genotype data; The eQTL site screening module is used to screen for a set of immune cell-specific eQTL sites based on genomic SNP genotype data and immune cell single-cell transcriptome sequencing data. The feature integration module is used to integrate multiple types of genetic features to generate an initial set of genetic features. The model building module is used to perform feature selection and weight coefficient optimization on the initial genetic feature set to build a GRS prediction model. The subject risk calculation module is used to input the genotype data of the subjects to be tested into the GRS prediction model and calculate the genetic risk score; The risk outcome output module is used to output the risk assessment results of the subject's risk of developing type 1 diabetes and / or early-onset type 1 diabetes.

[0009] As can be seen from the above description, the present invention does indeed have the following advantages: 1. This invention uses a Chinese population cohort for site screening and model training, and incorporates Chinese population-specific HLA risk alleles to overcome the shortcomings of insufficient predictive efficacy of GRS models from other regions of the world in the Chinese population, making it more suitable for the genetic background of Chinese and East Asian populations; 2. By using single-cell transcriptome data to screen for immune cell-specific eQTL sites, we can capture the regulatory effects of genetic variations on immune cell expression at the functional level, filter out irrelevant tissue-derived noise, and enhance the biological explanatory power of genetic characteristics. 3. Employing four-digit high-resolution HLA allele typing fully explores the detailed genetic risk information of HLA regions and avoids information loss caused by low-resolution typing; 4. It integrates three types of genetic characteristics from different sources, taking into account functional variations, HLA risk alleles, and publicly available clinical GWAS markers, and makes full use of multi-dimensional genetic information; 5. The backward elimination iterative screening algorithm is used to simplify features. The weights are optimized by combining the objective function with regularization term and L-BFGS-B quasi-Newton algorithm. At the same time, weight numerical constraints are set to alleviate the overfitting problem caused by high-dimensional features and take into account the risk assessment of both ordinary type 1 diabetes and early-onset type 1 diabetes. 6. Construct a corresponding prediction system to achieve modularization of the process and facilitate automated analysis of batches of clinical samples. Attached Figure Description

[0010] Figure 1 This is a schematic diagram of a preferred embodiment of the type 1 diabetes risk prediction method of the present invention; Figure 2 In a preferred embodiment, the ROC curves of the model of the present invention and the existing comparative model were found in the queue; Figure 3 The ROC curves of the present invention model and the existing comparative model in the verification queue are shown in the preferred embodiment. Detailed Implementation

[0011] To provide a clearer understanding of the technical features, objectives, and effects of the present invention, specific embodiments of the present invention will now be described with reference to the accompanying drawings.

[0012] It should be understood that the specific embodiments described herein are for illustrative purposes only and are not intended to limit the scope of protection of the invention. In each embodiment of this section, a general overview of the corresponding technical solution is first provided, followed by a detailed explanation of the parameters, calculations, and operational procedures for each step, in conjunction with specific implementation examples. The sample size, software version, threshold, queue division, etc., involved are all exemplary specific implementation examples, and those skilled in the art can make reasonable adjustments based on actual research conditions. Each embodiment concludes with a corresponding explanation, illustrating how this embodiment solves the technical problems raised in the background art and achieves the corresponding technical effects described in the invention, thus realizing the correspondence between technical problems, technical means, and technical effects.

[0013] Implementation Method 1: Basic Process for Type 1 Diabetes Risk Prediction Based on eQTL and HLA Typing This technical solution comprises seven core steps: cohort-level data acquisition, screening of immune cell-specific eQTL loci, integration of multiple genetic features, feature screening and weight optimization to construct a GRS prediction model, genotype detection of test subjects, calculation of genetic risk scores, and risk assessment. First, high-resolution HLA typing and genomic SNP genotypes are obtained using the research cohort, and immune cell-specific eQTLs are screened using single-cell transcriptome data. Then, eQTL loci, risky HLA alleles, and publicly available clinical phenotype SNPs are integrated to form an initial set of genetic features. After feature screening and weight optimization, a GRS prediction model is obtained. Subsequently, test subjects only need to have their genotypes from the corresponding feature set tested and input into the model to obtain a genetic risk score, thereby assessing the subject's risk of developing type 1 diabetes and / or early-onset type 1 diabetes.

[0014] The following is an example of a specific implementation: A research cohort was constructed, comprising individuals with type 1 diabetes and healthy controls. Biological samples were collected from this cohort, and high-resolution HLA typing and genomic SNP genotyping data were obtained. This cohort serves as the data foundation for subsequent eQTL site screening, HLA allele screening, and model training and optimization. As a specific implementation example, the research cohort can be divided into a discovery cohort, a single-cell sequencing subset, a validation cohort, and a flow cytometry validation cohort. The discovery cohort included 568 clinically diagnosed T1D patients under 18 years of age and 1543 healthy controls without a family history of autoimmune diseases or diabetes. The single-cell sequencing subset included 109 T1D patients (diagnosed before age 18, disease duration ≤ 5 years) and 56 age- and sex-matched healthy controls. The validation cohort included 141 T1D patients and 386 healthy controls. The flow cytometry validation cohort included 17 T1D patients and 20 healthy controls. Peripheral blood samples were collected from all the aforementioned Chinese participants.

[0015] Based on genomic SNP genotype data from a cohort population and combined with single-cell transcriptome sequencing data from immune cells, we screened for immune cell-specific eQTL sites to obtain a set of immune cell-specific eQTL sites. Unlike conventional batch tissue eQTL analysis, this protocol distinguishes different immune cell subsets based on single-cell sequencing data and screens for variant sites that play a regulatory role only in specific immune cell subsets. These variants are more likely to be involved in the autoimmune pathogenesis of type 1 diabetes.

[0016] After obtaining the set of immune cell-specific eQTL loci, the set of loci, HLA allele characteristics associated with the risk of type 1 diabetes, and publicly available SNP characteristics of clinical phenotypes of type 1 diabetes are integrated to form an initial set of genetic characteristics, thereby achieving the fusion of multi-source genetic information.

[0017] The initial genetic feature set is subjected to feature screening and weight coefficient optimization. Redundant features and features that do not contribute to the model's predictive performance are removed to obtain a simplified final feature set. At the same time, the risk weight coefficient corresponding to each feature is solved. Based on the above features and weights, a genetic risk score (GRS) prediction model is constructed.

[0018] After the model is built, the risk assessment phase for the test subjects begins. Biological samples are collected from the test subjects to obtain all genotype data for their final characteristic set. This genotype data is then input into the trained GRS prediction model to calculate the genetic risk score for each subject. Finally, based on the obtained genetic risk score, the risk assessment for the subject's development of type 1 diabetes and / or early-onset type 1 diabetes is completed. This protocol can be used for risk assessment of common type 1 diabetes as well as risk prediction for early-onset type 1 diabetes with a lower age of onset. It is suitable for screening the general population and risk assessment of high-risk family members.

[0019] Most existing GRS models are based on European and American population cohorts. Directly applying them to the Chinese population can lead to distorted predictions due to differences in allele frequencies and HLA haplotype distributions. This implementation method, however, uses a case-control cohort within the Chinese population for site selection, feature integration, and model training. This allows it to capture genetic risk variants unique to the Chinese population, particularly high-risk HLA-DR9-related alleles. Furthermore, the methodology incorporates functional genetic information at the immune cell level, rather than solely relying on statistically correlated SNPs. Subsequent feature selection and weight optimization result in risk assessments directly applicable to Chinese subjects, overcoming the poor racial suitability of traditional models and enabling accurate prediction of genetic risk for type 1 diabetes and early-onset type 1 diabetes in the Chinese and East Asian populations.

[0020] Implementation Method 2: Sample and HLA typing parameter limitations This implementation method specifies limitations on the biological sample type, HLA genotyping resolution, and SNP detection and imputation techniques. Peripheral blood genomic DNA samples are used as biological samples; HLA genotyping employs four-digit resolution allelic genotyping; after obtaining SNP genotypes, quality control and genotype imputation are performed to obtain a high-quality SNP dataset suitable for subsequent analysis.

[0021] The following is an example of a specific implementation: For both the cohort population and subsequent test subjects, peripheral blood-derived genomic DNA samples were used as biological samples for testing. Peripheral blood sampling is convenient in clinical settings, is a minimally invasive sampling method, and is easy to use for collecting samples from large cohorts. It does not require obtaining difficult-to-obtain human tissues such as the pancreas or immune organs, making it convenient for clinical application.

[0022] As a specific implementation example, genomic DNA can be purified from peripheral blood samples using the QIAGEN DNeasy Blood & Tissue Kit or QIAamp DNA Blood Mini Kit to obtain high-quality genomic DNA. HLA genotyping was performed using the GenDx NGSgo kit to specifically amplify target regions of the HLA gene. The kit covers the 5'UTR to 3'UTR regions of HLA-A, -B, -C, -DPA1, -DQA1, and -DPB1, as well as the exon 2-exon 4 regions of HLA-DRB1, -DQB1, and -DRB3 / 4 / 5. After constructing the sequencing library, sequencing was performed on the Illumina MiSeq platform, and the sequencing results were analyzed using NGSengine software to obtain four-digit resolution HLA alleles. Rare variants with minor allele frequencies (MAF) < 0.01 were filtered out to obtain the HLA genotyping dataset usable for subsequent analysis.

[0023] HLA genotyping uses four-digit resolution HLA allele typing. For example, HLA-DRB1 03:01, HLA - DRB1 04:05 belongs to the four-digit resolution allele; the two-digit genotype can only identify... 03、 The four main categories do not differentiate between subtypes. Within the same main category, different four-digit subtypes show significant differences in the risk of type 1 diabetes, for example, in the high-prevalence subtypes in the Chinese population. This gene belongs to a high-risk allele, while some subtypes within the same group have very low risk levels. Using low-resolution genotyping would confuse high-risk and low-risk subtypes, introducing assessment bias. Using four-digit high-resolution genotyping preserves detailed risk information within the HLA region, ensuring the accuracy of subsequent models.

[0024] As a specific implementation example, SNP genotyping was performed using the Illumina AsianScreeningArray chip for genotyping. This chip is specifically designed for genetic polymorphism in Asian populations and is more suitable for the Chinese population. The original idat intensity file was analyzed using IAAP software. The quality control conditions were: sample SNP detection rate > 97%, and invalid allele loci, multiple allele loci, single morphologies, and insertion / deletion variants were removed. The original chip yielded 686,592 high-quality SNPs, which were converted to hg38 genomic coordinates. Genotyping: For non-HLA regions, the ChinaMAP Chinese population reference panel was used, with phase reconstruction performed using Eagle2 and genotyping performed using Minimac4. For HLA regions, SNP2HLA was used in conjunction with the pan-Asian reference panel for genotyping. Post-genotyping quality control: r² < 0.8, MAF < 0.05, Hardy-Weinberg equilibrium P < 1 × 10⁻⁶. -6 The sites were removed; ultimately, 4,646,233 high-quality SNPs were obtained for eQTL analysis.

[0025] Existing techniques often use low-resolution two-digit HLA typing and frequently employ European population reference panels for genotyping, resulting in the loss of risk information for HLA subtypes specific to the Chinese population and making the imputation results prone to bias. This implementation method, on the one hand, uses high-resolution four-digit HLA typing, which can accurately distinguish key risk alleles in the Chinese population, such as... On the other hand, SNP chips designed for Asian populations and the ChinaMAP reference panel for the Chinese population were used to complete genotyping, improving the accuracy of genotyping detection and completion of Chinese population samples. This reduced systematic errors from the upstream sample and data generation stages, laying a solid data foundation for the subsequent construction of models adapted to the Chinese population.

[0026] Implementation Method 3: eQTL Types and Analysis Window Parameters This embodiment defines the types of immune cell-specific eQTLs and the parameters of the cis-eQTL analysis window. This invention selects cis-eQTL loci for subsequent analysis, and the cis-eQTL analysis window is set to ±100kb upstream and downstream of the target gene. This window size can be adjusted within a reasonable range based on the cohort size.

[0027] The following is an example of a specific implementation: The immune cell-specific eQTL sites identified in this study are cis-eQTL sites. Cis-eQTLs represent the regulation of expression of nearby target genes on chromosomes by SNPs; trans-eQTL variants are farther from the target genes, have more complex regulatory relationships, and have a higher false positive rate. Therefore, this invention prioritizes cis-eQTLs for subsequent analysis.

[0028] Traditional GRS models directly use GWAS to statistically analyze significant SNPs without distinguishing between cis and trans eQTLs, which introduces a large number of false positives and non-pathogenic genetic variations. This implementation method limits the use to cis eQTLs only and sets a reasonable analysis window to prioritize capturing real expression regulatory variations in genomic neighborhoods and filter out interference from trans eQTLs with high false positive rates. Based on cis eQTL results obtained from a Chinese population cohort, this method reflects the real regulatory patterns within immune cells in the Chinese population, distinguishing it from eQTL datasets from public databases in Europe and the United States, and further improving the effectiveness of the feature set in the Chinese population.

[0029] Implementation Method 4: Complete Screening Sub-process for Immune Cell-Specific eQTL Sites This implementation describes a complete screening process for a set of immune cell-specific eQTL sites, which consists of five sub-steps: single-cell transcriptome preprocessing and annotation, pseudo-batch matrix construction, cis-eQTL association analysis, conditional eQTL delinking, and cross-matching with GWAS datasets. Through multi-step processing, a set of immune cell-specific eQTL sites associated with the risk of type 1 diabetes is obtained.

[0030] The following is an example of a specific implementation: S2-1: Perform quality control, standardization, dimensionality reduction and clustering on the raw single-cell transcriptome sequencing data to obtain annotations for various immune cell types.

[0031] As a specific implementation example, the single-cell sequencing subset samples (109 T1D cases and 56 healthy controls, all from the Chinese population) were used for PBMC library construction using 10X Chromium and sequencing using Novaseq 6000; the original fastq was aligned to the hg38 reference genome using cellrangerv7.1.0; souporcell v2.0 was used to remove sample cross-contamination; Scanpy v1.9.3 was used for data analysis; and DoubletDetection v4.2 was used to filter double-cell events; ultimately, 1,149,633 high-quality single-cell transcriptomes were obtained.

[0032] The quality control process employed clear filtering criteria: cells with fewer than 200 or more than 6000 gene counts were removed; cells with mitochondrial genes accounting for more than 10% of total expression were removed; and low-quality cells with an excessively low total number of sequencing reads were removed. After cell filtering, gene expression was normalized to a total count of 1,000,000 and transformed using log2(x+1); 1514 interfering genes related to mitochondria, heat shock proteins, ribosomes, and cell dissociation were removed; the hypervariable gene HVG was identified; principal component analysis (PCA) was performed; batch effects were corrected using Harmony v0.0.6; UMAP dimensionality reduction was performed based on the corrected principal components; a nearest neighbor graph was constructed using the neighbors function (n_neighbors=35); and Leiden clustering (resolution=1.3) was used to complete the clustering.

[0033] Customized clustering parameters were used for different immune cell lineages: the top 20 principal components were used for T / NK cells, n_neighbors=25, resolution=1.3; the top 20 principal components were used for B cells, n_neighbors=30, resolution=1.3; and the top 20 principal components were used for monocytes / dendritic cells, n_neighbors=30, resolution=1.3.

[0034] Cell annotation was performed based on classical marker genes: CD3E, CD40LG, CD8A, NKG7, CD79A, and MZB1 were used for lymphoid lineages; CD14, FCGR3A, CD1C, and LILRA4 were used for myeloid lineages; and 10 types of immune cells were finally identified: CD14+ monocytes, CD16+ monocytes, conventional dendritic cells (cDC), plasmacytoid dendritic cells (pDC), CD4+ T cells, CD8+ T cells, proliferating T cells, NK cells, B cells, and plasma cells.

[0035] S2-2: Generate a pseudo-batch gene expression matrix for each type of immune cell.

[0036] Pseudo-batch matrix construction: Gene expression counts of all cells belonging to the same immune cell type are summed and merged; the matrix rows are genes, and the columns are sample numbers of the research cohort; one pseudo-batch expression matrix corresponds to one immune cell subpopulation; if a sample has no cells in a certain cell subpopulation, a deletion marker is added to the corresponding position in the matrix, and the sample is excluded in subsequent analysis. For each cell type, genes with an expression percentage of less than 10% in the cohort are filtered out; the expression values ​​of the remaining genes are logarithmically transformed using log2(x+1).

[0037] S2-3: For each immune cell type, a linear regression model was used to conduct cis-eQTL association analysis. Covariates were included to correct for confounding factors. The false discovery rate (FDR) was used as the significance screening threshold to obtain preliminary eQTL association results.

[0038] This step uses the MatrixEQTL v2.3 R package to perform cis-eQTL calculations.

[0039] For separate modeling of T1D patients and healthy controls, the complete model formula is as follows: ;

[0040] In the formula: : Expression level of gene g in sample i after log2(x+1) transformation; : Genotype encoding of sample i locus s (0, 1, 2); The first three principal components of the genotype; Individual age; Individual gender; Duration of illness; Indicator function: T1D patients = 1, healthy controls = 0; Sequencing batches; Residual.

[0041] For each target gene, all SNPs within a ±100kb window upstream and downstream of the target gene are extracted, and SNP-gene pairing regression is performed one by one. The regression coefficient, standard error, and p-value of each pairing are output.

[0042] Benjamini-Hochberg multiple test correction is performed on the set of p-values ​​obtained from all tests. The FDR calculation formula is as follows: ; m represents the total number of SNP-gene tests, and i is the sequence number after sorting the P-values ​​from smallest to largest; the FDR threshold is preset to 0.05, and associations with FDR < 0.05 are retained to obtain preliminary eQTL association results.

[0043] This specific implementation example identified 21,276 eGenes; The following was used:

[0044] Cell specificity was calculated, with 49.15% of eGene in the T1D group and 56.27% in the control group.

[0045] S2-4: Conditional eQTL analysis was used to remove linkage disequilibrium interference, and a maximum of 5 independent eSNPs were retained for each target gene.

[0046] Conditional eQTLs perform cyclic processing, processing each target gene individually in a loop: 1. Read the set of SNPs for all preliminarily significant eQTLs of the target gene; 2. Select the SNP with the smallest P-value as the primary condition SNP; 3. Using the genotype of the primary condition SNP as a new covariate, refit the aforementioned linear regression model, and re-perform association tests on all other SNPs under this gene window; 4. SNPs that still achieve FDR < 0.05 after retesting are retained as candidate secondary eSNPs; 5. Select the SNP with the smallest p-value from the candidate secondary eSNPs, add conditional covariates, and perform another cyclic regression test; 6. Iterate until no new SNPs can meet the FDR significance threshold, then terminate the cycle for the target gene; 7. For each target gene, a maximum of 5 independent eSNPs are retained, which constitutes the set of mutually independent eSNPs for that target gene. Through the above iterative approach, false positive association signals caused by genomic linkage disequilibrium are eliminated.

[0047] S2-5: Cross-matching was performed using publicly available GWAS risk variant datasets for type 1 diabetes to screen for a set of immune cell-specific eQTL loci associated with type 1 diabetes risk. T1D-GWAS dataset screening criteria: genome-wide significance P < 5 × 10⁻⁶. -8 MAF ≥ 0.01.

[0048] Matching rules: Cross-matching is performed based on the SNP's chromosome number, genomic physical location, and SNP identifier. Only independent eSNPs that appear simultaneously in the GWAS risk variant dataset and have passed the aforementioned eQTL screening are retained for inclusion in the immune cell-specific eQTL locus set; all other loci are discarded. In this specific implementation example, cross-matching yielded 120 risk eGenes, which were further screened to obtain 100 risk eSNPs.

[0049] For example, the 120 risk eGenes obtained from cross-matching are: A2M、A2M-AS1、AC002378.1、AC004832.5、AC007100.1、AC007216.2、AC007381.1、AC009163.7、AC010175.1、AC018521.5、AC025162.1、AC034102.3、AC068254.1、AC092821.3、AC093388.1、AC099508.2、AC108105.1、AC133065.1、AC145285.3、AGXT2、AL158136.1、AL356417.1、ALOX12-AS1、AMPD1、AP003774.2、AP005482.1、BATF3、BCAR1、C11orf21、CABP7、CACNG2、CCDC63、CCR3、CD226、CEP76、CLEC2D、CLECL1、CTRB1、CTSH、CUX2、DDX43、DDX59-AS1、DGKQ、DOK6、FHL3、FIGNL1、FLRT1、FLVCR1-DT、GCA、GCAT、GPR18、GSDMA、GSDMB、HLA-A、HLA-B、HLA-C、HLA-DMA、HLA-DMB、HLA-DOA、HLA-DOB、HLA-DPA1、HLA-DPB1、HLA-DQA1、HLA-DQA2、HLA-DQB1、HLA-DQB1-AS1、HLA-DQB2、HLA-DRB1、HLA-DRB5、HLA-F、HLA-F-AS1、HLA-G、HORMAD2、HSCB、IL15RA、INAFM1、INO80B、ITGB7、JAZF1-AS1、KCNQ1DN、KLRG1、LINC00612、LINC00706、LINC00707、LINC00987、LOXL3、M6PR、MAS1、MEIS3、MTFP1、NOTCH2、NPIPB6、NSL1、OLFML3、ORMDL3、PFKFB3、PHETA1、PRKCQ、PRKCQ-AS1、PZP、RGS1、RMI2、RPS26、SIRPB1、SKAP2、SLC22A18AS、SLX1B、STIP1、SUOX、SYT6、TAFA3、TCN2、TMEM116、TMEM170A、TMPRSS3、TNP2、TRIOBP、TUFM、UBASH3A、ZPBP2。

[0050] Example of a small computational case: A total of 220 SNPs under a certain target gene were tested. After BH-FDR correction, 12 significant eQTLs were obtained. After multiple rounds of conditional eQTL cycle analysis, non-independent signals caused by linkage were removed, and finally 3 independent eSNPs were retained. Then, compared with the T1DM-GWAS dataset, 2 sites appeared in the GWAS risk variant set. These 2 sites were then included in the set of immune cell-specific eQTL sites.

[0051] Most existing technologies directly reference bulk eQTLs from public databases in Europe and the United States, which cannot reflect the true expression regulation within immune cells in the Chinese population, nor can they distinguish the specific effects of different immune cell subsets. This implementation method, based on the single-cell sequencing subset of the Chinese population, completes a full set of single-cell quality control, pseudo-batch construction, group modeling, and conditional eQTL delinking to screen for Chinese population-specific eSNPs for immune cells. By separately modeling the patient and control groups, it can capture expression regulatory signals specific to the disease state. Conditional eQTL removes linkage disequilibrium, avoiding redundant noise. Cross-matching with GWAS risk variants ensures the association between the loci and the risk of type 1 diabetes. The resulting set of loci possesses a biological mechanism explanation, distinct from purely statistically associated SNPs, providing high-quality functional genetic features for subsequent construction of high-performance GRS models adapted to the Chinese population.

[0052] Implementation Method 5: Covariate Setting for Immune Cell Type and eQTL Analysis This implementation specifies the types of immune cells involved in eQTL analysis and the covariates used in the regression analysis. Multiple immune cells involved in the pathological process of type 1 diabetes are selected, and corresponding confounding covariates are included in the regression analysis to reduce false positive associations.

[0053] The following is an example of a specific implementation: The immune cells used in the analysis were selected from at least one of the following: CD14+ monocytes, CD16+ monocytes, conventional dendritic cells, plasmacytoid dendritic cells, CD4+ T cells, CD8+ T cells, proliferating T cells, NK cells, B cells, and plasma cells. All of these immune cells are directly involved in the autoimmune damage process in type 1 diabetes. Monocytes and dendritic cells are responsible for antigen presentation, T cells mediate pancreatic β-cell killing, B cells and plasma cells participate in autoantibody production, and NK cells participate in immune regulation. Prioritizing the analysis of all of these cell subpopulations ensures the full capture of eQTL regulatory signals in different immune subpopulations. As a specific implementation example, in scenarios with limited sample size, a selection of specific cell types can also be used for screening.

[0054] The covariates included in the cis-eQTL association analysis were genotype principal component, age, sex, and experimental batch. For type 1 diabetes patients within the cohort, disease duration was additionally included as a covariate. Genotype principal component analysis corrected for population stratification; age and sex eliminated individual differences in basic biological characteristics; experimental batch analysis eliminated batch effects from sequencing experiments; and disease duration was used to correct for the impact of patient disease course on immune cell expression profiles, reducing false positive associations caused by confounding factors.

[0055] Incomplete selection of immune cell types and missing covariate settings can easily introduce confounding factors such as population stratification, age and sex, sequencing batch, and disease course, resulting in a large number of false-positive eQTL signals. This implementation method covers all key immune cell types in the pathogenesis of type 1 diabetes and fully sets a complete set of covariates adapted to the Chinese population cohort; the genotype principal component can correct for population stratification caused by different regional substructures within the Chinese population; the additional inclusion of the patient's disease duration corrects for the interference of disease course on immune expression; thereby reducing false positives caused by confounding factors, ensuring the authenticity and reliability of the selected eQTL loci, and further improving the predictive stability of the subsequent model in the Chinese population.

[0056] Implementation Method 6: Screening Methods Based on HLA Allele Characteristics This implementation describes the screening techniques for HLA alleles associated with the risk of type 1 diabetes. A logistic regression model adjusted for sex and geographic confounding factors was used, combined with maximum likelihood estimation, Newton-Raphson iteration, and multiple positive control to screen for significant risk and protective HLA alleles.

[0057] The following is an example of a specific implementation: After obtaining four-digit HLA typing data from the cohort population (568 T1D cases and 1543 healthy controls, all of whom were Chinese), a logistic regression model was used to adjust for sex and geographic confounding factors, with case-control status as the dependent variable and each HLA allele as the independent variable.

[0058]

[0059] In the formula: T1D=1 represents a patient with T1D, and T1D=0 represents a healthy control. The value is the HLA allele dose, ranging from 0 to 2; gender is the sex; and region is the geographical region (southern / northern).

[0060] The regression coefficients are solved using maximum likelihood estimation (MLE). Likelihood function: ; ; Log-likelihood: ; Parameter iteration updates use Newton-Raphson iteration: ; The first derivative score function; The second derivative matrix is ​​the Hessian matrix; the change in the parameters of the iterative convergence criterion is less than... .

[0061] After iterative convergence, a Wald test is performed to determine the significance of the coefficients; multiple tests are performed using Bonferroni correction, with a significance threshold set. .

[0062] Odds ratio (OR) and 95% confidence interval calculation: ; ; Judgment rule: Simultaneously satisfying and It was determined to be a risk allele; it meets the following criteria. and The alleles were identified as protective. In this specific implementation case, 35 significant HLA allele characteristics were identified, including those specific to the Chinese population. High-risk alleles.

[0063] In this specific implementation example, the 35 HLA alleles significantly associated with T1D obtained by logistic regression screening and their data are shown in Table 1.

[0064] Table 1: Summary table of associations between HLA alleles and type 1 diabetes.

[0065]

[0066] Example illustration: For a specific four-digit HLA allele, regression yielded P = 2.3 × 10⁻⁶. -5 less than the threshold If the allele is positive, it is retained as an HLA risk trait; if P=0.02 is greater than the threshold, it is removed.

[0067] Existing technologies directly use HLA risk alleles reported in European and American literature, which are not suitable for the Chinese population and do not correct for the regional demographic structure between northern and southern China. This implementation method directly performs logistic regression analysis within a case-control cohort of the Chinese population, specifically correcting for sex and regional confounding factors within the Chinese population. This identifies the true risk and protective HLA alleles in the Chinese population, rather than directly applying HLA risk lists from European and American populations. Through rigorous multiple testing to control false positives, 35 HLA characteristics containing risk subtypes specific to China were selected, providing the model with HLA genetic risk information suitable for the Chinese population.

[0068] Implementation Method 7: Characteristic Composition of the Initial Genetic Characteristic Set This implementation describes the composition and integration process of the initial genetic feature set. The initial genetic feature set consists of three types of features: immune cell-specific risk eSNP features, type 1 diabetes risk-related HLA allele features, and type 1 diabetes clinical phenotype SNP features from publicly available literature. Deduplication is performed during the integration process.

[0069] The following is an example of a specific implementation: The initial set of genetic features obtained after feature integration contains genetic features from three different sources.

[0070] ; : 35 HLA allele characteristics; : 100 immune cell-specific risk eSNP features; : 81 clinical phenotype-related SNP features of type 1 diabetes reported in publicly available literature.

[0071] Screening rules for 81 disease classification SNPs: taken from publicly available literature related to T1D subtypes and clinical phenotypes; screening criteria association P < 0.05, effect beta has biological significance; original effect coefficients and standard errors were extracted; this set of features is used to assist in the assessment of clinical phenotypes such as T1D subtypes and age of onset.

[0072] For example, the 81 SNPs (species variants) representing T1D disease classifications derived from publicly available literature are as follows: 10:6087680:C:T、2:203866366:G:A、22:30108663:A:G、AA_DRB1_11_32660115_D、6:33081532:A:G、12:111569952:C:T、AA_DRB1_11_32660115_P、6:90248512:C:G、4:122408207:G:A、16:75213347:T:G、AA_A_9_30018537_S、9:4291747:C:G、6:32427740:T:C、22:37195278:G:A、2:162267541:C:T、12:24428822:T:C、19:46705224:T:C、6:31439802:G:A、22:30185733:C:T、15:78944951:C:T、AA_A_9_30018537_T、16:11086016:A:G、AA_C_275_31345841_E、13:99429512:C:T、AA_C_275_31345841_K、12:56088396:T:G、10:88291560:A:G、AA_DRB1_11_32660115_L、10:88275897:C:A、AA_DRB1_74_32659926_E、6:159044945:T:C、10:6055320:T:C、2:203874196:G:A、AA_DRB1_74_32659926_L、2:100208905:C:T、19:10381598:C:T、AA_DRB1_11_32660115_V、18:69864406:T:C、7:50959497:A:C、1:206766559:G:A、11:2159830:T:G、14:98032614:A:G、AA_C_275_31345841_G、AA_A_9_30018537_F、11:2160994:A:T、15:38554821:C:T、12:55991020:G:A、10:6080046:T:A、4:122456825:A:C、6:29972123:C:T、6:31376405:T:., 6:32658495:C:A, AA_DRB1_11_32660115_S, 6:29868825:C:T, AA_DRB1_74 _32659926_A, AA_DRB1_74_32659926_R, 9:4290823:A:T, 16:28528527:C:T , AA_A_9_30018537_Y, 21:42421219:C:A, 15:78943104:T:C, 14:100840110 :T:C, 6:90292775:A:G, AA_DRB1_11_32660115_G, 17:40614034:T:C, 1:113 834946:A:G、15:38610275:A:C、6:32479411:T:A、2:110857502:T:G、18:69 859408:T:C、2:162254026:A:G、12:56080696:A:G、AA_DRB1_74_32659926_ Q、20:1629905:T:C、21:42405613:G:C、14:68796882:T:C、12:9733403:G:A , 22:37191071:A:C, 18:12809341:A:G, 16:28580209:G:A, 10:8060309:G:A. .

[0073] The complete feature integration process is as follows: First, collect all immune cell-specific risk eSNPs. For the same SNP that appears repeatedly in different immune cell subsets, retain only one feature record. Second, load the screened HLA risk allele feature set. Third, load the clinical phenotype SNP feature set of type 1 diabetes from publicly available literature. Fourth, deduplicate all feature sets; if the same variant belongs to multiple categories, retain only one feature. Merge all deduplicated features to obtain the initial genetic feature set, which contains 216 features in this specific implementation example. These three types of features provide genetic information from three dimensions: immune function regulation, major HLA risk regions, and publicly available disease-associated markers, complementing each other and improving the overall information capacity of the model.

[0074] Traditional GRS models mostly use single-source features, either only GWAS-SNPs or only HLA, losing multi-dimensional genetic information. This implementation integrates three types of features: eQTL features from immune cells derived from the Chinese population, HLA risk alleles from the Chinese population screened in this cohort, and SNPs related to clinical phenotypes. This integrates both molecular functional variations and population-specific HLA risks, while also incorporating markers related to disease subtypes and early onset. Feature deduplication eliminates redundant sites. The initial feature set takes into account both disease risk and early-onset clinical phenotypes, providing comprehensive feature material for subsequent training of a GRS model adapted to the Chinese population and capable of assessing early-onset risk.

[0075] Implementation Method 8: Backward Elimination Iterative Feature Screening This implementation describes how to perform feature selection on the initial genetic feature set. A backward elimination iterative selection algorithm is used, with the AUC obtained from five-fold cross-validation as the evaluation metric, to iteratively eliminate redundant features, resulting in a streamlined final feature subset.

[0076] The following is an example of a specific implementation: The backward elimination iterative screening algorithm is used to screen features, and the AUC value obtained by K-fold cross-validation is used as the evaluation index for feature selection. As a specific implementation example, this scheme prefers 5-fold cross-validation.

[0077] set up This is the feature set after the (i-1)th iteration; for Five-fold cross-validation average AUC; The average AUC of the feature set after removing features is obtained through five-fold cross-validation.

[0078] initialization: ; ; Iterative loop: for i=1 to N: for each x j ∈Sᵢ₋1: (Remove feature x) j ) (50% cross-validation AUC) Find (Select the feature that improves AUC the most)

[0079] (Remove this feature) (Updated AUC) else break.

[0080] The 5-fold cross-validation process is as follows: The cohort case-control dataset is randomly and non-repeatingly divided into 5 subsets. In each round, 4 subsets are used as the training subset, and the remaining subset is used as the independent test subset. The training subset is used to train the model corresponding to the current feature subset, and the test subset is used to calculate the ROC curve (Receiver Operating Characteristic) and AUC. The process iterates through all 5 partitioning methods, taking the arithmetic mean of the 5 test set AUCs as the performance evaluation metric for the current feature subset. The ROC (Receiver Operating Characteristic) curve represents the receiver operating characteristic curve. The Y-axis (vertical axis) represents the true positive rate, the proportion of people who are actually infected but are correctly predicted by the model; the X-axis (horizontal axis) represents the false positive rate, the proportion of people who are not actually infected but are misjudged as high-risk by the model; the dashed line in the graph represents a baseline for random guessing. If the model coincides with this dashed line, it means the prediction is completely ineffective. AUC, or the area under the ROC curve, ranges from 0 to 1. AUC closer to 1 indicates a stronger model's discriminative ability; an AUC of 0.5 is equivalent to random guessing and has no predictive value. (See reference for more information.) Figure 2 and Figure 3 As shown.

[0081] Iteration termination condition: When removing any feature fails to improve the average AUC of cross-validation, the iteration stops; output a simplified feature subset, which in this specific implementation example yields 190 final features after iteration.

[0082] This screening method can efficiently remove redundant and noisy features, reduce model dimensionality, and alleviate overfitting.

[0083] The initial feature set has high dimensionality, contains some redundancy, and features that offer no predictive gain for the Chinese population. If all of these features are directly used in modeling, overfitting is highly likely, leading to performance degradation on independent samples, especially external validation samples from China. This implementation uses a backward elimination iterative selection method, relying on the AUC of the five-fold cross-validation in this cohort as an objective evaluation metric, to iteratively eliminate features with no predictive gain, resulting in a streamlined final feature subset of 190. By reducing the feature dimensionality while retaining genetic markers that are discriminative against the Chinese population, the risk of overfitting caused by high dimensionality is mitigated, ensuring that the subsequent model has good generalization ability on independent Chinese validation cohorts.

[0084] Implementation Method 9: Weight Coefficient Optimization and GRS Calculation Formula This implementation method introduces the weighting coefficient optimization method, constraints, GRS score calculation formula, and risk level classification. A constrained L-BFGS-B quasi-Newton algorithm is used to solve for the optimized weights, and risk levels are classified based on the distribution of the cohort samples to calculate the risk scores for disease incidence and early onset.

[0085] The following is an example of a specific implementation: After completing feature filtering and obtaining a simplified feature subset, coefficient optimization is performed.

[0086] The first step is to screen the features to be optimized: from the 190 final features, features with a relevance score not lower than the median are selected to participate in the weight optimization.

[0087] ; This is the original effect coefficient. This is the p-value for the chi-square test.

[0088] Construct the optimization objective function: ; The regularization parameter λ = 0.001. These are the optimized weighting coefficients.

[0089] Set coefficient boundary constraints: ; Queue partitioning: As a specific implementation example, it was found that the queue was randomly divided into 80% training set and 20% validation set.

[0090] The L-BFGS-B (quasi-Newton algorithm with boundary constraints) method is used for iterative solution; initialization... ; Calculate the objective function and its gradient: ; The algorithm automatically determines the step size α and updates the coefficients, and checks the boundary constraints in each round; convergence condition: the objective function changes by less than 10. -6 The iteration terminates, and the optimized weight coefficients of all 190 features are output. .

[0091] In this specific implementation example, after feature simplification and weight optimization calculation, the original effect coefficients and optimized risk weight coefficients of 190 features used to construct the genetic risk scoring model were obtained, and the data are shown in Table 2.

[0092] Table 2: Characteristic and original effect coefficients and optimized weight coefficients in the genetic risk scoring formula.

[0093]

[0094]

[0095]

[0096]

[0097]

[0098]

[0099]

[0100]

[0101] The GRS prediction model of this invention is a weighted summation model, and the genetic risk score is calculated according to the following formula: ; In the formula, n represents the total number of features within the final feature set; This represents the risk weight coefficient obtained after optimizing the j-th feature; This represents the genotype code value for the j-th genetic trait locus of subject i. Genotype coding strictly follows the additive inheritance model: the reference homozygous genotype is coded as 0, the heterozygous genotype as 1, and the risk allele homozygous as 2.

[0102] Example simulation calculation: Assume the final feature set has 3 features with weights of 0.32, 0.75, and 0.41 respectively; and a subject's genotype codes are 1, 2, and 0 respectively. ; Furthermore, risk levels were determined based on the quantile distribution of GRS scores from the entire Chinese population sample in the training cohort: the 20th and 80th quantiles of the cohort were used as cutoff thresholds; a GRS score below the 20th quantile was considered low risk; the 20%-80th quantile range was considered medium risk; and a score above the 80th quantile was considered high risk. A higher GRS score indicates a higher risk of genetic disease in the subject. This scoring system can also be used for risk stratification of early-onset type 1 diabetes.

[0103] Model evaluation metrics: ; T represents the total number of samples in the queue; as a specific implementation example, it was found that queue T = 163 (107 T1D + 56 HC); independent validation queue T = 527 (141 T1D + 386 HC).

[0104] Simultaneously assess the Spearman correlation between GRS score and age at diagnosis of T1D; this specific implementation case yielded... , The higher the GRS, the earlier the patient tends to be diagnosed, confirming that the model also has the ability to predict early-onset T1D.

[0105] Traditional GRS directly adopts the original effect coefficients from literature, most of which are derived from GWAS in European and American populations. Directly applying these coefficients can lead to weight mismatches in the Chinese population. This implementation method does not directly use publicly available original effects. Instead, it re-optimizes the weights based on Chinese population cohort data using the L-BFGS-B algorithm with L2 regularization and hard constraints on upper and lower weight limits. The regularization term suppresses excessive weight amplification, and the weight constraints prevent the optimization results from completely deviating from biological meaning. The risk grading threshold is calculated based on the GRS distribution of the Chinese population cohort in this invention, rather than simply adopting foreign thresholds. The model can not only assess disease risk but also reflect the tendency for early onset of disease. The resulting GRS score and risk grading are more closely aligned with the genetic background of the Chinese population, while also mitigating model overfitting.

[0106] Implementation Method 10: Risk Prediction System Example This embodiment describes a prediction system for executing all the aforementioned methodologies. The system consists of multiple functional modules, which work together to complete the entire process from sample data acquisition, site selection, feature integration, model construction, risk calculation, and result output.

[0107] The following is an example of a specific implementation: A type 1 diabetes risk prediction system based on eQTL and HLA typing includes a sample genotype acquisition module, an eQTL locus screening module, a feature integration module, a model construction module, a subject risk calculation module, and a risk result output module.

[0108] The sample genotype acquisition module is used to acquire biological samples from the cohort population and the subjects to be tested and obtain the corresponding genotype data, and to receive HLA typing results and SNP genotype detection results. The eQTL site screening module is used to screen for immune cell-specific eQTL site sets based on genomic SNP genotype data and immune cell single-cell transcriptome sequencing data. The module integrates all the analysis logic of quality control filtering, pseudo-batch matrix generation, eQTL linear regression association analysis, BH-FDR multiple testing, multi-round cyclic conditional eQTL analysis, and GWAS cross-matching. The feature integration module is used to integrate immune cell-specific eQTL sites, HLA allele features, and publicly available SNP features, complete feature deduplication, and generate an initial set of genetic features. The model building module is used to perform backward elimination iterative feature screening on the initial genetic feature set, weight optimization of the L-BFGS-B constrained regular objective function, complete the construction of the GRS prediction model, and store the final feature set and risk weights of each feature, as well as the risk quantile threshold calculated based on the Chinese cohort. The subject risk calculation module is used to read the genotype data of the subjects to be tested, call the pre-trained GRS prediction model, and substitute it into the weighted summation formula to calculate the genetic risk score; The risk outcome output module is used to output the risk assessment results of the subject's risk of developing type 1 diabetes and / or early-onset type 1 diabetes. The output content may include the original GRS score, risk level, and risk warning. The output format may be at least one of text, table, and visualization report.

[0109] As a concrete example, the system can be deployed on a local server or packaged as a web-based tool to achieve automated analysis of batch samples.

[0110] Existing technologies lack a complete and modular predictive software system tailored to the Chinese population. Researchers need to piece together multiple third-party software programs to complete the analysis, resulting in a high operational threshold and making it difficult to use for large-scale clinical sample screening. The system implemented in this way solidifies all modules of the analysis process adapted to the Chinese population; the model building module stores features, weights, and risk stratification thresholds obtained from training on the Chinese cohort; it only requires the genotype of the subject to be tested to directly output a risk report for the Chinese population; it lowers the threshold for clinical use and is suitable for large-scale genetic risk screening of type 1 diabetes in domestic hospitals and testing institutions.

[0111] Implementation Method 11: Model Performance Comparison Example This implementation method introduces a model performance verification and evaluation scheme, setting up multiple sets of comparison models and conducting ROC-AUC evaluations in the discovery queue and independent validation queue, combined with... Figure 2 , Figure 3 The predictive performance of different models was compared to verify the technical advantages of this solution. The cohort sample size and AUC value are specific implementation examples. ROC, or Receiver Operating Characteristic curve, is used to evaluate the discriminative power of binary classification prediction models; AUC is the area under the ROC curve, and the closer the AUC value is to 1, the stronger the model's ability to distinguish between diseased and healthy samples.

[0112] The following is a specific implementation example: Multiple sets of comparative GRS models were set up, and ROC-AUC performance evaluations were conducted in the discovery cohort and independent validation cohort (both containing Chinese population samples) to illustrate the predictive performance advantages of the present invention's technical solution compared to existing technologies. Performance results are available in [link to relevant documentation]. Figure 2 as well as Figure 3 . Figure 2 and Figure 3 In the diagram, the outermost curve at the top is the ROC curve of this study, and its AUC value is higher than that of the others.

[0113] The comparison model is set as follows: Comparison Model 1: The traditional GRS model only uses SNP loci from publicly available GWAS sources and directly uses the original GWAS effect values ​​as weights (mainly based on European and American population data). Comparison Model 2: GRS model constructed using only HLA allele features; The model of this invention is a GRS model obtained by fusing immune cell-specific eQTL sites, high-resolution HLA alleles, and publicly known clinical phenotype SNPs, and by backward elimination feature screening and constrained L-BFGS-B weight optimization.

[0114] Figure 2 To compare the ROC curves of the cohort (107 cases of T1D and 56 healthy controls, Chinese population); Figure 3 This is a comparison of ROC curves for an independent validation cohort (141 cases of T1D and 386 healthy controls, Chinese population). The independent validation cohort was not involved in model training and was used to objectively evaluate the model's generalization ability.

[0115] ROC generation principle: 1. Calculate the GRS score for each sample; sort all samples from lowest to highest GRS score. 2. Iterate through all samples and take the GRS value as the judgment threshold T; ; ; 3. With FPR as the horizontal axis and TPR as the vertical axis, connect all the points corresponding to the thresholds (FPR, TPR) to obtain the ROC curve; 4. AUC is the area under the ROC curve: ; An AUC of 0.5 represents random guessing, while an AUC of 1 represents perfect discrimination.

[0116] Combination Figure 2 It can be seen that: in the discovery queue, the ROC curve of the GRS model of this invention is located above the curves of the other two sets of comparative models; the AUC of the model of this invention is 0.92, 95% CI [0.87-0.96]; the highest AUC of the comparative model is 0.87; the confidence interval does not overlap with the highest AUC of the existing model.

[0117] Combination Figure 3As can be seen, in the independent validation cohort, the AUC of the GRS model of this invention is 0.89, with a 95% CI [0.85-0.92], which is still higher than the comparison model; this proves that the model of this invention still has better discrimination ability on independent Chinese population samples and has not experienced serious overfitting.

[0118] Simultaneously, in the evaluation of the early-onset type 1 diabetes subgroup, the model of this invention also showed superior AUC performance compared to the control model; the GRS was correlated with the age of diagnosis by Spearman (r = -0.34, P = 2.91 × 10⁻⁶). -4 This confirms its predictive value for early-onset T1D.

[0119] The highest AUC among GRS samples from Europe and America in a Chinese independent cohort is only 0.87. (This is achieved through...) Figure 2 , Figure 3 As can be seen from the comparison, the entire scheme of this invention achieves an AUC of 0.92 in the discovery cohort and maintains an AUC of 0.89 in the independent Chinese validation cohort, which is significantly better than the comparison model. At the same time, it can effectively predict the risk of early-onset type 1 diabetes. The experimental results directly verify that this invention solves the technical defects of insufficient predictive efficacy of traditional models in the Chinese population, and confirms the actual performance gain brought about by the integration of immune cell eQTL, Chinese population-specific HLA features, and feature screening-weight optimization techniques.

[0120] The present invention will now be described with reference to a more comprehensive specific embodiment. This comprehensive embodiment fully demonstrates a complete project implementation process. The number of queue samples, reagents, and software parameters described below are all specific implementation examples, and those skilled in the art can make adaptive adjustments according to actual circumstances.

[0121] The study collected data from 568 patients with type 1 diabetes and 1200 healthy controls in the discovery cohort; a single-cell subset included 109 patients with type 1 diabetes and 56 healthy controls; and a validation cohort included 141 patients with type 1 diabetes and 386 healthy controls. All participants were from the Chinese population, and peripheral blood genomic DNA samples were collected from all participants. Illumina Asian Screening Array (ISA) chip typing, quality control, and genotyping were performed. Simultaneously, 10X single-cell transcriptome sequencing was conducted on the subset samples. Four-digit HLA high-resolution typing was performed using NGS-MiSeq.

[0122] A full set of quality control filters was performed on the single-cell transcriptome data to remove low-quality cells; standardization, PCA, Harmony batch correction, and Leiden clustering were performed to complete the annotation of 10 immune cell classes; pseudo-batch expression matrices for each immune subset were generated; corresponding covariates were included, and cis-eQTL linear regression analysis was performed with a window of ±100kb and FDR < 0.05 for screening; multiple rounds of conditional eQTL were performed to remove linkage disequilibrium interference; cross-matching with the T1D-GWAS dataset was performed to obtain 100 risk eSNPs.

[0123] We used sex- and region-corrected logistic regression (Newton-Raphson iteration, Bonferroni correction) to screen for 35 HLA alleles significantly associated with T1D; we imported 81 clinical phenotype SNPs of T1D from literature; and combined the three to obtain an initial set of 216 genetic features after deduplication.

[0124] A backward elimination iterative algorithm was employed, with five-fold cross-validation AUC as the evaluation metric. Redundant features were iteratively removed to obtain 190 final features. Feature relevance scores were calculated to filter features for optimization. An objective function with L2 regularization and upper and lower bound constraints on coefficients was constructed, and the L-BFGS-B algorithm was used iteratively to obtain the optimized β weights of each feature. A weighted summation GRS model was then constructed. Model evaluation was performed in the discovery queue and independent validation queue, respectively. Figure 2 , Figure 3 The ROC-AUC results were obtained; the Spearman correlation between GRS and age of diagnosis was calculated; and the thresholds for low, medium, and high risk classification were determined based on the GRS quantiles of the Chinese population cohort sample.

[0125] One hundred Chinese subjects to be screened were selected, peripheral blood DNA was collected, and only the genotypes of the loci corresponding to the final feature set were detected; the DNA was imported into the prediction system, the GRS score was automatically calculated, and the risk level and risk warning were output.

[0126] A further embodiment: a combined risk assessment module incorporating clinical phenotypic indicators. Based on the GRS genetic risk score obtained from the aforementioned implementation method, this system can further add a combined risk assessment module to achieve a fusion assessment of the genetic risk score and clinical phenotypic indicators. This module is optional and can be configured according to actual usage needs, without affecting the basic function of the aforementioned GRS model to independently complete genetic risk assessment.

[0127] Specifically, the input data of the joint risk assessment module includes two parts: the first part is the subject's genetic risk score output by the GRS prediction model of this invention; the second part is the clinical phenotypic indicators collected from the subject, which can be at least one of the following: type 1 diabetes-related autoantibody test results and diabetes family history information; autoantibodies can include classic type 1 diabetes screening antibody indicators such as islet cell antibodies, glutamic acid decarboxylase antibodies, and insulin autoantibodies; family history indicators are used to identify whether the subject's first-degree relatives have a history of type 1 diabetes.

[0128] In terms of module execution logic, a genetic-clinical joint risk assessment reference benchmark is established in advance based on the Chinese population cohort of this invention: the GRS genetic risk score of each subject in the cohort, along with their corresponding autoantibody positive / negative results and family history information, are used as inputs to train the joint risk assessment criteria. During training, the optimized GRS model weights and final feature set of this invention are retained, and the internal parameters of the already trained genetic model are not modified. Clinical phenotypic information is only superimposed at the model output level to avoid damaging the performance of the original genetic model and ensure the stability and reliability of the genetic prediction results.

[0129] The output of combined risk can be implemented using a tiered correction logic: the genetic risk level given by the GRS is used as a baseline, and then clinical indicators are added to adjust the risk level. An exemplary criterion is as follows: if a subject's GRS indicates medium risk, and they also have at least one positive type 1 diabetes-related autoantibody, or a family history of type 1 diabetes, then the subject's overall risk is upgraded to high risk; if a subject's GRS indicates low risk, but they also have multiple positive autoantibodies, then enhanced clinical follow-up monitoring is required; if a subject's GRS indicates high risk, and they also have positive autoantibodies and a positive family history, then the highest level of combined risk warning is output.

[0130] In addition to retaining the original GRS score and genetic risk level, this module outputs a comprehensive risk level and explains the influencing factors of clinical indicators for clinicians' reference. It should be noted that this combined assessment module is an optional extension; the core genetic risk assessment still relies entirely on the GRS model obtained through feature screening and constraint optimization as described above. Clinical indicators only perform risk correction at the outcome level and do not participate in the feature training and weight optimization process within the GRS model, thus ensuring the independence and reproducibility of the genetic model of this invention.

[0131] With the addition of this module, the system can output both genetic risk assessment reports independently and comprehensive assessment reports combining genetic and clinical information, adapting to different clinical application scenarios: when only gene testing data is available, the genetic risk results can be used directly; when both gene and clinical testing data are available, the combined assessment mode can be activated to further improve the overall screening efficiency.

[0132] Relying solely on the genetic GRS score can only reflect an individual's genetic susceptibility. When an individual's genetic risk falls within the middle range, relying solely on the genetic score can lead to blurred distinctions. This implementation method, without altering the already trained GRS model, integrates readily available clinical phenotypic indicators such as autoantibodies and family history into the output layer to reasonably correct the genetic risk results. This further narrows the interpretation ambiguity range for the middle-risk population, enhancing the overall practical value of clinical screening. Simultaneously, it maintains the genetic model parameters unchanged, ensuring the stability and reproducibility of the core solution of this invention.

[0133] The above description is merely an illustrative embodiment of the present invention and is not intended to limit the scope of the invention. Any equivalent changes and modifications made by those skilled in the art without departing from the concept and principles of the present invention should fall within the scope of protection of the present invention.

Claims

1. A method for predicting the risk of type 1 diabetes based on eQTL and HLA typing, characterized in that, Includes the following steps: S1. Obtain biological samples from the cohort population to obtain high-resolution HLA typing data and genomic SNP genotype data of the cohort population; S2. Based on the genomic SNP genotype data and immune cell single-cell transcriptome sequencing data of the cohort population, a set of immune cell-specific eQTL sites is screened. S3. Integrate the set of immune cell-specific eQTL sites, HLA allele characteristics associated with the risk of type 1 diabetes, and SNP characteristics of publicly disclosed clinical phenotypes of type 1 diabetes to obtain an initial set of genetic characteristics. S4. Perform feature screening and weight coefficient optimization on the initial genetic feature set to obtain the final feature set and the risk weight coefficients corresponding to each feature, and construct a genetic risk score (GRS) prediction model for assessing the risk of developing type 1 diabetes. S5. Obtain biological samples from the subject to be tested, and detect the genotype data of the subject to be tested corresponding to the final feature set; S6. Input the genotype data of the subject to be tested into the GRS prediction model to calculate the genetic risk score of the subject to be tested; S7. Based on the genetic risk score, assess the risk of the subject developing type 1 diabetes and / or early-onset type 1 diabetes.

2. The method for predicting the risk of type 1 diabetes based on eQTL and HLA typing according to claim 1, characterized in that, In step S1, the biological samples used by the cohort population and the subjects to be tested are genomic DNA samples derived from peripheral blood; the high-resolution HLA typing is HLA allele typing with a four-digit resolution.

3. The method for predicting the risk of type 1 diabetes based on eQTL and HLA typing according to claim 1, characterized in that, In step S2, the immune cell-specific eQTL site is a cis-eQTL site; the analysis window for cis-eQTL is set to ±100kb upstream and downstream of the corresponding target gene.

4. The method for predicting the risk of type 1 diabetes based on eQTL and HLA typing according to claim 1, characterized in that, Step S2 includes the following sub-steps: S2-1: Perform quality control, standardization, dimensionality reduction and clustering on single-cell transcriptome sequencing data to annotate and obtain various immune cell types; S2-2: Generate a pseudo-batch gene expression matrix for each type of immune cell; S2-3: A linear regression model was used to conduct cis-eQTL association analysis. Covariates were included to correct for confounding factors. The false discovery rate (FDR) was used as the significance screening threshold to obtain preliminary eQTL association results. S2-4: Conditional eQTL analysis was used to remove linkage disequilibrium interference, and several independent eSNPs were retained for each target gene; S2-5: Cross-matching was performed using publicly available GWAS risk variant datasets for type 1 diabetes to screen for a set of immune cell-specific eQTL sites associated with type 1 diabetes risk.

5. The method for predicting the risk of type 1 diabetes based on eQTL and HLA typing according to claim 4, characterized in that, The immune cells are selected from at least one of CD14+ monocytes, CD16+ monocytes, conventional dendritic cells, plasmacytoid dendritic cells, CD4+ T cells, CD8+ T cells, proliferating T cells, NK cells, B cells, and plasma cells; the covariates include genotype principal component, age, sex, experimental batch, and, for type 1 diabetes patients, disease duration as a covariate.

6. The method for predicting the risk of type 1 diabetes based on eQTL and HLA typing according to claim 1, characterized in that, In step S3, HLA allele characteristics that are significantly associated with the risk of type 1 diabetes are obtained by using a logistic regression model that adjusts for gender and geographical confounding factors, and then filtering according to a preset statistical significance threshold.

7. The method for predicting the risk of type 1 diabetes based on eQTL and HLA typing according to claim 1, characterized in that, The initial set of genetic features obtained in step S3 includes three types of features: immune cell-specific risk eSNP features, HLA allele features related to type 1 diabetes risk obtained through screening, and SNP features related to type 1 diabetes clinical phenotype obtained from published literature.

8. The method for predicting the risk of type 1 diabetes based on eQTL and HLA typing according to claim 1, characterized in that, In step S4, the feature selection adopts the backward elimination iterative selection algorithm; the AUC value obtained by cross-validation is used as the evaluation index for feature selection, and features that have no benefit to the model prediction performance are iteratively removed to obtain a simplified feature subset.

9. The method for predicting the risk of type 1 diabetes based on eQTL and HLA typing according to claim 8, characterized in that, In step S4, the weight coefficient optimization uses an objective function with a regularization term for iterative solution. Numerical constraints are set on the optimized weight coefficients relative to the original effect coefficients, and a quasi-Newton optimization algorithm is used to solve the problem. After convergence, the optimized risk weight coefficients corresponding to each feature are obtained. The GRS prediction model is a weighted summation model, and the formula for calculating the genetic risk score is: ； Where n is the total number of features in the final feature set; The optimized risk weight coefficient for the j-th feature; is the genotype encoding value of the j-th genetic characteristic locus of the subject.

10. A type 1 diabetes risk prediction system based on eQTL and HLA typing, characterized in that, The method described in any one of claims 1-9 comprises: The sample genotype acquisition module is used to acquire biological samples from the cohort population and the subjects to be tested and obtain the corresponding genotype data; The eQTL site screening module is used to screen for a set of immune cell-specific eQTL sites based on genomic SNP genotype data and immune cell single-cell transcriptome sequencing data. The feature integration module is used to integrate multiple types of genetic features to generate an initial set of genetic features; The model building module is used to perform feature selection and weight coefficient optimization on the initial genetic feature set to build a GRS prediction model. The subject risk calculation module is used to input the genotype data of the subjects to be tested into the GRS prediction model and calculate the genetic risk score; The risk outcome output module is used to output the risk assessment results of the subject's risk of developing type 1 diabetes and / or early-onset type 1 diabetes.

Images