A genotyping method based on DNA methylation chip and its application

Abstract

The present invention relates to the field of epigenomic association research, and relates to a genotyping method based on DNA methylation chip and its application. The present invention has developed a method for detecting genotypes and inferring kinship between samples based on Infinium methylation chip data. The method uses chip raw data (.IDAT format) or processed data (β or M value matrix) as input, detects genotypes for SNP probes, type I probes and type II probes respectively, and generates standard VCF files. The present invention realizes for the first time the accurate detection of the genotypes of thousands of SNPs from methylation chip data, and realizes accurate inference of population structure and kinship based on these genotypes. If this method is used to exclude samples with kinship before conducting epigenomic association analysis (EWAS), false positive results caused by sample correlation can be reduced, the accuracy of research results can be improved, and the study of disease mechanisms can be promoted.

Description

Technical Field

[0001] The present invention relates to the field of epigenomic association research, and more specifically, to a genotyping method and application based on DNA methylation chip, and in particular to a method for detecting genotypes based on DNA methylation chip data, inferring population structure, and identifying kinship between samples. Background Art

[0002] DNA methylation primarily occurs at cytosine-guanine dinucleotides (CpG) and is the process by which a methyl group is attached to the fifth carbon atom of cytosine by the action of DNA methyltransferases. DNA methylation changes dynamically under environmental influences and plays a crucial role in gene regulation. With the rapid development of high-throughput methylation array technologies, DNA methylation has become the most widely studied epigenetic modification. Mainstream methylation array technologies include Infinium Human Methylation 450 (450K) and Infinium Human Methylation EPIC (EPIC). These technologies can simultaneously detect hundreds of thousands of CpG sites, enabling the widespread application of epigenomic association studies (EWAS). These technologies have identified a large number of CpG sites associated with a variety of complex diseases and environmental exposures, deepening our understanding of disease mechanisms (Rakyan et al. Nat Rev Genet 2011, Wei et al. Adv Sci (Weinh) 2021, Fraszczyk et al. Diabetologia 2022).

[0003] When conducting EWAS studies, the handling of related samples is particularly crucial. Because related samples are exposed to similar environments over long periods of time, their methylation levels are more similar. Including these samples in EWAS studies can easily lead to false positive signals (Rakyan et al. Nat Rev Genet 2011, Gross et al. BMC Genet 2017, Campagna et al. Clin Epigenetics 2021). However, despite the importance of confounding effects from relatedness, most EWAS studies fail to account for this effect due to the lack of methods to calculate relatedness from methylation data. Even in cohorts that generate both methylation and genotype array data, direct use of relatedness inferred from genotype array data for EWAS studies can be problematic due to differences in the quality control processes used for the two datasets and the resulting differences in the samples that pass quality control. Therefore, it is crucial to develop methods for inferring relatedness from methylation data to better facilitate EWAS studies.

[0004] If sufficient single nucleotide polymorphism (SNP) data are available, we can very accurately infer population structure and inter-sample kinship (Purcell et al. Am J Hum Genet 2007, Manichaikul et al. Bioinformatics 2010, Thornton et al. Am J Hum Genet 2012, Conomos et al. Am J Hum Genet 2016, Dou et al. PLoS Genet 2017). Although the Infinium methylation array has dozens of SNP probes designed to detect sample mix-up (for example, the EPIC array has 59 SNP probes, 6 of which are located on chromosome X) (Assenov et al. Nat Methods 2014, Heiss et al. Clin Epigenetics 2018, Muller et al. Genome Biol 2019), the number of these SNPs is far too small to be used to accurately infer kinship. On the other hand, when performing methylation data quality control, hundreds of probes located near common SNPs (minor allele frequency MAF>0.01) are deleted to avoid the influence of SNPs on methylation signal detection results (McCartney et al. GenomData 2016, Pidsley et al. Genome Biol 2016, Zhou et al. Nucleic Acids Res 2017). If the SNP is located at the single-base extension (SBE) position of the probe, the methylation signal of these probes usually varies with the SNP genotype, thus showing a multimodal distribution (usually trimodal) (Daca-Roszak et al. BMC Genom 2015, Andrews et al. Epigenetics Chromatin 2016, LaBarre et al. Epigenetics Chromatin 2019). Therefore, we believe that the methylation signal intensity of these probes can be used to infer genotype and further infer population structure and kinship.

