Method for constructing gene regulatory network by utilizing eQTL and application

Abstract

The invention discloses a method for constructing a gene regulatory network by utilizing eQTL and application, and the method specifically comprises the following steps: utilizing eQTL data, calculating linkage unbalance degree between eSNP in each eQTL and all other SNPs in a 1 million base pair range of a chromosome where the eSNP is located through an eQTL2GRN program, classifying according to the linkage unbalance degree, and utilizing physical coordinate positions of the eSNP and eGene to construct a gene regulatory network. Judging that the regulation and control type is cis or trans; and finally, constructing a gene regulation and control network, and carrying out visualization. By means of the method, the regulation and control relation between the genes can be recognized and analyzed, and then how genetic variation affects gene expression and phenotypic characteristics can be deeply understood.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics and molecular breeding, and more specifically, to a method and application of constructing a gene regulatory network using quantitative trait expression sites. Background Art

[0002] Gene regulatory network (GRNs) is an intrinsic regulatory structure formed by the interaction between transcription factors and genes in cells or specific genomes, and plays an important role in the expression regulation mechanism in cells. GRNs are complex dynamic systems that include interactions between DNA, RNA, proteins and other small molecules, reflecting the interaction of various biochemical substances in the gene regulation process at the molecular level, thereby regulating gene expression (Huang CD, Cao JD, Xiao M. Hybrid control on bifurcation for a delayed fractional gene regulatory network [J]. Chaos, Solitons and Fractals. 2016 Feb; 87: 19-29). These regulatory networks provide key clues for the formation of important agronomic traits of crops.

[0003] With the development of high-throughput transcriptome sequencing technology, a large amount of gene expression data can be obtained from different germplasm materials with high throughput (Delgado FM, Gómez-Vela F. Computational methods for Gene Regulatory Networks reconstruction and analysis: A review. Artificial Intelligence in Medicine. 2019

[0004] Apr; 95: 133-145). Expression Quantitative Trait Loci (eQTL) uses gene expression as a molecular phenotype and reveals the regulatory relationship between gene expression and genetic variation through association analysis. At present, eQTL has become a powerful tool for studying the association between gene expression variation and genetic structure, revealing the molecular connection between genetic variation and phenotypic diversity (Hormozdiari F, Gazal S, van de Geijn B et al.. Leveraging molecular quantitative trait loci to understand the genetic architecture of diseases and complex traits. Nature Genetics. 2018 Jul; 50 (7): 1041-1047). In crops such as rice, cotton, and tomatoes, eQTL positioning has made significant progress. Liu et al. (2022) sequenced the leaf transcriptomes of 287 rice varieties from all over the world, and then conducted a whole-genome association analysis study, and found 44,354 eQTLs that regulate the expression of 13,201 genes (Liu C, Zhu XY, Zhang J et al..eQTLs play critical roles in regulating gene expression and identifying key regulators in rice. Plant Biotechnology Journal. 2022Dec;20(12):2357-2371). Ma et al. (2021) sequenced the anther transcriptomes of 218 cotton varieties grown under high temperature stress, and identified 32,396 eQTLs through association analysis, which were involved in the regulation of 8,975 protein-coding genes and 4,157 long non-coding RNAs, covering approximately 14% of the genes in the cotton genome (Ma YZ, Min L, Wang JD. A combination of genome-wide and transcriptome-wide association studies reveals genetic elements leading tomale sterility during high temperature stress in cotton. New Phytologist. 2021 Jul; 231(1): 165-181).Ran et al. (2017) used introduced tomato populations developed from wild tomatoes and domesticated tomatoes in desert environments to analyze the impact of each introduction on the transcriptome, identifying approximately 7,200 eQTLs and regulating 5,300 genes (Ranjan A, Budke JM, Rowland SD et al..eQTL regulating transcript levels associated with diverse biological processes into mato. Plant Physiology. 2016 Sep; 172(1): 328-340).

