A low-depth sequencing population genotype filling calculation memory optimization method

Abstract

The application discloses a low-depth sequencing population genotype filling calculation memory optimization method, which comprises the following steps: carrying out block division on a single nucleotide polymorphism site (SNP), and setting a checkpoint value according to the block division; and carrying out calculation on a forward auxiliary variable and a backward auxiliary variable according to the set checkpoint value, and calculating a haplotype observation conditional probability according to the block division. The application has the advantages of simple principle, effective reduction of calculation memory consumption, improvement of genotype filling efficiency and the like.

Description

Technical Field

[0001] This invention mainly relates to the fields of gene sequencing and bioinformatics analysis technology, specifically a method for optimizing computational memory for filling genotypes in low-depth sequencing populations. Background Technology

[0002] With the decreasing cost of sequencing, large-scale population low-depth whole-genome sequencing has gradually become an indispensable technique in population genetics research. By performing low-depth sequencing and genotyping on a large number of samples, denser molecular markers can be obtained at the same or even lower cost than traditional microarray genotyping and simplified genome sequencing, thereby improving the efficiency of GWAS localization. However, due to the very low sequencing depth, this method inevitably results in the loss of a large number of single nucleotide polymorphism (SNP) genotypes, which in turn affects the accuracy of downstream genetic map construction, QTL localization, and population GWAS studies. Therefore, improving the integrity of genotypes at the population level is of great significance.

[0003] The basic principle of population genotype imputation is to infer the missing genotypes in the samples using statistical methods by using known haplotype reference datasets or shared haplotype information among a large number of samples in a population.

[0004] The existing steps for genotyping using known haplotype reference data are as follows:

[0005] 1) Perform haplotype phasing on the target sample to obtain haplotype information within the chromosomal region;

[0006] 2) Fill in the missing genotypes by referencing the haplotype set.

[0007] The methods described above require a pre-built reference haplotype dataset, which is often unavailable for species with limited haplotype mapping, thus rendering the method unusable. Furthermore, the accuracy of genotype imputation is affected when the target population differs from the population from which the reference haplotype was constructed. In 2016, researchers such as Davies developed the STITCH method, which models ancestral haplotypes using shared read information from sequencing data of a large number of individuals without relying on a reference haplotype. This method enables the inference of missing genotypes for each individual, solving the genotype imputation problem for species without reference haplotype data.

[0008] Another method, STITCH, models ancestral haplotypes in a population using a first-order Hidden Markov Model (HMM), treating chromosomes as chimeras containing K unknown ancestral haplotypes (for outcrossed populations, these haplotype sets can be considered to capture complete information on all different haplotypes within a region, so K may be large). This method iteratively updates the HMM model parameters using the EM algorithm. The E-step calculates the observed probability of an ancestral haplotype for each sample, and the M-step updates the ancestral haplotype frequency and other HMM parameters based on the observed probabilities of haplotypes from multiple samples. Computationally, the time complexity of a single iteration of this algorithm is related to the number of samples N, the number of SNPs T, and the number of ancestral haplotypes K. 2 The relationship is linear, and the computational bottleneck lies in the E-step, which calculates the observed probability of the ancestral haplotype for each sample.

[0009] The E-step calculation of the ancestral haplotype observation probability for each sample can be categorized as the evaluation problem among the three fundamental problems of the HMM model, namely, given the model λ = (Λ, B, π) and the observation sequence O = o1, o2, ..., o T Calculate the conditional probability p{O|λ}. If we directly perform probability calculations on this problem, the time complexity is {K}. 2} T This is clearly not feasible. A common solution to this problem is to use the forward-backward algorithm. This forward-backward algorithm effectively reduces the time complexity by introducing auxiliary variables. Its calculation steps are as follows:

[0010] 1. Calculate the forward auxiliary variable α t (i):

[0011] Forward auxiliary variable α t (i) is defined as a partial observation sequence o1, o2, ... o t The probability of termination in state i is mathematically represented as α. t (i) = p{o1, o2, ... o t q t =i|λ}, then we have:

[0012]

