Mitochondrial whole genome genotype detection method and application thereof

Abstract

The invention provides a mitochondrial whole genome genotype detection method based on multiple PCR amplification sequencing and application of the mitochondrial whole genome genotype detection method. Specifically, the method comprises the following steps: constructing a reference sequence set which comprises a mitochondrial genome reference sequence and a nuclear mitochondrial sequence which are subjected to extension treatment; comparing the sequencing result with the reference sequence set, evaluating according to the comparison condition to obtain a comparison result, and taking open gap as an optimal comparison result when the score penalty setting principle of evaluation is that mismatch and open gap exist at the same time; multiple PCR primer treatment is carried out, data filtering is carried out on the compared reads on the basis of the treated primers, and the data filtering comprises non-specific amplification reads filtering; and carrying out genotype analysis according to the filtered reads, wherein the genotype analysis comprises right alignment on insertion and deletion in the reads, primer sequence removal and primer specificity correction of site genotypes. According to the detection method, through a series of treatment and filtration, the accuracy of mitochondrial whole genome genotype detection can be effectively improved.

Description

Technical Field

[0001] The present invention relates to the field of mitochondrial genome analysis, and in particular to a mitochondrial whole genome genotype detection method based on multiplex PCR amplification sequencing and related applications thereof. Background Art

[0002] At present, there are five main types of mitochondrial whole genome variation detection technologies: Sanger sequencing based on PCR amplification of mitochondrial DNA products, mitochondrial whole genome second generation (Next-Generation Sequencing, NGS) sequencing, mitochondrial whole genome sequencing based on multiple PCR amplification, mitochondrial whole genome sequencing based on probe capture, and mitochondrial whole genome sequencing based on third-generation sequencing technology. Specifically, Sanger sequencing based on PCR amplification of mitochondrial DNA products: first extract total DNA, design primers that can cover the entire mitochondrial genome, then use PCR technology to amplify the mitochondrial genome, and after purification, use the Sanger sequencing method to obtain the mitochondrial sequence, and finally compare and analyze with the mitochondrial reference sequence to obtain the genotype results of the mitochondrial whole genome site. Mitochondrial whole genome second generation sequencing: first extract total DNA, directly amplify the mitochondrial genome together with nuclear DNA or use long-distance PCR technology, randomly break it into small fragments, construct a small fragment sequencing library and perform high-throughput sequencing, and finally align the obtained sequencing reads to the mitochondrial reference sequence to obtain the genotype results of the mitochondrial whole genome site. Mitochondrial whole genome sequencing based on multiplex PCR amplification: It is an improved method based on the second-generation sequencing method. It uses multiplex PCR amplification technology to enrich mitochondrial DNA and then perform high-throughput sequencing by connecting sequencing adapters. Finally, the obtained sequencing reads are aligned to the mitochondrial reference sequence to obtain the genotype results of the mitochondrial whole genome sites. Mitochondrial whole genome sequencing based on probe capture: Probes are designed for the whole genome or target region, and liquid phase hybridization is performed with mitochondrial genomic DNA to capture and enrich mitochondrial DNA. Then, high-throughput sequencing is performed with the help of a sequencing platform to obtain the genotype results of the mitochondrial whole genome sites. Mitochondrial whole genome sequencing based on third-generation sequencing technology: The long read length of the third-generation sequencing technology can directly span the entire mitochondria for sequencing, so as to obtain the genotype results of the mitochondrial whole genome sites.