[0005] In addition to mining candidate genes co-localized with associated sites, eQTL analysis can also be used to construct gene regulatory networks. eQTL hotspots are the manifestation of gene networks (Wang J, Yu HH, Weng XY et al.. An expression quantitative trait loci-guided co-expression analysis for constructing regulatory network using a rice recombinant inbred line population. Journal of Experimental Botany. 2014 Mar; 65 (4): 1069-1079). When genetic variation occurs in regulatory regions, such as enhancers or transcription factors, it will act on a large number of downstream target genes in a trans-mediated manner. Pierce et al. (2014) demonstrated that Trans-eQTL can first affect the expression of Cis-gene, and the expression product of Cis-gene then affects the expression of distal Trans-gene. This regulatory mode is called eQTL mediation effect (Mediation Effect) (Pierce BL, Tong L, Chen LS et al..Mediation analysis demonstrates that trans-eQTLs are often explained by cis-mediation: Agenome-wide analysis among 1,800 South Asians. PLOS Genetics. 2014 Dec; 10(12): e1004818.).Li et al. (2020) sequenced the fiber transcriptomes of 251 cotton varieties, identified 15,330 eQTL sites, and used 216 eQTL hotspots to construct a genome-wide genetic network that regulates the expression of 962 genes, which mainly regulates the synthesis of secondary cell wall of cotton fibers and contributes to the formation of cotton fibers (Li ZH, Wang PC, You CY et al.. Combined GWAS and eQTL analysis uncovers a genetic regulatory network orchestrating the initiation of secondary cell wall development in cotton. New Phytologist. 2020 Jun; 226(6): 1738-1752).Currently, the commonly used software and tools for eQTL analysis include Matrix eQTL (Shabalin, AA (2012). Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics, 28 (10), 1353-1358), FastQTL (Ongen, H., Brown, AA, Delaneau, O., Panousis, NI, Nica, AC, & Dermitzakis, ET (2016). Fastand efficient QTL mapper for thousands of molecular phenotypes. Bioinformatics, 32 (10), 1553-1555), and PLINK (Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, MA, Bender, D., ... & Sham, PC (2007). "PLINK: a toolset for whole-genome association and population-based linkage analyses." American Journal of Human Genetics, 81(3), 559-575.) and QTLtools (Delaneau, O., Ongen, H., Brown, AA, Fort, A., Panousis, NI, & Dermitzakis, ET (2016). "A complete tool set for molecular QTL discovery and analysis." Nature Communications, 7, 12474.). However, there is still a lack of a one-step program for constructing gene regulatory networks based on eQTL. Summary of the invention

[0006] The purpose of the present invention is to address the deficiencies of the prior art and propose a method and application of constructing a gene regulatory network using eQTL. Through this method, the regulatory relationship between genes can be identified and analyzed, thereby gaining a deeper understanding of how genetic variation affects gene expression and phenotypic characteristics.

[0007] To achieve the above object, the technical solution of the present invention is: the present invention provides a method for constructing a gene regulatory network using eQTL, the method comprising the following steps:

[0008] (1) Obtaining eQTL data, including associated genes, eSNPs, and significant P values ​​of the association between the associated genes and the regulated eSNPs;

[0009] (2) Plink software was used to calculate the degree of linkage disequilibrium between the eSNP in each eQTL and all other SNPs within 1 million base pairs of the chromosome where it was located. The genes regulated by eSNPs with a linkage disequilibrium > 0.1 were grouped into one cluster. The physical coordinate positions of the eSNP and eGene were used to determine whether the regulation type was cis or trans. The gene regulatory network was constructed.

[0010] Furthermore, the acquisition of eQTL data is specifically as follows: based on the genotype resequencing data, the samples are typed using the GATK software, and the obtained genotypes are stored in the VCF format; the gene expression levels of the samples are obtained using the hisat2 software, and stored in the gene expression matrix format, with each row representing a gene and each column representing a sample; and the eQTL list is generated by the matrixeQTL software, including associated genes, eSNPs, and the significant P values ​​of the association between the associated genes and their regulated eSNPs.

[0011] Furthermore, step (2) is achieved by inputting the eQTL list, the gene physical coordinate file stored in BED format, and the genotype file stored in VCF format into the eQTL2GRN program (Python environment) to generate a gene regulatory network file and a network visualization file. Each line of the gene regulatory network file represents a regulatory network and the genes it contains, and indicates the eSNPs and regulatory types that regulate the genes; the network visualization file contains each pair of genes, regulated eSNPs, regulatory types, regulatory networks to which they belong, and the significant P value of the association between the genes and their regulated eSNPs.

[0012] Furthermore, it also includes using Cytoscape software to visualize the network of the network visualization file.

[0013] Furthermore, the degree of linkage disequilibrium between the eSNP in each eQTL and all other SNPs within 1 million base pairs of the chromosome where it is located is calculated, specifically: for two genetic loci A and B on the chromosome, if there are alleles A1, A2 and B1, B2 respectively, the frequency of each allele and the frequency of the related genotype are calculated to obtain the difference between the actual frequency of the allele combination and the expected independent frequency, and evaluate the degree of linkage disequilibrium:

[0014] D = P(A1B1) - P(A1)P(B1);

[0015]

[0016] Among them, r2 reflects the statistical association strength value between two loci, which is used as a measure of the degree of linkage disequilibrium and ranges from 0 to 1. The genes regulated by eQTL are then classified according to the degree of linkage disequilibrium.

