Single-cell gene expression quantitative trait site identification method and system based on high-throughput multi-omics sequencing

Abstract

The invention relates to a single-cell gene expression quantitative trait site identification method and equipment based on high-throughput multi-omics sequencing. The method comprises the following steps: S1, constructing a data storage structure object; S2, standardizing a gene expression matrix; s3, filtering the data; s4, identifying single cell gene expression quantitative character sites; and S5, result visualization is carried out. The method provided by the invention realizes conjoint analysis of single cell transcriptomics data and genomics data, provides selection of a plurality of models for identifying gene expression quantitative trait sites for users, and can analyze data on a single cell level, thereby revealing association between genotypes and gene expression. And a new perspective is provided for understanding the molecular mechanism of complex diseases.

Description

Technical Field

[0001] The present invention belongs to the field of transcriptomic sequencing and genomic sequencing in high-throughput multi-omics sequencing, and specifically relates to a method and system for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing. Background Art

[0002] An expression quantitative trait locus refers to a genetic locus that can affect the amount of gene expression, where genetic variation is related to the expression level of certain genes. When a certain variation site affects gene expression, it is called a gene expression quantitative trait locus. Studies have shown that gene expression quantitative trait loci are population, tissue, and cell specific, and the gene expression quantitative trait loci in different tissues or cells are not the same. The development of single-cell RNA sequencing (scRNA-Seq) technology has greatly promoted our research on single-cell gene expression quantitative trait loci, giving us the opportunity to reveal the association between genotype and gene expression at the single-cell level. With the continuous development of the field of bioinformatics, some researchers have developed technologies for identifying gene expression quantitative trait loci, but these technologies generally have the disadvantages of using a single data model to identify gene expression quantitative trait loci and not highlighting the shortcomings of studying biological data at the single-cell level. Summary of the invention

[0003] In order to overcome the deficiencies in the prior art, the present invention aims to provide a single-cell gene expression quantitative trait locus identification technology based on high-throughput multi-omics sequencing.

[0004] The technical solution of the present invention is:

[0005] In a first aspect, the present invention provides a method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing, comprising the following steps:

[0006] S1. Construct data storage structure object: input transcriptomics data and genomics data, and construct data storage structure object;

[0007] S2. Standardized gene expression matrix: Standardize the transcriptome data to obtain a standardized gene expression matrix;

[0008] S3. Filtering data: filtering the data, wherein the data includes the variable site matrix and the normalized gene expression matrix data;

[0009] S4. Identify single-cell gene expression quantitative trait loci: Based on the filtered data, input the specified genome range, screen out the genes and variant sites within the range, pair them for data modeling, and obtain the pairing results of each variant site gene pair;

[0010] S5. Visualization of results.

[0011] Furthermore, in step S2, the gene expression data is normalized using a log normalization method. Single-cell data used for analysis often have extremely high sparsity and variability, and log normalization can help reduce such sparsity and variability, making the data easier to analyze.

[0012] Furthermore, in step S3, the filtering condition is: only consider genes whose expression level exceeds zero in cells that account for at least 10% of the total number of cells, and mutation sites that account for at least 10% of the total number of cells in each different genotype; this can effectively eliminate noise signals that may come from accidental phenomena or technical errors in a very small number of cells, ensuring that the analyzed genes and mutation sites have significant biological significance. At the same time, it also enhances the sensitivity to low-abundance genes, ensuring that even at low expression levels, genes with potential importance can be effectively captured, laying a solid foundation for subsequent in-depth analysis.

[0013] Furthermore, in step S4, before data modeling, cells are classified according to classification information of cell types or cell states stored in the constructed data storage structure object.

[0014] Furthermore, in step S4, a zero-inflated negative binomial model is used for modeling, and the formula is as follows:

[0015] F(n;μ,θ,π)=π·δ 0(n) +(1-π)·f {NB}(n；μ,θ)

[0016] Among them, μ is the mean parameter of the negative binomial distribution, θ is the shape parameter of the negative binomial distribution; π is the probability of a gene zero value. By introducing an additional zero-inflation parameter π, we can effectively distinguish between zero expression caused by technical noise or random factors and biologically true zero expression, thereby more accurately capturing the relationship between gene expression and genetic variation.

