CN115631789B - A Pan-Genome-Based Population Joint Variation Detection Method - Google Patents

Abstract

A pan-genome-based population joint variation detection method relates to the technical field of gene biology. The invention aims to solve the problems of low comparison efficiency, weak variation detection sensitivity and inaccurate detection results in the current population joint variation detection method. The present invention includes: obtaining a reference genome and a tested sample genome, using the reference genome and the tested sample genome to identify candidate variation sites on the sample genome, and classifying the candidate variation sites; Use candidate mutation sites to construct local haplotype sequences, and integrate the constructed local haplotype sequences among samples to obtain the integrated haplotype sequence set; use the integrated haplotype sequences to detect the genome of the tested sample The genotypes of the variant sites; traverse all the tested sample genomes in the population, and repeat steps 1 to 3 to obtain all the variant genotypes in the population. The invention is used for joint variation detection of population genome.

Description

A pan-genome-based population joint variation detection method

Technical Field

The present invention relates to the field of gene biology technology, and in particular to a method for detecting population joint variation based on a pan-genome.

Background Art

Variant detection is the detection of nucleotide differences between individual samples and reference genomes. It is the basis of many genomic studies and can discover the association between genomic variation and important phenotypes and diseases. In the genome, there are a large number of variations with small base sequence changes, such as SNVs and Indels. They are the most common types of variation in human heritable variation, accounting for more than 90% of known variations. The genome sequences of different individuals have 99.99% homology, and SNV/Indel plays a decisive role in the personalized characteristics and phenotypes of organisms.

With the advancement of high-throughput sequencing technology, more and more bioinformatics tools have been developed to detect genomic variations. Among them, a series of variation detection methods represented by GATK have been widely used and have good performance in detecting single-sample variations. Different groups have certain differences in inheritance and variation, thus showing different genetic characteristics. In order to characterize the differences between different groups and the frequency of variation within a group, joint variation detection based on a large population is required. The essence of large-scale genome projects and personal genome sequencing tasks is to use a series of high-reliability and high-precision data analysis algorithms to analyze the genetic and variation information contained in individual and group genomes, characterize the differences between different groups, discover pathogenic genes or susceptibility genes for various diseases, reveal the molecular mechanisms of disease occurrence and development, and achieve accurate diagnosis and treatment of diseases.

At present, the joint variation detection method is mainly used in large-scale genomic data analysis. The current mainstream algorithm is the BWA+GATK detection workflow developed by scientific research institutions in the United States, the United Kingdom and other countries, which has formed a closed algorithm ecosystem. However, this method still has problems in comparison efficiency, variation detection sensitivity and accuracy, and it is becoming increasingly difficult to adapt to the ever-increasing scale of large-scale genome projects.

Summary of the invention

The purpose of the present invention is to solve the problems of low comparison efficiency, weak variation detection sensitivity and inaccurate detection results in the current population joint variation detection methods, and propose a population joint variation detection method based on pan-genome.

A method for joint variation detection of a population based on a pan-genome. The specific process is as follows:

Step 1: Obtain the reference genome and the genome of the sample to be tested, and use the reference genome and the genome of the sample to be tested to identify candidate mutation sites on the genome of the sample to be tested, and classify the candidate mutation sites:

First, the reference genome is divided into interval blocks, and the reference genome after the interval blocks are compared with the genome of the sample to be tested to obtain a comparison information set; then the comparison information is stacked according to the coordinate order of the comparison information in the reference genome to obtain a stacking result; then the stacking results are counted, and the categories of candidate variant sites on the genome of the sample to be tested are identified using the stacking results;

The categories of candidate variation sites on the genome of the sample being tested include: high-confidence candidate variation sites, low-confidence candidate variation sites, and low-complexity region candidate variation sites;

Step 2: constructing local haplotype sequences using candidate variant sites according to the types of candidate variant sites, and integrating the constructed local haplotype sequences among samples to obtain an integrated haplotype sequence set;

Step 3: Use the integrated haplotype sequence to detect the genotype of the variant site in the genome of the sample being tested:

First, non-gap alignment is performed on the haploid sequence set and the sequencing fragments of the alignment information set to obtain the sum of the cumulative mismatch base masses of each sequencing fragment non-gap alignment as S, and then the sum of the cumulative mismatch base masses of each sequencing fragment non-gap alignment S is used to calculate the genotype of the variant site;

Step 4: Traverse the genomes of all tested samples in the population and repeat steps 1 to 3 to obtain all variant genotypes in the population.

Furthermore, the step 1 of obtaining a reference genome and a test sample genome, and using the reference genome and the test sample genome to identify candidate variant sites on the test sample genome, and classifying the candidate variant sites, includes the following steps:

Step 11: Obtain a reference genome, and divide the reference genome into interval blocks to obtain a divided reference genome B=B ₁ ,B ₂ ,…,B _n ;

Wherein, B ₁ , B ₂ ,…, B _n represent each divided interval block, and n is the total number of interval blocks;

The interval block division is performed according to the format of "chromosome: starting position - ending position";

There is no overlap between the interval blocks;

Step 1 and 2: obtain the genome of the sample to be tested, compare the genome of the sample to be tested with the divided reference genome, and obtain a comparison information set in which the genome of the sample to be tested falls within _Bi ;

Step 13: stack all the alignment information falling within _Bi obtained in step 12 according to the coordinate order of the alignment information in the reference genome to obtain a stacking result;

Step 14: Statistics are performed on the stacking results obtained in step 13 in each interval block of the divided reference genome, and the sequence information of the stacking results in each interval block is recorded through a dictionary variable. The candidate mutation sites of the tested sample are identified based on the stacking result sequence information recorded in the dictionary, and the candidate mutation sites are classified.

Furthermore, the step 13 stacks all the alignment information obtained in step 12 and falling within _Bi according to the coordinate order of the alignment information in the reference genome, including the following steps:

First, the alignment information of unaligned, suboptimal alignment and split alignment is filtered out according to the value of the alignment information ‘flag’ field in the alignment information set; then the alignment information with a value less than 10 in the ‘mapq’ field is filtered out to obtain a filtered alignment information set;

