Brassica juncea DNA-BSA (deoxyribonucleic acid-bovine serum albumin) sequencing analysis method

Abstract

The invention discloses a DNA-BSA (deoxyribonucleic acid-bovine serum albumin) sequencing analysis method for brassica juncea, and relates to the technical field of biological information, and the DNA-BSA sequencing analysis method for brassica juncea is technically characterized in that after DNA information of brassica juncea is screened and filtered, candidate genes associated with target traits can be directly screened, and subsequent research is facilitated.

Description

Technical Field

[0001] The present invention relates to the field of bioinformation technology, and in particular to a mustard rapeseed DNA-BSA sequencing analysis method. Background Art

[0002] DNA sequencing of Brassica juncea is the process of sequencing the genome of the plant, aiming to analyze its genetic information, gene structure, and function, providing fundamental data for breeding and improvement. The process includes sample collection and DNA extraction, DNA library construction, sequencing, and data analysis. Commonly used high-throughput sequencing technologies, such as Illumina and PacBio, can provide detailed genomic sequence data.

[0003] However, existing technologies also have some shortcomings. They take a long time to construct experimental populations, require large sample sizes, are expensive to develop, are insensitive to major genes, and have a large mapping interval that makes it difficult to quickly extract candidate genes associated with the target trait.

[0004] To this end, the present invention aims to provide a DNA-BSA sequencing analysis method for Brassica juncea to solve the above problems. Summary of the Invention

[0005] The purpose of the present invention is to solve the above problems and provide a DNA-BSA sequencing analysis method for Brassica juncea, which can screen out candidate genes associated with target traits.

[0006] In order to achieve the above object, the technical solution of the present invention is as follows:

[0007] The present invention provides a DNA-BSA sequencing analysis method for Brassica juncea, comprising the following steps:

[0008] S1. Extract the genetic information of the samples and resequence them; complete the resequencing of 4 samples, with the data volume meeting the contract standard and ensuring that Q30 reaches 80%;

[0009] S2. Perform data analysis and evaluation on the resequencing data and obtain CleanReads after data quality control;

[0010] S3, align CleanReads with the genome;

[0011] S4, variant detection and annotation;

[0012] S5. Association analysis: Determine the region associated with the target trait by calculating the genotype frequency of alleles between two mixed pools;

[0013] S6. Candidate SNP and InDel annotation: Annotate SNPs and InDels within the association region;

[0014] S7. Candidate gene annotation: Genes in the associated region were annotated using the GO, KEGG, COG, NR, and SwissProt databases.

[0015] In step S2, the analysis and data evaluation include statistics of sequencing data quantity, sequencing data quality and GC content.

[0016] In step S3, the comparison content includes comparison efficiency, genome sequencing depth, and genome coverage statistics.

[0017] In step S4, the variation detection is the detection of SNPs and InDels.

[0018] In step S6, the annotation content includes position information, non-synonymous mutation information and frameshift mutation information.

[0019] Compared with the existing technology, this solution has the following beneficial effects:

[0020] The present invention provides a DNA-BSA sequencing and analysis method for Brassica juncea. This method does not require the construction of a complex population. By screening, filtering, and analyzing the DNA information of Brassica juncea, candidate genes with a strong correlation with target traits can be quickly discovered. The method also has low R&D costs and facilitates subsequent research. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 It is an experimental flow chart in an embodiment of the present invention;

[0022] Figure 2 is a flow chart of bioinformatics analysis of resequencing BSA in an embodiment of the present invention;

[0023] Figure 3 is a distribution diagram of sample base error rates in an embodiment of the present invention;

[0024] Figure 4 This is a distribution diagram of the ratio of each base in the sample according to the embodiment of the present invention;

[0025] Figure 5 This is a distribution diagram of chromosome coverage depth of samples in an embodiment of the present invention;

[0026] Figure 6 This is a distribution diagram of inserted fragments in an embodiment of the present invention;

[0027] Figure 7 is a depth distribution diagram of a sample in an embodiment of the present invention;

[0028] Figure 8 is a SNP quality distribution diagram in an embodiment of the present invention;

[0029] Figure 9This is a distribution diagram of SNP mutation types in an embodiment of the present invention;