[0013] Where α1(i)=π i b i (o1), 1≤i≤K 2 b i (o t+1 Let α be the emission probability. ij π represents the hidden state transition probability. i Let α be the initial probability of state i, j represent the j-th ancestral haplotype combination, and t represent the t-th SNP locus. Using this recursive formula, α can be calculated. T(i), 1≤i≤K 2 ;

[0014] 2. Calculate the backward auxiliary variable β t (i):

[0015] backward auxiliary variable β t (i) is defined as a partial observation sequence o t+1 o t+2 ,...o T The probability of starting from state i is mathematically expressed as p{o t+1 o t+2 ,...o T |q t =i, λ}, similarly calculate β t The recursive formula for (i) is:

[0016]

[0017] Where β T (i)=α T (i), 1≤i≤K 2 b j (o t+1 Let α be the emission probability. ij Let be the hidden state transition probability, j represent the j-th ancestral haplotype combination, and t represent the t-th SNP site.

[0018] Then we have α t (i)β t (i)=p{O,q t =i|λ}, 1≤i≤K 2 Since 1 ≤ t ≤ T, we have:

[0019]

[0020] 3. The auxiliary variable α introduced in the above probability calculation process. t (i) and β t (i) While it can effectively reduce computational time complexity, it places very high demands on computer memory. Simple calculations show that for a chromosome containing 3 million SNP loci and 30 ancestral haplotypes, a single α... t (i) The storage of the matrix will consume 20G of memory, and the memory consumed by a single sample in the E-step will reach more than 60GB. This makes it impossible for this method to process multiple samples in parallel due to memory constraints in practical applications, and thus it cannot make full use of the CPU computing resources of many-core servers.

[0021] As can be seen from the above background, the memory consumption of the auxiliary variable matrix for a single sample is related to the number of SNPs T and the number of ancestral haplotypes K. 2The product is proportional. Considering the locality of haplotype regions, some practitioners have proposed a method to reduce memory usage. This method involves dividing the SNP set into N consecutive, non-contiguous subsets, each containing the same number of SNPs. Then, the algorithm is applied to each SNP subset. However, this method does not take into account the data dependencies between columns when calculating the auxiliary variable matrix, which affects the accuracy of genotype filling.

[0022] Another improvement proposed by practitioners is to allow adjacent sets to overlap by a fixed-length chromosome interval when dividing the SNP set. Although this improves the shortcomings of the previous method to some extent, the degree of linkage disequilibrium varies in different regions of the chromosome and different populations, and the fixed-length interval selected by humans will still inevitably lead to a loss of filling accuracy. Summary of the Invention

[0023] The technical problem to be solved by this invention is: in view of the technical problems existing in the prior art, this invention provides a simple method for optimizing the computational memory of genotype filling in low-depth sequencing populations, which can effectively reduce computational memory consumption and improve genotype filling efficiency.

[0024] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0025] A memory optimization method for genotype filling in low-depth sequencing populations, comprising:

[0026] Single nucleotide polymorphism (SNP) sites are segmented into blocks, and checkpoint values ​​are set based on these blocks.

[0027] Based on the set checkpoint values, forward and backward auxiliary variables are calculated in blocks, and the haplotype observation conditional probability is calculated.

[0028] As a further improvement to the method of the present invention: the checkpoint is determined by using the first round of forward auxiliary variable and backward auxiliary variable algorithm; based on the checkpoint value, the forward auxiliary variable and backward auxiliary variable within the block are calculated by using the second round of forward auxiliary variable and backward auxiliary variable algorithm.

[0029] As a further improvement to the method of the present invention, a block-parallel processing method is adopted in the process of calculating the forward auxiliary variable and the backward auxiliary variable in blocks according to the set checkpoint value.

[0030] As a further improvement to the method of the present invention, the process of segmenting single nucleotide polymorphism (SNP) sites includes:

[0031] Let T be the number of single nucleotide polymorphism (SNP) sites and K be the number of ancestral haplotypes. Then the size of a single forward auxiliary variable matrix is ​​T*K. 2 ;

[0032] If the block size is M, then the number of blocks is The required vector size for storing checkpoints is The size of the forward auxiliary variable for a single block is M*K 2 .

[0033] As a further improvement to the method of this invention: when all checkpoint auxiliary variable values ​​are stored in memory, the memory consumption of the forward auxiliary variable for a single block is...

[0034] As a further improvement to the method of the present invention: when At that time, peak memory is at its minimum, where M * To achieve the optimal block size, when the auxiliary variable checkpoint vector is stored on disk, the SNP block size is adjusted according to the host memory to adapt to the host memory capacity.

[0035] As a further improvement to the method of the present invention: the process of determining the checkpoint value of the forward auxiliary variable using the first round algorithm includes:

[0036] Recursively calculate the forward auxiliary variable α t (i);

[0037] Forward auxiliary variable α t (i) is defined as a partial observation sequence o1, o2, ... o t The probability of termination in state i is mathematically represented as α. t (i) = p{o1, o2, ... o t q t =i|λ}, then we have:

[0038]

[0039] Where α1(i)=π i b i (o1), 1≤i≤K 2 b i (o t+1 Let α be the emission probability. ij π represents the hidden state transition probability. i Let α be the initial probability of state i, j represent the j-th ancestral haplotype combination, and t represent the t-th SNP locus. Using this recursive formula, α can be calculated. T (i), 1≤i≤K 2 ;

[0040] Get checkpoint (0, M) * 2M * The auxiliary variable value at ...) And record it.

[0041] As a further improvement to the method of the present invention: for each block, only the forward auxiliary variable value of its first column is recorded.

[0042] As a further improvement to the method of the present invention: the process of determining the checkpoint value of the backward auxiliary variable using the first round algorithm includes:

[0043] Recursively calculate the backward auxiliary variable β t (i);

[0044] backward auxiliary variable β t (i) is defined as a partial observation sequence o t+1 o t+2 ,...o T The probability of starting from state i is mathematically expressed as p{o t+1 o t+2 ,...o T |q t =i, λ}, similarly calculate β t The recursive formula for (i) is:

[0045]

[0046] Where β T (i)=α T (i), 1≤i≤K 2 b j (o t+1 Let α be the emission probability. ij Let be the hidden state transition probability, j represent the j-th ancestral haplotype combination, and t represent the t-th SNP site.

[0047] Get checkpoints (T, ..., 2M) * M * Auxiliary variable at ) And record it.

[0048] As a further improvement to the method of the present invention: when the block is the last block, the backward auxiliary variable value of the last column of the block is recorded; otherwise, the backward auxiliary variable value of the first column of the next block is recorded.

[0049] As a further improvement to the method of the present invention: all single nucleotide polymorphism (SNP) sites are divided into blocks [0, M] in parallel. * ), [M * 2M * ), [2M * 3M * ... respectively apply the second round of forward auxiliary variable and backward auxiliary variable algorithms to calculate the haplotype observation conditional probability.

[0050] As a further improvement to the method of the present invention: for any interval [jM] * jM * +M * The calculation is as follows:

[0051] initialization For the checkpoint value, the forward auxiliary variable α is calculated. t (i) iteratively calculates the result sequentially.

[0052] initialization For the checkpoint value, apply the calculation of the backward auxiliary variable α. t (i) iteratively calculates the result sequentially.

[0053] Calculate in sequence

[0054] Compared with the prior art, the advantages of the present invention are as follows:

[0055] 1. The low-depth sequencing population genotype filling computation memory optimization method of the present invention is a low-depth sequencing population genotype filling computation memory optimization method without precision loss. By setting checkpoints and two rounds of forward-backward algorithm, it solves the precision loss caused by the traditional direct block method and realizes genotype filling without precision loss.

[0056] 2. The low-depth sequencing population genotype filling computation memory optimization method of this invention effectively reduces the memory consumption of intermediate auxiliary variables in the forward-backward algorithm through appropriate block partitioning and two rounds of forward-backward algorithm. The memory consumption of a single sample auxiliary variable is only a fraction of that before optimization.

[0057] 3. The low-depth sequencing population genotype filling computation memory optimization method of the present invention eliminates the data dependency between different SNP blocks by setting checkpoints, so that different SNP blocks of a single sample can be processed in parallel, which can make full use of the CPU computing resources of many-core servers and effectively improve the efficiency of genotype filling. Attached Figure Description

[0058] Figure 1 This is a flowchart illustrating the method of the present invention.

[0059] Figure 2 This is a schematic diagram comparing single-threaded memory consumption in a specific application example of the present invention.

[0060] Figure 3 This is a schematic diagram comparing the genotype filling results of samples after using different methods in a specific application example of the present invention. Detailed Implementation

[0061] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0062] like Figure 1 As shown, the low-depth sequencing population genotype filling computational memory optimization method of the present invention includes:

[0063] Single nucleotide polymorphism (SNP) sites are segmented into blocks, and checkpoint values ​​are set based on these blocks.

[0064] Based on the set checkpoint values, forward and backward auxiliary variables are calculated in blocks, and the haplotype observation conditional probability is calculated.

[0065] In a specific application example, the first round of forward and backward auxiliary variable algorithms is used to determine the checkpoints. Based on the checkpoint values, the second round of forward and backward auxiliary variable algorithms is used to calculate the forward and backward auxiliary variables within the block. Thus, through appropriate block partitioning and two rounds of forward-backward algorithms, the memory consumption of intermediate auxiliary variables in the forward-backward algorithm is effectively reduced, with the memory consumption of a single sample auxiliary variable being only a fraction of that before optimization.

[0066] In specific application examples, a block-parallel processing method is adopted in the process of calculating forward auxiliary variables and backward auxiliary variables in blocks according to the set checkpoint values.

[0067] By adopting the above-mentioned method of the present invention, the data dependency between different SNP blocks is eliminated by setting checkpoint values, which enables parallel processing of different SNP blocks of a single sample. This fully utilizes the CPU computing resources of many-core servers, effectively improves the efficiency of genotype filling, and ultimately optimizes the computing memory.

[0068] In specific application examples, the process of segmenting single nucleotide polymorphism (SNP) sites includes:

[0069] Let T be the number of single nucleotide polymorphism (SNP) sites and K be the number of ancestral haplotypes. Then the size of a single forward auxiliary variable matrix is ​​T*K. 2 In this embodiment, T = 1300000 and K = 24, so the memory size occupied by a single forward auxiliary variable matrix is ​​5.57 GB;

[0070] If the block size is M, then the number of blocks is The required vector size for storing checkpoints is The size of the forward auxiliary variable for a single block is M*K 2 ;

[0071] In a specific application example, if all checkpoint auxiliary variable values ​​are stored in memory, the memory required for the forward auxiliary variable of a single block will be...

[0072] In specific application examples, when At this time, peak memory is minimized; in this embodiment, the optimal block size is 1140, and the memory consumption for storing a single block matrix of forward auxiliary variables and the checkpoint matrix is ​​only 10MB. The single-threaded memory consumption of a single auxiliary variable matrix before and after optimization is compared as follows: Figure 2 As shown; if the auxiliary variable checkpoint vector is stored on disk, the SNP block size can be flexibly adjusted according to the host memory to adapt to the host memory capacity, which is highly flexible.

[0073] In a specific application example, the process of determining the checkpoint value of the forward auxiliary variable using the first round algorithm includes:

[0074] Recursively calculate the forward auxiliary variable α t (i);

[0075] Forward auxiliary variable α t (i) is defined as a partial observation sequence o1, o2, ... o t The probability of termination in state i is mathematically represented as α. t (i) = p{o1, o2, ... o t q t =i|λ}, then we have:

[0076]

[0077] Where α1(i)=π i b i (o1), 1≤i≤K 2 b i (o t+1 Let α be the emission probability. ij Given the hidden state transition probability, this recursive formula can be used to calculate α. T (i), 1≤i≤K 2 ;

[0078] Get checkpoint (0, M) * 2M * The auxiliary variable value at ...) And record it.

[0079] In a specific application example, during the above process, for each block, only the value of the forward auxiliary variable in its first column is recorded.

[0080] In a specific application example, the process of determining the checkpoint value of the backward auxiliary variable using the first round algorithm includes:

[0081] Recursively calculate the backward auxiliary variable β t(i);

[0082] backward auxiliary variable β t (i) is defined as a partial observation sequence o t+1 o t+2 ,...o T The probability of starting from state i is mathematically expressed as p{o t+1 o t+2 ,...o T |q t =i, λ}, similarly calculate β t The recursive formula for (i) is:

[0083]

[0084] Where β T (i)=α T (i), 1≤i≤K 2 .

[0085] Get checkpoints (T, ..., 2M) * M * Auxiliary variable at ) And record it.

[0086] In a specific application example, during the above process, when the block is the last block, the value of the backward auxiliary variable in the last column of the block is recorded; otherwise, the value of the backward auxiliary variable in the first column of the next block is recorded.

