CN117095747A - Method for detecting group inversion or transposon endpoint genotype based on linear ubiquitin genome and artificial intelligence model - Google Patents

Abstract

The invention discloses a method for detecting group inversion or transposon endpoint genotypes based on a linear ubiquitin genome and an artificial intelligent model, which comprises the steps of mounting group high-depth sequencing data on the linear ubiquitin genome; detecting the condition that each window is completely covered by the sequencing sequence according to the coverage condition of the second generation sequencing sequence of each sample in the population, and recording window position information and the number of reads completely covering the window; constructing an artificial intelligent model, and training the model through simulation data to judge whether a continuous window coverage area contains inversion or transposon endpoints; scanning all areas on one chromosome by using the model to obtain all inversion or transposon endpoint information on the chromosome, scanning a plurality of chromosomes in sequence, and summarizing the inversion or transposon endpoint information on all the chromosomes to form one sample; population-level inversion or transposon endpoint genotype matrices are sorted and screened based on the different samples. The invention realizes the detection of transposon and inversion end point genotypes by using the second generation sequencing data.

Description

A linear pan-genome and artificial intelligence model to detect population inversions or transpositions subendpoint genotype method

Technical field

The present invention relates to the technical field of natural population genotype detection, and more specifically to a method for detecting population inversion or transposon endpoint genotypes based on linear pan-genome and artificial intelligence models.

Background technique

In the field of population genetic research, the acquisition of population genotypes is the cornerstone of population research. Traditional population genotype analysis includes single nucleotide polymorphism (SNP) analysis and insertion and deletion (InDel) analysis. Studies have shown that most phenotypic variation is caused by structural variation. In recent years, structural variation (SV) analysis Some progress has also been made. Structural variations include inversions, translocations, duplications, and large insertions or deletions. At present, the price of third-generation sequencing is still relatively high, which greatly limits the mining of structural variations. The price of second-generation sequencing continues to decrease, but the read length of second-generation sequencing is short, which brings great difficulties to the analysis of structural variations, especially translocations and duplications caused by inversions and transposon jumping.

When a certain sequence of the genotype is inverted, the sequence is still included in the genome. Therefore, it is difficult to detect the existence of the inversion by using second-generation sequencing technology to mount short sequences to the reference genome. Currently, the construction of inversion maps mainly uses multiple assembled high-quality genomes to compare the genomes to find the inversion position. However, obtaining high-quality genomes requires huge costs, so the detected inversion maps are limited to Certain popular research species and limited individuals. Transposons can actively jump within the genome and may contain multiple repeats. Therefore, the short next-generation sequencing sequences within the transposon may be mounted to multiple locations at the same time, but a certain sample may only contain one repeat.

At present, there are massive public sequencing data for many species. However, limited by the shortcomings of second-generation sequencing read length, research on population transposons and inversions has been slow. Therefore, how to provide a method for detecting population transposons and inversions, and then analyze the genotype for subsequent research, is an urgent problem that those skilled in the art need to solve.

Contents of the invention

In view of this, the present invention provides a method for detecting population inversions or transposon endpoint genotypes based on linear pan-genome and artificial intelligence models to solve the problems in the above background technology and realize the use of second-generation sequencing data to detect transformations. Transposon and inversion endpoint genotypes.

In order to achieve the above objects, the present invention adopts the following technical solutions:

The present invention provides a method for detecting population inversions or transposon endpoint genotypes based on linear pan-genome and artificial intelligence models, which includes the following steps:

Step 1: Use comparison software to mount the population's high-depth sequencing data onto the linear pan-genome;

Step 2: Based on the second-generation sequencing sequence (reads) coverage of each sample in the population, use the windowing method to detect whether each window is completely covered by the sequencing sequence, record the window position information and the number of reads that completely cover the window, Generate data frame file;

Step 3: Construct an artificial intelligence model (artificial neural network model), and train the model through simulated data so that it can determine the area (continuous window coverage area) based on the fixed number, continuous window position and sequence number information obtained in step 2. whether it contains inversion or transposon endpoints;