[0017] Furthermore, in step S4, for the pairing results output by the model, the likelihood ratio test is first used as a preliminary screening tool, and then further verified by the chi-square test to obtain the pairing result of each variant site gene pair - the p-value; and the p-value is adjusted using the Bonferroni correction method. Only when the corrected p-value is less than 0.05, the variant site gene pair is considered to have a statistically significant association.

[0018] Furthermore, in step S5, the results are presented using a comprehensive view that integrates the box plot, violin plot, and scatter plot, which not only improves the aesthetics of the chart, but also significantly enhances the intuitiveness and clarity of the results.

[0019] In a second aspect, the present invention provides a single-cell gene expression quantitative trait locus identification system based on high-throughput multi-omics sequencing, comprising:

[0020] Data input module, used to input transcriptomics data and genomics data and construct data storage structure objects;

[0021] The normalization module is used to normalize the transcriptome data and obtain a standardized gene expression matrix;

[0022] Quality control module, used to filter variant site matrix and standardize gene expression matrix data;

[0023] The site identification module, based on the filtered data, inputs the specified genome range, screens out the genes and variant sites within the range, pairs them for data modeling, and obtains the pairing results of each variant site gene pair;

[0024] Visualization module, which visualizes the pairing results.

[0025] In a third aspect, the present invention further provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the above-mentioned method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing.

[0026] In a fourth aspect of the present invention, a computer device is provided, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing is implemented.

[0027] Compared with the prior art, the present invention has the following beneficial effects:

[0028] The method of the present invention realizes the joint analysis of single-cell transcriptomics data and genomics data, provides users with a variety of model options for identifying quantitative trait loci of gene expression and can analyze data at the single-cell level, thereby revealing the relationship between genotype and gene expression, providing a new perspective for understanding the molecular mechanisms of complex diseases. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 is a technical flow chart of the identification system of the present invention;

[0030] Figure 2 This is a schematic diagram of the partial operation results of the modeling function of Example 2 of the present invention;

[0031] Figure 3 This is an example of a gene expression level graph of Example 2 of the present invention. DETAILED DESCRIPTION

[0032] The present invention provides a method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing, comprising the following steps:

[0033] S1. Construct data storage structure object: input transcriptomics data and genomics data, and construct data storage structure object.

[0034] In order to support subsequent analysis work and effectively organize and process relevant data of single-cell gene expression quantitative trait site analysis, the present invention first creates a data storage structure to store user-input transcriptomics data and genomics data and subsequent operation results in the structure object, so that the data and operations can be managed and processed as a whole to facilitate subsequent analysis.

[0035] Furthermore, by constructing the transcriptomics data input by the user as a data storage structure object rather than a gene expression matrix, the cell type information or cell state information stored in the object can be extracted, making it easier to classify cells for subsequent analysis.

[0036] S2. Normalized gene expression matrix.

[0037] By standardizing the transcriptome data and obtaining a standardized gene expression matrix, the comparability and reliability of the data can be improved, making subsequent analysis results more accurate.

[0038] Preferably, the log normalization method is used to standardize the gene expression matrix data. Log normalization is a method for reducing data variability by taking the logarithm of gene expression. It is commonly used to compress a large range of data into a smaller interval, making low-expression and high-expression genes more visually balanced. The single-cell data used for analysis often have extremely high sparsity and variability. Log normalization can help reduce this sparsity and variability, making the data easier to analyze. The standardized results are stored in the data storage structure object created in the first step.

[0039] S3. Filter data.

[0040] By filtering, low-quality variant site matrices and standardized gene expression matrix data are filtered out, providing a basis for subsequent faster and more effective identification of gene expression quantitative trait sites.

[0041] Preferably, the filtering condition is: only consider genes whose expression level exceeds zero in cells that account for at least 10% of the total number of cells and variant sites that account for at least 10% of the total number of cells in each different genotype. Selecting 10% as the threshold can effectively eliminate noise signals that may come from accidental phenomena or technical errors in a very small number of cells, ensuring that the analyzed genes and variant sites have significant biological significance. At the same time, this strategy also enhances the sensitivity to low-abundance genes, ensuring that potentially important genes can be effectively captured even at low expression levels, laying a solid foundation for subsequent in-depth analysis.

[0042] S4. Identification of single-cell gene expression quantitative trait loci.

[0043] Based on the filtered data, the specified genome range is input, the genes and variant sites within the range are screened out, and the pairs are used for data modeling to obtain the pairing results of each variant site gene pair.