[0030] Figure 10 is a Venn diagram of SNP statistics between samples in an embodiment of the present invention;

[0031] Figure 11 It is a puzzle of the SNP annotation results in the embodiment of the present invention;

[0032] Figure 12 is a distribution diagram of InDel lengths in the whole genome and coding regions in an embodiment of the present invention;

[0033] Figure 13 INDEL statistics Venn diagram of samples in the embodiment of the present invention;

[0034] Figure 14 is the distribution of ED association values on chromosomes in the embodiment of the present invention Figure 1 ;

[0035] Figure 15 is a distribution diagram of SNP-index association values on chromosomes in an embodiment of the present invention;

[0036] Figure 16 is the distribution of ED association values on chromosomes in the embodiment of the present invention Figure 2 ;

[0037] Figure 17 is a distribution diagram of the InDel-index association value on the chromosome in an embodiment of the present invention;

[0038] Figure 18 This is a GO annotation cluster diagram of genes in the candidate region in an embodiment of the present invention;

[0039] Figure 19 is a pathway distribution map of genes in the candidate region in an embodiment of the present invention;

[0040] Figure 20 This is the gene COG annotation classification diagram of the SNPs in the candidate region in the embodiment of the present invention

[0041] Figure 21 : is a distribution diagram of the results of the samples in the embodiment of the present invention visualized on the chromosome, where A is the SNP density distribution and B is the InDel density distribution. DETAILED DESCRIPTION

[0042] In order to enable those skilled in the art to better understand the present invention, the technical solution of the present invention will be further described in detail below in conjunction with the embodiments of the present invention and the accompanying drawings. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

[0043] It should be noted that, in the absence of conflict, the embodiments and features of the embodiments of the present invention can be combined with each other. The present invention will be described in detail below with reference to the embodiments.

[0044] Example:

[0045] 1. Sample information:

[0046] Table 1 Sample information

[0047] Analyzed sample name fq1 fq2 F01 Unknown_BM035-02R0001_good_1.fq.gz Unknown_BM035-02R0001_good_2.fq.gz H03 Unknown_BM035-01R0001_good_1.fq.gz Unknown_BM035-01R0001_good_2.fq.gz L04 Unknown_BM035-01R0002_good_1.fq.gz Unknown_BM035-01R0002_good_2.fq.gz M02 Unknown_BM035-02R0002_good_1.fq.gz Unknown_BM035-02R0002_good_2.fq.gz

[0048] Table 2 Genome information

[0049]

[0050] 2. Experimental Procedure

[0051] 2.1 Experimental Procedure

[0052] The experimental process was carried out according to the standard protocol provided by Illumina, including sample quality testing, library construction, library quality testing and library sequencing. The specific process is as follows Figure 1 shown.

[0053] After the sample genomic DNA is tested and qualified, the DNA is fragmented by mechanical shearing (ultrasound), and then the fragmented DNA is purified, end-repaired, 3'-end A added, and sequencing adapters are connected. Then, agarose gel electrophoresis is used to select the fragment size, and PCR amplification is performed to form a sequencing library. The constructed library is first subjected to library quality inspection, and the library that passes the quality inspection is sequenced using Illumina.

[0054] 2.2 Information Analysis Process

[0055] The content of information analysis includes: data quality control (removal of adapters and low-quality data), comparison with the reference genome, variant detection and annotation (SNP, INDEL), association analysis, and annotation of candidate SNPs and candidate genes. The specific process of resequencing BSA bioinformatics analysis is as follows: Figure 2 shown.

[0056] 3. Bioinformatics Analysis and Results

[0057] 3.1 Introduction to Sequencing Data

[0058] The raw image data files obtained by high-throughput sequencing are converted into raw sequencing sequences (SequencedReads) through BaseCalling analysis, called RawData or RawReads. The results are stored in the FASTQ (abbreviated as fq) file format, which contains the sequence information of the sequencing sequence (Reads) and its corresponding sequencing quality information.

[0059] The detailed information of Illumina Sequence Identifiers is as follows:

[0060] Table 3 Illumina sequencing identification details

[0061]