[0003] Analysis of the above technologies revealed that mitochondrial genotype detection based on Sanger sequencing has the advantages of high accuracy and low false positive rate, and is often used as a reference standard for high-throughput sequencing detection. The disadvantage is that the throughput is low, which leads to a relatively high cost relative to high-throughput sequencing, and generally only sequencing is performed on mitochondrial hypervariable regions. In addition, due to the limitations of the technology itself, it is unable to detect low-frequency (<10%) heteroplasmic sites. There are many mitochondrial whole genome analysis software based on second-generation sequencing (such as MToolBox, MitoSeek, mtDNA-Server, MitoHPC, Mutect2, etc.). These software are mainly aimed at double-end sequencing after random interruption. For the data of this library sequencing strategy, false positive sites can be filtered out well, but it is not suitable for analyzing data for PCR amplification or probe capture and then sequencing. Mitochondrial whole genome analysis software based on multiple PCR amplification or probe capture are currently some commercial software (such as SoftGenetics's GeneMarkerHTS, Verogen's Universal Analysis Software, Thermo Fisher's Converge Forensic Analysis Software, etc.). Although some analysis software based on second-generation sequencing (such as Mutect2, bcftools call, etc.) can also perform mitochondrial whole genome analysis of multiple PCR amplification, these software are not designed specifically for mitochondria, so these software do not make corresponding processing for the circular structure of mitochondria themselves and mutations introduced by multiple PCR primer sequences, resulting in unsatisfactory genotype detection results. Although one read of mitochondrial whole genome sequencing based on third-generation sequencing technology can better cover some complex regions on mitochondria, the sequencing quality of third-generation sequencing is not as good as that of second-generation sequencing (high-throughput sequencing). It is necessary to correct its own sequencing errors by increasing the sequencing depth, so the sequencing cost is more expensive than that of second-generation sequencing (high-throughput sequencing). At present, there are few mitochondrial analysis software based on third-generation sequencing.

[0004] In summary, the current mitochondrial whole genome analysis software for multiplex PCR amplification or probe capture are all commercial software, and general analysis software cannot perform mitochondrial whole genome genotype analysis well. Therefore, it is necessary to develop dedicated software for mitochondrial whole genome high-throughput sequencing data based on multiplex PCR amplification sequencing in order to better solve the mitochondrial whole genome genotype detection and analysis. Summary of the invention

[0005] In view of the above problems, the first aspect of the present invention provides a mitochondrial whole genome genotype detection method based on multiplex PCR amplification sequencing. The detection method comprises: constructing a reference sequence set, the reference sequence set comprising a mitochondrial genome reference sequence and a nuclear mitochondrial sequence that has been extended. The sequencing result is compared with the reference sequence set, and an evaluation is performed according to the comparison to obtain a comparison result, wherein the penalty setting principle for the evaluation is that when mismatches and gaps exist at the same time, gaps are used as the optimal comparison result. The multiplex PCR primers are processed, and the reads after the comparison are filtered based on the processed primers (i.e., the multiplex PCR primers originally designed for amplification, including the primer length and the position information of the primers on the reference mitochondrial genome when the primers are designed), and the data filtering includes filtering of non-specific amplification reads. Genotype analysis is performed based on the filtered reads, and the genotype analysis includes right alignment of the insertions and deletions in the filtered reads and analysis after primer sequence removal, and primer-specific filtering is performed on the filtered reads.

[0006] The extension process includes: adding the starting end sequence of the mitochondrial genome reference sequence to the end thereof for sequence extension, and the extension length refers to the sequencing read length.

[0007] The penalty mechanism in the above evaluation includes: adjusting the penalty rules of the existing comparison software - by reducing the penalty for gap opening and gap extension, and increasing the penalty for mismatch - to ensure that gap opening is given priority when there is a mismatch and gap opening.

[0008] The above-mentioned multiplex PCR primer processing includes analyzing the position information of the multiplex PCR primers on the mitochondrial genome, and for primer pairs with overlapping amplified fragments on the mitochondrial genome, extending the 3' ends of the upstream primer and the downstream primer closest to the overlapping region by several bases.

[0009] When the above data processing is for single-end sequencing, it can further include alignment strand filtering. Alignment strand filtering is mainly for single-end sequencing. For reads from the positive strand of the mitochondrial genome sequence, if they are aligned to the antisense strand of the mitochondrial genome sequence, such reads need to be filtered out, and vice versa. Paired-end sequencing is not affected.

[0010] The above-mentioned non-specific amplification reads filtering includes: 1) judging the threshold value of the number of mismatches of the corresponding primer sequences in the reads, that is, first calculating the number of mismatches of the corresponding primer sequences in the reads, and performing a first threshold judgment on the number of mismatches of the reads and the sequences corresponding to the primers after comparison. If the number of mismatches is higher than the first threshold, the reads are filtered; 2) for the reads whose mismatch number judgment meets the threshold, the position information of the primer corresponding to the sequence aligned to the reference sequence is judged, and the actual position information of the primer corresponding to the reads on the mitochondria is calculated based on the position information of the reads on the mitochondria after comparison. When the calculated primer position differs from the position when the primer is designed by more than a second threshold, the reads are filtered.