[0044] According to the genome range input by the user, the Ensembl database can be connected to obtain position information, and genes and variant sites within the base range defined by the start and end positions of the specified gene or the start and end positions of all genes can be matched.

[0045] Preferably, before data modeling, cells are classified according to classification information of cell types or cell states stored in the constructed data storage structure object.

[0046] Preferably, this step uses a zero-inflated negative binomial model for data modeling. The formula is as follows:

[0047] F(n;μ,θ,π)=π·δ 0(n) +(1-π)·f {NB}(n；μ,θ)

[0048] Where μ is the mean parameter of the negative binomial distribution, θ is the shape parameter of the negative binomial distribution, and π is the probability of the gene being zero.

[0049] Since single-cell RNA sequencing (scRNA-seq) data are generally highly sparse, traditional statistical models often face challenges in processing such data. The zero-inflated negative binomial model (ZINB) effectively distinguishes zero expression caused by technical noise or random factors from biologically true zero expression by introducing an additional zero-inflated parameter π, and can more accurately capture the relationship between gene expression and genetic variation.

[0050] Preferably, for the pairing results output by the model, the likelihood ratio test is first used as a preliminary screening tool, and then further verified by the chi-square test to obtain the pairing result of each variant site gene pair - the p-value; and the Bonferroni correction method is used to adjust the p-value. Only when the corrected p-value is less than a certain threshold, the variant site gene pair is considered to have a statistically significant association.

[0051] The running results of this step include the p-value, corrected p-value, beta value and other information of the variant site gene pair. If the corrected p-value obtained by model fitting for a variant site gene pair is closer to zero, it means that the effect of the variant site on the expression of the gene is more significant, and this variant site is very likely to be a gene expression quantitative trait site.

[0052] Preferably, only when the corrected p-value is less than 0.05 is the variant site gene pair considered to have a statistically significant association.

[0053] For example, in a specific embodiment, cells can be first divided into three groups according to genotype: AA, Aa, and aa. Next, the zero-inflated negative binomial (ZINB) model is applied to fit the gene expression data of these three groups of cells, and the model parameters unique to each group are accurately extracted and derived. Subsequently, the likelihood ratio test (LRT) is used as a preliminary screening tool, and further verified by the subsequent chi-square test. Based on these carefully calculated parameters, we obtain the p-value of each variant site gene pair. Finally, in order to control the false positive risk caused by multiple hypothesis testing, the Bonferroni correction method is used to adjust the p-value. Only when the corrected p-value is less than 0.05, the variant site gene pair is considered to have a statistically significant association.

[0054] S5. Visualization of results.

[0055] The user can select the possible variant site gene pairs of the gene expression quantitative trait loci for result display based on the results of identifying the gene expression quantitative trait loci in the previous step. If the cells are classified into cell types or cell states in the previous step, the result graph of this step can output the result graphs of different cell types or cell states, so as to compare the differences in different cell types or cell states.

[0056] Preferably, using a comprehensive view that integrates box plots, violin plots, and scatter plots to present the results can not only improve the aesthetics of the chart, but also significantly enhance the intuitiveness and clarity of the results.

[0057] In addition, refer to Figure 1 The present invention also provides a single-cell gene expression quantitative trait locus identification system based on high-throughput multi-omics sequencing, comprising:

[0058] Data input module, used to input transcriptomics data and genomics data and construct data storage structure objects;

[0059] The normalization module is used to normalize the transcriptome data and obtain a standardized gene expression matrix;

[0060] Quality control module, used to filter variant site matrix and standardize gene expression matrix data;

[0061] The site identification module, based on the filtered data, inputs the specified genome range, screens out the genes and variant sites within the range, pairs them for data modeling, and obtains the pairing results of each variant site gene pair;

[0062] Visualization module, which visualizes the pairing results.

[0063] The system can execute the above-mentioned identification method to identify quantitative trait loci of single-cell gene expression.

[0064] The present invention also provides a computer-readable storage medium, which stores a computer program. When the computer program is executed by a processor, it implements the above-mentioned single-cell gene expression quantitative trait site identification method based on high-throughput multi-omics sequencing.

[0065] The present invention also provides a computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing is implemented.