Among them, the value of the ‘flag’ field of the alignment information of the unaligned alignment is 4; the value of the ‘flag’ field of the alignment information of the suboptimal alignment is 256; the value of the ‘flag’ field of the alignment information of the split alignment is 2048;

Then, the filtered alignment information is stacked according to the coordinate order of the alignment information in the reference genome using an alignment result analysis tool to obtain a stacking result;

The stacking result includes: reference sequence name, genome coordinates, reference base, number of sequencing fragments covering the current position, base sequence of the current position, and base quality ASCII code sequence of the current position.

Furthermore, the step 14 performs statistics on the stacking results obtained in step 13 in each interval block of the divided reference genome, and records the sequence information of the stacking results in each interval block through a dictionary variable, and identifies the candidate variant sites of the tested sample according to the stacking result sequence information recorded in the dictionary, and classifies the candidate variant sites, including the following steps:

Step 141: Perform statistics on the stacking results obtained in step 13 in each interval block of the divided reference genome, and record the bases, inserted and deleted sequence information, and the frequency information of each base sequence at each genomic position in each interval block through dictionary variables:

First, the base sequence information of bases, insertions, and deletions occurring at each genomic position within the interval block is recorded through dictionary variables;

The dictionary variables include: pileup_dict;

The pileup_dict is used to obtain the number of occurrences of base sequences, inserted and deleted base sequences;

Then, the frequency of occurrence of each base is obtained based on the number of occurrences of each base sequence recorded in the pileup_dict dictionary variable:

P _j1∈keys ＝C _j1 /T

Where T is the total number of occurrences of all bases, keys = {A, C, G, T, I, D} is a set of bases, j1 is a base in keys, and C _j1 represents the number of occurrences of base j1;

Step 142: Use step 141 to obtain the frequency of occurrence of each base to obtain candidate variant sites:

If there is at least one non-reference sequence base at a genomic position with a frequency P _j1∈keys >T _alt , then the current position is a candidate variant site;

Where, T _alt is the frequency threshold;

The candidate variant site is represented by a 4-tuple, including: genomic position, reference base, non-reference base, and frequency of occurrence of the non-reference base sequence;

Step 143: Classify the candidate variant sites obtained in step 142.

Furthermore, the candidate variant sites obtained in step 142 are classified in step 143, specifically:

The candidate variant sites are divided into the following three categories: high confidence candidate variant sites HCs, low confidence candidate variant sites LCs, and low complexity region candidate variant sites LCCs;

Candidate variant sites are HCs when they meet the following three characteristics:

(1) The occurrence frequency P _HCs of a non-reference base is greater than the second frequency threshold T _HCs , that is, P _HCs ＞T _HCs ;

(2) No other candidate variant sites are included in the window extending w-bp upstream and downstream of the current candidate variant site on the genome of the sample being tested;

(3) The location of the current candidate variant site on the genome of the sample being tested does not belong to the low-complexity region, including: the GAP region after genome splicing constitutes a heterochromatin domain; segmental duplication; pseudo-autosomal region of sex chromosomes; short tandem repeats;

Candidate variant sites are LCs when they meet one of the following four characteristics:

(1) When the occurrence frequency P _LCs of a non-reference base is between the second frequency threshold T _HCs set by the user and the third frequency threshold T _LCs preset by the user, that is, T _LCs < P _LCs < T _HCs ;

(2) The location of the current candidate variation site is a candidate multi-allelic variation site;

(3) Other candidate variant sites are included in the window extending w-bp upstream and downstream of the current candidate variant site on the genome of the sample being tested;

(4) The current location belongs to a high complexity area;

When a candidate mutation site does not meet the judgment criteria of HCs and LCs, it is classified as LCCs.

Furthermore, the step 2 of constructing a local haplotype sequence using the candidate variation site according to the type of the candidate variation site, and integrating the constructed local haplotype sequence between samples to obtain an integrated haplotype sequence set includes the following steps:

Step 21: Divide each interval block of the divided reference genome into windows of fixed length, traverse all candidate mutation sites in the interval block, and assign all candidate mutation sites in the interval block to the optimal window:

First, each interval block of the divided reference genome is divided into multiple windows of length W, and the overlap length between adjacent windows is set;

Then, all candidate variant sites are traversed in each interval block of the divided reference genome, and the candidate variant sites are assigned to the optimal window:

When the candidate variant site is offset from the start of the window When , the current window is the optimal window for the current candidate mutation site;

Step 22: constructing a haplotype sequence of the sample being tested in each window according to the type of candidate variant site assigned to the optimal window;

Step 2: Integrate the haplotype sequences of the tested samples in the window to obtain an integrated haplotype sequence set:

First, within a window, all haplotype sequences of the tested samples are traversed, and duplicate sequences are removed according to the variants they contain, and finally a set of non-duplicate haplotype sequences of the samples within the window is generated;

Then, all windows are traversed, and a haplotype sequence set h = {h ₁ , h ₂ , ..., h _n } of each window is generated in turn.

Furthermore, in step 22, constructing the haplotype sequence of the sample to be tested in each window according to the type of candidate variant site assigned to the optimal window comprises the following steps:

Step 221: for each window, determine whether each window contains a candidate variation site. If it does not contain a candidate variation site, the reference genome sequence is the haplotype sequence of the current window. If it contains a candidate variation site, determine whether the window contains a low-confidence candidate variation site. If the window contains a low-confidence candidate variation site, execute step 222 to obtain the constructed haplotype sequence of the sample being tested. If the window contains only a high-confidence candidate variation site, execute step 223 to obtain the constructed haplotype sequence of the sample being tested.

Step 222: Use the de Bruijn graph-based splicing method within the window to obtain the haplotype sequence of the sample being tested:

Extract all sequencing fragments within the window and filter out the sequencing fragments with an overlap of less than 10% with the window;

Set the starting k-mer to 41bp, use the remaining sequenced fragments to construct a local deBruijn graph, use the depth-first search algorithm to identify all possible consensus sequences from the local deBruijn graph, and sort the sequenced fragments according to their support for all possible consensus sequences, and select the contigs with the highest support;