Step 4: Use the model in Step 3 to scan all regions on a chromosome to obtain all inversion or transposon endpoint position information on the chromosome, scan multiple chromosomes in sequence, and summarize the inversions or transposons on all chromosomes of a sample. Socket endpoint information;

Step 5: Based on the inversion or transposon endpoint information of different samples, integrate the population-level inversion or transposon endpoint genotype matrix.

As the preferred technical solution for the above technical solution, in step one, the sequencing depth of the population high-depth sequencing data should be greater than 20×.

As the preferred technical solution for the above technical solution, in step 2, the windowing method is used to detect whether each window is completely covered by the sequencing sequence. The window refers to a 39bp region on the genome. The position of the window is named after the chromosome and the window. The middle position (the chromosome position is as if it were placed on a number axis, with the smaller on the left and larger on the right) combined to name, and the step size of the window is 20bp.

As the preferred technical solution for the above technical solution, in step 2, being completely covered means that the leftmost position of the sequencing sequence should be on the left side of the window (the left position of the sequencing number is smaller than the left position of the window), and the left position of the sequencing sequence plus the sequence matches the pan The sum of the number of bases on the genome and the number of detected missing bases should be greater than the right position of the window.

As the preferred technical solution for the above technical solution, in step 2, the window position information and the number of sequencing sequences that completely cover the window are recorded, and a data frame containing two columns of data can be obtained. The first column of the data frame is the window position, and the second column is the window position. Column is the number of sequences that completely cover the window, and is saved as a csv file; since all samples are aligned to the same linear pan-genome, the window is based on the linear pan-genome, and the size and step size of the window are consistent, so the data frame csv of all samples The number of file lines is the same, the information in the first column is exactly the same, and the sequence number information in the second column is different.

As the preferred technical solution for the above technical solution, the artificial neural network model in step three is a 4-layer fully connected layer network, the last layer activation function is the sigmoid function, and the other activation functions are ReLU functions; the input of the model is a column of length 25 Array, derived from the data frame csv file in step 2, perform windowing operation on the second column in the csv file, the window size is 25, the step size is 22, and the corresponding genome sequence length of each window is 39bp+(25-1) *20=519bp; the output of the model is a value used to determine whether the 519bp region contains an inversion or transposon endpoint.

As the preferred technical solution for the above technical solution, the data frame csv file in step 2 needs to be preprocessed before being input to the model. The number in the second column of the csv file that is less than or equal to 3 is defined as 0, indicating no sequence complete coverage. Numbers greater than 3 are defined as 1, indicating complete coverage by the sequence.

As the preferred technical solution for the above technical solution, in step three, the model is trained through simulated data. The simulated data is generated by artificial simulation and includes situations where there are obvious inversions or transposon endpoints (for example, 10 consecutive 1s, and the 11th bit is 0 , the 12th to 25th bits are all 1) and there is obviously no inversion or transposon endpoint (for example, the first 15 0s and the last 10 1s), a data set corresponding to the feature value and label is made for model training. .

As the preferred technical solution of the above technical solution, in step four, the chromosomes are processed separately and the samples are processed separately, because the data calculations are independent of each other and can be processed in parallel;

It can be seen from the above technical solutions that compared with the existing technology, the technical effects achieved are:

1) The present invention innovatively proposes to use the number of sequencing sequences (reads) that completely cover a certain region to identify coverage. In the past, to identify the coverage of a region, the number of single base sites covered by reads was first calculated, and then By averaging the number of reads covering multiple single bases. Previous methods were unable to detect inversions or transposon endpoints because the inversion or transposon endpoint is actually a site. After the sequencing data is mounted to the reference genome, the alignment software must include fault-tolerance parameters. Therefore, the mounted sequence will span 1-10 bp across the endpoint of the inversion or transposon, and the endpoint and nearby base positions will be covered by multiple reads. The average coverage of a region containing the endpoint is also non-zero, and The coverage average results of other continuous coverage areas that do not include endpoints are similar, so previous methods cannot detect endpoints. The present invention proposes to completely cover the number of reads of 39 bp, and use a step size of 20 to divide the window, eliminating the factor of mismatching of the mounting sequence. If there is an endpoint, it can be detected that the number of reads completely covering the 39 bp region including the endpoint is 0.