[0062] By using the ASCII value corresponding to each character in the fourth row to calculate, the sequencing quality value corresponding to the base in the second row is obtained. If the sequencing error rate is represented by e and the base quality value of Illumina is represented by Qphred, the following relationship exists:

[0063] Q Phred = -10log10(e)

[0064] 3.1.2 Data Quality Statistics

[0065] The original sequencing sequences (SequencedReads) or RawReads obtained by sequencing contain low-quality reads with adapters. To ensure the quality of information analysis, RawReads are filtered to obtain CleanReads for subsequent information analysis. The main steps of data filtering are as follows:

[0066] (1) Remove reads with adapters;

[0067] (2) filtering reads with N content exceeding 10%;

[0068] (3) Reads with bases with a quality value lower than 10 accounting for more than 50% were removed.

[0069] Table 4 Sample sequencing data evaluation statistics

[0070]

[0071]

[0072] 3.1.3 Base Sequencing Quality Distribution

[0073] The sequencing error rate of each base is obtained by converting the sequencing Phred value (Phredscore, Qphred) through Formula 1. The Phred value is calculated during the base calling process using a model that predicts the probability of base call errors. The corresponding relationship is shown in the following table:

[0074] Table 5 Predicted probability of base call error

[0075]

[0076] When sequencing on the Illumina sequencing system, the library is first prepared for chip preparation, with the aim of fixing the library DNA template on the chip. During the process of fixing the DNA template, each DNA molecule will form a cluster, and a cluster is a sequencing site. During the fixation process, a very small number of clusters will overlap physically. During sequencing, the sequencing software analyzes and identifies these overlapping points through the first 4 bases, separates the positions of these overlapping points, and ensures that each point measures a DNA molecule. Therefore, the error rate of the first few bases at the 5′ end of the sequencing sequence is relatively high. In addition, the sequencing error rate will increase with the increase in the length of the sequencing sequence (SequencedReads), which is caused by the consumption of chemical reagents during the sequencing process. Therefore, when performing base sequencing quality distribution analysis, the quality values of the sample's base quality distribution in the first 4 bases and the last dozen bases will be lower than the middle sequencing bases, but their quality values are all higher than Q30. According to the relationship between quality value and error rate, the quality value is converted into error rate, and the error rate distribution graph is drawn as shown below. Figure 3 shown.

[0077] 3.1.4 Base type distribution

[0078] The base type distribution check is used to detect the presence of AT and GC separation, which may be caused by sequencing or library construction and will affect subsequent analysis. The sequences measured by high-throughput are DNA fragments after random interruption of the genome. Since the distribution of sites on the genome is approximately uniform, the G / C and A / T contents are also approximately uniform. Therefore, according to the law of large numbers, in each sequencing cycle, the GC and AT contents should be equal, respectively, and equal to the GC and AT contents of the genome. Similarly, due to the relationship between overlapping clusters, the AT and GC of the first few bases of the sample will fluctuate greatly, which is higher than other sequencing segments, while the GC and AT contents of other segments are equal and evenly distributed without separation. Figure 4 shown.

[0079] 3.2 Comparison statistics with reference genome

[0080] CleanReads obtained from sequencing need to be remapped to a reference genome before subsequent variant analysis can be performed. bwa-mem2 (v2.2) software is primarily used to align short sequences obtained from second-generation high-throughput sequencing (such as those from Illumina sequencing platforms) with a reference genome. CleanReads are aligned to the reference genome using "bwa-mem2mem-t4-M." The alignment results are sorted using samtools (v1.9) sort. Based on the sorted results, statistics are generated for each sample, including sequencing depth and genome coverage.

[0081] 3.2.1 Comparison results statistics

[0082] Compare to the reference genome and use commands such as samtoolsflagstat / depth to calculate sample alignment efficiency, sequencing depth, and coverage. If the reference genome is selected appropriately and there is no contamination during the experimental process, the sequencing read rate will generally be higher than 70%. In addition, the alignment rate is affected by the close relationship between the sequenced species and the reference genome, the quality of the reference genome assembly, and the quality of the reads. The closer the species, the more complete the reference genome assembly, and the higher the quality of the sequencing reads, the more reads can be mapped to the reference genome, and the higher the alignment rate. The sample alignment results are shown in the following table:

[0083] Table 6 Comparison results statistics