Add all contigs generated when k-mer = 41 to the sequenced fragment, set k-mer = 46, reconstruct the local deBruijn graph, and regenerate contigs;

Set the step length to 5bp, end k-mer = 75bp, and repeat the above operations until the last step generates contigs, i.e. the haplotype sequence of the sample being tested;

Step 223: Use the alignment analysis tool to extract all sequencing fragments covering the current window from the alignment information set, and obtain the sequencing fragment ID that supports each candidate variation site in the alignment information set; take the intersection of the sequencing fragment IDs that support any two candidate variation sites. If the number of sequencing fragments contained in the intersection is not less than the quantity threshold, it is considered that the mutations corresponding to the two candidate variation sites appear in the same haplotype sequence, and the bases at the corresponding positions on the reference sequence are replaced with the variant bases to construct a haplotype sequence to obtain a haplotype sequence; if the two candidate variation sites do not have any support from the same sequencing fragments, they are considered to come from different haplotype sequences, and haplotype sequences are constructed separately according to the candidate mutations and reference genome sequences contained in different haplotypes.

Furthermore, the step 3 of using the integrated haplotype sequence to detect the genotype of the variable site of the genome of the sample being tested comprises the following steps:

Step 3.1. Perform non-gap alignment of the haploid sequence set obtained in step 2 with the sequencing fragments of the alignment information set to obtain the sum of the cumulative mismatch base masses S of the non-gap alignment of each sequencing fragment:

Step 3: Group the candidate variant sites and haplotype sequences in each window:

First, all haplotype sequences containing candidate mutation sites are selected from the haplotype sequence set, which are denoted as H; the remaining haplotype sequences do not contain candidate mutations, which are denoted as R;

Then, the candidate variant sites contained in the tested samples are grouped: if the distance between two candidate variant sites is less than the distance threshold, the two candidate variant sites are regarded as a group;

Step 3-2: Use the candidate variant sites and the grouped haplotype sequences in each window to obtain the sum of the cumulative mismatch base masses S of the non-gap alignment of each sequencing fragment:

First, for each group of candidate variant sites obtained in step 3, all haplotype sequences containing the group of variants are selected from H and recorded as h, and all sequencing fragments within the window are extracted from the alignment information set;

Then, the complete matches between all sequencing fragments and haplotype sequences and the reference sequence are mapped respectively, the longest complete match block is obtained, and the starting alignment position is determined as p;

Then, starting from the alignment start position, single base extension is performed to the left and right sides respectively. When a base mismatch is encountered during the extension process, the base mass of the mismatched base is accumulated; when the extension reaches both ends of the sequence or the sum of the accumulated mismatched base masses is greater than the base mass threshold, the base extension is stopped;

Finally, select non-repeating haplotype sequences from R and record them as r, perform non-gap base extension on each sequenced fragment, and obtain the sum of the cumulative mismatch base masses of non-gap alignment of each sequencing fragment and record them as S;

Step 32: Use the sum S of the cumulative mismatch base masses of non-gap alignments of each sequencing fragment obtained in step 31 to obtain the genotype of the variant site.

Furthermore, the step 32 uses the sum of the cumulative mismatch base masses S of the non-gap alignment of each sequencing fragment obtained in step 31 to obtain the genotype of the variant site, including the following steps:

Step 321: Use the sum of the cumulative mismatch base masses S of the non-gap alignment of each sequencing fragment obtained in step 31 to obtain the optimal alignment probability set for each sequencing fragment after all sequencing fragments are aligned with the haplotype sequences in h and r respectively.