2) The present invention innovatively studies inversions or transposons at the population level, using second-generation sequencing data combined with linear pan-genomes to find the endpoints of inversions or transposons (also called breakpoints in the present invention), and can obtain richer In the past, genomes were assembled through third-generation sequencing and high-quality genomes were compared to obtain the positions of inversions or transposons. The cost was huge. Due to the limited number of high-quality genomes, the number of inversions or transposons obtained was limited. . The present invention can utilize the massive amount of published second-generation sequencing data and has obvious advantages.

3) The present invention innovatively applies an artificial intelligence model to identify inversions or transposon endpoints (breakpoints). The artificial intelligence model can extract abstract features, handle complex situations, and accurately identify endpoints. When the material sequencing data containing the endpoints are mounted on the linear pan-genome, the comparison results have obvious characteristics: the endpoints (the 39bp interval including the endpoints) are not completely covered by reads, and both the left and right sides of the endpoints can be fully covered by multiple reads. The traditional method is easier to judge a certain site or a certain area (a value). The judgment of the endpoint requires the combination of multiple values, and the order of the numerical values has an impact on the judgment result. The traditional method is difficult to handle. Artificial intelligence models are good at handling such nonlinear and feature extraction problems.

4) The calculations within the window of the present invention are independent of each other, and multi-thread parallel acceleration can be used to reduce program running time.

Description of the drawings

In order to explain the embodiments of the present invention or the technical solutions in the prior art more clearly, the drawings needed to be used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings in the following description are only These are embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on the provided drawings without exerting creative efforts.

Figure 1 is a schematic flow chart of the method for detecting population inversions or transposon endpoint genotypes based on linear pan-genome and artificial intelligence models according to the present invention;

Figure 2 is a schematic diagram showing that the reads of the present invention completely cover a certain area;

Figure 3 is a schematic diagram of the windowing method of the present invention for detecting genome level coverage (complete coverage);

Figure 4 is a structural diagram of the training data set generation and artificial intelligence model of the present invention;

Figure 5 is a schematic diagram of the windowing method of the present invention converting data structure for use in artificial intelligence models.

Detailed ways

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only some of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts fall within the scope of protection of the present invention.

Example 1

Please refer to Figures 1-4. This embodiment provides a method for detecting population inversions or transposon endpoint genotypes based on linear pan-genome and artificial intelligence models, including:

A. Statistics on the coverage of next-generation sequencing data at the genome level

A1. Extract DNA from population samples.

A2. Use ultrasonic waves to randomly fragment the DNA and build a library, and then use second-generation DNA sequencing technology to perform high-depth sequencing (greater than or equal to 20x). The current read length of second-generation sequencing data is 150bp.

A3. Use alignment software to align the high-depth sequencing of all samples to the linear pan-genome. The comparison software can be BWA or Samtools, both of which have no Chinese names. They are directly expressed in English in the industry and are used to extract target sequences.

A4. Take 39bp as a window and count the number of reads that are completely covered in this window. Among them, reads refers to the sequencing sequence generated during the second-generation sequencing process. It is a 150bp base sequence and can be directly expressed in English in the industry.

As shown in Figure 2, if the leftmost end of a read is on the left side of the window and the rightmost end of the read is on the right side of the window, it means that the read completely covers (crosses) the window, and the number of reads that the window is completely covered +1. In the actual operation process, the rightmost end of the read should be the sum of the leftmost end plus the number of bases matched to the linear pan-genome, plus the number of detected missing bases.

A5. Use the windowing method to perform windowing operations on the linear pan-genome. The window size is 39bp and the step size is 20bp. Analyze whether 39bp of each window is completely covered by reads. Merge the position information of the window and the number of fully covered reads into a corresponding data frame with 2 columns. In order to reduce subsequent computing costs, compare the number of fully covered reads with a threshold (3 in this example, if the sequencing depth is deeper, it can be increased appropriately). If it is greater than 3, use 1 to indicate that the 39bp window is completely covered by reads; if it is less than Equal to 3, 0 indicates that the window is not completely covered. Combine the window position information and the 0 and 1 numbers indicating the coverage to form a data frame containing two columns of numbers (Figure 3).