[0084]

[0085] The coverage depth of each chromosome site is plotted. Due to the limited image space, a maximum of 20 chromosomes / scaffolds can be displayed. If the coverage depth is evenly distributed on the chromosome, it can be considered that the sequencing randomness is good. The chromosome coverage depth distribution of the sample is as follows: Figure 5 shown.

[0086] Depend on Figure 5 The genome is evenly covered, indicating good sequencing randomness. Uneven depth areas on the graph may be caused by repetitive sequences or PCR bias. If all samples have no read coverage at the same position, it may indicate a large gap in the genome assembly.

[0087] 3.2.2 Insert fragment distribution statistics

[0088] By detecting the start and end positions of paired-end sequences on the reference genome, we can determine the actual size of the sequenced fragments generated after fragmentation of the sample DNA, known as the insert size (InsertSize). This is a crucial parameter for information analysis. The distribution of insert sizes generally follows a normal distribution with a single peak. The InsertSize distribution plot can show the distribution of insert lengths for each sample. Insert size analysis of each sample's sequencing data was performed using the CollectInsertSizeMetric.jar software in the Picard software toolkit.

[0089] Depend on Figure 6 It can be seen that the distribution of insert fragment length conforms to the normal distribution, indicating that there is no abnormality in the construction of the sequencing data library.

[0090] 3.2.3 Depth distribution statistics

[0091] After the reads are located on the reference genome, the coverage of the bases on the reference genome can be counted. The percentage of bases covered by reads on the reference genome is called genome coverage; the number of reads covering the bases is the coverage depth. Genome coverage can reflect the completeness of variation detection on the reference genome. The more areas covered, the more variation sites can be detected. Coverage is mainly affected by the sequencing depth and the closeness of the relationship between the sample and the reference genome. The coverage depth of the genome will affect the accuracy of variation detection. In areas with higher coverage depth (non-repetitive sequence areas), the accuracy of variation detection is also higher. In addition, if the coverage depth of the bases on the genome is more evenly distributed, it also means that the sequencing randomness is better. The base coverage depth distribution curve and coverage distribution curve of the sample are shown in Figure 7 .

[0092] 3.3 Mutation Detection and Annotation

[0093] 3.3.1 Mutation Detection Tools and Methods

[0094] SNP (Single Nucleotide Polymorphism) and small INDEL (Small Insertion and Deletion) detection is primarily performed using the GATK software toolkit. Based on the CleanReads mapping results on the reference genome, samtools (v1.9) is used to filter redundant reads to ensure the accuracy of the detection results. SNP and INDEL variant detection is then performed using the HaplotypeCaller (local haplotype assembly) algorithm of GATK (v3.8). A gVCF is generated for each sample, followed by population joint-genotyping. Finally, a "hard" filtering process is performed to obtain the final set of variant sites, which are typically stored in VCF format.

[0095] 3.3.2 Quality Control of Variant Detection Results

[0096] The mutation results are strictly filtered to ensure the reliability of the mutation results. The main filtering parameters are as follows:

[0097] (1) SNPs within 5 bp of the InDel and adjacent InDels within 10 bp were filtered out based on the subroutine vcfutils.pl (varFilter-w5-W10) in bcftools;

[0098] (2) clusterSize2clusterWindowSize5, indicating that the number of variants within a 5 bp window should not exceed 2;

[0099] (3) QUAL < 30, the quality value in Phred format, indicates the possibility of variant variation at the site.

[0100] Those with a quality value lower than 30 were filtered out;

[0101] (4) QD < 2.0, the ratio of the variant quality value (Quality) divided by the coverage depth (Depth), the coverage depth is the sum of the coverage depths of all samples containing variant bases at this site. Samples with a QD lower than 2.0 are filtered out;

[0102] (5) MQ<40, the root mean square of the alignment quality of all reads aligned to this site. Reads with an MQ lower than 40 were filtered out;

[0103] (6) FS>60.0, a value converted from the p-value of the Fisher test, describes whether there is obvious positive and negative strand specificity for reads containing only variants and reads containing only reference sequence bases during sequencing or alignment. In other words, if there is no strand-specific alignment result, the FS should be close to zero. FS above 60 is filtered out;