[0011] The above genotype analysis includes: firstly, right-aligning the insertions and deletions in the filtered reads and removing the primer sequences, and then counting the bases detected at each site on the mitochondrial genome and the corresponding depth information. At the same time, the primer-specific filtering in the genotype analysis includes the Fisher test for sites covered by both pairs of primers to test whether there are significant differences in the bases sequenced by the two pairs of primers at the sites, and selecting primers for sites with significant differences, and the selection of primers is based on the site being closer to the center of the primer amplification fragment.

[0012] In a second aspect, the present invention provides a detection device based on the detection method of the first aspect, the detection device comprising a memory and a processor, wherein the memory is used to store programs; the processor is used to implement the detection method of the first aspect by executing the programs stored in the memory.

[0013] In a third aspect, the present invention provides a computer-readable storage medium based on the analysis method of the first aspect, wherein a program is stored on the medium, and the program can be executed by a processor to implement the detection method of the first aspect.

[0014] The beneficial effect of the present invention is that the detection method provided by the present invention can effectively improve the accuracy of mitochondrial whole genome genotype detection through a series of processing and filtering. At the same time, the detection method based on the present invention has improved running time and reduced memory compared with other methods, can increase concurrent tasks, and greatly save analysis costs. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 This is a schematic diagram of the mitochondrial whole genome genotype detection process in a specific embodiment of the present invention;

[0016] Figure 2Schematic diagram of the extension process of the mitochondrial reference sequence in the present invention (a sequence of the mitochondrial starting part is added to the mitochondrial ending position, i.e., the bold sequence in the figure, the length of which is determined according to the length of the sequenced reads);

[0017] Figure 3 It is a schematic diagram for explaining the BWA comparison penalty in the prior art in the present invention;

[0018] Figure 4 An example diagram of adjusting the positions of overlapping multiplex PCR primers on mitochondria in multiplex PCR primer processing (F1 becomes F1', R2 becomes R2');

[0019] Figure 5 IGV software display in the embodiment of the present invention Figure 1 ; The IGV displays the filtering of non-specific amplified sequences caused by alignment errors (the upper part is the IGV software display of reads alignment before filtering, and the lower part is the IGV software display of reads alignment after filtering by the software of the present invention; each gray line represents a reads alignment result, and by comparing the results before and after filtering, it is found that the filtered reads are reads whose alignment positions are inconsistent with the preset primer positions);

[0020] Figure 6 IGV software display in the embodiment of the present invention Figure 2 , the IGV displays the filtering of non-specific amplification sequences caused by mismatches (the upper part is the IGV software display diagram of reads alignment before filtering, and the lower part is the IGV software display diagram of reads alignment after filtering by the software of the present invention; each gray line represents a reads alignment result, and by comparing the results before and after filtering, it is found that the filtered reads are all reads whose mismatch number in the primer binding region exceeds the set threshold). DETAILED DESCRIPTION

[0021] With the application of high-throughput sequencing technology, mitochondrial whole genome sequence analysis has gradually become the main means of mitochondrial DNA detection. Forensic workers can easily obtain comprehensive mitochondrial genome information through a single test, while improving the recognition rate of mitochondrial DNA detection, they can obtain more sample source information, which is very suitable for the inspection of trace and severely degraded samples. Mitochondrial whole genome sequencing technology based on multiplex PCR amplification is currently the main technical means for mitochondrial analysis in the field of forensic medicine. However, the complexity and heterogeneity of mitochondrial DNA sequences have always affected the accuracy of mitochondrial DNA genotype detection. The present invention conducts in-depth analysis and research on the problems currently existing in the mitochondrial whole genome genotype analysis based on multiplex PCR amplification, and finds the following problems existing in the prior art:

