Sequencing data generation method based on Tn5 transposase

Abstract

The invention discloses a sequencing data generation method based on Tn5 transposase, and belongs to the field of bioinformatics. The method comprises the following steps: for a reference genome, generating a variation haplotype sequence with a special code based on a preset variation rule or a variation file; constructing a basic sequencing feature model based on the real sequencing data; generating a fragment sequence after insertion of the Tn5 transposase based on a Tn5 transposase feature model, wherein the Tn5 transposase feature model comprises low complexity of 7-9bp overlapping of adjacent fragments and base site fragment types; generating a double-end sequencing result with a joint based on the basic sequencing feature model; sequencing errors are introduced based on a high-precision sequencing error model. According to the method, the technical defect that the Tn5 transposase database building data features are not simulated in an existing sequencing data generation method is overcome, the key features of real sequencing data are reserved, and data support is provided for the expression of an evaluation analysis tool on a Tn5 type data set.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics technology and relates to a method for generating sequencing data based on Tn5 transposase. Background Technology

[0002] Cells, as the basic structural and functional units of life, are the cornerstone of organisms. Their state directly affects organ function and overall health, making the deciphering of cellular molecular mechanisms crucial for revealing disease mechanisms and developing novel therapies. With the development of Next-Generation Sequencing (NGS) technology, its high-throughput and low-cost characteristics have revolutionized the research paradigm at the cellular and molecular levels. NGS, through massively parallel sequencing and modular library construction technologies, has achieved precise analysis from population-based cell sequencing to single-cell resolution, such as single-cell transcriptome sequencing (scRNA-seq), single-cell chromatin accessibility sequencing (scATAC-seq), single-cell whole-genome sequencing, and multi-omics.

[0003] Before sequencing, sequencing libraries need to be constructed. In recent years, many library construction methods have utilized Tn5 transposase. As a core functional enzyme in next-generation sequencing library preparation, Tn5 transposase can simultaneously achieve random DNA fragmentation and adapter sequence integration, greatly simplifying the library construction process and significantly improving efficiency. Furthermore, Tn5 transposase requires less sample volume, enhancing the sensitivity of the detection system. However, the sequencing data obtained after Tn5 transposase library construction differs significantly from conventional library construction methods, mainly in the following two aspects: Firstly, adjacent fragments have a 7-9 bp overlap; secondly, unlike DNA polymerases with strand substitution activity such as φ29, the start and end positions of the fragmented product of the same DNA molecule after Tn5 transposase fragmentation are fixed, and read stacking does not occur after alignment and deduplication.

[0004] After obtaining sequencing data, it is necessary to develop or use analysis tools and algorithms to extract useful information from the data. However, the characteristics of real sequencing data are unknown, making it difficult to determine whether tools or algorithms can correctly process the data and obtain correct results. In the bioinformatics field, sequencing data generation methods based on real data characteristics have become a key platform for evaluating the performance of analysis tools. These methods generate synthetic data containing preset benchmark values, providing a benchmark dataset for algorithm validation. However, existing methods generally neglect the characteristics of Tn5 transposases, and the fragmentation patterns modeled in Tn5 transposases are homogeneous with traditional randomized library construction methods, leading to systematic biases in the evaluation of Tn5-related data analysis algorithms. Summary of the Invention

[0005] To address the shortcomings of existing technologies, the present invention aims to provide a sequencing data generation method based on Tn5 transposase. This method integrates the fragmentation characteristics of Tn5 transposase with the basic characteristics of real sequencing data to generate high-fidelity Tn5 transposase library construction and sequencing simulation data. This data can be used as a benchmark dataset to evaluate the performance of existing tools, algorithms, and analysis procedures on Tn5 type data, and can also help in the development of new tools for Tn5 type data.

[0006] The objective of this invention is achieved through the following technical solution:

[0007] Firstly, a method for generating sequencing data based on Tn5 transposase is provided. The method includes: introducing variations into a reference genome to obtain a variant haplotype sequence whose positional mapping relationship is known; constructing a basic sequencing feature model based on real sequencing data, including insert length distribution, coverage distribution, and quality value distribution; for the variant haplotype sequence, sequentially extracting fragments along the haplotype sequence according to the insert length distribution of the real sequencing data, setting adjacent fragments to have 7-9 bp overlap, and generating fragment sequences after Tn5 transposase insertion; dividing the reference genome into continuous non-overlapping regions, counting the number of fragments in each region after deduplication of the real sequencing data as the number of fragment noise models, and adjusting the number of fragments extracted from each region according to the insert length distribution to equal the number of fragment noise models; generating paired-end sequencing results with adapters based on the basic sequencing feature model; extracting error features from the real sequencing data and introducing sequencing errors into the paired-end sequencing results.

[0008] Furthermore, the output variant haplotype sequence introduces variations based on the original reference genome. These variations are specified by single nucleotide variant files, copy number variant files, or generated by preset variation rules. The single nucleotide variant file is a standard VCF (Variant Call Format) file. The copy number variant file is a variant of the standard BED (BrowserExtensible Data) file, consisting of four columns: chromosome name, start position, end position, and copy number. The position coordinates follow a 0-base start and a 1-base end. The preset variation rules include mutation rate, the proportion of insertion / deletion variations among all variations, the probability of extension of insertion / deletion variations, the maximum length of the insertion / deletion portion, the ratio of transitions and transversions in single nucleotide variants, and the proportion of homozygous variations.

[0009] Furthermore, the basic sequencing feature model consists of an insert length distribution model, a coverage model, and a quality value model. Before generating the feature model, the sequencing data needs to be preprocessed, including adapter trimming, alignment with the reference genome, and alignment result filtering, to finally obtain a high-quality BAM (Binary Alignment / Map format) file. This file is a compressed version of the SAM (Sequence Alignment / Map) file, which has smaller storage space and faster processing speed.

[0010] Furthermore, the inserted fragment length distribution model is constructed as follows: the number of reads corresponding to different inserted fragment lengths in the BAM file is counted, specifically represented by the values ​​greater than or equal to 0 in the ninth column. Then, the frequency of different inserted fragment lengths is calculated, and the relationship between the frequency and the inserted fragment length is fitted using a Gaussian kernel.

[0011] Furthermore, the coverage model is constructed as follows: the reference genome is divided into continuous non-overlapping intervals of 1000 bp in length, the start and end positions of each pair of reads in the bam file on the reference genome are obtained, and the total number of bases in each interval is counted based on the start and end positions.

[0012] Furthermore, the quality value model consists of two parts: the probability distribution of each quality value at each position of R1 and R2 is calculated from the fastq file of the real sequencing data, and the conditional probability distribution of the current base quality value under the first few base sequence combinations.

[0013] Furthermore, the positional mapping between the introduced haplotype sequence and the reference genome is achieved through the following special encoding: the A, C, T, and G bases of the inserted portion are encoded as H, I, J, and K, respectively; the bases of the deleted portion are encoded as Z; and the bases of regions with smaller copy numbers are encoded as Y.

[0014] Furthermore, the generated insert fragments are adjusted using a mathematical model and a fragment noise model, and their distribution exhibits characteristics of real Tn5 transposase insert fragments. The mathematical model stipulates that adjacent fragments generated from the same haplotype file have 7-9 bp overlap, and after alignment with the reference genome, each site is covered by only one fragment, except for the overlapping regions of adjacent fragments. The fragment noise model divides the reference genome into non-overlapping regions, counts the number of fragments in each region after deduplication of the real data, and the region in which a fragment is located is determined by the midpoint of the fragment. The fragment noise model adjusts the fragments generated by the mathematical model using the following two strategies: if the number of fragments in the mathematical model is greater than that in the fragment noise model for a certain region, the number of fragments in the fragment noise model is downsampled; if the number of fragments in the mathematical model is less than that in the fragment noise model for a certain region, the missing fragment is generated by decoding the variant haplotype sequence of that region and combining it with the insert fragment length distribution model.