[0066] In order to better demonstrate the technical scheme, purpose and advantages of the present invention, the present invention is further described in detail in conjunction with the following specific examples and accompanying drawings. The process, conditions, experimental methods, etc. for implementing the present invention, except for the contents specifically mentioned below, are all common knowledge and common common sense in the art, and the present invention has no special restrictions. The biochemical reagents used in the following examples are all commercially available.

[0067] Example 1 Comparison of different modeling models

[0068] This example uses genomic and transcriptomic data from acute myeloid leukemia (AML) patient samples to perform single-cell eQTL site identification analysis and compare the zero-inflated negative binomial model (ZINB) with the Poisson model and linear model. The specific steps are as follows:

[0069] 1. First, input transcriptomics data and genomics data. The transcriptomics data covers the expression of 2705 cells and their corresponding 200 genes. The genomics data contains the information of 10,000 variant sites in 2705 cells. The species information is specified as "human". Create a data storage structure named eQTLObject.

[0070] 2. Then the original gene expression matrix in the eQTLObject is standardized using the log normalization method, and the standardized gene expression matrix is ​​stored in the eQTLObject.

[0071] 3. Then, the genes and variant sites were filtered, and only genes whose expression levels exceeded zero in at least 10% of the total cells and variant sites whose expression levels in each genotype accounted for at least 10% of the total cells were considered.

[0072] 4. Then, the zero-inflated negative binomial model, Poisson model and linear model were used to pair the filtered genes and variant sites. The likelihood ratio test (LRT) was used as a preliminary screening tool, and the subsequent chi-square test was used to further verify the p-value of each variant site gene pair. The p-value was adjusted using the Bonferroni correction method, and finally a data frame containing the p-value of the pairing result of the variant site gene pair and the corrected p-value was returned.

[0073] The zero-inflated negative binomial model formula is as follows:

[0074] F(n;μ,θ,π)=π·δ 0(n) +(1-π)·f {NB}(n；μ，θ) , where μ is the mean parameter of the negative binomial distribution, θ is the shape parameter of the negative binomial distribution, and π is the probability of the gene being zero.

[0075] The results used the corrected p-value less than 0.05 as the screening criterion. When the corrected p-value of a site meets this condition, we can reasonably infer that the site is very likely to be a quantitative trait locus (eQTL) that affects gene expression. Based on this criterion, the number of single-cell gene expression quantitative trait loci of different modeling models was counted and obtained, as shown in the following table:

[0076] Table 1 Comparison of the results of three modeling methods

[0077] Zero-inflated negative binomial model Poisson Model Linear Model Gene expression quantitative trait loci (number) 828 567 376

[0078] By comparison, it can be seen that the zero-inflated negative binomial model is more suitable for single-cell identification of gene expression quantitative trait sites. Although the Poisson model and the linear model are classic tools for identifying gene expression quantitative trait sites (eQTLs), their effectiveness appears relatively limited in the face of the unique challenges of single-cell data-high sparsity and a wide range of zero-value distributions. In contrast, the zero-inflated negative binomial model stands out with its excellent adaptability. The model not only effectively distinguishes zero expression from technical noise or random fluctuations from true zero expression in a biological sense by cleverly introducing the zero-inflated parameter π, but also greatly enhances the ability to capture the subtle relationship between gene expression and genetic variation in complex single-cell data. Therefore, when the zero-inflated negative binomial model is used for analysis, a larger number of gene expression quantitative trait sites can be identified, which not only verifies the superiority of the model when processing single-cell data, but also further emphasizes its indispensable role in accurately analyzing the gene expression regulatory mechanism. Therefore, the zero-inflated negative binomial model is selected in the embodiment of the present invention for modeling.

[0079] Example 2

[0080] Based on the results of Example 1, the genomics and transcriptomics data of acute myeloid leukemia (AML) patient samples were used to perform single-cell eQTL site identification analysis using a zero-inflated negative binomial model. The specific steps are as follows:

[0081] 1. First, input transcriptomics data and genomics data, specify the species information as "human", define the cell grouping information as "celltype", and create a data storage structure named eQTLObject.

[0082] 2. Then the original gene expression matrix in the eQTLObject is standardized using the log normalization method, and the standardized gene expression matrix is ​​stored in the eQTLObject.

[0083] 3. Then, the genes and variant sites were filtered, and only genes whose expression levels exceeded 0 in at least 10% of the total cells and variant sites whose expression levels in each genotype accounted for at least 10% of the total cells were considered.