B. Artificial intelligence model training to identify inversions or transposon endpoints (breakpoints)

B1. Create a training set. The coverage of sequencing reads at the endpoint (breakpoint) of an inversion or transposon has obvious characteristics: the area before the breakpoint is fully covered by a large number of reads, and the area behind the breakpoint is also fully covered by a large number of reads. However, the 39bp window containing the breakpoint is not. Completely covered by reads. Therefore, in the data frame generated by A5, if a region (519bp, including 25 small windows) appears with multiple consecutive 1s, only one 0 appears, followed by multiple 1s, indicating that there is a breakpoint in this region. Through the visualization of actual genome third-generation sequencing data, we found that most breakpoints are one base position, but they may also be several consecutive bases or a small sequence. Therefore, after multiple 1s appear in succession, 1-10 0s may appear, followed by all 1s, which also indicates that the area contains a breakpoint. Using the windowing method, the middle (one or more) 0s may appear in different positions of the window, so the situation is complicated and difficult to judge by traditional algorithms. Create all data sets containing breakpoints (Figure 4): the number of windows is set to 25, and the feature values are arrays containing 25 numbers. Simulate the existence of 1 0, 0 may appear at any position in the array (25 numbers); simulate the existence of 2 0, 0 may appear at any position in the middle of the array... Create all data sets that do not contain breakpoints. When one or more 0s appear at the far left or right end of the array, it means that the window is the boundary of large fragment insertion and deletion; if there are more than 10 0s in the array, it means that the window Contains an indel. Set 25 consecutive 0 or 1 data as feature values, and whether there is a breakpoint (1 or 0) as a label value. Correspond the feature values and label values to form a data set.

B2. Create an artificial intelligence model and construct an artificial neural network model containing 3 hidden layers and 1 output layer. The number of contacts in the hidden layer is 64, 16, and 4 respectively. The activation function is ReLU. The output layer uses the sigmoid function to output a value. The model input is an array of 25 feature values, matched to the training set.

B3. Use the training set to train the artificial intelligence model, input the data set into the artificial intelligence model, perform calculations on the feature values and model parameters, output the results through the multi-layer network, compare the output results with the real label values, and use the gradient descent algorithm Improve model parameters to continuously narrow the gap between output results and real label values. Input multiple sets of data and train 4,000 times so that the artificial intelligence model can more accurately (99.9%) identify the existence of breakpoints.

C. Use artificial intelligence models to identify population inversions or transposon endpoints (breakpoints)

C1. For each sample, a data matrix representing the whole genome coverage is generated after step A5 (the first column is the window position information, and the second column is whether the window is completely covered). Use the windowing method to perform windowing operation on the data matrix, with 25 values as a window and 22 as the step size (Figure 4). The data matrix is converted into a matrix containing 25 eigenvalues (the column is 25, and the row names are renamed after the interval positions containing these 25 data).

Among them, when drawing windows, first determine that every 25 small windows are on the same chromosome. If the last window on a chromosome cannot contain 25 small windows (for example, when the first chromosome is windowed to the end, only the last 10 numbers are left). ), discard the window.

C2. 25 numbers in each row of the data matrix are input to the artificial intelligence model and output whether the area (the interval where the 25 small windows are located) contains a breakpoint.

C3. Organize the artificial intelligence output results so that position information and breakpoint information are output correspondingly. Organize the results according to chromosomes and output them to a file. Each sample outputs one file result.

D. Merging and filtering breakpoint data

D1. Since all samples use the same linear pan-genome as the reference genome, the first column of the output result file of step A5 is exactly the same and contains the position information of the window on the chromosome. After the C1 windowing method, the output result files of each sample still have the same number of lines, and the line names of each line are the same.

D2. Vertically merge the breakpoint files of all samples.

