CN108004302A - A kind of association analysis method of transcript profile reference and its application - Google Patents

Abstract

The present invention relates to a method for association analysis of transcriptome reference, comprising: 1). Establishing a reference population and measuring the transcriptome sequences of all individuals in the reference population to obtain reference transcripts; 2). Performing high-level analysis on the predicted population Throughput RNA sequencing, and compare and splicing the obtained sequencing data with the reference transcripts to obtain the population transcriptome sequencing results; 3). Analyze the population SNP data from the population transcriptome sequencing results, using The association analysis between the population SNP and the trait of interest is carried out to determine the region where the target gene is located; then the gene expression variation in the target region is used to correlate with the trait to obtain the candidate gene of the target trait. Compared with GWAS, transcriptome sequencing is cheap, fast and not limited by genome complexity. This technology can theoretically carry out trait association analysis in all species, regardless of whether the species has a reference genome or whether the genome is complex.

Description

A method for association analysis of transcriptome reference and its application

technical field

The present invention relates to the field of biotechnology, in particular to a transcriptome reference association analysis method and its application.

Background technique

With the development of sequencing technology and the popularization of high-density single nucleotide polymorphism (single nucleotide polymorphism, SNP) chips, genome-wide association study (GWAS) has increasingly become the focus of human disease research and animal breeding. A powerful tool. GWAS is based on the principle of linkage disequilibrium, by analyzing the correlation between tens of thousands of SNP markers in the genome and target traits in natural populations to achieve the positioning of target trait genes.

With the in-depth research on GWAS, it has gradually exposed the following defects, specifically:

First, in order to realize the identification of SNP markers in the whole genome, it is necessary to resequence the samples and compare them with the reference genome. However, only a few species have completed whole-genome sequencing and have the reference genomes required for GWAS. Thus, a major disadvantage of GWAS methods is that the application of the technique relies on the reference genome of the corresponding species. Although the recent GWAS method based on the simplified genome does not rely on the reference genome of the species, it is difficult to detect the candidate gene of the target trait because the gene information near the genetic marker is unknown, which is of little significance in biological research.

Second, GWAS technology can only determine the target gene in a genomic region, and it is difficult to directly determine the target trait candidate genes.

Third, in order to avoid statistical errors in gene mapping, GWAS analysis usually requires a larger population sample.

Fourth, GWAS technology is difficult to determine the relationship between candidate genes of target traits.

In view of this, the present invention is proposed.

Contents of the invention

The purpose of the present invention is to provide a transcriptome reference association analysis method (transcriptome-referenced association study, TRAS), in order to solve the above problems.

The present invention relates to a method for association analysis of transcriptome reference, comprising:

1). Establish a reference population and measure the transcriptome sequences of all individuals in the reference population to obtain reference transcripts;

2). Perform high-throughput RNA sequencing on the predicted population, and compare and splice the obtained sequencing data with the reference transcript to obtain the population transcriptome sequencing result;

3). Analyze the population SNP data from the population transcriptome sequencing results, use the population SNP and traits of interest to perform association analysis, and determine the region where the target gene is located; then use the gene expression variation and traits in the target region to perform Correlation to get the candidate genes of the target traits.

Unlike genome analysis, transcriptome sequencing is cheap, fast, and not limited by genome complexity. In the present invention, the full-length transcriptome of the corresponding research species is used as a reference sequence to replace the genome to carry out SNP identification and gene expression quantification of the research population. Therefore, this technology can theoretically carry out trait association analysis in all species, regardless of whether the species has a reference genome or whether the genome is very complex will not affect the implementation of the technology.

Detailed ways

The present invention relates to a method for association analysis of transcriptome reference, comprising:

1). Establish a reference population and measure the transcriptome sequences of all individuals in the reference population to obtain reference transcripts;

2). Perform high-throughput RNA sequencing on the predicted population, and compare and splice the obtained sequencing data with the reference transcript to obtain the population transcriptome sequencing result;

3). Analyze the population SNP data from the population transcriptome sequencing results, use the population SNP and traits of interest to perform association analysis, and determine the region where the target gene is located; then use the gene expression variation and traits in the target region to perform Correlation to get the candidate genes of the target traits.

In the present invention, the TRAS technology is carried out in two steps. First, it is the same as the GWAS technology, using SNP and traits for association analysis to determine the region where the target gene is located. After that, TRAS uses the gene expression of the target region to correlate with the trait, and then determines the candidate gene of the target trait.

Preferably, the above-mentioned association analysis method of transcriptome reference, said method also includes:

Determine the eQTL of the candidate gene, and then find whether the remaining target trait candidate genes are contained in the eQTL.

eQTL (Expression quantitative trait loci) takes the expression data of each genotype from a segregated population as a quantitative trait, and uses traditional QTL analysis methods for analysis. In the present invention, in order to determine the interaction relationship between candidate genes, TRAS technology uses the expression level of candidate genes (named A gene) as phenotypic data, and carries out association analysis with SNPs in the population, thereby determining the eQTL of these candidate genes. After that, find out whether the candidate gene of the target trait is contained in these eQTLs, and if so (named as B gene). Because gene B is not only a target trait candidate gene, but also an eQTL of another target trait candidate gene, it can be determined that gene A and B gene have an interactive relationship in the regulation of target traits.

Preferably, in the transcriptome-referenced association analysis method described above, the gene expression variation includes gene sequence variation and gene expression variation.

GWAS only uses one parameter (sequence variation in the population—SNP) for association analysis with traits, so in order to reduce errors, the analysis requires as large a population sample as possible. The TRAS technology not only carries out the association analysis between SNP and traits in the population, but also carries out the association analysis between gene expression variation and traits in the population. It has two parameters of sequence variation and gene expression variation, so the sample size of the analysis population does not need to be as large as that of GWAS. big.

Preferably, the method for association analysis of transcriptome reference as described above, in step 1), after obtaining the reference transcripts, further includes: performing functional annotation on the reference transcripts.

Preferably, in the above-mentioned association analysis method of transcriptome reference, in step 1), the number of species in the reference population is ≥ 1.

In the present invention, the role of the reference population is to serve as a reference for the data processing of the population transcriptome. Therefore, considering the cost, the number of reference populations should not be set too much.

Preferably, in the above-mentioned transcriptome reference association analysis method, in step 1), the transcriptome sequence determination method is single-molecule real-time sequencing.

Preferably, in the above-mentioned transcriptome reference association analysis method, in step 3), before using the population SNP to perform association analysis with the trait of interest, first perform co-expression analysis on the population transcriptome sequencing results, Transcripts exhibiting co-expression in a population are assigned to a module, and the association of said module with a trait of interest is analyzed to identify co-expressed modules associated with said trait of interest.

Preferably, in the above-mentioned transcriptome-referenced association analysis method, in step 3), when using the population SNP to perform association analysis with the trait of interest, the association analysis model used is a general linear model or a mixed linear model.

Preferably, in the above-mentioned transcriptome reference association analysis method, in step 3), when using the population SNP to perform association analysis with the trait of interest, the quality score of the SNP is ≥ 40, and the coverage depth is ≥ 2, Minor allele frequency ≥ 0.05, deletion genotype frequency ≤ 0.5.

The application of the above-mentioned method in the screening and identification of main effect genes, and the correlation analysis of gene functions.

Embodiments of the present invention will be described in detail below in conjunction with examples, but those skilled in the art will understand that the following examples are only for illustrating the present invention, and should not be considered as limiting the scope of the present invention. Those who do not indicate the specific conditions in the examples are carried out according to the conventional conditions or the conditions suggested by the manufacturer. The reagents or instruments used were not indicated by the manufacturer, and they were all commercially available conventional products.

Genome-wide association study (GWAS) is a powerful tool for identifying genes for complex traits, but its application is limited by the need to study species reference genome sequences. The present invention proposes a method, the transcriptome-referenced association study (TRAS), which uses the transcriptome generated by single-molecule real-time sequencing as a reference sequence to evaluate the population variation of gene sequence and gene expression. Candidate genes were identified when both scores were associated with a trait, and their potential interactions were determined by expression quantitative trait locus analysis. In the following example, by applying this approach to characterize garlic bulb traits in 102 landraces, we identified 23 candidate transcripts, most of which showed extensive interactions. Thirteen transcripts were lncRNAs, and the others were proteins mainly involved in carbohydrate metabolism and protein degradation. TRAS serves as an efficient tool for association studies independent of reference genomes, extending the applicability of association studies to a wide range of species.

Example

1. Method

1. Plant Material, Experimental Design and Phenotypic Measurements

We collected 92 Chinese local varieties and 10 overseas varieties. A total of 102 garlic local varieties were planted in the experimental fields of the Institute of Hemp Research, Chinese Academy of Agricultural Sciences in September 2014. The field experiment adopts a random complete block design with two repetitions in each block. Each replicate planted 36 cloves of garlic of a variety in three rows with an interval of 10 cm between each row and an interval of 20 cm between two replicates. Harvest mature garlic in 2015, and take 30 plants in the middle of the three rows. Then, air dry the surfaces and measure their length, width, thickness with calipers.