[0015] Furthermore, the generation of paired-end sequencing results with adapters based on the basic sequencing feature model includes: allocating the number of reads per chromosome according to chromosome coverage and the proportion of total coverage; using the ratio of the normalized weight of the fragment start position coverage to the number of fragments in the interval as the sampling probability of the fragment; generating read pairs with different alignment directions; and completing read pairs with shorter insertion fragment lengths (less than 150bp) by using adapter sequences and N-dominant random sequences.

[0016] Furthermore, the extraction of error features from real sequencing data includes: based on real sequencing data, firstly, obtaining the overlapping region of R1 and R2, identifying the different positions of R1 and R2 bases in the overlapping region, statistically analyzing the frequency of the bases at that position in all data, and identifying those with fewer occurrences as sequencing errors, and then calculating the probability of sequencing errors at each position of R1 and R2, the probability of sequencing errors under the combination of the first three and the last base sequence, the type of error and the probability of its corresponding quality value when sequencing errors occur.

[0017] Secondly, a method for evaluating the performance of a variant detection tool is provided, wherein benchmark data is constructed using the method described in the claims, detection results are generated using the variant detection tool, and variants with known mapping relationships to reference genomic locations are compared to evaluate the performance of the variant detection tool.

[0018] The beneficial effects of this invention are as follows:

[0019] 1. By constructing a hybrid generative model through multi-source feature extraction and adopting a model collaborative adjustment strategy, the synthetic data maintains a high degree of consistency with real experimental data in terms of statistical distribution and biological characteristics, and can be used as a benchmark dataset.

[0020] 2. An innovative mathematical model of the Tn5 transposase insertion fragment was established and coupled with the fragment noise characteristics modeled based on real data, so that the generated benchmark dataset accurately reproduces the characteristics of the Tn5 transposase library construction and sequencing results.

[0021] 3. By identifying overlapping regions in paired-end sequencing, sequencing errors in these regions are detected and statistically analyzed, resulting in a high-precision sequencing error model that allows the generated benchmark dataset to retain systematic errors introduced by the sequencer.

[0022] 4. The modular architecture design supports compatibility between the generation system and third-party components such as haplotype simulators, sequencing quality models, and sequencing error models, ensuring the method's portability and universality in different application scenarios.

[0023] 5. It can flexibly generate complete datasets containing raw sequencing reads and quality control indicators, which are suitable for benchmark data construction and related bioinformatics tool evaluation in Tn5 transposase application scenarios such as whole genome sequencing, chromatin accessibility analysis, and transcriptome sequencing. Attached Figure Description

[0024] Figure 1 This is a flowchart of the sequencing data generation method based on Tn5 transposase provided by the present invention;

[0025] Figure 2 This is a visualization of the generated data and the original data IGV provided by the present invention;

[0026] Figure 3 This is a comparison chart of the length distribution of the generated data and the original data inserted segments provided by the present invention;

[0027] Figure 4 This is a comparison chart of the coverage distribution between the generated data and the original data provided by this invention; wherein, Figure 4 In the diagram, A represents the coverage data. Figure 4 In the diagram, B represents a scatter plot showing the consistency between the generated data and the original data coverage.

[0028] Figure 5 This is a comparison chart of the average sequencing quality of the generated data and the original data provided by this invention; wherein, Figure 5 In the graph, A represents a comparison of the mean quality value at each position of R1 between the generated data and the original data. Figure 5 B in the graph is a comparison chart of the mean mass values ​​at each position of R2;

[0029] Figure 6 This is a comparison chart between the sequencing error model constructed from the generated data provided by this invention and the actual generated data; wherein, Figure 6 In the diagram, A represents a comparison of R1 error cases. Figure 6 B in the graph represents a comparison of R2 error cases;

[0030] Figure 7 This is a comparison chart of the sequencing error models of the generated data and the original data provided by this invention; wherein, Figure 7 In the diagram, A represents a comparison of R1 error cases. Figure 7 B in the graph represents a comparison of R2 error cases;