D3. The merged file contains a matrix with the row name of the breakpoint position information on the chromosome (a 519bp region) and the column name of each sample. This matrix represents preliminary breakpoint genotypes. Each row of the sample was screened based on the minimum allele frequency (MAF), and rows with MAF less than 0.05 were removed. Get the filtered breakpoint genotypes and output them to a file.

The application principle of the present invention will be described in detail below with reference to the accompanying drawings:

Example 2 Obtaining genome-level coverage of sequencing reads of a single sample based on high-depth sequencing

The chemical decontamination method (CTAB method, Cetyltrimethylammonium Bromide) was used to extract DNA from the leaves of rice population materials, and the IlluminaNovaSeq6000 sequencing platform was used to complete high-depth sequencing of the DNA. Each sample generated approximately 20Gb of data. The rice genome size is 400Mb and the sample depth is 50X. This example takes IRRI2K_91 as an example. After compression, the paired-end sequencing files of IRRI2K_91 sequencing data are 6.8Gb (R1.fq.gz) and 7.0Gb (R2.fq.gz) in size.

Use BWA software to align the sample paired-end sequencing data to the rice linear reference genome to obtain an alignment file in bam format. Use picard software to sort the bam files, mark duplicate positions and establish index operations. At this time, the bam size is approximately 8.1Gb.

Run the self-written program a0_read_one_bam_absolute_depth1.09.py and set the window size parameter for full coverage. In this example, it is set to 39bp and the step size parameter is set to 20bp. This program takes the bam file as the input file and outputs a csv file. The csv file contains two columns of data, separated by commas. The first column records the position of each window on the chromosome, and the second column records the window position corresponding to the first column. Completely covering the number of reads, the size of the output file IRRI2K_91absolute_depth1.09.csv is 347M, with a total of 21,771,459 lines.

Example 3: Training an artificial intelligence model to identify inversions or transposon endpoints (breakpoints)

Run the self-written program a1_create_training_data.py, set the window number parameter to 25, and set the breakpoint to contain a maximum of 10 consecutive 0s. The program will generate a training set, including multi-line feature values and corresponding label values, and output it as a csv file. .

The self-written program will use Python to call the tensorflow module to set the network structure of the artificial intelligence model, which includes 3 hidden layers, with the number of contacts being 64, 16, and 4 respectively. The activation function is ReLU, and the output layer uses the sigmoid function to output a value. The model input is an array of 25 feature values, matched to the procedurally generated training set. The self-written program inputs the training set into the model to train the model, sets the number of cycles, and after training 4,000 times, the model's prediction accuracy for the training set breakpoint reaches more than 99.9%, and the trained model is saved.

Example 4 Artificial Intelligence Model Detection of Inversions or Transposon Endpoints (Breakpoints)

Use the self-written program a2_depth_shape_cast.py to process the IRRI2K_91absolute_depth1.09.csv file, and set the minimum coverage threshold to 3. If the second column of IRRI2K_91absolute_depth1.09.csv is less than or equal to 3, it is converted to 0, and if it is greater than 3, it is converted to 1. . The programmed window size parameter is 25. The self-written program will use the windowing method to convert the two columns of data into matrix data containing 25 eigenvalues (Figure 4), and output it to the file IRRI2K_91absolute_depth_shape25.csv using the self-written program The program a3_detecting_breakpoint.py loads the trained model, detects breakpoints on each line in the IRRI2K_91absolute_depth_shape25.csv file, and outputs the position information and prediction results to an IRRI2K_91_breakpoint.csv file.

Perform the same steps sequentially on multiple samples in the population to obtain the IRRI2K_**_breakpoint.csv of multiple samples (** represents other numbers in the population). Since the execution of different samples does not depend on other samples, this step can use multi-process parallel operations to speed up the analysis process; since multiple samples are based on the same linear pan-genome as the reference genome, the final IRRI2K_**_breakpoint.csv file The location information is consistent. Using the pandas module in python, read all IRRI2K_**_breakpoint.csv and merge them into a genotype matrix, which contains position information as a column and the name of each sample as a column, and outputs it to a file. Use the self-written program a4_filter_population_breakpoint.py to screen the genotype matrix for minimum allele frequency (MAF), and obtain the screened population inversion or transposon endpoint genotypes.