[0017] The beneficial effects of the present invention are as follows:

[0018] This invention provides a program based on Python programming language to analyze eQTL and construct gene regulatory network. This program reduces the operational and technical barriers faced by professionals engaged in genome analysis when processing eQTL data, and also simplifies the implementation process of eQTL in a wider range of application fields. Through this innovative technology, users can more efficiently and accurately analyze complex gene regulatory relationships, which in turn helps to deeply understand gene functions and biological processes. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 It is a workflow diagram for constructing gene regulatory networks using eQTL;

[0020] Figure 2 This is a diagram illustrating the parameters of the eQTL2GRN program; in the figure:

[0021] The eQTL2GRN program parameters include:

[0022] -h, --help: Displays program parameter description information, which is not required for running;

[0023] -v, --vcf: add the input genotype VCF file, run the required parameters;

[0024] -e, --eQTL_file: add the input eQTL file, necessary parameters for running;

[0025] -b, --bed_file: add the input gene annotation BED file, necessary parameters for running;

[0026] -o1, --output_file1: output gene regulatory network result file, not a required parameter for running, the default is GRN.txt;

[0027] -o2, --output_file2: Output network visualization file, not a required parameter for running, the default is Visualize.txt.

[0028] Figure 3 is the output gene regulatory network file; in the figure:

[0029] The first column is the cluster number of the classification, the second column is the number of genes included in the cluster, and the third column is the genes included in the cluster, where the genes and each eSNP (expression single nucleotide polymorphism site) regulated by them are listed one by one, and the type of regulation (cis or trans) is marked;

[0030] Figure 4 This is the output network visualization file format diagram; in the figure:

[0031] The first column is the classification cluster number, the second column is the associated gene, the third column is the eSNP (expressed single nucleotide polymorphism site), the third column is the type of eSNP-regulated gene, which is divided into cis and trans; the third column is the significant P value (P-value) of the association between the gene and eSNP. DETAILED DESCRIPTION

[0032] Reference will now be made in detail to exemplary embodiments, examples of which are illustrated in the accompanying drawings.

[0033] Embodiment 1:

[0034] Using eQTL information, the linkage disequilibrium (LD) coefficient between eQTL sites was calculated by writing the program eQTL2GRN to construct a gene regulatory network. The program achieves this goal based on the mutual coordinate relationship between eSNP (expression-related single nucleotide polymorphism) and eGene (gene whose expression is affected by variation) and the calculation of linkage disequilibrium coefficient between eSNPs.

[0035] 1. Acquisition of cotton 1-DPA fiber eQTL:

[0036] eQTL (Expression Quantitative Trait Loci) research is mainly used to identify genetic variations that affect gene expression levels, helping researchers understand how genes affect phenotypes. eQTL analysis results usually include eSNPs (expression-related single nucleotide polymorphisms) and eGenes (genes whose expression is affected by the variation), as well as their associated statistical test values.

[0037] The initial number of cotton fibers is significantly correlated with fiber yield, and 1 day after flowering (1DPA) is the critical period for fiber initiation. In this study, we used the core germplasm of upland cotton as the research object and conducted field statistics on fiber yield. At the same time, transcriptome sequencing and genotype sequencing analysis were performed on the fibers of the corresponding germplasm at the 1DPA stage. First, based on the genotype resequencing data, we used GATK (Genome Analysis Toolkit) software to type the genotype data of each sample and stored these data in VCF (Variant Call Format) format. This format facilitates subsequent genetic variation analysis and association studies. For the analysis of gene expression, we used hisat2 software to compare and quantitatively analyze the transcriptome data to obtain the expression level of each gene in each germplasm material at the 1DPA ovule stage. All gene expression data are stored in the gene expression matrix (expression.csv) format, in which each row represents a gene, each column represents a sample, and the corresponding value is the expression level of the gene.

[0038] Next, we converted these gene expression levels into molecular phenotypes and associated them with the genetic variation of the samples. By using the matrixeQTL software (https: / / github.com / andreyshabalin / MatrixEQTL), we were able to generate a list of eQTL (expression quantitative trait loci), such as Figure 2 The format of the list is as follows:

[0039] First column: associated genes;

[0040] Second column: eSNP (expressed single nucleotide polymorphism site);

[0041] The third column: P-value of significance; used to quantify the significance of the association between eSNP and gene expression level.

[0042] This analytical approach enabled us to reveal the association between genetic variation and gene expression, providing important clues for understanding the molecular mechanisms of cotton fiber initiation and yield, however, information on gene regulatory networks remains unknown.

[0043] 2. Use the software eQTL2GRN to build the network:

[0044] Run the program eQTL2GRN with the eQTL file, gene physical coordinate file (stored in BED format) and genotype VCF file as input. The program eQTL2GRN can be found at https: / / github.com / sxppppp / eQTL2GRN Get.

[0045] The main function of the program is to calculate the degree of linkage disequilibrium between the eSNP in each eQTL and all other SNPs within 1 million base pairs (1Mb) of the chromosome where it is located. Specifically, for two genetic loci A and B on the chromosome, if there are alleles A1, A2 and B1, B2 respectively, the program will calculate the frequency of each allele and the frequency of the related genotype. Through these calculations, the difference between the actual frequency of the allele combination and the expected independent frequency can be obtained, thereby evaluating the degree of linkage disequilibrium. The calculation method is as follows:

[0046] D=P(A1B1)-P(A1)P(B1),

[0047]

[0048] Its r2 value ranges from 0 to 1. r2 reflects the statistical association strength between two loci. In the analysis, the genes regulated by eQTL can be classified into the same network according to the degree of linkage disequilibrium.

[0049] After executing eQTL2GRN in the Python environment, eQTL2GRN can extract the eSNP list contained in the eQTL file from the genotype VCF file, calculate the LD value between eSNPs using Plink software, and group the genes regulated by eSNPs with LD>0.1 into a cluster. At the same time, eQTL2GRN will also use the physical coordinates of eSNPs and genes to determine whether the type of regulation is cis or trans, and finally generate a comprehensive result file together with the significant P value.

[0050] The generated result files include:

[0051] 1. Gene regulatory network files: such as Figure 3 As shown in the figure, each line in the file represents a regulatory network and the genes it contains, and each eSNP regulated by the gene indicates whether it is in cis or trans mode.

[0052] 2. Network visualization files: such as Figure 4 As shown, the file contains each pair of genes, regulated eSNPs, regulation types, regulatory networks, and significant P-values ​​(P-values) associated with the genes and their regulated eSNPs. The network can be visualized using Cytoscape software (https: / / cytoscape.org / ).

[0053] Through these steps, users can effectively construct and visualize gene regulatory networks.

[0054] It should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Without making creative work, any modification, equivalent replacement, improvement, etc., shall be included in the protection scope of the present invention.

Claims

1. A method for constructing a gene regulatory network using eQTL, characterized in that: The method comprises the following steps: (1) Obtaining eQTL data, including associated genes, eSNPs, and significant P values ​​of the association between the associated genes and the regulated eSNPs; (2) Plink software was used to calculate the degree of linkage disequilibrium between the eSNP in each eQTL and all other SNPs within 1 million base pairs of the chromosome where it was located. The genes regulated by eSNPs with a linkage disequilibrium > 0.1 were grouped into one cluster. The physical coordinate positions of the eSNP and eGene were used to determine whether the regulation type was cis or trans. The gene regulatory network was constructed.

2. The method according to claim 1, characterized in that The eQTL data acquisition is specifically as follows: based on the genotype resequencing data, the samples are typed using the GATK software, and the obtained genotypes are stored in the VCF format; the gene expression levels of the samples are obtained using the hisat2 software, and stored in the gene expression matrix format, where each row represents a gene and each column represents a sample; and the eQTL list is generated using the matrixeQTL software, including associated genes, eSNPs, and the significant P values ​​of the association between the associated genes and their regulated eSNPs.

3. The method according to claim 1 or 2, characterized in that: The step (2) is achieved by inputting the eQTL list, the gene physical coordinate file stored in the BED format, and the genotype file stored in the VCF format into the eQTL2GRN program to generate a gene regulatory network file and a network visualization file. Each line of the gene regulatory network file represents a regulatory network and the genes it contains, and indicates the eSNPs and regulatory types that regulate the genes; the network visualization file contains each pair of genes, regulated eSNPs, regulatory types, regulatory networks to which they belong, and the significant P value of the association between the genes and their regulated eSNPs.

4. The method according to claim 3, characterized in that: The method also includes using Cytoscape software to visualize the network on the network visualization file.

5. The method according to claim 1, characterized in that The degree of linkage disequilibrium between the eSNP in each eQTL and all other SNPs within 1 million base pairs of the chromosome where it is located is calculated as follows: for two genetic loci A and B on the chromosome, if alleles A1, A2 and B1, B2 exist respectively, the frequency of each allele and the frequency of the related genotype are calculated to obtain the difference between the actual frequency of the allele combination and the expected independent frequency, and to evaluate the degree of linkage disequilibrium: D = P(A1B1) - P(A1)P(B1); Among them, r2 reflects the statistical association strength value between two loci, which is used as a measure of the degree of linkage disequilibrium and ranges from 0 to 1. The genes regulated by eQTL are then classified according to the degree of linkage disequilibrium.

Images