[0031] Figure 8 This is a performance statistics chart of different comparison tools provided by the present invention in generating data; wherein, Figure 8 In the graph, A represents the accuracy statistics. Figure 8 A statistical graph showing the predominant chromosome percentage of B chromosomes. Figure 8 C in the graph represents the runtime statistics.

[0032] Figure 9 This invention describes the performance of different SNV detection tools in generating data. Detailed Implementation

[0033] To enable those skilled in the art to better understand the technical solutions in the embodiments of this application, the technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art should fall within the protection scope of the embodiments of this application.

[0034] Example 1

[0035] In this specific embodiment, single-cell whole-genome data from the publicly available database META-CS is used. The META-CS method fragments genomic DNA from single-cell lysates based on the Tn5 transposase. Double-stranded DNA is labeled and complementary strands are amplified using a single-tube reaction, resulting in DNA fragments with specific transposon adapters (containing the ME core sequence) at both ends. Figure 1 The benchmark dataset for simulating this sequencing data includes the following steps:

[0036] Dataset Acquisition. The dataset used in this embodiment comes from the NCBI public database, numbered SRR8981979, and its sequencing method is 2*150bp. This data is used to generate the benchmark dataset in this embodiment.

[0037] Data preprocessing. The `cutadapt` tool was used to trim the adapters from the raw FastQ files. The 5' adapter sequences trimmed from R1 and R2 were both AGATGTGTATAAGAGACAG (ME sequence, SEQ ID NO.1), and the 3' adapter sequences were both CTGTCTCTTATACACATCT (rME sequence, SEQ ID NO.2). The shortest adapter matching length was 3 bp, and the acceptable error rate was 10% (rounded down). The `bwa-mem2` tool was used to align the trimmed FastQ files with the reference genome hg38. The alignment results were then piped to the `samtools` tool, and the output was converted to a binary BAM file using the `samtools view` function with the parameter `-Sb`. The `sambamba sort` tool was used to sort the aligned BAM files by alignment position. Use `samtools view` to filter the alignment results, retaining reads that correctly match the reference genome at both ends and have an alignment quality of 60 or higher, excluding secondary and supplementary alignments, to obtain a high-quality, un-duplicated BAM file, avoiding interference from duplicate or low-quality data. The specific parameters are `-f 2 -F 256 -F 2048 -q 60`. Obtain the high-quality, deduplicated BAM file, and use `sambamba markdup` to perform initial deduplication on the pre-quality control BAM file. Then, use `samtools view` to filter it, and use the smaller, filtered BAM file for secondary deduplication. The secondary deduplication filtering conditions are: retain only one fragment with the same alignment start and end positions, and retain only reads whose R1 and R2 values ​​are not filtered.

[0038] Generate a uniquely encoded haplotype sequence based on preset mutation rules or mutation files. If introducing mutations via mutation files, the file introducing single nucleotide variants must be a standard vcf (Variant Call Format) file, containing a header file, chromosome name, variant location, reference sequence, and variant sequence. The file introducing copy number variants must be a variant of the standard bed (Browser Extensible Data) file, consisting of four columns: chromosome name, start position, end position, and copy number, with position coordinates following a 0-base start and a 1-base end. If introducing mutations based on preset mutation rules, the mutation-related parameters need to be modified according to the actual situation, including mutation rate, the proportion of insertion / deletion variants among all variants, the probability of insertion / deletion variant extension, the maximum length of the insertion / deletion portion, the ratio of transitions and transversions in single nucleotide variants, and the proportion of homozygous variants. In this embodiment, mutations are introduced through the 00-common_all.vcf.gz file in the dbSNP database. Note that the chromosome name in the vcf file must be consistent with the reference genome. First, a conflict check is performed on the variants in the file, retaining only those without positional conflicts. The proportion of homozygous variants specified by the parameters then determines whether each variant is applied to all haplotypes or only one. In this embodiment, the copy number on all chromosomes is set to 2, generating diploid data that does not include copy number variations. The output files in this section include different haplotype files and the introduced variants and their homozygous information. In the output haplotype sequences, the A, C, T, and G bases of the insertion portion of the insertion variant are encoded as H, I, J, and K, respectively; the bases of the deletion portion of the deletion variant are encoded as Z; and the bases in regions with lower copy numbers are encoded as Y.