[0022] 1) Reference sequence problem: Since mitochondria are circular DNA, the actual reference sequences for comparison are generally linear. Therefore, for sequences near the head-to-tail connection point, there are generally poor comparison results or no comparison. Nuclear Mitochondrial DNA (NUMTs) contamination problem: NUMTs have a high homology with mitochondrial sequences. During PCR amplification, the NUMTs sequence is amplified.

[0023] 2) Problem of mismatching: Due to the existence of some complex regions in mitochondria (for example, the sequence ACCCCCCCTCCCCCG from 302 to 316), mismatching is prone to occur under the existing penalty rules, thus affecting the detection of genotypes.

[0024] 3) Genotype alignment problem: The current software outputs the detected insertions and deletions (Indel) in a left-aligned manner, while GA / T 1979-2022 "Forensic Science Mitochondrial DNA Next-Generation Sequencing Technical Specifications" requires that insertions and deletions be right-aligned.

[0025] 4) Problems caused by non-specific amplification: Due to the differences between primers and templates, dimers formed by primer pairing and errors introduced during primer amplification, non-specific amplification reads will occur.

[0026] 5) Genotype problems introduced by multiplex PCR amplification primers: Since the base at a certain position of the primer is inconsistent with the base of the template DNA, the base will continue to accumulate during the PCR amplification process and be sequenced, resulting in the detection of an incorrect mitochondrial genotype.

[0027] Based on the above problems, the present invention provides a series of filtering and processing methods, and provides a mitochondrial whole genome genotype detection method based on multiplex PCR amplification sequencing. For ease of description, taking a specific embodiment of the present invention as an example, the specific process of the detection method is as follows: Figure 1 As shown, including:

[0028] 1) Construct a reference sequence set.

[0029] 2) Compare the sequencing results with the reference sequence set, and evaluate the comparison to obtain the comparison results.

[0030] 3) Perform multiplex PCR amplification primer processing, and filter the aligned reads based on the processed primers.

[0031] 4) Perform genotype analysis based on the filtered reads.

[0032] The above detection method is described in detail as follows.

[0033] To construct a reference sequence set, different from the processing of the reference gene set in the prior art, the present invention needs to first extend the mitochondrial genome. Specifically, based on the circular DNA structure of mitochondria, a sequence of the beginning of the mitochondrial reference sequence (such as Figure 2 The black bold square shown in the figure) ensures that the reads that cross the head-tail connection point can be correctly aligned. It should be understood that in actual operation, for the reads that are aligned to this part, the mitochondrial coordinate information can be automatically converted internally. In addition, when constructing the reference sequence set, nuclear mitochondrial sequences (NUMTs) need to be added to exclude contamination from nuclear genome sequences. For reads that are only aligned to NUMTs, it is determined that the reads are derived from nuclear mitochondrial sequences and are filtered out.

[0034] In the comparison of the sequencing results with the reference sequence set, regarding the reads comparison processing, the bwa comparison is adopted in the specific embodiment of the present invention. The following is a specific description using the bwa comparison as an example, and those skilled in the art can make appropriate adjustments based on the actual selected comparison method to meet the comparison judgment principle in the present invention. In the bwa comparison of the prior art, when the sequencing reads are compared to the constructed reference gene set, there may be a variety of comparison situations, such as mismatching, gap opening, gap extension and other different processing for sites that fail to fully match, thereby generating a variety of comparison results. In order to judge the best comparison result of the reads under different comparison conditions, a set of penalty rules will be defined first (for example: complete match +2 points, mismatch -3 points, gap opening -8 points, gap extension -4 points, etc.), and then the results of different comparison situations of each read will be scored according to this set of rules, and the highest score is considered to be the best comparison. Since there are some complex regions in the mitochondrial sequence (for example, continuous C bases in the sequence), the default penalty rule in the prior art will output a comparison result that does not conform to the actual situation, resulting in misjudgment. Figure 3As shown, according to the penalty rule in the existing software, -15 (preferential mismatch)>-16 (preferential gap opening), therefore, -15 will be retained as the optimal alignment result, but this alignment result is not actually in line with the actual alignment result, so it is necessary to adjust the penalty for gap opening, so as to output an alignment result that is suitable for the mitochondrial genome. Therefore, through repeated research on the data, the present invention adjusts the penalty for mismatch, gap and gap extension on the basis of the existing penalty to balance the scores of mismatch and gap opening. Specifically, in the penalty rule of bwa ​​alignment, the penalty for gap opening and gap extension can be reduced, while the penalty for mismatch can be increased to ensure that gap opening is given priority in the case of mismatch and gap.