[0005] The Infinium methylation array is designed with two types of probes: Type I and Type II. Type I probes feature two probes at each CpG site to be tested, pairing with the methylated and unmethylated sequences, respectively. Successful pairing results in an extension of one base (SBE) and a fluorescent signal. If the SBE is A or T, red fluorescence is generated; if it is G or C, green fluorescence is generated. Mutations other than A / T and G / C at the SBE position result in a shift in the color channel (CCS). By comparing the signal intensities of the different color channels, the genotype at the SBE position can be inferred (Zhou et al. Nucleic Acids Res 2017). Type II probes feature one probe at each CpG site to be tested, with the SBE corresponding to the C site to be tested. After sulfite conversion, unmethylated Cs are converted to Ts, which emit red fluorescence upon probe binding; methylated Cs remain unchanged and emit green fluorescence upon probe binding. In the presence of a SNP, if the C site being tested mutates to an A or T, it will appear red; if it mutates to a G, it will appear green. Because both methylation and SNPs affect the color channel, inferring genotypes for type II probes is more challenging. Although the majority of probes on the Infinium array are type II (Pidsley et al. Genome Biol 2016), there is currently no method for inferring genotypes for type II probes. A method that can accurately infer genotypes and identify kinship relationships based on DNA methylation arrays is needed to improve the accuracy of EWAS results. Summary of the Invention

[0006] This paper has developed a method for detecting genotypes and inferring kinship between samples based on Infinium methylation array data (EPIC or 450K arrays). This method uses raw array data (.IDAT format) or processed data (beta or M-value matrix) as input, detects genotypes for SNP probes, type I probes, and type II probes, and infers population structure and kinship based on these genotypes.

[0007] According to a first aspect of the present invention, a genotyping method based on a DNA methylation chip is provided, comprising the following steps:

[0008] (1) Calculate the signal ratio RAI supporting the mutant allele for each candidate SNP probe, candidate type I probe, and candidate type II probe, and generate a RAI matrix, in which each probe corresponds to a RAI value in each sample;

[0009] (2) constructing a mixed model of RAI distribution for each of the three probes, wherein the mixed model includes three beta distributions and one uniform distribution, wherein the three beta distributions represent the three genotypes of reference homozygote, heterozygote, and mutant homozygote, respectively, and the uniform distribution represents background noise;

[0010] (3) solving the parameters of the mixed model described in step (2) and calculating the background probability of each probe in each sample and the three genotype probabilities of reference homozygote, heterozygote and mutant homozygote;

[0011] (4) Based on the background probability and the three genotype probabilities obtained in step (3), calculate the genotype of each probe in each sample.

[0012] Preferably, in step (1), the candidate SNP probe is a SNP probe designed on a methylation chip; the candidate type I probe and the candidate type II probe are probes on a methylation chip, and a common SNP exists at the base extension position thereof, and the common SNP has a minor allele frequency MAF>0.01; the RAI calculation formulas for the three probes, namely, the SNP probe, the type I probe, and the type II probe, are as follows:

[0013] For SNP probes, let S(p REF ) and S(p ALT ) correspond to the signal intensity of the cytosine allele and the signal intensity of the mutant allele, respectively. The RAI calculation formula is:

[0014]

[0015] For type I probes, the RAI calculation formula is:

[0016]

[0017] where p M and p U Refers to methylated and unmethylated probes, respectively. oob Indicates the out-of-band signal strength value, S ib Indicates the in-band signal strength value;

[0018] For type II probes, het It is defined as the β value corresponding to the position of the middle peak, which corresponds to the heterozygous genotype, and the proportion of methylated C bases in each CpG unmutated C is calculated as p M Perform the calculation:

[0019]

[0020] Then, calculate RAI as follows:

[0021]

[0022] Preferably, in step (2), the hybrid model is constructed as follows: suppose the RAI value is an m×n matrix X, where m and n represent the number of probes and the number of samples, respectively; assume that X obeys a mixed distribution of three beta distributions and one uniform distribution, where the three beta distributions correspond to the three genotypes respectively, and the uniform distribution represents background noise:

[0023]

[0024]

[0025] where X ij represents the RAI value of probe i in sample j, k represents the genotype corresponding to the reference homozygote, heterozygote and mutant homozygote, and the values ​​are 0, 1 and 2 respectively. Beta(α k ,β k ) represents beta distribution, U(0,1) represents uniform distribution, λ represents the probability that RAI comes from background noise, (1-λ)w ik represents the probability of RAI corresponding to genotype k, w ik represents the weight determined by the allele frequency AF, φ i represents the allele frequency AF of the SNP corresponding to probe i.

[0026] Preferably, in step (3), the expectation maximization algorithm is used to solve the model parameters, and the specific steps are as follows:

[0027] Let B ijk =Beta(X ij ; α k ,β k ) is X ij According to Beta(α k ,β k ) distribution, assuming that X ij Assuming that all probes and samples are independent of each other, the log-likelihood function can be written as:

[0028]

[0029] The expectation maximization algorithm is divided into two steps: the first step is to calculate the expectation of the hidden variable; the second step is to calculate the value of the parameter by maximizing the likelihood function; the parameters are solved by repeating these two steps;