[0104] (7) Other variant filtering parameters are processed using the default values officially specified by GATK.

[0105] 3.3.3 High-quality mutation result display

[0106] SNP variations are categorized into transitions and transversions. Mutations between bases of the same type are called transitions, such as purine-to-purine and pyrimidine-to-pyrimidine variations. Mutations between bases of different types are called transversions, such as purine-to-pyrimidine variations. Generally speaking, transitions are more likely to occur than transversions, so the transition / transversion (Ti / Tv) ratio is typically greater than 1. The specific value depends on the species being tested. For diploid or polyploid species, if all SNPs on homologous chromosomes contain the same base, the SNP is considered homozygous. If the SNPs on homologous chromosomes contain different bases, the SNP is considered heterozygous. A greater number of homozygous SNPs indicates a greater divergence between the sample and the reference genome, while a greater number of heterozygous SNPs indicates a higher degree of heterozygosity. The specific results depend on the sample material used. The SNP detection results for the sample and the reference genome are shown in the table below.

[0107] Table 7 Statistics of SNPs detected

[0108] BMKID SNP number Transition Transversion Ti / Tv Heterozygosity Homozygosity Het-ratio F01 1,411,018 819,440 591,578 1.38 620,222 790,796 43.95％ H03 1,993,879 1,157,875 836,004 1.38 1,594,465 399,414 79.96％ L04 1,976,107 1,147,213 828,894 1.38 1,571,858 404,249 79.54％ M02 1,274,443 739,664 534,779 1.38 288,427 986,016 22.63％ Total 2,203,490 1,280,994 922,496 1.38 - - -

[0109] 3.3.4 Detection of SNPs between samples (multiple samples)

[0110] Based on the comparison results between the sample and the reference genome, all the different variant sites between the samples are summarized. The SNP list file format between the samples is as follows:

[0111] Table 8 Schematic diagram of SNP list between samples

[0112]

[0113]

[0114] The coding of SNP genotypes uses standard nucleotide symbols, and the symbol table is as follows:

[0115] Table 9 SNP genotype coding

[0116] Nucleotide code significance Nucleotide code significance A Adenosine M AC(aMinogroup) C Cytosine S GC (Strong interaction) G Guanine W AT (Weak interaction) T Thymidine B GTC(notA)(BcomesafterA) U Uracil D GAT(notC)(DcomesafterC) R GA(puRine) H ACT(notG)(HcomesafterG) Y TC(pYrimidine) V GCA(notT, notU)(VcomesafterU) K GT(Ketone) N AGCT(aNy)

[0117] Whole genome SNP mutations can be divided into 6 categories. Taking T:A>C:G as an example, this type of SNP mutation includes T>C and A>G. Since sequencing data can be aligned to both the positive and negative strands of the reference genome, when T>C type mutations appear on the positive strand of the reference genome, A>G type mutations are at the same position on the negative strand of the reference genome, so T>C and A>G are divided into one category. The statistical results of SNPs between samples are as follows: Figure 10 shown.

[0118] The detailed SNP annotation classification of each sample (All represents the whole) is shown in the following table:

[0119] Table 10 SNP annotation results statistics

[0120]

[0121]

[0122] Based on the whole genome SNP annotation classification results of all samples, a pie chart is drawn as follows Figure 11 shown.

[0123] The SNP annotation types are described in the following table

[0124] Table 11 Functional area meaning

[0125]

[0126] 3.3.6 Detection of SmallInDels between Samples and Reference Genomes

[0127] Based on the mapping results of the sample's CleanReads on the reference genome, we detect whether there are small insertions and deletions (SmallInDels) between the sample and the reference genome. Sample insertions and deletions are detected using GATK. SmallInDel variations are generally less common than SNP variations and also reflect differences between the sample and the reference genome. InDels in coding regions can cause frameshift mutations, leading to changes in gene function. Specific data are shown in the table below.

[0128] Table 12 InDel statistics for the whole genome and coding regions

[0129]

[0130]

[0131] Statistics were made based on the length of InDels in the CDS region and the whole genome. The length distribution diagram is shown in Figure 12 (If there are many samples, 50 will be displayed by default).

[0132] 3.3.7 Detection of SmallInDels Between Samples (Multiple Samples)