[0035] After the alignment, the multiplex PCR primers need to be processed. Since the error rate of the primers at the start of amplification is relatively high, after the design and optimization of the present invention, it is necessary to find the primers that overlap at the mitochondrial genome position (i.e., the corresponding amplified fragment) based on the position information of the multiplex PCR primers on the mitochondrial genome, and process these primers. It should be understood that when the fragment amplified by the primer does not contain the above-mentioned overlapping fragment, the primer does not need to be extended, and its sequence and corresponding position information are retained in the state at the time of design and used for subsequent data filtering. When the above-mentioned overlapping fragment appears, the present invention will automatically adjust the position of the overlapping primers on the mitochondria (such as Figure 4 As shown in FIG. 1 ), the 3' ends of the upstream primer (F1) and the downstream primer (R2) closest to the overlapping region are extended by several bases (for example, 5 bp), that is, the upstream primer (as shown in FIG. Figure 4 F1 in becomes F1') and the downstream primer (such as Figure 4 In the above analysis, the length of the primers is changed from R2 to R2'), thereby shortening the length of the overlapping amplified fragments. At the same time, this process also adjusts the specific position information of each primer pair with the overlap on the genome, which will affect the statistical results of the corresponding reads in the subsequent reads data filtering. Through this adjustment, the first few bases amplified by the primers can be avoided as much as possible during analysis. This avoids the errors introduced by sequencing or amplification affecting the subsequent mitochondrial genotype detection.

[0036] Based on the primers after the above treatment, the reads after the comparison are filtered. In the present invention, the data filtering needs to be combined with the primers to filter the non-specific amplification reads. The non-specific amplification reads filtering includes:

[0037] a. Since the mismatch between the primer and the mitochondrial binding region directly affects the binding efficiency of the primer and the mitochondrial binding region, the present invention first counts the number of mismatches between the primer and the mitochondrial binding region. If the number of mismatches is greater than a set first threshold (e.g., 4 mismatches), it is considered that the primer and the mitochondrial binding region have undergone non-specific amplification, thereby filtering out the reads (a specific example is as follows Figure 6 shown).

[0038] b. For the unfiltered reads in a, the position information of the primers on the mitochondria can be further calculated based on the position information of the reads on the mitochondria, and then compared with the pre-set position information of the primers on the mitochondria. If the two position information differ greatly (for example: 5bp), it is considered that the reads are caused by non-specific amplification and need to be filtered. This difference is the second threshold (a specific example is Figure 5 shown).

[0039] It should be noted that in the above-mentioned non-specific amplification reads filtering, if the primers corresponding to the reads belong to the extended primers in the multiplex PCR amplification primer processing, the relevant mismatch number should be calculated and processed with reference to the extended primers.

[0040] In addition to the above non-specific amplification reads filtering, the data filtering in the present invention also includes conventional filtering that general software has, such as but not limited to alignment chain filtering, alignment quality value filtering, low base quality value filtering, minimum reads length filtering, and reads overall alignment mismatch number filtering. Among them, the alignment chain filtering analysis is mainly for single-end (SE) sequencing. If it is known that the reads come from the sense chain or antisense chain of the genome sequence, the reads with incorrect alignment chain information are automatically filtered out by setting parameters (fwd or rev). The above filtering can be performed before non-specific reads filtering.

[0041] Regarding genotype analysis, the present invention first automatically adjusts the insertion or deletion in the reads to the rightmost end in accordance with the requirements of GA / T1979-2022 "Technical Specifications for Second Generation Sequencing of Mitochondrial DNA in Forensic Science". In addition, since the base at a certain position of the primer is inconsistent with the base of the mitochondrial DNA detection, the base is sequenced and the wrong genotype is detected. The present invention will detect whether the sequenced reads include the position area of ​​the primer according to the pre-set position information of the primer in the mitochondria. If it does, the primer area will be automatically filtered out (i.e., the corresponding primer sequence on the reads will be removed). The bases detected at each site on the mitochondrial genome and the corresponding depth information are counted, and the Fisher test is performed on the sites covered by both pairs of primers to check whether there is a significant difference in the bases sequenced by the two pairs of primers (for example: p-value<0.01, it is considered that there is a significant difference), and the primers are selected for the sites with significant differences. The judgment method of primer selection is based on the site being closer to the center of the primer amplification fragment. If the distance to the center of the primer amplification fragment is the same, the one with the higher average base quality value on different primers shall prevail. Finally, according to the ratio and depth threshold set for variant detection, the genotype and corresponding related information of each site on the mitochondrial genome are output, including but not limited to the total sequencing depth, sequencing depth of each base, detection quality value, average base quality value of the detected base, average base quality value of the upstream and downstream of the detected base, depth of each base on different primers, and other information.

[0042] The above description of the detection method mainly provides a detailed explanation of the special processing steps in the present invention that are different from the conventional processing methods in the prior art. Based on the principles of the detection method described above, technical personnel in this field can appropriately adjust the specific parameters and steps in combination with actual detection needs, and perform detection processing in combination with other necessary or preferred processing steps in the existing conventional processing methods.

[0043] The present invention will be further described in detail below through specific implementation modes in conjunction with the accompanying drawings.

[0044] Example

[0045] The present invention performs human mitochondrial whole genome genotype detection and analysis based on multiplex PCR high-throughput sequencing data. The main input files required include mitochondrial genome sequence (e.g., GeneBank accession number: NC_012920.1), nuclear mitochondrial DNA sequence NUMTs (downloadable from UCSC database), high-throughput sequencing data (Fastq format) and position information of multiplex PCR primers on the mitochondrial genome. The program automatically processes the mitochondrial genome sequence and multiplex PCR primers according to the sequence length of reads. The processed mitochondrial genome sequence and nuclear mitochondrial DNA sequence are used to construct a reference sequence set, and then the sequencing data is aligned to the reference sequence set. After filtering the alignment results according to the sequencing type, base statistics are performed for each site of the mitochondrial genome to obtain the final result.

[0046] Example 1: perl / path / to / mtVar.pl-ref / path / to / mito.fa-numts / path / to / numts.fa-fastq / path / to / test.fastq.gz-type fwd-rlen 100-align bwa-opt'-O1-B1-at 4'-extend 5-primer / path / to / primer.list-mlen 40-mis 4-qual 10-depth 100-rate 0.15

[0047] -sample test-outdir / path / to / outdir

[0048] Example 2: perl / path / to / mtVar.pl-ref / path / to / mito.fa-numts / path / to / numts.fa-fq1 / path / to / test_1.fq.gz-fq2 / path / to / test_2.fq.gz-type unknown-rlen100-align bwa-opt'-O1-B1-at 4'-extend 5-primer / path / to / primer.list-mlen 40-mis 4-qual10-depth 100-rate 0.15-sample test-outdir / path / to / outdir

[0049] Parameter Description:

[0050] -ref: mitochondrial genome reference sequence (fasta format) file

[0051] -numts: nuclear mitochondrial sequence (fasta format) file

[0052] -fastq, -fq1, -fq2: Input sequencing fastq (gz format) files, single-end sequencing (-fastq), double-end sequencing (-fq1 and -fq2). Double-end sequencing will be treated as single-end sequencing data, and it is best to splice before input.

[0053] -type sequencing type (three options: fwd, rev, unknown), fwd means that the reads come from the forward sequence of the reference sequence, rev means that the reads come from the reverse complementary sequence of the reference sequence, and unknown means that the source of the reads is unknown and no filtering is performed.

[0054] -rlen: Sequencing reads length parameter, the program will process the mitochondrial genome reference sequence according to this reads length. (Default: 100)

[0055] -align: Alignment software parameters, software that aligns reads to the reference genome sequence (default: bwa)

[0056] -opt: Alignment software parameters, by adjusting the gap opening and mismatching penalties, in the case of mismatches and gap openings, gap opening is prioritized. (For example: bwa defaults to: '-O 1-B 1-at 4', where the -O parameter indicates the gap opening penalty and the -B parameter indicates the mismatching penalty)

[0057] -extend: Maximum extension length of primers for overlapping amplified fragments. (Default: 5)

[0058] -primer: The location information of the primer pair on the mitochondrial genome. The program will find the primer pairs with overlapping amplified fragments and extend the corresponding primer pairs based on the location information of the primers on the mitochondrial genome. The maximum extension length is controlled by the -extend parameter.

[0059] -mlen: minimum read length, if the read length is less than this length, it will be filtered out. (Default: 40)

[0060] -mis: The maximum number of mismatches between reads and primers. If the maximum number of mismatches between primers and reads is greater than this value, the reads are considered to be non-specific amplification products and need to be filtered out. (Default: 4)

[0061] -qual: minimum quality value for alignment. If the alignment quality value of reads is less than this value, the reads will be filtered out.

[0062] (Default: 0)

[0063] -depth: The minimum sequencing depth for genotype analysis. If the sequencing depth of a site is less than this value, it is considered that the depth of the site is insufficient and genotype analysis cannot be performed. (Default: 100)

[0064] -rate: Pure heterozygous genotype determination threshold. If the genotype ratio of a site is greater than the threshold and less than 1-threshold, it is determined to be heterozygous. In other cases, it is homozygous. (Default: 0.15)

[0065] -sample: Sample name when outputting results

[0066] -outdir: output result path or directory

[0067] The results obtained after analysis include: VCF format mitochondrial genome site genotype (sampleID.vcf.gz), the number of reads for each pair of primers in the alignment (sampleID.pread.xls), mitochondrial variant site information file (sampleID.mutation.xls), the depth file of each site in the mitochondrial genome (sample.depth.xls.gz) and the relevant statistical result file of reads alignment (sampleID.map_stat.xls). The reads alignment file (before filtering) and the alignment file (BAM format) after filtering by the present invention are displayed using IGV (Integrative Genomics Viewer), as shown in Figure 1. Figure 5 and Figure 6 As shown. In the IGV diagram, the upper part represents the location information of the genome, and the lower part represents the alignment of the reads in this genomic interval. Under the default settings, each gray line in IGV represents the alignment of a read. Insertion mutations will be displayed as a purple "I", and deletion mutations will be displayed as a black "—". Point mutations will be colored according to the actual base, and the specific base will be displayed after zooming in. Green is A, red is T, blue is C, and yellow is G. Compared with the alignment results after filtering, there are more reads alignment results in the alignment results before filtering. The amplified fragments of these alignments are obviously shorter, and there are more mutations (color) and insertions and deletions (black). These reads are reads generated by non-specific amplification, which affect the normal typing results and need to be filtered out before analysis (such as Figure 5 shown).

[0068] By using the sequencing data of the commonly used standards in forensic science (including 2800M, 9947A, 9948, etc.) on the MGI G99 sequencing platform (single-end sequencing, read length: 100bp), the genotype accuracy of the three standards (three replicates) is 100% under the condition of the detection threshold of 15% (each base depth / total depth of the site) commonly used in the forensic field using the software developed by the present invention. Except for the primer sequence removal step (the primer sequence removal of the present invention is processed inside the software, and other software needs to remove the primer sequence by other software), other conditions are the same, and other software tested (such as: Mutect2 (version number: gatk-4.1.9.0), mutserve (version number: v2.0.0) and MitoHPC (version number: 20240222)) will have 1 to 2 sites with incorrect genotype detection, which may be due to the fact that these software are not specifically for multiplex PCR sequencing, and lack necessary processing in removing errors introduced by primers, resulting in incorrect genotype detection of some sites.

[0069] Those skilled in the art will appreciate that all or part of the functions of the various methods in the above-mentioned embodiments can be implemented by hardware or by computer programs. When all or part of the functions in the above-mentioned embodiments are implemented by computer programs, the program can be stored in a computer-readable storage medium, and the storage medium can include: read-only memory, random access memory, disk, optical disk, hard disk, etc., and the program is executed by a computer to implement the above-mentioned functions. For example, the program is stored in the memory of the device, and when the program in the memory is executed by the processor, all or part of the above-mentioned functions can be implemented. In addition, when all or part of the functions in the above-mentioned embodiments are implemented by computer programs, the program can also be stored in a storage medium such as a server, another computer, disk, optical disk, flash disk or mobile hard disk, and can be downloaded or copied and saved in the memory of the local device, or the system of the local device is updated, and when the program in the memory is executed by the processor, all or part of the functions in the above-mentioned embodiments can be implemented.

[0070] The above specific examples are used to illustrate the present invention, which is only used to help understand the present invention and is not intended to limit the present invention. For those skilled in the art, according to the concept of the present invention, some simple deductions, modifications or substitutions can be made.

Claims

1. A method for detecting mitochondrial whole genome genotype based on multiplex PCR amplification and sequencing, characterized in that: The method comprises: Constructing a reference sequence set, the reference sequence set comprising a mitochondrial genome reference sequence and a nuclear mitochondrial sequence that have been extended; Comparing the sequencing results with the reference sequence set, and evaluating the comparison to obtain the comparison results, wherein the penalty setting principle for the evaluation is that when mismatches and open gaps exist at the same time, the open gap is taken as the optimal comparison result; Perform multiple PCR primer processing, and filter the aligned reads based on the processed primers, including filtering non-specific amplification reads; Genotype analysis is performed based on the filtered reads, wherein the genotype analysis includes right-aligning the insertions and deletions in the filtered reads and removing the primer sequences before analysis, and performing primer-specific filtering on the filtered reads.

2. The detection method according to claim 1, characterized in that The extension process includes: adding the starting end sequence of the mitochondrial genome reference sequence to the end thereof to extend the sequence, and the extension length refers to the sequencing read length.

3. The detection method according to claim 1, characterized in that The penalty mechanism in the evaluation is to reduce the penalty for gap opening and gap extension, while increasing the penalty for mismatch.

4. The detection method according to claim 1, characterized in that The multiplex PCR primer processing includes analyzing the position information of the multiplex PCR primers on the mitochondrial genome, and for the primer pairs with overlapping amplified fragments on the mitochondrial genome, extending the 3' ends of the upstream primer and the downstream primer closest to the overlapping region by several bases; Preferably, the number of the several bases is 5 bp.

5. The detection method according to claim 1, characterized in that The non-specific amplification reads filtering includes: 1) judging the mismatch number threshold of the corresponding primer sequence in the reads, performing a first threshold judgment on the mismatch number of the reads and the corresponding sequence of the primer, and filtering the reads if the mismatch number is higher than the first threshold; 2) judging the position of the corresponding primer sequence for the reads that meet the mismatch number threshold judgment in 1), calculating the actual position information of the primer corresponding to the read on the mitochondria according to the position information of the reads on the mitochondria, and filtering the reads when the calculated primer position differs from the position when the primer is designed by more than a second threshold; Preferably, the first threshold is 4; Preferably, the second threshold is 5 bp.