[0030] In the first step, calculate X ij The probabilities of U(0,1) are respectively derived from and Beta(α k ,β k )

[0031]

[0032]

[0033] In the second step, the model parameters are reestimated using moment estimation:

[0034]

[0035]

[0036]

[0037]

[0038] in,

[0039]

[0040]

[0041] Repeat these two steps until the log-likelihood function converges to the maximum value and obtains the maximum likelihood estimate of the model parameters (α, β, φ, λ) and the genotype probability and background probability

[0042] Preferably, in step (4), the specific process of inferring genotypes is as follows: the genotype with a larger background probability is set as missing, and the genotype probability is set as Updated to To ensure that for any probe i and sample j For each genotype not set to missing, define is the genotype with the highest probability, and its values ​​are 0, 1, and 2 respectively. is the dosage genotype, with a value between 0 and 2; it is calculated using the following formula and

[0043]

[0044]

[0045] Where n is the sample size, for The variance in all samples; the genotype test results generated above include all candidate sites and samples to be studied genotype Dosage genotype Genotype probability and allele frequency AF.

[0046] Preferably, the genotype with a larger background probability is set to be missing, specifically: The genotype is set to missing.

[0047] According to another aspect of the present invention, there is provided an application of any one of the DNA methylation chip-based genotyping methods for inferring population structure.

[0048] Preferably, a reference ancestral space is constructed using a genetic dataset with a known population structure, and the research sample is projected into the reference space using a LASER algorithm to infer the population structure of the research sample.

[0049] According to another aspect of the present invention, there is provided an application of any one of the DNA methylation chip-based genotyping methods for inferring kinship.

[0050] Preferably, for samples from a single population and a mixed population, the SEEKIN-hom method and the SEEKIN-het method are used to infer kinship, respectively.

[0051] In general, the above technical solutions conceived by the present invention have the following technical advantages compared with the existing technology:

[0052] (1) The existing technology can only infer genotypes for hundreds of SNP sites in SNP probes and type I probes. The present invention innovatively uses RAI values ​​to represent the signal ratio supporting mutant alleles, and for the first time, it can genotype a large number of SNPs located on type II probes. In addition, when fitting the mixed distribution of RAI values, the present invention assigns different weights based on the AF of each SNP, further improving the accuracy of genotyping. The present invention is the first to accurately detect the genotypes of thousands of SNPs from methylation chip data (EPIC chip: ~4000, 450K chip: ~2000).

[0053] (2) Based on these genotyping results, the population structure, kinship, inbreeding coefficient, etc. can be accurately inferred. When dealing with samples with large population heterogeneity, adding population structure information can greatly improve the accuracy of kinship results.

[0054] (3) If this method is used to exclude samples with kinship before conducting EWAS, it can reduce false positive results caused by sample correlation, improve the accuracy of research results, and promote the study of disease mechanisms.

[0055] (4) The present invention integrates all the above processes into the MethylGenotyper software package based on R language, so that the entire analysis process from the raw data of the methylation chip to the generation of kinship and other results can be solved in one stop. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 It is a data processing flow chart involved in the present invention.

[0057] Figure 2 The distribution of RAI in the simulated data is shown in Figure 1. A is a type I probe, containing 400 SNPs and 3200 samples; B is a type II probe, containing 4000 SNPs and 3200 samples.

[0058] Figure 3 The results of testing MethylGenotyper using simulated data generated based on the characteristics of type I probes. Where AC is the error rate of predicting α; DF is the error rate of β; G is the error rate of λ; H is the error rate of φ; and I is the genotype consistency rate with the true value.

[0059] Figure 4 The results of testing MethylGenotyper using simulated data generated based on the characteristics of type II probes. Where AC is the error rate of predicting α; DF is the error rate of β; G is the error rate of λ; H is the error rate of φ; and I is the genotype consistency rate with the true value.

[0060] Figure 5 The performance of MethylGenotyper on the EPIC methylation data of the Dongfeng Tongji cohort is shown in Figure 2. AC represents the fit of the RAI distribution of SNP probes, type I probes, and type II probes to a mixed beta distribution; DF represents the consistency comparison between the AF obtained for these three probes and the genotyping array results.

[0061] Figure 6 Results of kinship inference using EPIC methylation data from the Dongfeng Tongji cohort. A represents the results based on SNP probes; B represents the results based on SNP probes and type I probes; C represents the results based on SNP probes, type I probes, and type II probes; and D represents the results based on genotyping array data, which served as the gold standard.

[0062] Figure 7 This is the result of inferring kinship by taking the genotypes of the common probes of the 450K chip from the Dongfeng Tongji cohort.