[0087] In a specific application example, all single nucleotide polymorphism (SNP) sites are segmented in parallel [0, M]. * ), [M * 2M * ), [2M * 3M * In the second round, the forward and backward auxiliary variable algorithms (forward-backward algorithm) are applied respectively to calculate the haplotype observation conditional probability.

[0088] Specifically, for any interval [jM] * jM * +M * The calculation is as follows:

[0089] initialization For the checkpoint value, the forward auxiliary variable α is calculated. t (i) iteratively calculates the result sequentially.

[0090] initialization For the checkpoint value, apply the calculation of the backward auxiliary variable α. t (i) iteratively calculates the result sequentially.

[0091] Calculate in sequence

[0092] In a specific application example of this invention, the method of this invention is applied to genotyping of low-depth mouse sequencing data to further verify the advantages of this invention compared to traditional methods:

[0093] (1) Download test data

[0094] This dataset contains 2073 samples and 1516 SNP loci. Among them, 4 samples underwent high-depth sequencing, and the genotypes of the corresponding SNP loci were known, which can be used to assess the accuracy of genotype filling.

[0095] (2) Select K=4, and use traditional method 1, traditional method 2 and the method of the present invention to fill the genotypes of the dataset. The specific steps are as follows:

[0096] a) Divide the SNP set into 15 parts, each containing 100 SNP loci (the last block contains less than 100), with no overlapping SNPs between adjacent blocks. Perform genotyping on each block sequentially and estimate the genotyping accuracy index R for 4 samples with known genotypes. 2 ;

[0097] b) Divide the SNP set into 15 parts, and add 50 SNP loci to the first and last parts of the partitioned set from step a), so that there are overlapping SNP loci between adjacent blocks. Perform genotyping on each block sequentially and estimate the genotyping accuracy index R of 4 samples with known genotypes. 2 ;

[0098] c) Set the block size M = 100, directly use the method of this invention to perform genotype filling, and estimate the genotype filling accuracy index R for 4 samples with known genotypes. 2 ;

[0099] d) The R 2 The calculation method is as follows: First, the sample genotype is encoded as {0, 1, 2}, where 0 indicates that the sample has no variation at the SNP site, 1 indicates that the sample is heterozygous, and 2 indicates that the sample is homozygous. Then, the square of the Pearson correlation coefficient between the genotype filling result and the high-depth sequencing gold standard result is calculated.

[0100] (3) Figure 3 As shown, the average R values ​​obtained by conventional method 1, conventional method 2, and the present invention are... 2 The genotyping results of four samples with scores of 0.66, 0.67, and 0.84 respectively indicate that the present invention has higher genotyping accuracy than the traditional direct segmentation method.

[0101] The above are merely preferred embodiments of the present invention. The scope of protection of the present invention is not limited to the above embodiments. All technical solutions falling within the scope of the present invention's concept are within the scope of protection of the present invention. It should be noted that for those skilled in the art, any improvements and modifications made without departing from the principles of the present invention should be considered within the scope of protection of the present invention.

Claims

1. A memory optimization method for genotype filling in low-depth sequencing populations, characterized in that, include: Single nucleotide polymorphism (SNP) sites are segmented into blocks, and checkpoint values ​​are set based on these blocks. Based on the set checkpoint values, forward and backward auxiliary variables are calculated in blocks, and the conditional probability of haplotype observations is calculated. The checkpoints are determined using the first-round forward and backward auxiliary variable algorithms. Based on the checkpoint values, the forward and backward auxiliary variables within the block are calculated using the second-round forward and backward auxiliary variable algorithms. A block-parallel processing method is adopted in the process of calculating forward and backward auxiliary variables in blocks based on the set checkpoint values. The process of segmenting single nucleotide polymorphism (SNP) sites includes: Let the number of single nucleotide polymorphism (SNP) sites be . The number of ancestral haplotypes is Then the size of a single forward auxiliary variable matrix is ; The block size is Then the number of blocks is The vector size required to store checkpoints is ; The size of the forward auxiliary variable for a single block is .

2. The memory optimization method for low-depth sequencing population genotype filling according to claim 1, characterized in that, When all checkpoint auxiliary variable values ​​are loaded into memory, the memory consumption of the forward auxiliary variable for a single block is... .

3. The computational memory optimization method for low-depth sequencing population genotype filling according to claim 1, characterized in that, when At that time, peak memory was at its minimum, of which To achieve the optimal block size, when the auxiliary variable checkpoint vector is stored on disk, the SNP block size is adjusted according to the host memory to adapt to the host memory capacity.

4. The computational memory optimization method for low-depth sequencing population genotype filling according to claim 1, characterized in that, For each block, only the value of the forward auxiliary variable in its first column is recorded.

5. The memory optimization method for low-depth sequencing population genotype filling according to claim 1, characterized in that, When the block is the last block, record the value of the backward auxiliary variable in the last column of the block; otherwise, record the value of the backward auxiliary variable in the first column of the next block.

Images