6. The detection method according to claim 1 or 5, characterized in that: The data filtering for single-end sequencing also includes alignment chain filtering analysis: when it is known that the reads come from the positive chain or antisense chain of the genome sequence, the reads with incorrect alignment chain information are automatically filtered out by setting parameters.

7. The detection method according to claim 1, 5 or 6, characterized in that: The data filtering also includes: alignment quality value filtering, base quality value filtering, minimum reads length filtering, and reads overall alignment mismatch number filtering.

8. The detection method according to claim 1, characterized in that The genotype analysis includes: The insertions and deletions in the filtered reads were right-aligned and primer sequences were removed, and then the bases detected at each site on the mitochondrial genome and the corresponding depth information were counted; Primer specificity filtering is performed, and the primer specificity filtering includes performing a Fisher test on the sites covered by both pairs of primers to test whether there are significant differences in the sequencing bases of the two pairs of primers at the sites, and selecting primers for sites with significant differences, and the selection of primers is based on the site being closer to the center of the primer amplification fragment.

9. A mitochondrial whole genome genotype detection device based on multiplex PCR amplification sequencing, comprising: Memory, used to store programs; A processor, configured to implement the method according to any one of claims 1 to 8 by executing the program stored in the memory.

10. A computer-readable storage medium, characterized in that: The medium stores a program, which can be executed by a processor to implement the method according to any one of claims 1 to 8.

Images