[0063] Figure 8 Comparison of AF based on AIBL methylation data with AF from the 1KGP European population. AC represents SNP probes, type I probes, and type II probes, respectively.

[0064] Figure 9Figure 2 shows the performance of MethylGenotyper on the AIBL dataset. A shows the results of population structure inference; B shows the comparison of kinship inferred from methylation data and genotyping data. DETAILED DESCRIPTION

[0065] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely for the purpose of explaining the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.

[0066] This paper has developed a method for detecting genotypes and inferring kinship between samples based on Infinium methylation array data (EPIC or 450K array). This method uses raw array data (.IDAT format) or processed data (β or M value matrix) as input, detects genotypes for SNP probes, type I probes, and type II probes, and infers population structure and kinship based on these genotypes. The specific steps are as follows: Figure 1 It is a data processing flow chart involved in the present invention.

[0067] In the first step, background correction was performed using the noob method (Triche et al. Nucleic Acids Res 2013, Zhou et al. Nucleic Acids Res 2018) and staining bias correction was performed using the RELIC method (Xu et al. BMC Genom 2017).

[0068] In the second step, we calculated the RAI for each of the candidate SNP probe, type I probe, and type II probe. We use RAI to represent the proportion of signal supporting the mutant allele. For each of the three probes, the calculation method is as follows:

[0069] For SNP probes, two probes were designed for each candidate SNP site, corresponding to the signal intensity of the cytosine allele and the mutant allele, respectively, and recorded as S(p REF ) and S(p ALT ). The RAI calculation formula is:

[0070]

[0071] For type I probes, we calculated RAI according to the algorithm of Zhou et al. (Zhou et al. Nucleic Acids Res 2017):

[0072]

[0073] where p M and p U Refers to methylated and unmethylated probes, respectively. oob Indicates the out-of-band signal strength value (the actual color is opposite to the design), S ib Indicates the in-band signal strength value (the actual color is the same as the design).

[0074] For type II probes, we used the multimode method (Ameijeiras-Alonso et al. J StatSoftw 2021) to examine the distribution of methylation levels β for each CpG, including the number, position, and height of peaks. We only retained probes with two or more peaks (height > 0.001) (bandwidth set to 0.04). If the number of peaks detected for a probe was greater than 3, we removed the lowest peaks until only three peaks remained. We will l het It is defined as the β value corresponding to the position of the middle peak (if there are only two peaks, it is the position of the peak closest to 0.5), corresponding to the heterozygous genotype, and the proportion of methylated C bases in the unmutated CpG is estimated for each CpG (p M ):

[0075]

[0076] Finally, we calculate RAI as follows:

[0077]

[0078] Limit the RAI range to 0.01 to 0.99: if it is greater than 0.99, it is forced to be set to 0.99; if it is less than 0.01, it is forced to be set to 0.01; if it is greater than or equal to 0.01 and less than or equal to 0.99, the actual calculated RAI value is used;

[0079] In the third step, we further filtered type I and type II probes based on their RAI distribution. We used the multimode method to examine the RAI distribution of each probe across all input samples, ignoring peaks with heights less than 0.001. We then retained type I probes with ≥2 peaks and type II probes with ≥3 peaks.

[0080] Step 4: Genotype detection and quality control. For each probe, let its RAI value be an m×n matrix X, where m and n represent the number of probes and samples, respectively. Assume that X follows a mixed distribution of three beta distributions (corresponding to the three genotypes) plus a uniform distribution (representing background noise):

[0081]

[0082]

[0083] where X ij represents the RAI value of probe i in sample j, k represents the genotype (with values ​​of 0, 1, and 2 corresponding to reference homozygote, heterozygote, and mutant homozygote, respectively), Beta (α k ,β k ) represents beta distribution, U(0,1) represents uniform distribution, λ represents the probability that RAI comes from background noise, (1-λ)w ik represents the probability of RAI corresponding to genotype k, w ik represents the weight determined by the allele frequency (AF), φ i represents the allele frequency (AF) of the SNP corresponding to probe i.

[0084] Let B ijk =Beta(X ij ; α k ,β k ) is X ij According to Beta(α k ,β k ) distribution, assuming that X ij Assuming that all probes and samples are independent of each other, the log-likelihood function can be written as:

[0085]

[0086] The parameters (α, β, φ, λ) of the mixture distribution are estimated using the expectation maximization (EM) algorithm. The initial values ​​of these parameters are set as follows: α0 = 5, β0 = 60, α1 = β1 = 30, α2 = 60, β2 = 5, λ = 0.01. The EM algorithm consists of two steps: the E step and the M step. The parameters are solved by iterating these two steps repeatedly.

[0087] In step E, we calculate X ij The probabilities of U(0,1) are respectively derived from and Beta(α k ,β k )