[0084] 4. Then, the zero-inflated negative binomial model was used to pair the filtered genes and variant sites, and the likelihood ratio test (LRT) was used as a preliminary screening tool, and further verified by the subsequent chi-square test to obtain the p-value of each variant site gene pair. The p-value was adjusted using the Bonferroni correction method, and finally a data frame containing the p-value of the variant site gene pair pairing result and the corrected p-value was returned, and a threshold of 0.05 was used to limit the corrected p-value in the data frame. Figure 2Some results of the zero-inflated negative binomial model are given in Figure 1. Only when the corrected p-value is less than 0.05, the variant site gene pair is considered to have a statistically significant association.

[0085] 5. Finally, based on the results in the data frame, select the variant site gene pairs of interest to visualize the gene expression levels. Figure 3 The gene expression level diagram of some of them is shown in FIG.

[0086] The technical features of the above-described embodiments may be arbitrarily combined. To make the description concise, not all possible combinations of the technical features in the above-described embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0087] The above-mentioned embodiments only express several implementation methods of the present invention, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the invention patent. It should be pointed out that, for ordinary technicians in this field, several variations and improvements can be made without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the patent of the present invention shall be subject to the attached claims.

Claims

1. A method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing, characterized in that: The following steps are involved: S1. Construct data storage structure object: input transcriptomics data and genomics data, and construct data storage structure object; S2. Standardized gene expression matrix: Standardize the transcriptome data to obtain a standardized gene expression matrix; S3. Filtering data: filtering the data, wherein the data includes the variable site matrix and the normalized gene expression matrix data; S4. Identify single-cell gene expression quantitative trait loci: Based on the filtered data, input the specified genome range, screen out the genes and variant loci within the range, pair them for data modeling, and obtain the pairing results of each variant locus gene pair; S5. Visualization of results.

2. The method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing according to claim 1, characterized in that: In step S2, the gene expression data is normalized using the log normalization method.

3. The method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing according to claim 1, characterized in that: In step S3, the filtering condition is: only consider genes whose expression levels exceed zero in cells that account for at least 10% of the total number of cells, and only consider variant sites that account for at least 10% of the total number of cells in each different genotype.

4. The method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing according to claim 1, characterized in that: In the step S4, before data modeling, cells are classified according to the classification information of cell types or cell states stored in the constructed data storage structure object.

5. The method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing according to claim 1, characterized in that: In step S4, a zero-inflated negative binomial model is used for modeling, and the formula is as follows: F(n;μ,θ,π)=π·δ 0(n) +(1-π)·f {NB}(n；μ，θ) Where μ is the mean parameter of the negative binomial distribution, θ is the shape parameter of the negative binomial distribution, and π is the probability of the gene being zero.

6. The method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing according to claim 5, characterized in that: In step S4, for the pairing results output by the model, the likelihood ratio test is first used as a preliminary screening tool, and then further verified by the chi-square test to obtain the pairing result of each variant site gene pair - the p-value; and the p-value is adjusted using the Bonferroni correction method. Only when the corrected p-value is less than 0.05, the variant site gene pair is considered to have a statistically significant association.

7. The method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing according to claim 1, characterized in that: In step S5, a comprehensive view of a fused box plot, a violin plot, and a scatter plot is used to display the results.

8. A single-cell gene expression quantitative trait locus identification system based on high-throughput multi-omics sequencing, characterized in that: include: Data input module, used to input transcriptomics data and genomics data and construct data storage structure objects; The normalization module is used to normalize the transcriptome data and obtain a standardized gene expression matrix; Quality control module, used to filter variant site matrix and standardize gene expression matrix data; The site identification module, based on the filtered data, inputs the specified genome range, screens out the genes and variant sites within the range, pairs them for data modeling, and obtains the pairing results of each variant site gene pair; Visualization module, which visualizes the pairing results.

9. A computer-readable storage medium, characterized in that: The storage medium stores a computer program, and when the computer program is executed by the processor, the method for identifying single-cell gene expression quantitative trait loci based on high-throughput multi-omics sequencing as described in any one of claims 1 to 7 above is implemented.

10. A computer device, characterized in that: The invention comprises a memory, a processor and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the method for identifying quantitative trait loci of single-cell gene expression based on high-throughput multi-omics sequencing as described in any one of claims 1 to 7 is implemented.

Images