[0039] A basic sequencing feature model was constructed based on real sequencing data. Using a high-quality, un-duplicated BAM file as input, an insert length distribution model and a coverage model were constructed. The insert length distribution model was constructed as follows: the number of reads corresponding to different insert lengths in the BAM file was counted, specifically represented by values ​​greater than or equal to 0 in the ninth column. Then, the frequency of different insert lengths was calculated, and the relationship between the frequency and the insert length was fitted using a Gaussian kernel. The coverage model was constructed as follows: the reference genome was divided into continuous non-overlapping intervals of 1000 bp in length. The start and end positions of each pair of reads in the BAM file on the reference genome were obtained, and the total number of bases in each interval was counted based on the start and end positions. Using the original FASTQ file as input, a quality value model was constructed. The quality value model consists of two parts: R1 and R2 probability distributions for each quality value at each position; and the conditional probability distribution of the current base quality value under different combinations of the first two base sequences.

[0040] The Tn5 transposase feature model was used to generate the inserted fragment sequence. The Tn5 transposase feature model includes low complexity of 7-9 bp overlap between adjacent fragments and the types of base sites, represented by a mathematical model and a fragment noise model. The mathematical model stipulates that the process of generating the insert fragment from the haplotype sequence is independent. Taking a haplotype sequence as an example: special characters such as H, I, J, K, Z, and Y in the sequence are decoded to obtain the positional mapping relationship between the haplotype sequence and the reference genome hg38; the region to be used to generate the insert fragment is selected according to the region specified in the input parameters. This option can exclude regions that do not need to be considered or simulate chromatin openness in the chromatin accessibility data generation; for each region, the fragment is truncated from the smallest position, and the fragment length is generated by the insert fragment length distribution model through random probability; adjacent fragments have 7-9 bp overlap. Ideally, after alignment and deduplication, diploid sequencing data should have two distinct fragments covering a specific base site in the reference genome (ignoring overlapping regions from adjacent fragments of the same haplotype). However, in real-world sequencing data, due to various factors such as PCR bias, sequencing coverage, and library preparation defects, each location may correspond to fewer or more fragments than the number of plotypes. The fragment noise model reflects this characteristic of real-world data. The fragment noise model divides the reference genome into non-overlapping regions; in this example, the region size is 300 bp, and it is recommended to slightly exceed the average length of the inserted fragments. The number of fragments in each region is counted in the high-quality deduplication BAM file, and the region in which a fragment is located is determined by its midpoint. The fragment noise model adjusts the inserted fragments using two strategies: if the mathematical model has more fragments than the fragment noise model in a given region, the fragment noise model's count is downsampled; if the mathematical model has fewer fragments than the fragment noise model in a given region, the missing fragment is generated by randomly truncating the relevant haplotype sequences after decoding, based on the inserted fragment length distribution model.