Obtained in the following manner: suppose the list of haplotype sequences containing variants is h'={h ₁ , h ₂ , ..., _hn' }, the probability set of sequencing fragments aligned to haplotypes is p={p ₁ , p ₂ , ..., _nn' }, select the maximum probability p _max as the probability of the current sequencing fragment aligned to the haplotype containing the variant, and the corresponding haplotype is h _max ;

Where k∈[1,m] is the probability index in the optimal alignment probability set for each sequencing fragment; m is the total number of probabilities obtained, and n’ is the number of haplotypes in the window;

pj is obtained by the following formula:

Among them, j∈[1,n′] is the label of the probability of the sequenced fragment being aligned to n haplotypes;

Step 322: Use the P and Pr obtained in step 321 to calculate the likelihood probabilities of different genotypes, and take the genotype with the largest likelihood function as the mutated genotype in the current window, as follows:

Furthermore, the step 4 traverses all the genomes of the tested samples in the population, and repeats steps 1 to 3 to obtain all the variant genotypes in the population, specifically:

Traverse the candidate mutation sites in all windows. For a certain mutation site, if all samples have no mutation signal at this site, or have mutation signal but the genotype of the detected mutation is 0/0, it is considered that all individuals in the current population do not have this mutation and will not be output; if there is a sample with a genotype of 0/1 or 1/1 at this site, traverse all samples in turn and output the genotype of each tested sample at this site. If a site is multi-allelic, split it into multiple double alleles and output them in turn.

The beneficial effects of the present invention are:

The present invention utilizes a large amount of local haplotype information in the population genome to re-align sequencing fragments, and utilizes the population information to perform feature statistical analysis, effectively improving the sensitivity and accuracy of variation detection; in addition, the present invention utilizes the high-speed alignment characteristics of local haplotypes to significantly improve the efficiency of variation detection.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram of the present invention.

DETAILED DESCRIPTION

Specific implementation method 1: As shown in FIG1 , the specific process of a population joint variation detection method based on a pan-genome in this implementation method is as follows:

Step 1: Obtain a reference genome and a test sample genome, and use the reference genome and the test sample genome to identify candidate mutation sites on the test sample genome, and classify the candidate mutation sites, including the following steps:

Step 11: Obtain a reference genome and divide the reference genome into interval blocks to obtain a divided reference genome B = {B ₁ , B ₂ , …, B _n }:

The length of the interval block division of the reference genome is 10 Mbps;

The interval block division is performed according to the format of "chromosome: starting position - ending position";

There is no overlap between the interval blocks;

Wherein, B ₁ , B ₂ ,…, B _n represent each divided interval block respectively, n is the total number of interval blocks, i=1,..,n;

Step 1 and 2: obtain the genome of the sample to be tested, compare the genome of the sample to be tested with the divided reference genome, obtain the comparison position where the genome of the sample to be tested falls within _Bi , and thus obtain the comparison information set where the genome of the sample to be tested falls within _Bi ;

Among them, the comparison information is stored in the following three file formats: SAM, BAM, and CRAM;

Step 13: stack all the alignment information obtained in step 1 and 2 that falls within _Bi according to the coordinate order of the alignment information in the reference genome to obtain the stacking result:

First, before stacking (projection), the alignment information of unaligned (flag = 4), suboptimal alignment (flag = 256) and split alignment (flag = 2048) is filtered out according to the value of the ‘flag’ field in the alignment information; the alignment information with a value less than 10 in the ‘mapq’ field is filtered out to obtain a filtered alignment information set;

Then, the filtered alignment information set is input into the mpileup function module in the existing alignment result analysis tool SAMtools to pile up the coordinates in the reference genome, that is, to project the spatial position and obtain the stacking result;

The stacking result includes the base coverage of each position in the reference genome and the corresponding base quality, etc. There are 6 columns in total, namely: reference sequence name, genome coordinates, reference base, number of sequencing fragments covering the current position, base sequence of the current position, and ASCII code sequence of the base quality of the current position.

The sample genome with candidate mutations is mainly reflected in the base sequence at the current position, which contains relevant information such as matches, mismatches, insertions, deletions, alignment strands and alignment quality, and the structure is relatively complex. The explanation is as follows: 1) ‘.’ represents a match with the positive strand of the reference sequence; 2) ‘,’ represents a match with the negative strand of the reference sequence; 3) ‘ATCGN’ represents a mismatch on the positive strand; 4) ‘atcgn’ represents a mismatch on the negative strand; 5) ‘*’ represents an ambiguous base; 6) ‘^’ represents that the matched base is the beginning of a read; the ascii code immediately following ‘^’ minus 33 represents the alignment quality; these two symbols modify the following bases, and the base immediately following them (ATCGatc gNn) represents the first base of the read; 7) ‘$’ represents the end of a read, and this symbol modifies the previous base; 8) The regular expression ‘\+[0-9]+[ACGTNacgtn]+’ represents the base inserted after the site, for example, ‘+3agg’ represents the insertion of 3 bases here; 9) The regular expression ‘-[0-9]+[ACGTNacgtn]+’ represents the base deleted after the site, for example, ‘-4CTGA’ represents the deletion of 4 bases here;

Note that after deleting the "^any character", "$", "-[0-9]+[ACGTNacgtn]+", and "+[0-9]+[ACGTNacgtn]+" in the fifth column, the number of bases is consistent with the number of alignment quality ASCII characters in the sixth column.

Step 14: statistics the stacking results obtained in step 13 in each interval block of the divided reference genome, and record the sequence information of the stacking results in each interval block through a dictionary variable, and identify the candidate variant sites of the tested sample according to the stacking result sequence information recorded in the dictionary, and classify the candidate variant sites;

Step 141: Count the stacking results obtained in step 13 in each interval block of the divided reference genome, and record the bases, inserted and deleted sequence information and the frequency information of each base sequence at each genome position in each interval block through two dictionary variables;

First, the base sequence information of the bases, insertions, and deletions that appear at each genomic position within the interval block is recorded through a dictionary variable:

The dictionary variable is pileup_dict, which is used to obtain the number of occurrences of base sequences, inserted and deleted base sequences, specifically: Key is the genome position, Value is all possible base sequences covering the current position; all possible base sequences covering the current position include matches, mismatches ("ACGTN"), insertions ("I") and deletions ("D"), and the number of occurrences of each base sequence is recorded. For example: pileup_dict = {"A":Ca,"C":Cc,"G":Cg,"T":Ct,"N":Cn,"I":Ci,"D":Cd}, where A, C, G, T represent four bases, I and D represent insertion and deletion marks, and N represents an unknown base. Ca, Cc, Cg, Ct, Ci, Cd, Cn respectively represent the number of times the above-mentioned base sequences appear at the current position.

Then, the frequency of occurrence of each base is obtained based on the number of occurrences of each base sequence recorded in the pileup_dict dictionary variable:

Where, T = ∑ _j1∈keys C _j1 is the total number of occurrences of all bases, C _j1 represents the number of occurrences of base j1, keys = {A, C, G, T, I, D}, j1 is a base in keys;

Step 142: Obtain candidate variant sites using the frequency of occurrence of each base obtained in step 141:

If there is at least one non-reference sequence base at a genomic position that satisfies its frequency of occurrence P _j1∈keys ＞T _alt (default value: 0.2), where T _alt is the frequency threshold preset by the user, the current position is called a candidate variant site and is represented by a 4-tuple, which records the genomic position (ctgpos), reference base (refB), non-reference base (altB), and frequency of occurrence of the non-reference base sequence (P _alt ).

The bases include: reference bases and non-reference bases, where non-reference bases represent possible variations;

It is worth noting that if there are multiple non-reference sequence bases whose frequencies are all greater than the frequency threshold T _alt , the current position is called a candidate multi-allelic site, and multiple possible non-reference base information of the site is recorded at the same time.

Step 143: Classify the candidate variant sites obtained in step 142:

The candidate variant sites are divided into the following three categories: high confidence candidates (HCs), low confidence candidates (LCs) and low-complexity candidates (LCCs);

When a candidate variant site has the following three characteristics, it is classified as HCs:

(1) The occurrence frequency P _HCs of the non-reference base sequence classified as HCs is greater than the second frequency threshold T _HCs preset by the user, that is, P _HCs ＞T _HCs (default: 0.5);

(2) In the window of w-bp (w = 100bp by default) extending upstream and downstream from the position of the current candidate variant site on the genome of the sample being tested, no other candidate variant sites are included. When at least two candidate variant sites appear in the same window, due to the systematic deviation of the alignment backbone or the deviation of the Smith-Waterman scoring system, the two candidate variant sites that are close to each other will have position deviation and variant sequence deviation.