[0133] Based on the test results of the sample and the reference genome, the comparison results of the sample sequencing data are shown in the table below. Figure 13 shown.

[0134] Table 13 Sample SmallInDel sequencing data statistics

[0135]

[0136] 3.3.8 Notes on SmallInDel

[0137] Based on the location of the SmallInDel site in the reference genome obtained from sample testing and comparison with the reference genome's gene and CDS location information (typically in the gff file), we can annotate whether the InDel site occurs in an intergenic region, a gene region, or a CDS region, and whether it is a frameshift mutation. SmallInDel annotation is performed using SnpEff software. InDels with frameshift mutations may alter gene function. Specific annotation results for each sample (All represents the overall data) are shown in the table below.

[0138] Table 14 InDel annotation results statistics

[0139]

[0140]

[0141] Based on the whole genome InDel annotation classification results of all samples, the InDel annotation types are described in the following table:

[0142] Table 15 Note type description is shown in the following table

[0143]

[0144]

[0145] 3.4 Association analysis (SNP)

[0146] 3.4.1 High-quality SNP screening

[0147] Before association analysis, SNPs were filtered according to the following criteria:

[0148] 1. Filter out SNP sites with multiple genotypes and retain only diallelic genotype sites; 2. Filter out SNP sites with mixed pool read support less than 4;

[0149] 3. Filter out SNP sites with homozygous and consistent genotypes between the mixed pools; 4. Filter out SNP sites with homozygous and consistent genotypes between the two parents;

[0150] 5. Filter out SNP sites whose recessive mixed pool genotypes do not come from the recessive parent, and filter out SNP sites whose dominant mixed pool genotypes do not come from the dominant parent; ultimately, 937,395 high-quality reliable SNP sites were obtained.

[0151] Table 16 SNP filtering statistics

[0152]

[0153] 3.4.2 ED method correlation results

[0154] The Euclidean distance (ED) algorithm uses sequencing data to identify markers with significant differences between pools and evaluate regions associated with traits. Theoretically, except for the target trait-associated loci, all other loci between two pools constructed using BSA tend to be consistent, so the ED values for non-target loci should approach 0. The ED calculation formula is shown below. A larger ED value indicates a greater difference in the marker between the two pools.

[0155]

[0156] A mut is the frequency of base A in the mutation pool, A wt is the frequency of base A in the wild-type pool; C mut is the frequency of C base in the mutation pool, C wt is the frequency of C base in the wild-type pool; G mut is the frequency of G base in the mutation pool, G wt is the frequency of G base in the wild-type pool; T mut is the frequency of T base in the mutation pool, T wt is the frequency of T base in the wild-type pool.

[0157] During the analysis, the SNP sites with different genotypes between the two pools were used to count the depth of each base in different pools and calculate the ED value of each site. In order to eliminate background noise, the original ED value was squared. The fifth power of the original ED was used as the correlation value to eliminate background noise. Then, the SNPNUM method was used to fit the ED value. The correlation value distribution is shown as follows: Figure 14 shown.

[0158] The median + 3SD of the fitted values of all sites was taken as the association threshold for the analysis, which was calculated to be 0.01. Based on the association threshold, a total of 41 regions with a total length of 32.40 Mb were obtained.

[0159] Table 17 Statistics of associated area information

[0160]

[0161]

[0162] 3.4.3 SNP-index method association results

[0163] SNP-index is a recently published method for marker association analysis that uses differences in genotype frequencies between pools. It primarily identifies significant differences in genotype frequencies between pools and uses Δ(SNP-index) as a statistical measure. The stronger the association between a marker SNP and a trait, the closer Δ(SNP-index) is to 1.

[0164] SNPindex(aa)=Maa / (Maa+Paa),

[0165] SNPindex(ab)=Mab / (Mab+Pab),

[0166] ΔSNPindex=SNPindex(aa)-SNPindex(ab),

[0167] Maa represents the depth of the aa pool (recessive mixed pool) derived from the recessive parent; Paa represents the depth of the aa pool (recessive mixed pool) derived from the dominant parent. Mab represents the depth of the ab pool (dominant mixed pool) derived from the recessive parent; Pab represents the depth of the ab pool (dominant mixed pool) derived from the dominant parent.