[0041] Paired-end sequencing results with adapters are generated based on the basic sequencing feature model. The input parameters specify the total number of reads generated, the paired-end adapter sequences for R1 and R2, and the sequencing read length. The total number of reads is the number of filtered reads from the actual data divided by 97%-99%. The 5' adapters for R1 and R2 are both specified as TCGGCACCAATAAGATGTGTATAAGAGACAG (SEQ ID NO.3), and the 3' adapter is specified as CTGTCTCTTATACACATCTTATTGGTGCCGA (SEQ ID NO.4). The sequencing read length is 2*150bp. The coverage sum of all bases on each chromosome is obtained from the coverage model. The proportion of each chromosome's coverage sum to the total coverage is calculated. Based on this proportion, the total number of reads is allocated to each chromosome (rounded down). Finally, any excess reads are allocated to the chromosome with the most reads. Each chromosome is processed using a separate thread. Fragments for generating sequencing reads are selected through random sampling. The sampling probability of each fragment is the coverage normalized weight of the fragment's starting position divided by the number of fragments in its interval. When generating reads, the 5' adapter sequence is introduced first, followed by the fragment sequence, until the read length reaches 150 bp. If the fragment sequence length is insufficient, the 3' adapter sequence is introduced. If still insufficient, a random sequence with N as the predominant byte (80% probability of N, 5% probability of A, C, T, G) is introduced. The fragment sequence has a 50% probability of being introduced into reads in the forward direction (consistent with the reference genome direction) and a 50% probability of being introduced into reads in the reverse complementary direction. Sequencing quality is generated based on the read sequence and quality value model. Quality values ​​are first generated based on base positions and contextual information. If the contextual information is not present in the model, quality values ​​are generated only based on base positions.

[0042] Sequencing errors are introduced based on a high-precision sequencing error model. The raw FastQ data files are trimmed with the 3' adapter CTGTCTCTTATACACATCT (SEQ ID NO.2) using cutadapt, and then aligned with the reference genome hg38. The alignment results are filtered to obtain a high-quality BAM file. The alignment results of the reads in the BAM file are mapped onto the reference genome to obtain the overlapping region of R1 and R2. Positions where R1 and R2 bases differ within the overlapping region are selected, preserving their positional information on the reference genome, their positional information on R1 and R2, the base information of R1 and R2 at that position, and the corresponding sequencing quality information. The number of occurrences of each base at the preserved positions in all reads is counted in the BAM file. Reads with significantly fewer occurrences are considered true sequencing errors, thus obtaining the position of the sequencing error in R1 or R2, the error pattern (from which base changed to which base), the read context information at the error location (the first three bases and the last base), and the quality value of the error location. If the frequency of each base at a certain position is approximately the same, the existing SNP database is used for identification; if identification is still impossible, the position is disregarded. R1 and R2 sequencing error models are constructed, and the probability of an error at each position (the number of errors at a given position divided by the number of times that position appears in all overlapping regions), the probability of an error under each context condition (the number of errors under a given context condition divided by the number of times that context sequence appears in all overlapping regions), the probability of different error patterns under different contexts, and the probability of different sequencing quality values ​​at different error positions under a certain error condition are calculated. Sequencing errors are then introduced into the generated data, and the quality values ​​are modified accordingly. The probability of an error at a particular base is obtained by weighting the positional error rate and the context error rate. The error pattern is determined through context analysis, further determining the sequencing quality value after the error occurs.

[0043] Example 2

[0044] This example uses the data described in Example 1 to compare the generated data with the original data in multiple dimensions, further illustrating the high fidelity of the sequencing data generation method based on Tn5 transposase of the present invention.

[0045] The generated data exhibits Tn5 transposase library construction characteristics and has identical copies to the original data. After deduplication and quality control by comparing the generated data with the original data, the resulting high-quality BAM file is displayed in IGV software, such as... Figure 2 The top and bottom images show the generated data and the real data, respectively. There is a 7-9bp overlap between adjacent segments. Both segments come from two copies. The segment connected by the red line in the image is one of the copies.

[0046] The generated data possesses the basic characteristics of the original data. The frequency of inserted fragment lengths was statistically analyzed in the high-quality BAM files obtained after quality control and comparison of the generated and original data. Figure 3 ; Collect coverage data with intervals of 1000bp from the high-quality BAM files obtained after quality control and comparison of the generated data and the original data, such as Figure 4 Calculate the mean mass value at each position of R1 and R2 from the generated data and the original data fastq files respectively, such as... Figure 5 The comparison results show that the generated data has a high correlation with the original data.