(3) The current genomic position does not belong to the low-complexity regions, including: the GAP region after genome splicing constitutes the heterochromatin domain (constitutive heterochromatin domains); segmental duplications (segmental duplications); the pseudo-autosomal regions of the sex chromosomes (the pseudo-autosomal regions of the sex chromosomes); short tandem repeats (short tandem repeats);

The reference genome is obtained by splicing sequencing fragments, but it is not a complete sequence. There are some positional regions in the middle, which are called GAPs.

The GAP region after the genome is spliced includes: centromeres and telomeres.

A candidate variant site is classified as LCs when it has one of the following four characteristics:

(1) When the non-reference base occurrence frequency P _LCs is between the second frequency threshold T _HCs set by the user and the third frequency threshold T _LCs preset by the user, that is, T _LCs < P _LCs < T _HCs , T _LCs defaults to 0.2.

(2) The location of the current candidate variation site is a candidate multi-allelic site.

(3) Other candidate variant sites are included in a window extending w-bp upstream and downstream (default w = 100bp) centered on the position of the current candidate variant site on the genome of the sample being tested.

(4) The current location does not belong to the low complexity area.

Otherwise, when a candidate variant site does not meet the judgment criteria of HCs and LCs, it is classified as LCCs.

Step 2: construct a local haplotype sequence according to the type of the candidate variant site, and integrate the constructed local haplotype sequence among samples to obtain an integrated haplotype sequence, including the following steps:

Step 21: Divide each interval block of the divided reference genome into windows of fixed length, traverse all candidate mutation sites in the interval block, and assign all candidate mutation sites in the interval block to the optimal window:

Since the reference genome sequence block length is large and contains genome structural variation, some regional variations are relatively complex and difficult to construct haplotype sequences. Therefore, the present invention divides the interval block into multiple windows of fixed length W=300bp, with an overlap length of 150bp between adjacent windows, and constructs local haplotype sequences within the W window for local genome variation detection.

Traverse all candidate mutation sites in the interval block and assign them to an optimal window, that is, the candidate site falls in the center of the window as much as possible. Specifically, when the offset of the candidate mutation site compared to the starting position of the window is When , the current window is considered to be the optimal window.

Step 22: constructing a haplotype sequence of the sample being tested in each window according to the type of candidate variant site assigned to the optimal window;

Step 221: for each window, determine whether each window contains a candidate variation site. If it does not contain a candidate variation site, the reference genome sequence in the window is the haplotype sequence of the current window; if it contains a candidate variation site, determine whether the window contains a low-confidence candidate variation site. If the window contains a low-confidence candidate variation site, execute step 222 to obtain the constructed haplotype sequence of the sample being tested. If the window contains only a high-confidence candidate variation site, execute step 223 to obtain the constructed haplotype sequence of the sample being tested.

Step 222: When the window contains a low-confidence candidate variant site, a deBruijn graph-based splicing method is used in the window to generate high-quality consensus sequences (contigs) of the sequencing sequences around the candidate variant site, which are possible haplotype sequences:

First, all sequencing fragments within the window are extracted, and sequencing fragments with an overlap with the window of less than 10% are filtered out.

Set the starting k-mer to 41bp, construct a local deBruijn graph, identify all possible consensus sequences from the graph using a depth-first search algorithm, sort the consensus sequences based on the support of the sequencing fragments, and select the contigs with the highest support (default 2);

Add all contigs generated when k-mer = 41 to the sequenced fragment, set k-mer = 46, reconstruct the local deBruijn graph, and regenerate contigs according to the above steps;

Set the step size to 5bp, end k-mer = 75bp, and repeat the above operations until the last step generates contigs.

Step 223: When the window contains only high-confidence candidate sites, extract all sequencing fragments covering the current window from the alignment information set through the view function module in the alignment result analysis tool SAMtools, and obtain the sequencing fragment ID supporting each candidate variation site in the alignment information set (when the alignment information set contains the same base sequence as the candidate site, it is considered that the sequencing fragment actually supports the candidate variation site). Take the intersection of the sequencing fragment IDs supporting any two candidate variation sites. If the number of sequencing fragments contained in the intersection is not less than the user-defined number threshold (the default is 2), it is considered that the mutations corresponding to the two candidate variation sites appear in the same haplotype sequence, and the bases at the corresponding positions on the reference sequence are replaced with variant bases to construct the haplotype sequence to obtain the haplotype sequence; otherwise, it is considered to come from different haplotype sequences, and the haplotype sequences are constructed separately according to the candidate mutations and reference genome sequences contained in different haplotypes.

Step 2: Integrate the haplotype sequences of the tested samples in the window to obtain an integrated haplotype sequence set:

First, in a window, all haplotype sequences of the tested samples are traversed, and duplicates are removed according to the candidate variants contained therein, and finally a set of non-duplicate haplotype sequences of the samples in the window is generated. All windows are traversed, and the haplotype sequence set h = {h ₁ ,h ₂ ,…,h _n } of each window is generated in turn.

The deduplication is specifically as follows: if the haplotype sequences of two samples contain the same variation, only one haplotype sequence is retained. In a window, different sites may contain different types of variation.

Step 3: Detecting the genotype of the variable site of the genome of the population sample according to the integrated haplotype sequence, including the following steps:

Step 31: Perform non-gap alignment (non-gap allowed alignment) on the haploid sequence set obtained in step 2 and the sequencing fragments of the alignment information set, and obtain the sum of the cumulative mismatch base masses of the non-gap alignment of each sequencing fragment as S:

Step 3: Group the candidate variant sites and haplotype sequences in each window:

First, within the window, the candidate variation site information of each tested sample is traversed in turn. According to the tested candidate variation site information, all haplotype sequences containing the candidate variation site are selected from the haplotype sequence set, which are denoted as H. The remaining haplotypes do not contain the candidate variation and are denoted as R.

Then, the candidate variant sites contained in the tested samples are grouped: if the distance between two candidate variant sites is less than the distance threshold set by the user (default: 5bp), the two candidate variant sites are considered as a group.