[0168] In order to eliminate false positive sites, the ΔSNP-index values of markers on the same chromosome can be fitted using the position of the markers on the genome. The present invention uses the SNPNUM method to fit the ΔSNP-index, and then selects the region above the threshold as the region associated with the trait based on the association threshold. The distribution of SNP-index and ΔSNP-index of the two mixed pools is shown in the figure below. Figure 15 shown.

[0169] According to the results of computer simulation experiments, when the confidence level is 95, a total of 16 regions are obtained with a total length of 0.27Mb.

[0170] Table 18 Statistics of associated area information

[0171]

[0172] Note: Chromosome_ID: chromosome number; Start: starting position of the association region; End: ending position of the association region; Size: size of the association region, in Mb.

[0173] 3.4.4 SNP mapping results

[0174] The intersection of the SNP association region results obtained by the two association analysis methods is shown in the following table:

[0175] Table 19 Statistics of associated area information

[0176] Chromosome_ID Start End Size(Mb) A08 11,932,637 11,932,698 0 A08 11,941,577 11,941,577 0 A08 11,942,104 11,962,246 0.02 A08 11,962,273 11,962,828 0 A09 6,295,197 6,300,554 0.01 A09 6,301,145 6,301,331 0 A09 6,362,126 6,362,456 0 A09 6,366,933 6,367,173 0 A09 6,524,542 6,600,453 0.08 A09 6,659,772 6,660,022 0 A09 6,664,613 6,743,254 0.08 A09 6,781,700 6,785,293 0 A09 6,952,448 6,952,806 0 A09 6,953,384 7,029,604 0.08 A09 7,060,238 7,067,731 0.01 A09 7,085,142 7,088,419 0 Total - - 0.28

[0177] 3.5 Association Analysis (InDel)

[0178] 3.5.1 High-quality InDel screening

[0179] Before association analysis, InDels were first filtered using the same filtering criteria as for SNP analysis, ultimately obtaining 288,380 high-quality, reliable InDel sites.

[0180] Table 20 InDel filtration statistics

[0181]

[0182] 3.5.2 ED method correlation results

[0183] The same analysis method (ED method) as that used for SNP association analysis was used for the analysis, and the association value distribution was as follows: Figure 16 shown.

[0184] The median + 3SD of the fitted values of all sites was taken as the association threshold for the analysis, which was calculated to be 0.01. Based on the association threshold, a total of 18 regions with a total length of 31.77 Mb were obtained.

[0185] Table 21 Statistics of associated area information

[0186]

[0187]

[0188] 3.5.3 InDel-index method association results

[0189] When performing the analysis, the same analysis method (SNP-index method) as that used for SNP association analysis was used, and the association value distribution was as follows Figure 17 shown.

[0190] According to the results of computer simulation experiments, when the confidence level is 90, a total of 24 regions are obtained with a total length of 6.41Mb.

[0191] Table 22 Statistics of associated area information

[0192]

[0193]

[0194] 3.5.4 InDel positioning results

[0195] The intersection of the InDel association region results obtained by the two association analysis methods is shown in the following table:

[0196] Table 23 Statistics of associated area information

[0197]

[0198] 3.6 Candidate region screening and functional annotation

[0199] 3.6.1 Candidate Region Screening

[0200] The intersection of the results of the association regions corresponding to SNPs and InDels is taken, and the obtained intersection is shown in the following table:

[0201] Table 24 Statistics of associated area information

[0202]

[0203]

[0204] 3.6.2 SNP and InDel Annotation in Candidate Regions

[0205] The SNP annotation results within the candidate region among the samples of the present invention are shown in the following table:

[0206] Table 25 Statistics of SNP annotation results within candidate regions

[0207]

[0208] According to statistics, there are 96 SNPs with non-synonymous mutations between parents and 5 SNPs with non-synonymous mutations between pools. These SNPs are likely to be directly related to the traits, and the genes where they are located are called non-synonymous mutation genes.

[0209] The annotation results of the samples within the candidate region are shown in the following table:

[0210] Table 26 Statistics of InDel annotation results within candidate regions

[0211]

[0212]