[0047] Sequencing errors in the generated data exhibit characteristics of the raw data. To verify that the sequencing error model constructed using the overlapping regions of paired-end sequencing reads can reflect the situation in all regions, a sequencing error model was constructed from the generated data and compared with the actual error situation in the generated data, such as... Figure 6 The two are highly correlated, meaning the sequencing error model can reflect the original characteristics. To verify that the sequencing errors in the generated data possess the characteristics of the original data, the error rates at each position in R1 and R2 of the input sequencing error model and the generated data were compared, such as... Figure 7 The two are highly consistent.

[0048] Example 3

[0049] This example uses the data described in Example 1, taking common alignment tools and SNV detection tools as examples, to further illustrate the practicality of the sequencing data generation method based on Tn5 transposase of the present invention.

[0050] Currently, there are four mainstream whole-genome alignment tools: bwa mem, bwa-mem2, bowtie2, and minimap2. These were evaluated based on accuracy, main chromosome alignment rate, and runtime. All tools utilize 40 CPU cores, pre-build indexes, and then use index-related files for alignment, outputting SAM files which are piped to samtools view (8 cores) for conversion into BAM files. Alignment correctness is determined based on information in the read names. Accuracy is obtained by dividing the number of reads with correct alignment start positions by the total number of reads. The main chromosome alignment rate is the number of reads aligned to all 24 autosomes divided by the total number of reads. Runtime is obtained using the `time` command in Linux. Figure 8 bwa-mem and bwa-mem2 have the highest accuracy and the highest main chromosome rate, with bwa-mem2 having a faster running speed. minimap2 has the fastest running speed, but its accuracy and main chromosome rate are the worst.

[0051] The performance of GATK4 and lianti pileup (used in META-CS) in detecting SNVs was evaluated. Four processing methods were used: GATK4 base quality correction, variant site identification, variant quality correction, and variant filtering for high-quality deduplicated BAM files; GATK4 base quality correction, variant site identification, variant quality correction, and variant filtering for high-quality BAM files; lianti pileup processing for high-quality deduplicated BAM files; and lianti pileup processing for high-quality BAM files. Precision, recall, and F1 score were calculated for variants with coverage greater than or equal to a certain value under each processing method. 85% of the variants in the SNP database were treated as germline variants, and bulk sequencing data (randomly selected fragment start positions, fragment lengths generated by an insert fragment length model) were generated from these variants. The remaining 15% of variants were treated as somatic variants. Figure 9 The horizontal axis represents the minimum coverage of the considered variants, with "yes" and "no" indicating whether variant quality control was performed. The vertical axis represents the SNV detection tool used, with "yes" and "no" indicating whether deduplication was performed. As shown in the graph, deduplication is not conducive to identifying SNVs in Tn5 type data. GATK4's filtering operation leads to the exclusion of a large number of true variants, making it unsuitable for Tn5 type data.

[0052] Other embodiments of this application will readily occur to those skilled in the art upon consideration of the specification and practice of the disclosure herein. This application is intended to cover any variations, uses, or adaptations of this application that follow the general principles of this application and include common knowledge or customary techniques in the art not disclosed herein. The specification and embodiments are to be considered exemplary only.

[0053] It should be understood that this application is not limited to the steps described above and shown in the accompanying drawings, and various modifications and changes can be made without departing from its scope.

Claims

1. A method for generating sequencing data based on Tn5 transposase, characterized in that, The method includes: Based on the introduction of variations into the reference genome, variant haplotype sequences with known positional mapping relationships are obtained; Basic sequencing feature models were constructed based on real sequencing data, including insert length distribution, coverage distribution, and quality value distribution. For variant haplotype sequences, fragments are sequentially extracted along the haplotype sequence based on the insertion fragment length distribution of the real sequencing data, with adjacent fragments having a 7-9 bp overlap, to generate fragment sequences after Tn5 transposase insertion; the reference genome is divided into continuous non-overlapping regions, and the number of fragments in each region after deduplication of the real sequencing data is counted as the number of fragment noise models, and the number of fragments extracted from each region according to the insertion fragment length distribution is adjusted to be equal to the number of fragment noise models; Generate paired-end sequencing results with adapters based on the basic sequencing feature model; Error features are extracted from real sequencing data, and sequencing errors are introduced into the paired-end sequencing results.