Step 3-2: Use the candidate variant sites and the grouped haplotype sequences in each window to obtain the sum of the cumulative mismatch base masses S of the non-gap alignment of each sequencing fragment:

First, for each group of candidate variant sites obtained in step 3, all haplotype sequences containing the group of variants are selected from H and recorded as h, and all sequencing fragments within the window are extracted from the alignment information set.

Then, by mapping the exact match (Exactmatch) between all sequencing fragments and haplotype sequences in the local area and the reference sequence, the maximum exact match, that is, the longest exact match block, is obtained, and the starting alignment position is determined to be p. Starting from the alignment starting position, single base extension is performed to the left and right sides respectively. When a base mismatch is encountered during the extension process, the base quality of the mismatched base is accumulated (ASCII value -33). When the extension reaches both ends of the sequence or the sum of the cumulative mismatched base quality is greater than the base quality threshold set by the user (default: 300), the base extension is stopped. Similarly, non-repeating haplotype sequences are selected from R and recorded as r, and non-gap base extension is performed separately. The sum of the cumulative mismatched base quality of each sequencing fragment non-gap alignment is obtained and recorded as S.

Step 3.2: Use the sum of the cumulative mismatch base masses S of the non-gap alignment of each sequencing fragment obtained in step 3.1 to obtain the genotype of the variant site:

Step 321: Use the sum of the cumulative mismatch base masses S of each sequencing fragment non-gap alignment to obtain the optimal alignment probability set for each sequencing fragment after all sequencing fragments are aligned with the haplotype sequences in h and r respectively.