[0088]

[0089]

[0090] In the M step, we re-estimate the model parameters using moment estimation:

[0091]

[0092]

[0093]

[0094]

[0095] in,

[0096]

[0097]

[0098] Repeat the iterative E and M steps until the log-likelihood function converges to the maximum value and obtains the maximum likelihood estimate of the model parameters (α, β, φ, λ) and the genotype probability and background probability

[0099] Will The genotype of Updated to To ensure that for any probe i and sample j For each genotype that is not set to missing, we define is the genotype with the highest probability (values ​​are 0, 1, and 2), is the dosage genotype (value ranges from 0 to 2). According to the method of Li et al. (Li et al. Genet Epidemiol 2010), we use the following formula to calculate and

[0100]

[0101]

[0102] Where n is the sample size, for The variance in all samples. We will be missing or The genotype ratio was defined as the missing rate at the SNP level. To ensure data quality, we excluded the samples with MAF < 0.01 and HWE p < 10 -6 、 or SNPs with missing rate > 0.1, and the dosage genotype matrix of the remaining SNPs For subsequent steps.

[0103] The fifth step is to infer the population structure and calculate the sample-specific AF. If the sample under study comes from a mixed population, we use the LASER method (Wang et al. Nat Genet 2014, Wang et al. Am J Hum Genet 2015) to infer the population structure. This method first performs principal component analysis (PCA) on a genetic dataset with known population structure (such as the 1000 Genomes Project 1KGP) and extracts the first four principal components (PCs, denoted as ) is defined as the reference ancestral space. Next, based on the dosage genotype matrix obtained in the previous step We use the trace program in the LASER method to project each sample to be studied into the reference space, and then we can obtain the PCA results of all the samples to be studied. Similarly, we use the first four PCs (denoted as ) is used to represent the population structure of the sample to be studied.

[0104] Based on the above reference ancestral space, we can obtain the genotype for each SNP in 1KGP With the first four principal components The linear relationship: in is the genotype of SNP i in sample j in 1KGP, is the PC coordinate of sample j, β i· is the regression coefficient of SNP i. Next, the specific in is the PC coordinate of the sample to be studied. We define The range is 0.001 to 0.999, and the value greater than 0.999 Will be forced to 0.999, less than 0.001 Will be forced to 0.001.

[0105] Step 6: Infer kinship and inbreeding coefficient. We use the SEEKIN method (Dou et al. PLoS Genet 2017) to calculate the kinship coefficient between samples. Considering the uncertainty of the inferred genotype, this method is based on the fourth step. Different weights are applied to each SNP. For samples from a single population and a mixed population, we use SEEKIN-hom and SEEKIN-het to estimate kinship, respectively. The calculation formula for SEEKIN-hom is:

[0106]

[0107] where φ xyis the kinship coefficient between sample x and sample y, and is the dosage genotype obtained in the fourth step, is the AF obtained in the fourth step, is the weight. The calculation formula of SEEKIN-het is

[0108]

[0109] in and is the sample-specific AF obtained in the fifth step.