2. The sequencing data generation method based on Tn5 transposase according to claim 1, characterized in that, The mutations are specified by single nucleotide variant files, copy number variant files, or generated by preset mutation rules; The single nucleotide variant file is a standard VCF file; The copy number variation file consists of four columns, including chromosome name, start position, end position, and copy number, where the position coordinates follow a 0-base start and a 1-base end. The preset mutation rules include mutation rate, the proportion of insertion / deletion mutations to all mutations, the probability of insertion / deletion mutation extension, the maximum length of the insertion / deletion part, the ratio of conversion and transversion in single nucleotide mutations, and the proportion of homozygous mutations.

3. The sequencing data generation method based on Tn5 transposase according to claim 1, characterized in that, The insert fragment length distribution was obtained through the following steps: statistically analyzing the insert fragment length frequency after alignment with real sequencing data, and using a Gaussian kernel to fit the relationship between the frequency and the insert fragment length. The coverage distribution was obtained through the following steps: the reference genome was divided into continuous non-overlapping intervals of 1000 bp in length, and the total number of bases in each interval was counted after the alignment of the real sequencing data. The quality value distribution is obtained through the following steps: statistically analyzing the probability of each quality value at each position in R1 and R2 from the fastq file of the real sequencing data, and the conditional probability distribution of the current base quality value under the first few base sequence combinations.

4. The sequencing data generation method based on Tn5 transposase according to claim 1, characterized in that, The mapping relationship between the introduced haplotype sequence and the reference genome is achieved through the following special encoding: the A, C, T, and G bases of the inserted part are encoded as H, I, J, and K, respectively; the bases of the deleted part are encoded as Z; and the bases of the region with a small copy number are encoded as Y.

5. The sequencing data generation method based on Tn5 transposase according to claim 1, characterized in that, Based on the open chromatin regions of real sequencing data or by manual setting, select regions to generate insert fragments. For each region, extract fragments sequentially along the haplotype sequence according to the insert fragment length distribution of real sequencing data. When counting the number of fragments in each interval of the actual sequencing data after deduplication, the interval in which the fragment is located is determined by the midpoint of the fragment; The adjustment process involves selecting the number of segments from each interval based on the inserted segment length distribution, which is equal to the number of segment noise models. When the number of segments selected from the inserted segment length distribution within an interval exceeds the number of segment noise models, the number of segment noise models is downsampled. When the number of segments selected from the inserted segment length distribution within an interval is less than the number of segment noise models, the missing segments are generated by decoding the mutated haplotype sequence of that interval and combining it with the inserted segment length distribution.

6. The sequencing data generation method based on Tn5 transposase according to claim 1, characterized in that, The generation of paired-end sequencing results with adapters based on the basic sequencing feature model includes: allocating the number of reads per chromosome according to chromosome coverage and the proportion of total coverage; using the ratio of the normalized weight of the fragment start position coverage to the number of fragments in the interval as the sampling probability of the fragment; generating read pairs with different alignment directions; and completing read pairs with insert fragment lengths of less than 150bp by using adapter sequences and N-dominant random sequences.

7. The sequencing data generation method based on Tn5 transposase according to claim 1, characterized in that, The extraction of error features from real sequencing data includes: Based on real sequencing data, the overlapping region of R1 and R2 is obtained, and the positions where R1 and R2 bases differ in the overlapping region are identified. The frequency of the bases at these positions is counted in all data. Those with fewer occurrences are sequencing errors. Then, the probability of sequencing errors at each position of R1 and R2, the probability of sequencing errors under the combination of the first three and the last base sequences, the type of error and the probability of its corresponding quality value are calculated.

8. A method for evaluating the performance of a variant detection tool, characterized in that, A baseline dataset is constructed using the method described in any one of claims 1-7, and a variant detection tool is used to generate detection results. The results are then compared with variants whose mapping relationship with the reference genome position is known, in order to evaluate the performance of the variant detection tool.

Images