[0213] According to statistics, there are 6 InDels with frameshift mutations between parents, and 0 InDels with frameshift mutations between pools. These InDels are likely to be directly related to the traits, and the genes where they are located are frameshift mutation genes.

[0214] 3.6.3 Gene Annotation in the Candidate Region

[0215] BLAST software was used to deeply annotate the coding genes within the candidate region using multiple databases (NR, SwissProt, GO, KEGG, COG). Through detailed annotation, candidate genes were quickly screened. A total of 39 genes were annotated within the candidate region, including 18 non-synonymous mutation genes and 5 frameshift mutation genes between the parents. The annotation results are shown in the table below:

[0216] Table 27 Statistics of gene function annotation results in SNPs and InDels in candidate regions

[0217]

[0218] 3.6.4 GO enrichment analysis of genes in the candidate region

[0219] The GO database is a structured, standardized biological annotation system that establishes a standardized vocabulary for the functions of genes and their products, applicable across species. The database is structured into multiple hierarchies, with nodes representing more specific functions at lower levels. GO analysis is used to categorize genes according to cellular component, molecular function, and biological process.

[0220] The corresponding gene GO classification statistics in the candidate regions are shown in Figure 18 .

[0221] The enrichment analysis results of the corresponding genes in the candidate regions are shown in the table below.

[0222] Table 28 Schematic diagram of topGO enrichment results of genes corresponding to SNPs in candidate regions

[0223]

[0224] 3.6.5 KEGG enrichment analysis of genes in the candidate region

[0225] In organisms, different genes coordinate with each other to perform biological functions. The same pathway between different genes is called a pathway. Pathway analysis helps further understand gene function. KEGG is the main public database for pathways.

[0226] The KEGG annotation results of the corresponding genes in the candidate regions are classified according to the pathway type, as shown in the classification diagram. Figure 19 shown.

[0227] The corresponding KEGG enrichment analysis results in the candidate regions are shown in the table below:

[0228] Table 29 KEGG enrichment results of genes corresponding to SNPs in candidate regions

[0229]

[0230] 3.6.6 COG classification statistics of corresponding genes in candidate regions

[0231] The COG database is constructed based on the phylogenetic relationship between bacteria, algae, and eukaryotes. The COG database can be used to classify gene products into orthologous groups. The COG classification statistics of genes in the associated region are shown in Figure 20 .

[0232] 3.7 Results Visualization

[0233] Combining the above content, the variation results of the samples and the BSA association analysis results were plotted using circos software. The distribution of the results between samples on the chromosome is visualized in Figure 21 .

[0234] The above specific embodiments are merely explanations of the present invention and are not limitations of the present invention. After reading this specification, those skilled in the art may make non-creative modifications to the embodiments as needed. However, as long as they are within the scope of the claims of the present invention, they are protected by patent law.

Claims

1. A method for DNA-BSA sequencing analysis of Brassica juncea, characterized by: The method comprises the following steps: S1. Extract the gene information of the sample and resequence the gene information of the sample; S2. Perform data analysis and evaluation on the resequencing data and obtain CleanReads after data quality control; S3, align CleanReads with the genome; S4, variant detection and annotation; S5. Association analysis: Determine the region associated with the target trait by calculating the genotype frequency of alleles between two mixed pools; S6. Candidate SNP and InDel annotation: Annotate SNPs and InDels within the association region; S7. Candidate gene annotation: Genes in the associated region were annotated using the GO, KEGG, COG, NR, and SwissProt databases.

2. The method for DNA-BSA sequencing analysis of Brassica juncea according to claim 1, wherein: In step S2, the analysis and data evaluation include statistics of sequencing data quantity, sequencing data quality and GC content.

3. The method for DNA-BSA sequencing analysis of Brassica juncea according to claim 1, wherein: In step S3, the comparison content includes comparison efficiency, genome sequencing depth, and genome coverage statistics.

4. The method for DNA-BSA sequencing analysis of Brassica juncea according to claim 1, wherein: In step S4, the variation detection is the detection of SNPs and InDels.

5. The method for DNA-BSA sequencing analysis of Brassica juncea according to claim 1, wherein: In step S6, the annotation content includes position information, non-synonymous mutation information and frameshift mutation information.

Images