[0110] Based on the above formula, the kinship coefficient (φ xy According to the study of Manichaikul et al. (Manichaikul et al. Bioinformatics 2010), we set the kinship coefficient at 2 -t-1.5 to 2 -t-0.5 The sample pairs between are defined as t-level kinship. The calculation formula of inbreeding coefficient is: f x =2φ xx -1, where φ xx This coefficient can be calculated using the SEEKIN-hom formula. It reflects whether the parents are closely related. A coefficient greater than 0.0156 indicates that the father and mother share a common great-grandparent.

[0111] This method preferably uses raw data as input, as it automatically corrects the data for background and staining bias. If a β-value matrix or M-value matrix is ​​used as input, ensure that the data has been corrected for background and staining bias and that no other corrections (such as BMIQ) have been applied. Furthermore, genotype inference for type I probes will not be possible if a β-value matrix or M-value matrix is ​​used as input.

[0112] Preferably, this method mainly analyzes common SNPs, and it is recommended to first remove probes with MAF < 0.01 in the corresponding population.

[0113] Preferably, if the sample comes from a mixed population, it is necessary to infer the population structure before inferring kinship. When inferring population structure, it is recommended to use the LASER method (Wang et al. Am J Hum Genet 2015) and combine it with a genetic dataset with known population structure (such as 1KGP, etc.) to construct a reference ancestral space, and obtain a unique AF for each SNP. These AFs are used as SEEKIN het Input to infer kinship.

[0114] Preferably, although the method can be used with a sample size of 100, it works best with a sample size greater than 800.

[0115] The above invention has been compiled into a software package (MethylGenotyper) using the R language for user convenience. This package can generate standard VCF files to store genotype information. In addition, the genotypes detected by the present invention can also be used to calculate the inbreeding coefficient, which can be used to determine whether a sample is contaminated.

[0116] The present invention is used to detect genotypes and infer kinship from methylation chip data. The specific operating steps are as follows:

[0117] Step 1: Import the MethylGenotyper package

[0118] library(MethylGenotyper)

[0119] Step 2: Read the data and perform noob and staining bias correction

[0120] This step is limited to inputting raw chip data in the ".IDAT" format. If you use a methylation β or M-value matrix as input, you can skip this step. The command is as follows:

[0121] rgData<-correct_noob_dye(target,cpu=3)

[0122] This command reads the sample information in target. target is a data frame containing two columns, "Sample_Name" and "Basename." "Basename" should contain the path information for each sample. Each row corresponds to a single sample. This command reads each sample individually and performs noob and staining bias correction. You can trigger a multi-process task by specifying the "cpu" parameter. After the task completes, the variable "rgData" will contain the corrected signal values ​​for all candidate probes.

[0123] Step 3: Genotype detection

[0124] This method can detect genotypes for SNP probes, type I probes, and type II probes respectively. The commands are as follows:

[0125] genotype_snp<-callGeno_snp(rgData,input="raw",vcf="TRUE,pop="EAS")

[0126] genotype_typeI<-callGeno_typeI(rgData,vcf=TRUE,pop="EAS")

[0127] genotype_typeII<-callGeno_typeII(rgData,input="raw",vcf="TRUE,pop="EAS")

[0128] These three functions take the variable "rgData" generated in the previous step as input and output standard VCF files after genotyping. For more accurate results, the user should specify the appropriate population (EAS, AMR, AFR, EUR, SAS, and ALL) based on the sample size. For data from multiple populations, it is recommended to specify the population with the largest sample size.

[0129] By specifying the parameter input="beta" or input="mval", this method can also directly use the methylation beta or M value matrix as input (except for type I probes). The commands for SNP probes and type II probes are as follows:

[0130] genotype_snp<-callGeno_snp(beta_matrix,input="beta",vcf="TRUE,pop="EAS")

[0131] genotype_typeII<-callGeno_typeII(beta_matrix,input="beta",vcf=TRUE,pop="EAS")

[0132] The input matrix for this command uses probes as rows and samples as columns. It is important to note that the input data must have been corrected for background and staining bias, and no other corrections (such as BMIQ) have been performed, as otherwise the accuracy of the results may be affected.

[0133] Next, the genotypes obtained from the three probes can be combined using the following command:

[0134] dosage<-rbind(genotype_snp$dosage,genotype_typeII$dosage)

[0135] Step 4: Inferring the population structure

[0136] The specific code is as follows:

[0137] pc<-projection(dosage,plotPCA=TRUE,cpu=3)

[0138] data(cpg2snp)

[0139] snpvec<-cpg2snp[c(

[0140] rownames(genotype_snp$genotypes$RAI),

[0141] rownames(genotype_typeI$genotypes$RAI),

[0142] rownames(genotype_typeII$genotypes$RAI) )]

[0144] indAF<-get_indAF(snpvec,pc$refPC,pc$studyPC)

[0145] Where pc$studyPC is the population structure of the sample to be studied (i.e., the first four PCs), and the variable "indAF" is the AF specific to each probe and sample in the sample to be studied.

[0146] Step 5: Inferring kinship and inbreeding coefficients

[0147] For samples from a single population, SEEKIN is recommended. hom To estimate the kinship coefficient, the code is as follows:

[0148] res<-getKinship(dosage)

[0149] For samples from mixed populations, SEEKIN is recommended. het To estimate the kinship coefficient, the code is as follows:

[0150] res<-getKinship_het(dosage,indAF)

[0151] Next, you can use the following command to retrieve the results of the kinship coefficient and inbreeding coefficient:

[0152] kinship<-res$kinship#kinship coefficients

[0153] inbreed<-res$inbreed#inbreeding coefficients

[0154] The following are specific embodiments

[0155] Example 1

[0156] The present invention tests MethylGenotyper by generating simulated data based on the parameters of type I and type II probes, respectively. The parameters of type I probe are: α0=3, β0=35, α1=β1=20, α2=65, β2=4, λ=0.025, and the parameters of type II probe are: α0=2, β0=20, α1=35, β1=40, α2=40, β2=3, λ=0.015. For both probes, the three Beta distributions estimated by MethylGenotpyer can perfectly fit the RAI distribution ( Figure 2 ). As the sample size and the number of SNPs increase, the estimated values ​​of these four parameters (α, β, φ, λ) are getting closer to the true values ​​( Figure 3 、 Figure 4 For type II probes, the error rate of parameter estimation stabilizes when the sample size reaches 800. For example, when the sample size is 3200 and the number of SNPs is 4000, the error rate of α is 0.0046 (95% confidence interval: 0.0036-0.0055); the error rate of β is 0.0041 (95% confidence interval: 0.0033-0.0049); the error rate of λ is 0.086 (95% confidence interval: 0.084-0.088); and the error rate of φ is 0.031 (95% confidence interval: 0.031-0.031). The genotype concordance rate of type II probes is approximately 98.4%. In contrast, the genotype concordance rate of type I probes is lower (97.7%), which may be due to the higher background noise of type I probes (λ is 0.025 for type I probes; λ is 0.015 for type II probes) and the smaller number of SNPs.

[0157] Example 2

[0158] Based on the EPIC chip methylation data of 4662 Chinese blood samples from the Dongfeng Tongji cohort, a total of 5271 probes were retained for genotyping, including 53 SNP probes, 168 type I probes, and 5050 type II probes. The RAI values ​​of these probes can be well fitted by three Beta distributions ( Figure 5 The background noise ratios of these three probes were 0.014, 0.024 and 0.017 respectively. The present invention achieved high-quality genotyping for 4319 SNPs (MAF≥0.01, HWE p≥10-6, The panel included 53 SNP probes, 111 type I probes, and 4155 type II probes. Compared with the genotyping chip results of the same sample, the genotype consistency of these SNPs was 98.26%, the heterozygote consistency was 96.64%, and the AF also had very high consistency ( Figure 5DF in ), indicating that the genotype detected by MethylGenotyper has a very high accuracy.

[0159] Based on the genotyping chip, a total of 123 pairs of duplicate samples, 110 pairs of first-level relatives, 22 pairs of second-level relatives, and 53 pairs of third-level relatives were detected from these 4662 samples. Figure 6 The present invention uses this as the gold standard, considers the second-level and higher-level kinship relationships as positive sets, and calculates the accuracy, sensitivity, and F1 value according to the following formula to evaluate the performance of MethylGenotyper:

[0160]

[0161]

[0162]

[0163] When only genotyping data from 53 SNP probes and 111 type I probes were used, it was impossible to distinguish samples with and without kinship, and the accuracy was almost zero ( Figure 6 When the genotype data of 4155 type II probes were additionally included, the variance of the obtained kinship coefficient was greatly reduced, and the kinship relationships at all levels could be clearly distinguished. The accuracy was 0.9659 and the sensitivity was 1 ( Figure 6 C in, F1=0.9827).

[0164] More than 90% of the probes on the 450K chip are included in the EPIC chip (Pidsley et al. Genome Biol 2016). This paper extracts these common probes from the Dongfeng Tongji EPIC chip data and evaluates the performance of MethylGenotpyer on the 450K chip. A total of 2212 SNPs were detected with high quality (MAF ≥ 0.01, HWE p ≥ 10-6, The chip has a missing rate of ≤0.1), including 53 SNP probes, 104 type I probes, and 2055 type II probes. Based on these SNPs, the inference of kinship achieved almost the same accuracy and sensitivity as the EPIC chip ( Figure 7 ).

[0165] Example 3

[0166] Based on the EPIC chip methylation data of 702 samples in the Australian AIBL cohort (Ellis et al. Int Psychogeriatr 2009, Fowler et al. J Alzheimers Dis Rep 2021), a total of 4217 SNPs were detected with high-quality genotypes (MAF ≥ 0.01, HWE p ≥ 10-6, The AF of these probes was highly consistent with that of the 1KGP European population, with the exception of a few SNPs ( Figure 8 ).

[0167] The reference ancestral space was constructed using 1KGP data, and the population structure of the AIBL cohort samples was inferred based on the LASER method. It was found that although most of the samples were European, a small number of samples were East Asian and South Asian ( Figure 9 Based on this, the present invention fully considers the population heterogeneity of the samples studied and uses SEEKIN-het to infer kinship. The study found that the kinship obtained based on methylation data is highly consistent with the kinship obtained based on genotyping chips ( Figure 9 (B) The present invention detected four pairs of samples with first-degree genetic relationships from the AIBL sample for the first time. These results fully demonstrate the superior performance of MethylGenotyper and its reliability in processing mixed population samples.

[0168] It will be easily understood by those skilled in the art that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A genotyping method based on DNA methylation chip, characterized in that: The following steps are involved: (1) Calculate the signal ratio RAI supporting the mutant allele for each candidate SNP probe, candidate type I probe, and candidate type II probe, and generate a RAI matrix, in which each probe corresponds to a RAI value in each sample; (2) constructing a mixed model of RAI distribution for each of the three probes, wherein the mixed model includes three beta distributions and one uniform distribution, wherein the three beta distributions represent the three genotypes of reference homozygote, heterozygote, and mutant homozygote, respectively, and the uniform distribution represents background noise; (3) solving the parameters of the mixed model described in step (2) and calculating the background probability of each probe in each sample and the three genotype probabilities of reference homozygote, heterozygote and mutant homozygote; (4) Based on the background probability and the three genotype probabilities obtained in step (3), calculate the genotype of each probe in each sample.

2. The genotyping method based on DNA methylation chip according to claim 1, wherein In step (1), the candidate SNP probe is a SNP probe designed on a methylation chip; the candidate type I probe and the candidate type II probe are probes on a DNA methylation chip, and a common SNP exists at the base extension position thereof, and the common SNP has a minor allele frequency MAF>0.01; the RAI calculation formulas for the three probes, namely, the SNP probe, the type I probe, and the type II probe, are as follows: For SNP probes, let S(p REF ) and S(p ALT ) correspond to the signal intensity of the cytosine allele and the signal intensity of the mutant allele, respectively. The RAI calculation formula is: For type I probes, the RAI calculation formula is: Among them S oob (p M ) and S oob (p U ) refer to the out-of-band signal intensity values ​​of methylated and non-methylated probes, respectively, S ib (p M ) and S ib (p U ) refer to the signal intensity values ​​of the methylated and unmethylated bands, respectively; For type II probes, het It is defined as the β value corresponding to the position of the middle peak, which corresponds to the heterozygous genotype, and the proportion of methylated C bases in each CpG unmutated C is calculated as p M Make an estimate: Then, calculate RAI as follows:

3. The genotyping method based on DNA methylation chip according to claim 1, wherein In step (2), the hybrid model is constructed as follows: Let the RAI value be an m×n matrix X, where m and n represent the number of probes and samples, respectively; assume that X obeys a mixed distribution of three beta distributions and one uniform distribution, where the three beta distributions correspond to the three genotypes respectively, and the uniform distribution represents the background noise: where X ij represents the RAI value of probe i in sample j, k represents the genotype corresponding to the reference homozygote, heterozygote and mutant homozygote, and the values ​​are 0, 1 and 2 respectively. Beta(α k ,β k ) represents beta distribution, U(0,1) represents uniform distribution, λ represents the probability that RAI comes from background noise, (1-λ)w ik represents the probability of RAI corresponding to genotype k, w ik represents the weight determined by the allele frequency AF, φ i represents the allele frequency AF of the SNP corresponding to probe i.

4. The genotyping method based on DNA methylation chip according to claim 3, wherein In step (3), the expectation maximization algorithm is used to solve the model parameters. The specific steps are as follows: Let B ijk =Beta(X ij ; α k ,β k ) is X ij According to Beta(α k ,β k ) distribution, assuming that X ij Assuming that all probes and samples are independent of each other, the log-likelihood function can be written as: The expectation maximization algorithm is divided into two steps: the first step is to calculate the expectation of the hidden variable; the second step is to calculate the value of the parameter by maximizing the likelihood function; the parameters are solved by repeating these two steps; In the first step, calculate X ij The probabilities of U(0,1) are respectively derived from and Beta(α k ,β k ) In the second step, the model parameters are reestimated using moment estimation: in, Repeat these two steps until the log-likelihood function converges to the maximum value and obtains the maximum likelihood estimate of the model parameters (α, β, φ, λ) and the genotype probability and background probability 5. The genotyping method based on DNA methylation chip according to claim 1, wherein In step (4), the specific process of inferring genotypes is as follows: the genotype with a larger background probability is set as missing, and the genotype probability is set as Updated to To ensure that for any probe i and sample j For each genotype not set to missing, define is the genotype with the highest probability, and its values ​​are 0, 1, and 2 respectively. is the dosage genotype, with a value between 0 and 2; it is calculated using the following formula and Where n is the sample size, for The variance in all samples; the genotype test results generated above include all candidate sites and samples to be studied genotype Dosage genotype Genotype probability and allele frequency AF.

6. The genotyping method based on DNA methylation chip according to claim 5, characterized in that: The genotype with a larger background probability is set to missing, specifically: The genotype is set to missing.

7. Use of the DNA methylation chip-based genotyping method according to any one of claims 1 to 6 for inferring population structure.

8. The use according to claim 7, characterized in that A reference ancestral space is constructed using a genetic dataset with known population structure, and the research samples are projected into the reference space using the LASER algorithm to infer the population structure of the research samples.

9. Use of the DNA methylation chip-based genotyping method according to any one of claims 1 to 6 for inferring kinship.

10. The use according to claim 9, characterized in that For samples from a single population and a mixed population, the SEEKIN-hom method and SEEKIN-het method were used to infer kinship, respectively.

Images