2. PacBio single-molecule real-time sequencing library construction and sequencing

The present invention selects a local variety——Yangxi from 102 local varieties of garlic, and performs single-molecule sequencing. About 1 μg of total RNA from developing bulbs was used for reverse transcription to synthesize full-length cDNA using the Clontech SMARTer cDNA Synthesis Kit (Clontech Laboratories, Mountain View, CA, USA) and oligo dT primers. The reaction was repeated three times. PCR products were purified using AMPure PB beads (Pacific Biosciences, Menlo Park, CA, USA). The obtained dsDNA was size selected and subsequently reamplified using the BluePippin system. Finally, the product was purified and subjected to Iso-Seq SMRTBell library prep (https://pacbio.secure.force.com/SamplePrep) to generate three libraries (1-2kb, 2-3kb and 3-6kb), using P6- The C4 chemistry kit sequenced the three libraries on the PacBio RSII platform.

3. PacBio sequencing data analysis

Sequencing data were processed using SMRT analysis software (http://pacificbiosciences.github.com/DevNet/). Circular sequencing reads generated from subread files were obtained using the following program: minimum length 300; max_drop_fraction, 0.8; min_passes, 1; min_predicted_accuracy, 0.8. According to whether there is a polyA tail and the minimum sequencing length of 300, the generated circular sequencing reads are divided into full-length and non-full-length. By default, both full-length and non-full-length circularized sequences belong to the same type of cluster. Transcripts are defined as those sequences that are non-redundant and have no ends. Because the PacBio reads have higher nucleotide errors compared with the next-generation sequencing of shorter sequencing reads, they were corrected using the proovread software based on the Illumina RNA sequencing data of Yangxi garlic. Redundant sequences are put into CD-HIT.

4. Annotation of Transcriptome

We used the Coding Potential Calculator to detect the coding potential of transcripts based on transcriptome quality, completeness, and open reading frame alignment with similar sequences in current databases. Those without protein-coding potential were categorized as long noncoding RNAs (lncRNAs), and the rest were annotated from seven public databases as in previous studies. The coding sequence (CDS) of each transcript was searched using BLAST in the NCBI non-redundant protein sequence database and the SwissProt protein database to predict function, and the ESTscan program was performed.

5. Illumina RNA sequencing library construction and sequencing

To characterize population SNPs (single nucleotide) and GE (gene expression) variation, all 102 cultivars were subjected to Illumina RNA sequencing. Use kit【 UltraTM RNALibrary Prep Kit for (New England BioLabs, Ipswich, MA, USA)] The total RNA of each individual was constructed into a mixed fragment library consisting of cDNA with a fragment length of 250bp (±25bp). Then, the Illumina sequencing platform (HiSeqTM2500) combined with HiSeq PE Cluster Kit v4 cBot was used for sequencing. Finally, filter the raw data (rawreads) to obtain the analyzable data (clean reads).

6. Expression Analysis

In order to quantify the expression of transcripts of 102 species, all clean Illuminareads of each species were aligned to the transcripts obtained by PacBio sequencing by Bowtie 2. Using RSEM, the expression level of each transcript was analyzed for each sample by evaluating the measured expected value of fragments per kilobase per million bases of transcript sequence. Coexpression analysis of transcripts in 102 genotypes was performed using Weighted Gene Coexpression Network Analysis (WGCNA) in the R language package, and transcripts showing coexpression in the population were assigned to a module. The minimum module size was set to 30 transcripts, and modules were merged if they shared no less than 25% similarity. Then, the eigenvalues of each module were estimated, and the co-expression modules associated with the trait were identified based on the Pearson correlation factor analysis of the module's correlation with the trait.

7. Population SNP detection and phylogenetic analysis

To identify SNPs from 102 landraces, all Illumina reads for each sample were aligned to the reference transcriptome by using Samtools in Picard. After removing reads without unique positions, the SNP calling of the population was performed using GATK software; this process requires SNP quality ≥ 40. To rule out SNP calling errors caused by incorrect alignments, only high-quality SNPs (coverage depth ≥ 2, minor allele frequency ≥ 0.05, and deletion genotype frequency ≤ 0.5) were retained. An individual-based neighbor-joining tree was constructed using TreeBest (http://treesoft.sourceforge.net/treebest.shtml) to clarify the phylogenetic relationships among 102 local populations, and in Figtree (http://tree.bio .ed.ac.uk/software/figtree/).

8. Identify Trait Candidate Transcripts

To mine candidate transcripts for CL (flap length), CW (flap width) and CT (flap thickness), we performed a transcriptome-wide analysis to detect possible loci associated with the trait. In a combined panel containing 102 samples, a mixed linear model was employed to perform association analysis for three CSTs (flap-shaped traits) using a total of 19,912 high-quality SNPs. A suggested (1/N) P-value threshold was set to control the error rate, resulting in a P-value threshold of 5×10 ⁻⁵ . To validate associations between proposed loci and traits, Pearson correlations between expression of transcripts associated with proposed loci and trait phenotypes were calculated; significant correlations were assumed at P<0.05. Transcripts associated with traits at both sequence and expression levels were defined as candidate transcripts for traits.

9. eQTL analysis to detect potential interactions

eQTLs, which link changes in gene expression levels to genotype, have been shown to detect gene interactions. To determine the interrelationships of candidate transcripts, eQTL analysis of 102 local cultivars from TRAS (Association Analysis of Transcriptome Reference) was performed using mixed linear models. In our eQTL analysis, 19,912 SNPs were defined as genotype and expression of candidate transcripts as phenotype. A significant (0.05/N) p-value threshold was set at 2.5 × 10 ^-6 to control the error rate. If one transcript maps to another transcript's eQTL, we think that these two transcripts may interact.

2. Results

1. Three-generation sequencing (Pacific Biosciences) to obtain the transcriptome data of the developing flap

The local variety Yangxi had a total of 36,321 transcripts totaling 54.48 million bases. Transcript length ranges from 120bp to 4803bp, with an average of 1500bp. Detecting PloyA and PloyT, 70% were found to be full-length transcripts. 31,125 transcripts had functional annotations, 287 were not annotated, and 4,909 were lncRNAs (long non-coding RNAs).

2. Population variation in traits, sequences, and gene expression

The variation of petal shape characteristics of 102 varieties was 2.3 times of petal length, 3.7 times of petal width and 3.6 times of petal thickness. The three traits have a significant correlation P<0.01, indicating that they influence each other.

To perform sequence and expression analysis of the genotypes of 102 landraces, we sequenced the transcriptome of developing bulbs in the population and obtained approximately 6.05 billion clean reads, with an average of 59.3 million reads per landrace. These sequences were aligned to the reference transcriptome (three-generation sequencing) with a coverage of 51.3%–80.2%. These aligned reads were further used to score variation in gene sequence (SNP) and gene expression (GE). 55,012 SNPs were identified from 8,245 transcripts, indicating that 77% of the transcripts were conserved among the 102 varieties. After filtering, 19,912 high-quality SNPs were obtained from 5,408 transcripts, half of which were located within the coding sequence region, and the remaining SNPs were located in the 3' and 5' non-coding regions. Phylogenetic analysis based on SNP genotypes divided 102 garlic landraces into three distinct groups. Interestingly, the landraces from China showed a clear geographical distribution. We also quantified GE and characterized the genome-wide expression profile of each landrace during the flap expansion stage. Based on the co-expression network analysis of weighted genes, 36,321 transcripts were involved in 46 co-expression models (CEMs), each containing genes ranging from 37-7,132.

3. Candidate transcripts related to flap length, flap width and flap thickness

We identified 20 SNPs from 18 transcripts associated with flap thickness, 40 SNPs from 24 transcripts associated with flap width, and 21 SNPs from 15 transcripts associated with flap length (P<5×10 ^–5 ). A total of 42 transcripts were associated with flap-type signatures (CSTs). 5 of them are associated with all three traits, and 5 are associated with two traits. Because linkage disequilibrium-based association analyzes can only identify trait-associated genomic regions that may contain many genes, it remains difficult to identify candidate genes for traits. To determine whether candidate transcripts are involved in CSTs, we performed correlation analysis on GE and corresponding trait phenotypes of 42 transcripts, and the results showed that the expression of 5, 17 and 5 transcripts were significantly correlated with CT, CW and CL ( P<0.05). Among them, 2, 16 and 4 transcripts were significantly negatively correlated with CT, CW and CL, respectively, suggesting that these transcripts negatively regulate garlic clove growth. Only 3, 1 and 1 transcripts were found to positively regulate CT, CW and CL, respectively. Because both gene sequences and GE levels were associated with traits, we concluded that the 5, 17 and 5 transcripts identified above were candidate transcripts for CT, CW and CL, respectively. Since four of these were pleiotropic, the total number of candidate transcripts was reduced to 23.

4. Potential interactions between flap shape-related candidate transcripts

To detect potential interactions among trait-associated transcripts, we performed eQTL analysis on 23 trait-associated transcripts. The results showed that some SNPs were associated with the expression of 23 transcripts. No SNPs were associated with the expression of ASTG3334, ASTG34822, ASTG35427, and ASTG9068 transcripts, suggesting that the expression of these transcripts is rarely regulated by other transcripts and that they may function upstream in pathways involved in trait regulation. The expression of ASTG35908, ASTG1416, ASTG4276 and ASTG4283 was associated with at least 40 SNPs (in 29, 26, 26 and 207 transcripts, respectively), suggesting that these four candidate transcripts are regulated by the expression of many genes and that they may play a role in trait regulation downstream of the pathway.

We define the transcript where the eQTL resides as the eQTL-transcript, and we hypothesize that if a transcript and its eQTL-transcript are associated with the same trait, then both transcripts control the trait through interaction. Next, we identified eQTL transcripts for each candidate trait-associated transcript. The results showed that among 5 CT and 5 CL candidate transcripts, only two CT and two CL transcripts showed potential interactions, while 14 of 17 CW transcripts formed a potential interaction network. CW candidate transcripts were divided into three groups. Six transcripts whose expression was not regulated by other transcripts were assigned to the first group and considered as upstream transcripts for trait regulation. Six transcripts that not only regulate the expression of downstream transcripts but also upstream transcripts were classified in the second group. The second set of six transcripts showed extensive interactions with each other, resulting in a subnetwork of interactions. The remaining two transcripts not involved in the regulation of expression of other transcripts were assigned to the third group. Thus, the results provide a fundamental insight into the interaction networks governing CSTs.

It should be noted that at last: above each embodiment is only in order to illustrate technical scheme of the present invention, and is not intended to limit; Although the present invention has been described in detail with reference to foregoing each embodiment, those of ordinary skill in the art should understand that: It is still possible to modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements for some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the technical solutions of the various embodiments of the present invention. range.

Claims (10)

1. A method for association analysis of transcriptome reference, characterized in that, comprising: 1). Establish a reference population and measure the transcriptome sequences of all individuals in the reference population to obtain reference transcripts; 2). Perform high-throughput RNA sequencing on the predicted population, and compare and splice the obtained sequencing data with the reference transcript to obtain the population transcriptome sequencing result; 3). Analyze the population SNP data from the population transcriptome sequencing results, use the population SNP and traits of interest to perform association analysis, and determine the region where the target gene is located; then use the gene expression variation and traits in the target region to perform Correlation to get the candidate genes of the target traits. 2. the association analysis method of transcriptome reference according to claim 1, is characterized in that, described method also comprises: Determine the eQTL of the candidate gene, and then find whether the remaining target trait candidate genes are contained in the eQTL. 3. The method for association analysis of transcriptome reference according to claim 1, wherein the gene expression variation includes gene sequence variation and gene expression variation. 4. The method for association analysis of transcriptome reference according to claim 1, characterized in that, in step 1), after obtaining the reference transcript, it further comprises: performing functional annotation on the reference transcript. 5. The method for association analysis of transcriptome reference according to claim 1 or 4, characterized in that, in step 1), the number of species in the reference population is ≥ 1. 6. The method for association analysis of transcriptome reference according to claim 1, characterized in that, in step 1), the method for measuring the transcriptome sequence is single-molecule real-time sequencing. 7. The method for association analysis of transcriptome reference according to claim 1, characterized in that, in step 3), before utilizing the population SNP to perform association analysis with the trait of interest, the population transcriptome is sequenced Results Co-expression analysis was performed to assign transcripts showing co-expression in the population to a module, and the association of said modules with a trait of interest was analyzed to identify co-expressed modules associated with said trait of interest. 8. The association analysis method of transcriptome reference according to claim 1, characterized in that, in step 3), when utilizing the population SNP and the trait of interest to carry out association analysis, the association analysis model used is a general linear model or a mixed linear model. 9. The method for association analysis of transcriptome reference according to claim 1, characterized in that, in step 3), when using the population SNP to perform association analysis with the trait of interest, the quality score of the SNP ≥ 40 , coverage depth ≥ 2, minor allele frequency ≥ 0.05, deletion genotype frequency ≤ 0.5. 10. The application of the method according to any one of claims 1 to 9 in the screening and identification of major genes and the correlation of gene functions.