Randomly select 30 positions where breakpoints were detected in the screened genotypes, and use IGV visualization software to conduct manual observation and detection of these 30 regions. After testing, these 30 positions are completely consistent with the inversion or transposon endpoints. Features, using the present invention can replace manual observation to detect breakpoints.

Each embodiment in this specification is described in a progressive manner. Each embodiment focuses on its differences from other embodiments. The same and similar parts between the various embodiments can be referred to each other. As for the device disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple. For relevant details, please refer to the description in the method section.

The above description of the disclosed embodiments enables those skilled in the art to make or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be practiced in other embodiments without departing from the spirit or scope of the invention. Thus, the present invention is not intended to be limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims (9)

1. A method for detecting population inversions or transposon endpoint genotypes based on linear pan-genome and artificial intelligence models, which is characterized by including the following steps: Step 1: Use comparison software to mount the population's high-depth sequencing data onto the linear pan-genome; Step 2: Based on the second-generation sequencing sequence coverage of each sample in the population, use the windowing method to detect whether each window is completely covered by the sequencing sequence, record the window position information and the number of reads that completely cover the window, and generate a data frame document; Step 3: Construct an artificial intelligence model and train the model through simulated data so that it can determine whether the continuous window coverage area contains inversion or transposon endpoints based on the fixed number, continuous window position and sequence number information obtained in step 2; Step 4: Use the model in Step 3 to scan all regions on a chromosome to obtain all inversion or transposon endpoint position information on the chromosome, scan multiple chromosomes in sequence, and summarize the inversions or transposons on all chromosomes of a sample. Socket endpoint information; Step 5: Based on the inversion or transposon endpoint information of different samples, organize and screen the population-level inversion or transposon endpoint genotype matrix. 2. The method according to claim 1, characterized in that in step one, the sequencing depth of the population high-depth sequencing data is greater than 20×. 3. The method according to claim 1, characterized in that: the window in step 2 refers to a 39 bp region on the genome, and the position of the window is named based on the combination of the chromosome and the middle position of the window; draw the window. The step size is 20bp. 4. The method according to claim 3, characterized in that: the complete coverage in step 2 means that the leftmost position of the sequencing sequence is on the left side of the window, and the left position of the sequencing sequence plus sequence matching is matched to the pan-genome. The sum of the number of bases and the number of detected missing bases is greater than the right position of the window. 5. The method according to claim 4, characterized in that: in step two, the window position information and the number of reads that completely cover the window are recorded to obtain a data frame containing two columns of data, the first of which is The first column is the window position, and the second column is the number of sequences that completely cover the window, which is saved in a csv file. 6. The method according to claim 5, characterized in that: the artificial intelligence model in step three is an artificial neural network model, which is a network containing 4 layers of fully connected layers, and the activation function of the last layer is a sigmoid function. , other activation functions are ReLU functions; the input of the model is an array of length 25, derived from the data frame csv file in step 2; Perform a windowing operation on the second column in the csv file. The window size is 25 and the step size is 22. The length of the genome sequence corresponding to each window is 39bp+(25-1)*20=519bp; the output of the model is a value, Used to determine whether the 519bp region contains inversions or transposon endpoints. 7. The method according to claim 6, characterized in that: the data frame csv file in step 2 needs to be preprocessed before being input into the model, and the number in the second column of the csv file that is less than or equal to 3 is defined as 0, which means No sequence complete coverage; a number greater than 3 is defined as 1, indicating complete sequence coverage. 8. The method according to claim 6, characterized in that: training the model through simulated data in step three: Among them, the simulation data is generated by artificial simulation, including situations where inversions or transposon endpoints obviously exist and situations where inversions or transposon endpoints obviously do not exist, and is produced into a data set corresponding to feature values and labels for use in the model. train. 9. The method according to claim 1, characterized in that in step four, chromosomes are processed separately and samples are processed separately, using independent parallel operations.