It is obtained in the following way: let h' = {h ₁ , h ₂ , ..., _hn' }, representing the list of haplotype sequences containing variants, the probability set of sequencing fragments aligning to n haplotypes is p = {p ₁ , p ₂ , ..., p _n' }, and the p _max with the largest probability is selected as the probability of the current sequencing fragment aligning to the haplotype containing the variant, and the corresponding haplotype is h _max .

Where k∈[1,m] is the probability index in the optimal alignment probability set for each sequencing fragment; m is the total number of probabilities obtained, and n’ is the number of haplotypes in the window;

p _j is obtained by the following formula:

Among them, j∈[1,n′] is the label of the probability of the sequenced fragment being aligned to n′ haplotypes;

Step 3: Calculate the likelihood of three possible genotypes (0/0, 0/1, 1/1) based on P and Pr, and take the genotype with the largest likelihood function as the mutated genotype in the current window. The formula is as follows:

If the genotype is 0/1 or 1/1, the optimal haplotype in h is inferred through each probability in P, and each variant genotype contained in the optimal haplotype is set to 0/1 or 1/1.

Step 4: Traverse the genomes of all tested samples in the population and repeat steps 1 to 3 to obtain all variant genotypes in the population:

Traverse the candidate mutation sites in all windows. For a certain mutation site, if all samples have no mutation signal at this site, or have mutation signal but the genotype of the detected mutation is 0/0, it is considered that all individuals in the current population do not have this mutation and will not be output; if there is a sample with a genotype of 0/1 or 1/1 at this site, traverse all samples in turn and output the genotype of each tested sample at this site. If a site is multi-allelic, split it into multiple double alleles and output them in turn.

Claims (8)

1. A population joint variation detection method based on a pan genome is characterized by comprising the following specific processes:

step one, obtaining a reference genome and a tested sample genome, identifying candidate mutation sites on the tested sample genome by using the reference genome and the tested sample genome, and classifying the candidate mutation sites:

Firstly, dividing a reference genome into interval blocks, and comparing the reference genome after dividing the interval blocks with a detected sample genome to obtain a comparison information set; then stacking the comparison information according to the coordinate sequence of the comparison information in the reference genome to obtain a stacking result; then counting the stacking results, and identifying the category of the candidate mutation site on the genome of the tested sample by using the stacking results;

the categories of candidate mutation sites on the genome of the sample to be tested include: high confidence candidate mutation sites, low complex region candidate mutation sites;

the first step comprises the following steps:

step one, obtaining a reference genome, and dividing the reference genome in interval blocks to obtain a divided reference genome B=B ₁ ,B ₂ ,…,B _n ；

Wherein B is ₁ ,B ₂ ,…,B _n Each divided interval block is respectively represented, and n is the total number of the interval blocks;

the interval block division is divided according to the following format: chromosome: start position-end position;

no overlap exists between the interval blocks;

step one, obtaining a genome of a tested sample, comparing the genome of the tested sample with a divided reference genome to obtain a genome of the tested sample falling on B _i A comparison information set in the database;

step one, the step one is that the product obtained in the step one falls on the B _i Stacking all the comparison information in the genome according to the coordinate sequence of the comparison information in the reference genome to obtain a stacking result;

step four, counting the stacking results obtained in the step three in each interval block of the divided reference genome, recording sequence information of the stacking results in each interval block through dictionary variables, identifying candidate variation sites of the tested sample according to the stacking result sequence information recorded by the dictionary variables, and classifying the candidate variation sites;

step two, constructing a local haplotype sequence by utilizing the candidate mutation sites according to the categories of the candidate mutation sites, and integrating the constructed local haplotype sequence among samples to obtain an integrated haplotype sequence set;

detecting the genotype of the mutation site of the genome of the tested sample by utilizing the integrated haplotype sequence set:

firstly, non-gap comparison is carried out on an integrated monomer type sequence set and sequencing fragments of a comparison information set to obtain the sum of accumulated mismatched base masses of the non-gap comparison of each sequencing fragment as S, and then the genotype of a variation site is obtained by utilizing the sum S of the accumulated mismatched base masses of the non-gap comparison of each sequencing fragment;

The genotype of the mutation site of the genome of the tested sample is detected by utilizing the integrated haplotype sequence set in the step three, and the method comprises the following steps:

step three, non-gap comparison is carried out on the integrated haplotype sequence set obtained in the step two and the sequencing fragments of the comparison information set, and the sum S of accumulated mismatched base masses of the non-gap comparison of each sequencing fragment is obtained:

step three, grouping candidate mutation sites and haplotype sequences of each window one by one:

firstly, selecting all haplotype sequences containing candidate mutation sites from an integrated haplotype sequence set, and marking the haplotype sequences as H; all the other haplotype sequences do not contain candidate variation and are marked as R;

then, candidate mutation sites contained in the tested sample are grouped: if the distance between the two candidate mutation sites is less than the distance threshold, treating the two candidate mutation sites as a group;

step three, two, the sum S of accumulated mismatched base masses of non-gap comparison of each sequencing fragment is obtained by utilizing the candidate mutation sites after each window grouping and the haplotype sequences after the grouping:

firstly, selecting all haplotype sequences containing the group of mutation from H as H for each group of candidate mutation sites obtained in the third step, and extracting all sequencing fragments in a window from a comparison information set;

Then, matching all sequencing fragments with a haplotype sequence, matching all sequencing fragments with a reference sequence to obtain a matching result, obtaining a complete matching block with the longest length from the matching result, and determining an initial comparison position as p;

then, starting from the initial comparison position, respectively performing single base extension to the left side and the right side, and accumulating the base quality of the mismatched base when base mismatch is encountered in the extension process; stopping base extension when the sum of the accumulated mismatched base masses or the base extension ends of the sequence is greater than a base mass threshold;

finally, selecting a non-repeated monomer sequence from R, marking the monomer sequence as R, respectively performing base extension of the non-gap, and marking the sum of accumulated mismatched base masses of the non-gap comparison of each sequencing fragment as S;

step three, obtaining the genotype of the mutation site by utilizing the sum S of the accumulated mismatched base quality of each sequencing fragment non-gap comparison obtained in the step three;

traversing all tested sample genomes in the population, and repeatedly executing the first to third steps to obtain genotypes of all mutation sites in the population.

2. The method for detecting the joint mutation of the population based on the pan genome according to claim 1, wherein the method comprises the following steps: the first step comprises the following steps:

Firstly, filtering comparison information of un-comparison, sub-optimal comparison and split comparison according to the value of the comparison information flag field in the comparison information set; filtering out the comparison information with the median value smaller than 10 in the mapq domain to obtain a filtered comparison information set;

wherein, the value of the non-aligned information flag field is 4; the value of the comparison information flag field of the sub-optimal comparison is 256; the value of the flag field of the comparison information of the split ratio pair is 2048;

then, stacking the filtered comparison information according to the coordinate sequence of the comparison information in the reference genome by using a comparison result analysis tool to obtain a stacking result;

the stacking result includes: reference sequence name, genome coordinates, reference base, number of sequenced fragments covering current position, base sequence of current position, base quality ASCII code sequence of current position.

3. The method for detecting the joint mutation of the population based on the pan genome according to claim 2, wherein the method comprises the following steps: in the fourth step, the stacking result obtained in the third step is counted in each interval block of the divided reference genome, the sequence information of the stacking result in each interval block is recorded through dictionary variables, candidate mutation sites of the tested sample are identified according to the stacking result sequence information recorded by the dictionary variables, and the candidate mutation sites are classified, and the method comprises the following steps:

Step one, counting the stacking result obtained in step one in each interval block of the divided reference genome, and recording the sequence information of bases, insertions and deletions appearing at each genome position in each interval block and the frequency information of each base sequence appearance through dictionary variables:

firstly, base sequence information of bases, insertions and deletions appearing at each genome position in a section block is recorded through dictionary variables;

the dictionary variables include: pileup_subject;

the pileup_subject is used for obtaining the occurrence times of the base sequence, the inserted base sequence and the deleted base sequence;

then, the frequency of occurrence of each base is obtained from the number of occurrences of each base sequence recorded in the pileup_subject dictionary variable:

P _j1∈keys ＝C _j1 /T

wherein T is the total number of occurrences of all bases, keys= { A, C, G, T, I, D } is the collection of bases, j1 is a base in keys, C _j1 Indicating the number of occurrences of base j 1;

step one, four and two, obtaining candidate mutation sites by utilizing the frequency of occurrence of each base obtained in step one, four and one:

if at least one base of the non-reference sequence is present at a genomic position, the frequency P of occurrence is _j1∈keys ＞T _alt The current position is the candidate mutation site;

Wherein T is _alt Is a frequency threshold;

the candidate mutation site is represented by a 4-tuple comprising: genomic position, reference base, non-reference base, frequency of occurrence of non-reference base sequence;

and step one, step four and step three, classifying the candidate mutation sites obtained in the step one, step four and step two.

4. A method for detecting group joint variation based on pan-genome according to claim 3, wherein: the step one, four and three classifies the candidate mutation sites obtained in the step one, four and two, specifically:

candidate mutation sites are classified into the following three classes: high confidence candidate mutation sites HCs, low confidence candidate mutation sites LCs, low complexity region candidate mutation sites LCCs;

HCs when the candidate mutation site satisfies the following three characteristics:

(1) Frequency of occurrence P of certain non-reference base _HCs Greater than a second frequency threshold T _HCs I.e. P _HCs ＞T _HCs ；

(2) The current candidate mutation site is positioned on the genome of the tested sample as a center, and other candidate mutation sites are not contained in windows extending w-bp at the upstream and downstream;

(3) The location of the current candidate mutation site on the genome of the sample to be tested does not belong to a low complexity region, comprising: GAP regions after genome splicing constitute heterochromatin domains; segment repetition; pseudo-autosomal regions of sex chromosomes; short tandem repeats;

LCs when the candidate mutation site satisfies one of four characteristics:

(1) Frequency of occurrence P of certain non-reference base _LCs At a second frequency threshold T set by the user _HCs Third frequency threshold T preset by user _LCs When in between, i.e. T _LCs ＜P _LCs ＜T _HCs ；

(2) The current candidate mutation site is located at a candidate multiallelic mutation site;

(3) The current candidate mutation sites are positioned on the genome of the tested sample as the center, and other candidate mutation sites are contained in windows extending w-bp at the upstream and downstream respectively;

(4) The current location does not belong to a low complexity region;

when none of the candidate mutation sites satisfies the judgment conditions of HCs and LCs, the candidate mutation sites are classified into LCCs.

5. The method for detecting the joint mutation of a population based on the pan genome according to claim 4, wherein the method comprises the following steps: in the second step, a local haplotype sequence is constructed by utilizing the candidate mutation sites according to the categories of the candidate mutation sites, and the constructed local haplotype sequence is integrated among samples to obtain an integrated haplotype sequence set, and the method comprises the following steps:

dividing each interval block of the divided reference genome into windows with fixed length, traversing all candidate mutation sites in the interval block, and distributing all candidate mutation sites in the interval block into an optimal window:

Firstly, dividing each interval block of a divided reference genome into a plurality of windows with the length W, and setting the overlapping length between adjacent windows;

then, traversing all candidate mutation sites in each interval block of the divided reference genome, and distributing the candidate mutation sites into an optimal window:

when the offset of the candidate mutation site compared with the window starting position is inWhen the current window is the optimal window of the current candidate mutation site;

step two, constructing a tested sample haplotype sequence in each window according to the candidate mutation site category distributed to the optimal window;

step two, integrating the haplotype sequences of the tested samples in the window to obtain an integrated haplotype sequence set:

firstly, traversing all tested sample haplotype sequences in a window, de-duplicating according to the included variation, and finally generating a non-repeated sample haplotype sequence set in the window;

then, traversing all windows, and sequentially generating a haplotype sequence set h= { h of each window ₁ ,h ₂ ,…,h _n }。

6. The method for detecting the joint mutation of a population based on the pan genome according to claim 5, wherein the method comprises the following steps: in the second step, according to the candidate mutation site category allocated to the optimal window, a detected sample haplotype sequence is constructed in each window, and the method comprises the following steps:

Step two, judging whether each window contains a candidate mutation site or not according to each window, and if the window does not contain the candidate mutation site, the reference genome sequence is the haplotype sequence of the current window; judging whether the window contains low-confidence candidate mutation sites or not if the window contains the low-confidence candidate mutation sites, executing the step two-two to obtain the constructed tested sample haplotype sequence if the window contains only the high-confidence candidate mutation sites, and executing the step two-two to three to obtain the constructed tested sample haplotype sequence if the window contains the low-confidence candidate mutation sites;

step two by two, a splicing method based on a de Bruijn graph is used in a window to obtain a haplotype sequence of a tested sample:

extracting all sequencing fragments in the window, and filtering out sequencing fragments with the overlapping degree of less than 10% with the window;

setting initial k-mer=41 bp, constructing a local deBruijn graph by using the rest sequencing fragments, identifying all possible consistent sequences from the local deBruijn graph by a depth-first search algorithm, sequencing the support degree of all the possible consistent sequences according to the sequencing fragments, and selecting a plurality of contigs with highest support degree;

Adding all the connigs generated at k-mer=41 to the sequencing fragment, setting k-mer=46, reconstructing the local deBruijn map, and reconstructing the connigs;

setting the step length to be 5bp, ending k-mer=75bp, and sequentially repeating the operation until the contigs is generated in the last step, namely the haplotype sequence of the tested sample;

step two and three, extracting all sequencing fragments covering the current window from the comparison information set by using a comparison result analysis tool, and obtaining sequencing fragment IDs supporting each candidate variation site in the comparison information set; taking intersection of sequencing fragment IDs supporting any two candidate mutation sites, if the number of sequencing fragments contained in the intersection is not less than a number threshold, considering that the mutation corresponding to the two candidate mutation sites appears on the same haplotype sequence, and replacing bases at corresponding positions on a reference sequence with mutation base construction haplotype sequences to obtain the haplotype sequence; if two candidate mutation sites are not supported by any identical sequencing fragment, they are considered to be from different haplotype sequences, and the haplotype sequences are constructed separately from the candidate mutation and reference genomic sequences contained on different haplotypes.

7. The method for detecting the joint mutation of the population based on the pan genome according to claim 6, wherein the method comprises the following steps: the genotype of the mutation site is obtained by utilizing the sum S of the accumulated mismatched base masses of each sequencing fragment non-gap alignment obtained in the step three, and the method comprises the following steps:

step III, obtaining the optimal comparison probability set of each sequencing fragment after comparing all sequencing fragments with haplotype sequences in h and r respectively by utilizing the sum S of accumulated mismatched base masses of non-gap comparison of each sequencing fragment obtained in step III

Obtained by: let a haplotype sequence list h' = { h containing variation ₁ ,h ₂ ,…,h _n’ Probability set of sequenced fragments aligned to haplotype sequences is p= { p ₁ ,p ₂ ,…,p _n ' selecting the highest probability p _max As the probability of the current sequencing fragment being aligned to the variant-containing haplotype sequence, the corresponding haplotype sequence is h' _max ；

Wherein k epsilon [1, m ] is the index of the probability in the optimal alignment probability set of each sequencing fragment; m is the total number of probabilities obtained and n' is the number of haplotype sequences that contain variation within the window;

p _j obtained by the following formula:

wherein j.epsilon.1, n' ] is the index of the probability of the sequenced fragment being aligned to n haplotype sequences;

Step III, two calculating likelihood probabilities of different genotypes by using P and Pr obtained in the third step and the first step, the genotype with the largest likelihood function as the variation in the current window is represented by the following formula:

8. the method for detecting the joint mutation of a population based on the pan genome according to claim 7, wherein the method comprises the following steps: traversing all tested sample genomes in the group in the fourth step, and repeatedly executing the first to third steps to obtain genotypes of all variation sites in the group, wherein the genotypes are specifically as follows:

traversing candidate mutation sites in all windows, aiming at a mutation site, if all samples have no mutation signal at the site or have mutation signals but the genotype for detecting mutation is 0/0, considering that all individuals in the current population do not have the mutation, and outputting the current population; if the genotype of the sample at the locus is 0/1 or 1/1, traversing all samples in turn, outputting the genotype of each tested sample at the locus, and if the locus is multiallelic, splitting into a plurality of bi-alleles and outputting the bi-alleles in turn.