CN106715711B - Method for determining probe sequence and method for detecting genome structure variation - Google Patents

Abstract

The present invention provides a method for determining a probe sequence based on a reference sequence and a method for detecting genomic structural variation. Wherein, the method for determining a probe sequence based on a reference sequence includes: constructing a first candidate probe set based on multiple discrete high-frequency SNP sites, the first candidate probe set is composed of multiple candidate probes, and the multiple candidate probes Each of the needles contains at least one discrete high-frequency SNP; the plurality of candidate probes in the first candidate probe set are aligned with the reference sequence to obtain an alignment result; based on the alignment result, the first candidate probe is aligned. The first screening of the needle set is performed to obtain a second candidate probe set; the reference sequence is divided into a plurality of windows with a predetermined length, and the plurality of candidate probes in the second candidate probe set are respectively allocated to the respective matching windows, so as to obtain a second candidate probe set. The respective position information of the plurality of candidate probes is determined; based on the position information and the allele frequencies of the discrete high-frequency SNPs, a second screening is performed on the second candidate probe set, so as to determine the probe sequence.

Description

Methods for determining probe sequences and methods for detecting genomic structural variation

priority information

none

technical field

The invention relates to the technical field of genomics and bioinformatics, in particular to a method for determining a probe sequence and a method for detecting genome structural variation.

Background technique

DNA copy number variation (CNV) and loss of heterozygosity (LOH) are different types of genomic variation. CNV is a common genomic structural variation, with fragments ranging from 1 kb to several Mb, mainly manifested as deletions and duplications at the submicroscopic level. LOH refers to the deletion of a gene on one chromosome of a pair of chromosomes, but it still exists on the paired chromosome, which shows that only homozygous SNPs appear in a long region of DNA. When there is no copy number change in LOH, that is, only two copies are inherited from one parent, it is called uniparental disomy (UPD). CNV, LOH, and UPD are associated with many common genetic disorders, cancer and other complex diseases. Establishing an accurate, comprehensive, efficient, rapid, simple and economical method for the detection of CNV, LOH and UPD is of great value for studying chromosomal mutation events, clarifying the etiology of related diseases and taking corresponding treatment plans.

At present, there are some inspection technologies, such as PCR technology, including real-time fluorescence quantitative PCR technology and Multiplex Ligation-dependent Probe Amplification (MLPA), real-time fluorescence PCR technology analyzes one or several targets at a time, MLPA It can analyze more than 40 sequences at a time, with high sensitivity, and the detection range is limited by the chromosomes and regions targeted by the probe; FISH technology, FISH technology is generally used to detect a few specific chromosomes, but cannot detect unknown regions; chip-based technology, Including chip-based comparative genomic hybridization (array-basedComparative Genomic Hybridization, aCGH) and SNP chip-based technology (SNP-array), aCGH can detect genome-wide CNV, can not detect polyploidy, the loss of small fragments and sequencing technology, based on whole genome sequencing (WGS) to detect genome-wide structural variation and target region-based sequencing to detect variation in target regions, there are four main methods to analyze CNV, including: paired Paired-end mapping, read-depth analysis, split-read strategies and sequence assembly comparisons.

With the development of sequencing technology, it is necessary to study the means to discover genomic structural abnormalities based on sequencing results, especially local region sequencing results, including the discovery of chromosomal aneuploidy, CNV, insertion-deletion (indel), LOH, UPD and the means of SNP.

SUMMARY OF THE INVENTION

One aspect of the present invention provides a method for determining a probe sequence based on a reference sequence, comprising the steps of: constructing a first candidate probe set based on a plurality of discrete high-frequency SNP sites, and the first candidate probe set consists of a plurality of candidate probes Probe composition, and each of the plurality of candidate probes contains at least one discrete high-frequency SNP; the plurality of candidate probes in the first candidate probe set are aligned with the reference sequence, so as to obtain an alignment result; based on Compare the results, perform a first screening on the first candidate probe set to obtain a second candidate probe set; divide the reference sequence into a plurality of windows with a predetermined length, and separate the plurality of candidate probes in the second candidate probe set respectively. needles are assigned to their respective matching windows to determine the respective position information of the plurality of candidate probes; based on said position information and the allele frequencies of the discrete high frequency SNPs, a second screening of the second candidate probe set is performed so as to The probe sequence is determined. Among them, the discrete high-frequency SNP site has an allele frequency greater than 10%, preferably not greater than 90%, and the physical distance from any other discrete high-frequency SNP site on the reference genome is not less than the length of the candidate probe. The probe length is 50-250mer.

The probe obtained by the method for determining the probe sequence of the present invention is used to hybridize the captured genome to obtain a plurality of local regions of the genome, and the captured local regions can represent the whole genome and reflect the variation information of the whole genome, and are used for discovering the whole genome. The occurrence of gene-wide structural variation.

Another aspect of the present invention provides a method for detecting genomic structural variation, suitable for detecting chromosomal aneuploidy, copy number variation and indels, comprising the steps of: sequencing the target sample genomic nucleic acid to obtain a genomic sequencing result , the genome sequencing result is composed of a plurality of reading segments, optionally, the sequencing includes using a probe to screen, wherein the probe is a method for determining a probe sequence based on a reference sequence provided by an aspect of the present invention. acquired. The genome sequencing results can be obtained by extracting genomic DNA, library construction and on-machine sequencing according to the existing high-throughput platform instruction manual; the genome sequencing results can also be obtained by capturing and sequencing the genome of the target sample with probes. Obtained by the method for determining a probe sequence based on a reference sequence provided in one aspect of the present invention; dividing the reference genome into m regions, and calculating the coverage depth TD _i of the target sample genome region i by using the reads that fall into the region i in the genome sequencing result , where m and i are natural numbers, 1≤i≤m, 10<m; based on the difference between the coverage depth of the target sample genome region i and the coverage depth of the k reference samples, the structural variation of the target sample region i is determined occurs, where k is a natural number, k≥2, and the method for obtaining the coverage depth of the region i of each reference sample can refer to the method for obtaining the coverage depth of the target sample region i. By merging adjacent regions with structural variation, it is further detected whether large structural variation occurs in the merged region, or whether the structural variation occurring in region i spans several regions.

Yet another aspect of the present invention provides a method suitable for detecting another kind of genome structural variation—loss of heterozygosity, comprising the following steps: obtaining a genome sequencing result of a target sample, optionally, the genome sequencing result is obtained by The probe is obtained by capturing the genome of the target sample and performing sequencing, and the probe is obtained according to the method for determining the probe sequence based on the reference sequence provided in one aspect of the present invention; the genome is divided into m' regions, and based on the genome sequencing results, the The reads in the region i and the population region i data, obtain the SNP set shared by the target sample genomic region i and the population region i, calculate the heterozygosity of each SNP site in the target sample and the population shared SNP set respectively, and obtain the target The heterozygosity set U i of the sample genome region _i , and the heterozygosity set U _0i of the population region i, compare the target sample U _i and the population U _0i to determine whether the loss of heterozygosity in the target sample region i occurs; The allele frequency of each SNP is greater than 0.1. The fragment of a SNP site in the common SNP set is bounded by the upstream and downstream two SNPs adjacent to the SNP. m' and i are natural numbers, and m'≥ i≥1, m'≥6. How many samples to draw can truly reflect the population can be determined according to the accuracy required for detection, statistical methods, and sample data distribution. The population data consists of multiple sample data of the same species, which can be obtained through whole genome sequencing, or by obtaining target samples based on data, or obtained from publicly available databases or websites, such as 1000 Genomes.

Yet another aspect of the present invention provides a computer-readable storage medium for storing a program for computer execution. Those of ordinary skill in the art can understand that when the program is executed, the above-mentioned detection of genomic structural variation can be accomplished by instructing relevant hardware. All or part of the steps of various methods. The so-called storage medium may include: read-only memory, random access memory, magnetic disk or optical disk, and the like.

According to a final aspect of the present invention, there is provided an apparatus for detecting genome structure variation, comprising: a data input unit for inputting data; a data output unit for outputting data; a storage unit for storing data, including an executable program; The processor is connected to the data input unit, the data output unit and the storage unit, and is used for executing the executable program stored in the storage unit. .

Using the probes obtained by the method for determining the probe sequence based on the reference sequence of the present invention, and using the probes or the solid-phase/liquid-phase chips containing these probes to capture and sequence the target region, the low sequencing cost can be achieved in the whole genome. Detect structural variation, including detection of CNV, LOH and UPD covering 23 pairs of human chromosomes, and the detection resolution can be adjusted according to needs by adjusting the average spacing distribution of probes, that is, increasing/decreasing SNP sites. The target region capture sequencing combined with the bioinformatics analysis method of the present invention realizes high-resolution, high-accuracy, high-throughput, and low-cost detection of CNV, LOH and UPD in the whole genome, and the genome structure variation of the present invention The detection method is also suitable for the detection of chromosomal aneuploidy variation, SNP and Indel, and for the detection of structural variation based on whole-genome sequencing data.

Additional aspects and advantages of the present invention will be set forth, in part, from the following description, and in part will be apparent from the following description, or may be learned by practice of the invention.

Description of drawings

The above and/or additional aspects and advantages of the present invention will be apparent and readily understood from the following description of embodiments in conjunction with the accompanying drawings, wherein:

Fig. 1 is a schematic diagram of the characteristics of the SeTR probe in one embodiment of the present invention on the whole genome, (A) the length distribution map of the SeTR probe sequence; (B) the physical distance distribution map of the two probes in the SETR probe .

Fig. 2 is a graph of the test results of the SeTR probe in one embodiment of the present invention, (A) a coverage depth distribution graph of the target region (B) a read distribution graph supporting ref base type and non-ref base type.

FIG. 3 is a schematic diagram of the detection flow of CNV, LOH and UPD in one embodiment of the present invention.

FIG. 4 is a graph of the R _i reference line in one embodiment of the present invention.

5 is a schematic diagram of the detected genomic structure variation of a sample (GM50275) according to an embodiment of the present invention, the circle is from the outside to the inside, followed by I) chromosome information, II) change of _ri value (wavy line) ; III) P value change corresponding to R _het , IV) R _het value change (dots).

Detailed description of the invention

Embodiments of the present invention are described in detail below. The embodiments described below with reference to the accompanying drawings are exemplary, only used to explain the present invention, and should not be construed as a limitation of the present invention.

It should be noted that the terms "first" and "second" are only used for descriptive purposes, and cannot be understood as indicating or implying relative importance or implying the number of indicated technical features. Thus, a feature defined as "first" or "second" may expressly or implicitly include one or more of that feature. Further, in the description of the present invention, unless otherwise specified, "plurality" means two or more.

According to one embodiment of the present invention, there is provided a method for determining a probe sequence based on a reference sequence, comprising the following steps:

Step 1: Construct the first candidate probe set

A first candidate probe set is constructed by using discrete high-frequency SNP sites distributed in the genome. Each candidate probe in the first candidate probe set contains at least one discrete high-frequency SNP site, and the discrete high-frequency SNP sites are alleles The gene frequency is greater than 10%, and the physical distance from any other discrete high-frequency SNP site on the reference genome is not less than the length of the candidate probe, and the length of the candidate probe is 50-250mer.

In a specific embodiment of the present invention, the discrete high-frequency SNPs are obtained from the 1000 Genomes data, and the discrete high-frequency SNPs with allele frequencies less than 90% can also be further selected from other published genome data or obtained. Point, determine the candidate probe length is 100mer.

In a specific embodiment of the present invention, each candidate probe contains a discrete high-frequency SNP site, and the discrete high-frequency SNP site is located in the middle of the candidate probe. In this way, each candidate probe contains only one high-frequency SNP site, and adjacent candidate probes may or may not overlap. The "middle section" here is relative to the "front section" and "rear section", which can be understood conventionally. For example, for a sequence, the upper and lower 1/3 of the sequence are designated as the "front section" and the "rear section", respectively. 1/3 of is the "middle segment"; further, the discrete high-frequency SNP sites are located at the midpoint of the candidate probe, where the "midpoint" position, for example, a sequence contains 2n+1 nucleotides, the middle The point is the position of the n+1th nucleotide, and when a sequence contains 2n nucleotides, the midpoint of the sequence is the position of the nth or n+1th nucleotide, which can enhance the probe to the target Capture efficiency of discrete high-frequency SNP sites.

In a specific embodiment of the present invention, the first candidate probe set is pre-screened based on the GC content and/or single base repeats of the candidate probe sequences in the first candidate probe set, and the first candidate probes are retained Candidate probes with a GC content of 35%-65% and/or a single base gravity of less than 7 were pooled. The single-base repetition degree refers to the number of consecutive occurrences of a base type in a sequence. For example, in TGAAAAAAAGC, where A appears 8 times in a row, the A base repetition degree of the sequence is 8. The high or low GC content and high heterozygosity of the sequence can easily affect the PCR or hybridization capture process of the sequence, resulting in GC bias (GC bias), etc., which reduces the capture specificity. The needle set will not hybridize to these sequences, thus avoiding the effect of GC bias or low specificity capture on the results.

Step 2: Align the first candidate probe set with the reference sequence in order to obtain the alignment result

The first candidate probe set is aligned with the reference sequence to obtain an alignment result, and the position information of the first candidate probe set on the reference sequence is obtained. The reference sequence used is a known sequence, which can be any reference template in the biological category to which the target sample is obtained in advance. For example, if the target sample is human, the reference sequence can be selected from HG18 or HG19 provided by the National Center for Biotechnology Information (NCBI). Selecting a closer reference sequence based on factors such as gender, race, and region of the target sample is conducive to obtaining more targeted probe sequences.

Step 3: Perform a first screening on the first candidate probe set to obtain a second candidate probe set

In a specific embodiment of the present invention, the candidate probes retained after the first screening need to meet either of the following two conditions: 1) The candidate probes in the first candidate probe set that are aligned to the unique position of the reference genome 2) The first candidate probe set is aligned to multiple positions in the reference sequence, and the mismatch ratio with at least two positions in the multiple positions of the reference sequence is less than 10%; for example, the length of the candidate probe is 100mer, 10 The base mismatch means that the mismatch ratio is 10%, and the mismatch ratio is low. When used for hybridization, it can be close to complete complementary pairing with the target region, with good capture effect and high specificity.

Step 4: Divide the reference sequence into multiple windows, and assign the second candidate probe set to each matching window

The reference sequence is divided into multiple windows with a predetermined length, and the multiple candidate probes in the second candidate probe set are allocated to the matching windows by alignment, and the position information of each candidate probe on the respective window is obtained.

The lengths of multiple windows of predetermined lengths may be consistent or inconsistent, and may or may not overlap. In a specific embodiment of the present invention, the reference sequence is a reference genome, and the reference genome is divided into multiple windows of consistent length, and the window length is 10Kb, and two adjacent windows are connected but do not overlap.

Step 5: Based on the location information and the allele frequencies of the discrete high-frequency SNPs, the second candidate probe Sets for second screening to determine probe sequences

In a specific embodiment of the present invention, performing the second screening includes two steps: (a) if there are multiple candidate probes located in the same window, determine the candidate probe with the highest allele frequency of the discrete high-frequency SNP (b) If there is only one candidate probe with the highest allele frequency of a discrete high-frequency SNP, select the candidate probe with the highest allele frequency of the discrete high-frequency SNP as the probe. If the candidate probe with the highest allele frequency of the high frequency SNP is selected, the candidate probe closest to the center of the window among the candidate probes with the highest allele frequency of the multiple discrete high frequency SNPs is selected as the probe. The distance of the candidate probe from the center of the window may be the distance between the midpoint of the candidate probe and the center of the window. The target position should be located in the center of the probe sequence as much as possible, which is beneficial to improve the capture efficiency.

In a specific embodiment of the present invention, after the second screening is performed on the second candidate probe set, when two adjacent candidate probes in the second candidate probe set that fall into two adjacent windows respectively are in the reference genome When the distance on the reference genome is greater than the length of any one of the two adjacent windows, optionally, further add the short tandem repeats or a part of the short tandem repeats between the two adjacent candidate probes on the reference genome. A probe sequence is formed together in the second candidate probe set after the second screening. In this way, when the probe sequences obtained by these designs are used to capture the whole genome, the spacing of the captured regions can be relatively uniformly distributed, and the combination of captured regions can better comprehensively reflect the entire genome information.

According to another embodiment of the present invention, there is provided a method for detecting genomic structural variation, wherein the genomic structural variation includes at least one of chromosomal aneuploidy, copy number variation and indel, comprising the following steps:

(1) Sequencing the genomic nucleic acid of the target sample in order to obtain the genome sequencing result. The genome sequencing result consists of multiple reads. The genome sequencing result can be obtained by whole-gene sequencing, for example, by extracting genomic DNA. Instruction manuals for throughput platforms, such as library construction and on-board sequencing to obtain reads using Illumina Hiseq2000/2500, Roche 454, Life technologies Ion Torrent, single-molecule or nanopore sequencing platforms, etc.; The genome of the target sample is obtained by sequencing, and the probe can be designed and determined by the probe determination method provided in one aspect of the present invention, and then synthesized or prepared according to the existing method.

(2) Divide the reference genome into m regions, and calculate the coverage depth TD _i of the target sample genome region i by using the reads that fall into the region i in the reads in the sequencing result, where m and i are natural numbers, and i is the region number, 1≤i≤m, 10<m.

或者

其中，i表示区域的编号。读段落到基因组上位置可以通过序列比对确定，比对可使用各种比对软件，例如SOAP(Short Oligonucleotide Analysis Package)，bwa(Burrows-Wheeler Aligner)，samtools，GATK(Genome Analysis Toolkit)等。In a specific embodiment of the present invention, the calculation formula of the coverage depth of the area i is:

or

Among them, i represents the number of the area. The position of the reads on the genome can be determined by sequence alignment, and the alignment can be performed using various alignment software, such as SOAP (Short Oligonucleotide Analysis Package), bwa (Burrows-Wheeler Aligner), samtools, GATK (Genome Analysis Toolkit) and the like.

(3) Determine the occurrence of structural variation in the target sample region i based on the difference between the coverage depth of the target sample genome region i and the coverage depth of the k reference sample regions i, where k is a natural number and k≥2.

其中，TD_j表示n个连续区域中的第j个区域的覆盖深度，j为自然数，1≤j≤n；(b)在获得区域i的第一校正覆盖深度TD_ai后，进一步对TD_ai进行均一化获得

进而获得

在本发明一个具体实施方式中，对区域i的第一校正覆盖深度TD_ai进行均一化获得的

在本发明的一个具体实施方式中，在获得目标样本的R_i后进一步包括对R_i进行第二校正以获得r_i，

其中，R_ai为k个参照样本基因组区域i的覆盖深度系数的平均值，

y为自然数表示参照样本编号，R_i，y表示参照样本y基因组区域i的覆盖深度系数。In a specific embodiment of the present invention, the comparison of the degree of difference between the coverage depth of the genomic region i of the target sample and the coverage depth of the region i of the k reference samples is performed by comparing the coverage depth of the genomic region i of the target sample and the reference sample. The determination of the coverage depth coefficient R _i of the target sample genome region i includes the following steps: (a) performing a first correction on TD _i to obtain TD _ai , and the first correction is performed by adjusting 2n including the region i It is realized by performing linear regression on the coverage depth values of the consecutive regions, where n is a natural number, and 10<n≤m/2.

Among them, TD _j represents the coverage depth of the jth area in the n continuous areas, j is a natural number, 1≤j≤n; (b) After obtaining the first corrected coverage depth TD _ai of the area i, further measure the TD _ai to obtain a homogenization

to obtain

In a specific embodiment of the present invention, the first corrected coverage depth TD _ai of region i is obtained by normalizing

In a specific embodiment of the present invention, after obtaining the R _i of the target sample, it further includes performing a second correction on R _i to obtain r _i ,

where R _ai is the average of the coverage depth coefficients of the k reference sample genome regions i,

y is a natural number representing the reference sample number, R _{i, y} represents the coverage depth coefficient of the genomic region i of the reference sample y.

其中，R_ai为k个参照样本和一个目标样本的基因组区域i的覆盖深度系数的平均值，

y为自然数表示参照样本编号，R_i，y表示参照样本y基因组区域i的覆盖深度系数。In another specific embodiment of the present invention, after obtaining the R _i of the target sample, it further includes performing a second correction on R _i to obtain r _i ,

where _Rai is the average of the coverage depth coefficients of the genomic region i of k reference samples and one target sample,

y is a natural number representing the reference sample number, R _{i, y} represents the coverage depth coefficient of the genomic region i of the reference sample y.

In the process of calculating and processing the coverage depth coefficient R _i of the target sample genomic region i, the correction and normalization of the intermediate values can reduce the errors caused by fluctuations in experimental conditions and differences between samples, so that the final The r _i can truly reflect R _i and the fluctuation range around ₁ is smaller than R _i , and the r _i of multiple samples conform to the normal distribution; The value is normalized, which is equivalent to the process of calculating the mean value twice, that is, before the average coverage depth of n continuous areas including area i is intended to represent the coverage depth of area i, the coverage depth value of each area in the n areas is calculated. The calculation is expressed by the average coverage depth of n consecutive regions with this region as the first region, which is equivalent to correcting TD _i by using the coverage depth values of 2n continuous regions including the target region i, which can make the continuous region The depth of coverage remains stable. It should be noted that those skilled in the art can use other correction or averaging process to keep the coverage depth value of several adjacent areas stable, for example, to correct the target area with the average coverage depth of several areas separated from the target area. The coverage depth belongs to the concept of the present invention. The calculation process of the coverage depth coefficient of the reference sample genome region i can refer to the calculation process of the coverage depth coefficient of the target sample genome region i, and the reference sample data can be pre-calculated and processed for use, or can be performed synchronously with the calculation process of the target sample. get.

其中，

表示k个参照样本的r_i,y的平均值，r_i,y为参照样本y基因组区域i的经第二校正的覆盖深度系数，

为k个参照样本标准差，

基于目标样本基因组区域i的t_i值，获得显著水平P_i，当P_i<0.05，判定所述区域i发生结构变异；反之，则判定所述区域i不发生结构变异。在本发明的另一个具体实施方式中，基于目标样本基因组区域i的t_i值和预先确定的显著水平P_i0，获得t_i理论值t_i0，当t_i≥t_i0，判定所述区域i发生结构变异，反之，则判定所述区域i不发生结构变异，预先确定的P_i0≤0.05。根据t检验的t值表，预定P_i0后可查得对应的t_i0。In a specific embodiment of the present invention, the degree of difference between the coverage depth of the genomic region i of the target sample and the coverage depth of the region i of the k reference samples is judged by t-testing whether the difference in the coverage depth coefficients of the two is significant. realized. In a specific embodiment of the present invention, the calculation formula of the t-test statistic of the genomic region i of the target sample is:

in,

represents the average value of ri _,y of k reference samples, ri _,y is the second corrected coverage depth coefficient of genomic region i of reference sample y,

is the standard deviation of k reference samples,

Based on the t _i value of the target sample genome region i, a significant level P _i is obtained. When P _i <0.05, it is determined that the region i has structural variation; otherwise, it is determined that the region i has no structural variation. In another specific embodiment of the present invention, based on the t _i value of the target sample genome region i and the predetermined significance level P _i0 , the theoretical value t _i0 of t _i is obtained, and when t _i ≥t _i0 , the region i is determined Structural variation occurs, otherwise, it is determined that the region i does not undergo structural variation, and the predetermined P _i0 ≤ 0.05. According to the t value table of the t test, the corresponding t _{i0 can be found after the predetermined P i0 .}

In one embodiment of the present invention, in order to detect larger CNVs or indels, after step (3), merge W consecutive regions in the same direction to obtain a first-level merged region, and merge two first-level merges Region When two first-level merged regions are in the same direction and the span between them does not exceed L regions, the second-level merged region is obtained, and the structural variation of the second-level merged region is detected; among them, the same-direction region refers to the coverage depth of the region t In the region where the statistics are all greater than 0 or less than 0, W and L are both natural numbers, W≥2, L-W≤1. To further detect larger structural variations, it can be deduced by analogy, such as further merging qualified secondary merged regions, the merge conditions can be similar to two secondary merged regions in the same direction and the distance between them on the reference genome does not exceed L regions or L secondary merged regions.

In a specific embodiment of the present invention, detecting the structural variation of the secondary merged region is based on the degree of difference between the coverage depth of the secondary merged region of the target sample genome and the coverage depth of the corresponding regions on the genomes of multiple reference samples , to judge whether structural variation occurs in the secondary merged region, or to judge whether the structural variation occurring in region i spans W regions. For the acquisition of the coverage depth of the corresponding secondary merged region on the reference sample genome, the calculation of the t-statistic of the coverage depth of the secondary merged region on the target sample genome, and the process of structural variation judgment, please refer to the structural variation of the relatively small region i above. Computational judgment process.

According to yet another embodiment of the present invention, there is provided a method suitable for detecting loss of heterozygosity in genomic structural variation, comprising the following steps:

(1) Sequence the genomic nucleic acid of the target sample in order to obtain the genomic sequencing result. The genomic sequencing result is composed of multiple reads. The genomic sequencing result can be obtained by whole-gene sequencing, for example, by extracting genomic DNA. Instruction manuals for throughput platforms, such as library construction and on-board sequencing to obtain reads using Illumina Hiseq2000/2500, Roche 454, Life technologies Ion Torrent, single-molecule or nanopore sequencing platforms, etc.; The genome of the target sample is obtained by sequencing, and the probe can be designed and determined by the probe determination method provided in one aspect of the present invention, and then synthesized or prepared according to the existing method.

(2) Divide the reference genome into m' regions, and obtain the SNP sets shared by the target sample genome region i and the population region i based on the read information in the reference genome region i and the population region i data in the sequencing results, and calculate separately The heterozygosity of each SNP site in the shared SNP set of the target sample and the population is obtained, and the heterozygosity set U _i of the target sample genome region i and the heterozygosity set U _0i of the population region i are obtained, and the target sample U _i and the population are compared. U _0i to determine whether the loss of heterozygosity in the target sample region i occurs; wherein, the allele frequency of each SNP in the shared SNP set is greater than 0.1, and the fragment where a SNP site in the shared SNP set is located is a The upstream and downstream two SNPs adjacent to the SNP are boundary points, m' and i are natural numbers, m'≥i≥1, m'≥6.

In a specific embodiment of the present invention, the heterozygosity of the fragment where a SNP locus is located is represented by the sub-allele frequency coefficient of the SNP locus, and the sub-allele frequency coefficient of the SNP locus R _het = MAF/(1-MAF), MAF is the minor allele frequency of the high frequency SNP.

和U_0i的方差

是否有显著性差异，，若U_i和U_0i的方差差异显著，则判定所述目标样本区域i存在杂合性丢失，反之，则判定所述目标样本区域i没有存在杂合性丢失。In a specific embodiment of the present invention, comparing the target sample U _i with the population U _0i to determine whether the loss of heterozygosity in the target sample region i occurs, includes using an F test to determine the variance of U _i

and the variance of U _0i

Whether there is a significant difference, if the variance of U _i and U _0i is significantly different, it is determined that the target sample region i has loss of heterozygosity; otherwise, it is determined that the target sample region i has no loss of heterozygosity.

和群体U_i0的方差

计算获得两个互为倒数的统计量F_upper和F_under，利用互为所说的倒数的统计量获得显著水平p_F，比较p_F与预定显著水平p_F0的大小，p_F≤p_F0说明两方差差异显著，反之则差异不显著，F检验包含计算公式，In a specific embodiment of the present invention, the F test includes calculating the variances of U _i and U _i0 respectively, and using the variance of the obtained target sample U _i

and the variance of the population U _i0

Calculate and obtain two mutually reciprocal statistics F _upper and F _under , use the mutually said reciprocal statistics to obtain the significant level p _F , compare the size of p _F and the predetermined significant level p _F0 , p _F ≤ p _F0 description The difference between the two variances is significant, otherwise the difference is not significant. The F test includes the calculation formula,

p_F＝p_upper+(1-p_under)，其中，v为目标样本基因组区域i和群体区域i共有SNP集中SNP的编号，q为目标样本基因组区域i和群体区域i共有SNP集中SNP的个数，R_het,i,v为目标样本基因组区域i的共有SNP集中的第v个SNP的次等位基因频率系数，

为目标样本基因组区域i的共有SNP集中的q个SNP的次等位基因频率系数的平均值，R_het,i0,v群体样本基因组区域i的共有SNP集中的第v个SNP的次等位基因频率系数，

为群体样本基因组区域i的共有SNP集中的q个SNP的次等位基因频率系数的平均值，p_upper和p_under分别根据F_upper和F_under获得，p_F0≤0.05。p_F0可以取通常设置的值、或者根据所掌握的已知信息、对检测准确性的要求等调整设置。

p _F =p _upper +(1-p _under ), where v is the number of the SNP in the SNP set shared by the target sample genome region i and the population region i, q is the number of SNPs in the SNP set shared by the target sample genome region i and the population region i number, R _{het, i, v} is the sub-allele frequency coefficient of the vth SNP in the shared SNP set of the target sample genome region i,

is the average of the minor allele frequency coefficients of q SNPs in the shared SNP set of the target sample genome region i, R _het,i0,v is the minor allele of the vth SNP in the shared SNP set of the i population sample genome region i frequency coefficient,

is the average value of the sub-allele frequency coefficients of q SNPs in the common SNP set in the genomic region i of the population sample, p _upper and p _under are obtained according to F _upper and F _under , respectively, p _F0 ≤ 0.05. p _F0 can take the value that is usually set, or adjust the setting according to the known information and the requirements for detection accuracy.

In one embodiment of the present invention, in order to detect a larger LOH, after step (2), the W' loss of heterozygosity and continuous regions are combined to obtain a three-level combined region, and two three-level combined regions are combined. When the span between the two third-level merged regions does not exceed L' regions, obtain a fourth-level merged region, respectively obtain the heterozygosity set of the fourth-level merged region of the target sample and the heterozygosity set of the same region of the population, and compare Two sets of heterozygosity are used to determine whether loss of heterozygosity occurs in the four-level merged region of the target sample, where W' and L' are both natural numbers, W'≥2, and W'/2≥L'. In a specific embodiment of the present invention, W'≥4. To detect LOH in a larger region, it can be deduced by analogy, for example, further merging eligible fourth-level merged regions, the merging conditions can be similar that the distance between two fourth-level merged regions on the reference genome does not exceed L' area or L' tertiary merged areas.

According to another embodiment of the present invention, a method for detecting uniparental diploidy is provided. When there is a loss of heterozygosity in a genomic region of a target sample, the copy number of this region is calculated. When the copy number of this region is the same as that of the normal genome of the same species When the copy number of the above region is the same, it is determined that there is UPD in this genomic region of the target sample; whether there is LOH in the genomic region can be determined by the LOH detection method disclosed in one aspect of the present invention.

Those of ordinary skill in the art can understand that all or part of the steps of the various methods in the above-described embodiments can be completed by program instructions related to hardware, and the program can be stored in a computer-readable storage medium, and the storage medium can include: read-only memory, Random access memory, magnetic disk or CD, etc.

According to the last embodiment of the present invention, there is also provided an apparatus for detecting genome structure variation, comprising: a data input unit for inputting data; a data output unit for outputting data; a storage unit for storing data, including An executable program; a processor, data-connected to the above-mentioned data input unit, data output unit and storage unit, for executing the executable program stored in the storage unit, the execution of the program includes completing all of the various methods in the above-mentioned embodiments or part of the steps.

The operation results of the specific probe design method and the structural variation detection method according to the present invention will be described in detail below in conjunction with specific target individuals. The name definitions or specific parameter settings involved in the following process are selected as:

1. The designed probe is called Selected Target Region Primers (SeTR);

2. The following "coverage depth", "sequencing depth" and "depth" can be used interchangeably; the following "area" and "target area" can be used interchangeably;

3. Library construction and sequencing are operated according to the operation instructions for small fragment library construction and on-machine sequencing instructions provided by the Hiseq 2000 platform. The size of the library is 300bp-350bp, pair-end sequencing (pair-end sequencing), and the read length is 91bp (sequencing). Type is PE91+8+91);

4. The reference genome or reference sequence selected for alignment is the human reference genome (hg19, Build 37).

If the specific technique or condition is not indicated in the embodiment, according to the technique or condition described in the literature in this area (for example, with reference to J. Sambrook etc., "Molecular Cloning Experiment Guide" translated by Huang Peitang etc., 3rd edition, Science Press) or follow the product instructions. The reagents or instruments used without specifying the manufacturer are conventional products that can be obtained commercially, for example, can be purchased from Illumina.

Example 1: Chip Design, Preparation, and Testing

Usually, high (>60%) or low (<35%) GC content and high heterozygosity easily bring adverse effects on the DNA fragments during PCR or probe capture. In order to avoid this phenomenon, we designed The special probe is called SeTR. When designing the SeTR probe, the following principles are followed: a) The uniqueness and stability of the probe sequence are high, requiring low heterozygosity and moderate GC (35%～60%) content, b) contains discrete high-frequency SNP markers (SNPmarker), the allele frequency of each SNP (allele frequency, 0.9>AF>0.1) to better detect the LOH of the whole genome, c ), the final target area exhibits a relatively uniform distribution.

The SeTR probe design or the selection process of the target region is as follows:

1) Based on the 1000 Genes Database (ftp://ftp.ncbi.nih.gov/1000genomes/ftp/release), select the candidate SNP set with an allele frequency (AF) of 10% to 90%, Then, one SNP whose physical distance between two SNPs is less than 100pb is removed from the SNP set, thereby forming the SNP maker1 set.

2) Take each SNP in the SNP maker1 set as the midpoint, cut off 50 pb of the reference genome sequence upstream and downstream of it to form a theoretical probe sequence set of 100 bp, and then compare this probe sequence set back to the reference sequence. If the best alignment of a certain probe sequence has no mismatches, and the next best alignment has only less than 5% mismatches, then its corresponding SNP is retained, thus forming the SNP maker2 set.

3) Based on the SNP maker2 set, we select the SNPmaker that is uniform in the physical position of the reference genome as the final SNP maker set. In our study, we selected the SNP maker set with a physical distance of about 10Kbp.

4) If the distance between two adjacent SNPs is greater than 10Kbp in the final SNP maker set, the short tandem repeat (STR) between them is selected to fill evenly.

After designing the SeTR probe, we commissioned Roche to produce the SeTR LC chip. The SeTR liquid chip contains 278,800 probes with a total size of 41,795,106 bp, which covers an effective genome-wide (2.89G) 1.45% region. The average length of the SeTR probes reached 149 bp, and the average physical distance between adjacent probes was 10.6 kbp, as shown in Table 1 and Figure 1.

Table 1 Distribution of SeTR probes on each chromosome

Using 3 DNA samples that have passed the quality inspection, YH (Yanhuang sample, Chinese genomic DNA), HG00537 (a sample in the Thousand Genomes Project) and GM50275 (obtained from the human adult of Coriell Institute for Medical Research) fibroblast samples) to test the availability of the SeTR chip to ensure that this probe chip can be used for subsequent detection studies. All three samples were sequenced using SeTR capture library construction to obtain sequencing sequences (reads). First, after removing the reads that are contaminated by adapters and low quality, such as the average quality value lower than 20, the remaining reads are called clean reads (clean reads), and the clean reads are aligned to the reference sequence hg19 to get 98.13% to 99.29 of the reads were aligned to the reference genome, of which 67.43% to 67.87% were aligned to the target region. In addition, 99.73% to 99.95 of the target region was covered by at least one read, and more than 99% The area was covered at least 10 times, as shown in Table 2, and these performances were superior to the same type of exome capture chip, such as the exome liquid phase chip produced by Roche Nimblegen. In addition, the depth distribution of the target region, as shown in Figure 4A, is similar to the Poisson distribution, and Figure 4B shows the non-reference sequence base type (the non-reference sequence base type of most of the high heterozygous sites in the target region) The number of reads supported by -reference allele is almost the same as the number of reads supported by the base type of the reference sequence (reference allele), that is, the number of positive and negative reads obtained from a high heterozygous site during alignment is equivalent (the positive and negative reads come from two sources, respectively). Homologous chromosomes), these all show that this probe has no obvious haplotype (commonly referred to as the reference sequence base type, ie ref type) capture bias, and the capture uniformity of the target region is better.

Table 2 Comparison results of three samples

Example 2: Target Region Library Construction and Sequencing

1. Test materials, reagents and instruments

Samples: 15 target gDNA samples (human genomic DNA, the sample number is shown in Table 3 below, "GM"" begins with human fibroblasts), 24 reference DNA samples.

Main reagent instruments: PCR instrument, pipette, centrifuge, comfortable constant temperature mixer, DNA fragmentation instrument, vortex oscillator, magnetic stand, electrophoresis instrument, Hiseq2000 sequencer, Nanodrop UV spectrophotometer, etc. Reagents used Or the instrument does not indicate the manufacturer, it is a conventional product that can be obtained through the market.

Probe design and synthesis: obtained in Example 1, in the human genome, a target region of about 41M was selected, and the NimbleGen SeqCap EZ liquid-phase probe was customized from Roche. The probe set can capture the corresponding probes. The designed target area.

2. Library construction

1) Genomic DNA extraction

Using the QIAGEN DNA Extraction Kit (DNA Mini Kit) and following the kit instructions, about 3-5 μg of genomic DNA was extracted from the target sample for subsequent experiments. Run the extracted DNA (30-100ng) for electrophoresis to determine whether it is complete and the degree of degradation.

2) Genomic DNA fragmentation and purification

Genomic DNA was disrupted using a covaris E210 instrument (refer to the instructions for use of the instrument). The DNA was broken into 200-250bp. Use the QIAquick PCR Purification kit (250) kit and operate according to the kit instructions to purify the fragmented DNA fragments, and electrophoresis to detect whether the size of the main band meets the requirements, that is, whether the main band size is 200-250bp.

3) End repair, end adding A, adding adapter, pre-amplification

According to the library construction requirements, according to the steps of the double-end tag library construction instructions and the listed reagents, reaction conditions, etc., the above-mentioned fragmented and purified DNA fragments are end-repaired and purified; add a base A to the end-repaired purification The two ends of the DNA fragment are purified by adding A product to the end; the two ends of the end adding A product are connected with sequencing adapters, and the DNA fragments with adapters are purified by magnetic beads that can be complementary to the sequencing adapter. Prepare a PCR reaction system, amplify the DNA fragments with adapters, purify the PCR products with magnetic beads, and detect by electrophoresis whether the size of the main band of the amplified products is 300-350 bp; use a Nanodrop UV spectrophotometer to detect the amount of DNA, and the total amount should be greater than 1.0 μg.

4) SeTR probe hybridization and elution, amplification

Follow the instructions of the commercially available NimbleGen SeqCap EZ hybridization and elution kit, and purchase or configure the relevant reagents for hybridization and elution in the kit instructions. Prepare a 1.5mL centrifuge tube, add Cot-1 DNA, universal blocking sequence (Block), tag blocking sequence (index N Block) and the DNA sample after step 3). Then centrifuged for 1 min, concentrated and dried in vacuum at 60 °C, then added hybridization buffer, etc., and centrifuged by shaking, put it in a metal dry bath at 95 °C for denaturation for 10 minutes, and centrifuged at high speed after shaking. Add 4.5ul probe to the centrifuge tube and hybridize on PCR machine (47°C, 64-72hours). Elution is performed after hybridization is complete. Then carry out PCR according to the last amplification step of the library construction instructions, prepare a PCR reaction system as required, and eluate the DNA obtained by hybridization, polymerase, substrate, PCR reaction buffer, Flowcell primers (according to the sequencing chip flowcell of the sequencer). The primers with the fixed sequence design) are mixed evenly. The PCR program was pre-denaturation at 94 °C for 2 min, denaturation at 94 °C for 15 s, annealing at 58 °C for 30 s, extension at 72 °C for 30 s, and after 15 cycles of reaction, extension at 72 °C for 5 min. After the PCR is completed, the PCR product is taken out, centrifuged, purified by magnetic beads, and the target region library is obtained. Measure the concentration of the library with Nanodrop UV spectrophotometer and prepare it for on-board sequencing.

3. Hiseq2000 high-throughput sequencing

Qualified DNA libraries were sequenced according to the Hiseq2000 operating instructions. The data volume of each sample is about 4.5G, and the average sequencing depth reaches 100X, but because the efficiency of the capture chip is difficult to reach 100%, through analysis, the effective sequencing depth of the final target region is 30X~45X.

Example 3: Detection of CNV, LOH and UPD

The overall process is shown in Figure 3. After the sequencing is completed, the off-machine data is in the fastq file format. Then, the high-quality reads obtained after filtering were compared with the reference genome (Hg19, Build 37) using BWA software to obtain an alignment file in SAM format, and then the SAM alignment file format was converted into a binary BAM file using samtools software, and The results are compared and deduplicated and sorted. Then, the samtools software will be used again to convert the BAM format to the PILEUP format. For details, please refer to the section on bioinformatics analysis strategies.

1. Sequencing data filtering and comparison

First, perform simple data filtering on the sequencing data of the illumina Hiseq2000 in the above example, and remove the reads that are contaminated by the adapter, contain more than 5% of N, and whose average quality value is lower than Q20. Then use the bwa alignment software to align the filtered data to the human reference genome (hg19, Build 37), and output the sequence alignment result that is the alignment file in SAM (sequence alignment/map) format (SAM file for short), and then Use Samtools software to convert SAM files into binary BAM files, remove PCR duplicates and perform sorting, and use GATK software to compare the results for multiple alignment and re-correction.

2. Calculate the second corrected coverage depth coefficient ri of the target area and the heterozygosity R _het of the _fragment

The ri and fragment heterozygosity _Rhet values of each target region are calculated according to the information contained in the probe region files obtained after the above filtering and _alignment . The t-test was used to predict CNV according to the _ri value, and the F-test was used to predict LOH and UPD according to R _Het .

3. Analysis to detect CNV, LOH and UPD

1. CNV detection

1.1 Calculate the depth coefficient (R _i ) of each target region

Calculate the depth of the target area and express it with TD _i (such as formula 1). In order to maintain the stability of several consecutive target areas TD, the method of formula 2 is used to correct TD _i , that is, use the n' behind the i-th area The depth of the region is used to correct TD _i to obtain TD _ai , and then the TD _ai is normalized by formulas 3 and 4, and the depth coefficient R _i of each target region is obtained at this time.

Formula 1: TD _i =T _i base/T _i len,

Formula 2:

Formula 3:

T_ibase：比对到目标区域i的碱基数；T_ilen：目标区域i的长度。Formula 4:

T _i base: the number of bases aligned to the target region i; T _i len: the length of the target region i.

1.2 Use multiple reference sample (k=24) data to create a baseline, correct R _i to obtain r _i

的分布如图波动很大，而R_i波动相对小些，R_i通过基准线的校正后得到r_i，其波动更小，更敏感更易检测CNV的发生。理论上，认为不同的样品中，不发生CNV的情况下，在同一个目标区域内，R_i值理论上是符合泊松分布的，并且都围绕各自特有的值相对稳定的上下波动，为了保持各自特有值的稳定性，通过调查多个样品的同一区域的R_i值，采用R_i平均值(mean R_i)来代替这个各自特有的值，为每个目标区域构建一个各自特有的基准线(robust baseline)。基于每个目标区域的R_i是围绕着mean R_i值上下波动的假设，我们将R_i除以mean Ri转化成r_i，进而使得r_i围绕1上下波动的正态分布。Due to the fluctuation of each experimental condition and the differences between samples, there is also a certain fluctuation in the efficiency of each capture, which in turn causes fluctuations in R _i , which easily lead to CNV false signals. Therefore, it is very beneficial to create a unified baseline based on multiple sample fluctuations. Figure 4 is a good example of creating a baseline for this assay, the precursor R _i (preR _i )

The distribution of , as shown in the figure, fluctuates greatly, while the fluctuation of R _i is relatively small. After R _i is corrected by the baseline, _ri is obtained, and its fluctuation is smaller, which is more sensitive and easier to detect the occurrence of CNV. Theoretically, it is considered that in different samples, if CNV does not occur, in the same target area, the R _i value theoretically conforms to Poisson distribution, and fluctuates relatively stably around their unique values. The stability of each unique value, by investigating the R _i value of the same area of multiple samples, using the average value of R _i (mean R _i ) to replace this unique value, and constructing a unique baseline for each target area. (robust baseline). Based on the assumption that the R _i of each target area fluctuates around the mean R _i value, we convert R _i by the mean Ri into r _i , which in turn makes r _i a normal distribution that fluctuates around 1.

1.3 Detecting CNVs in the target area

Theoretically, the ri values of the same target area from multiple samples should conform to the normal distribution, so when investigating the target area _i of a certain sample, you can compare the _ri values of this area of multiple samples and use the T test, The calculation formula of the t statistic is as follows, to detect the copy number change of the target region i of this sample.

表示n₁个待测样本的r_i的平均数，

为n₂个参照样本的r_i的平均数，μ₁为理论上所有待测样本的r_i平均数，μ₂理论上全部参照样本r_i2平均数,S₁和S₂分别为待测样本和参照样本的标准差，df为自由度，df＝n₁+n₂-2。1 in the subscript of each parameter in the formula represents the target sample, 2 represents multiple reference samples,

represents the mean of ri of _{n 1} samples to be tested,

is the average number of ri of _{n 2} reference samples, μ ₁ is the average number of ri of all samples to be tested in theory, _{μ 2} is the average number of _ri of all reference samples in theory, S ₁ and S ₂ are the samples to be tested respectively and the standard deviation of the reference sample, df is the degree of freedom, df=n ₁ +n ₂ -2.

When the sample to be tested is 1, that is, n ₁ =1, the theoretical mean of the sample to be tested and the reference sample are the same, and the above formula is simplified as:

Through the above simplified formula, each target area corresponds to a t value of detectable CNV, and then the P value (confidence) is obtained. When the P value of a certain area is less than 0.05, this area is an area where CNV occurs.

1.4 Detection of large CNVs

Based on the p-value of the t-test of a single region, a pseudo-signal value is attached to each region to indicate whether it is considered by the next step of CNV region connection, and then along the chromosome, the target regions that may have the same CNV are connected into blocks, so as to determine The final size and copy number of the CNV.

The labeling rule for false signal values is that when the measured values of at least four consecutive target areas are in the same direction (t values are greater than or less than 0 at the same time), that is, when they deviate from the corresponding areas of the reference sample, if there are three areas whose P values are less than the first. Threshold (such as 0.05, commonly used significant level threshold), and the fourth does not exceed the second threshold (0.2, four times the first threshold), then the four regions are marked as deviating directions (for example, large areas are marked as +, partial The small mark is -), which is merged into a block; here, the number of continuous and same-direction regions and the first and second threshold values are adjustable. If the distance between one block and another block does not exceed the span of 5 areas, the two blocks will be merged into one large block, and so on, and finally a block will be obtained; referring to the method formula in the previous 1.3, this block The r value of is represented by the average value of r _i of all the areas it contains, and the r value of the block area of the sample to be tested and the reference sample is tested by t test, and the P value of the block is calculated. When the P<0.05 of the block, CNV occurred in this block, so the boundary and size of the block were determined, and the boundary and size of the large CNV were obtained.

Through the analysis of the target 15 samples, the CNV results we obtained were highly consistent with the known validation results (SNP-array results), and there were no false positives and false negatives, as shown in Table 3. Furthermore, we simulated 8 30X genome-wide data, including 5 normal samples and 3 CNV-containing samples, and compared the currently reported CNV predictions in the exome region by performing CNV detection and analysis on these 8 simulated data. The software CONTRA (Li J, Lupat R, et al, CONTRA: copynumber analysis for targeted resequencing, Bioinformatics. 2012 May 15; 28(10): 1307-13), the results show that the sensitivity and specificity of our method have reached 100%, and their respective copy numbers are also accurately detected, the detection accuracy of CNV can reach 500Kb and can be accurately located, while the sensitivity of CONTRA is 88.9%, the specificity is only 66.7%, and the copy number is not given. As shown in Table 4.

table 3

Table 4

2. LOH detection

2.1 Detection of heterozygous status in each region of the whole genome

In a certain region of the genome of the sample to be tested, find the SNP loci with gene frequency (AF) of 0.1 to 0.9 in the thousand-person data, and calculate the location of these SNP loci in the thousand-person neutralization and the sample to be tested according to the following formula The R _Het value of the area. When the region i in the sample to be tested is in an absolute heterozygous state, R _het =1; otherwise, when it is an absolute homozygous state, R _het =0.

R _het =MAF/(1-MAF), MAF (minor allele frequency) is the minor allele frequency.

和千人样本的同样区域的杂合度集的方差

以及待测样本该区域杂合度集Sm的p值。In the sample to be tested, any SNP site m in a certain region is taken as the starting point, and n consecutive SNP sites are taken backward as the heterozygosity set Sm in the region, that is, S _m ={R _het,sm ,R _het,s(m+1) ,...,R _het,s(m+n) }, in the same way, in the thousand-person database, take the SNP sites at the same position to form the heterozygosity set Pm, that is P _m ={R _het,phetm ,R _het,p(m+1) ,...,R _het,p(m+n) }, F tests whether the variances of the two heterozygosity sets are equal, specifically, respectively Calculate the variance of the heterozygosity set for this region of the sample to be tested

The variance of the heterozygosity set in the same area as the thousand-person sample

and the p-value of the heterozygosity set Sm in this region of the sample to be tested.

H ₀ :σ _s =σ _p

H _A :σ _s ≠σ _p

p=p _upper +(1-p _under )

When p≤0.01, we accept _HA and judge that the heterozygosity set Sm loses the heterozygosity in the population, that is, LOH occurs in the area where the set Sm is located.

2.2 Detection of large LOH

Combined with the results of 2.1, the method of step 1.4 for detecting large CNVs was used to record the consecutive 4 subsets of loss of heterozygosity as a minimum unit. If there are no more than 2 subset spans between the two units, the two units will be merged into a larger unit, and so on, and finally connected into a block. Perform F test on the R _Het value of , and calculate the p value of the block. When p≦0.01, we consider this block to have LOH, otherwise it is a non-LOH block.

Alternatively, the merging conditions can be set more stringently for more accurate detection. For example, in order to avoid false positives caused by some random errors, a region larger than 5M can be defined as a real LOH. On this basis, under the condition that the block fault tolerance is set to 1 (that is, the p value of 1 subset in the block is allowed to be greater than 0.01), the continuous subset with p≤0.01 in the vicinity of the subset with p≤0.01 in 2.1 will be merge with it. Finally, an F test is performed on the RHet in the merged area. If the p-value is less than 0.01, the block is considered to be a real LOH.

3. UPD detection

Combined with the above genome-wide CNV and LOH detection results, according to the Mendelian law of inheritance, UPD detection was performed.

If a DNA region shows a heterozygous state in the thousand-person data, that is, R _Het = 1, but in the actual detection, its heterozygous state disappears, that is, R _Het approaches 0, then it is determined that this region has LOH. , and if there is a CNV in this region with two copies at the same time (CN=2), that is, the copy number does not change (the sample in this example is a diploid sample, and each region of the normal diploid sample genome has two copies. copies), it was determined that uniparental diploidy (UPD) occurred in this region.

In 13 of the 15 samples, 10 LOH greater than 5M and 4 UPD greater than 5M were detected. The results are shown in Table 5. The detection of LOH and UPD was performed in the absence of paired samples (usually the use of own lesions). Compared with normal tissue, this is a paired sample, that is, a sample that has a certain correlation, and this embodiment detects LOH and UPD by comparing the target sample with multiple reference sample sets. The target sample and the reference sample are compared. The sample set has no correlation, so it is not a paired sample), the LOH test result of ≥ 5M is consistent with the CNV result of CN=1 (the CNV test result can be used to verify the accuracy of the LOH test result), the solution of the present invention detects the accuracy of LOH and UPD High performance, and can reach 5M level of accuracy.

The Circos plot (Fig. 5) comprehensively shows the CNV, LOH and UPD detection results of GM50275 samples.

table 5

Industrial Applicability

The method for determining the probe sequence based on the reference sequence of the present invention can be effectively used to determine the probe sequence, and the obtained probe is used to hybridize the captured genome to obtain a plurality of local regions of the genome, and the captured local regions can represent the whole Genome, which can reflect genome-wide variation information, is used to discover the occurrence of genome-wide structural variation.

Although specific embodiments of the present invention have been described in detail, those skilled in the art will understand. Various modifications and substitutions of those details may be made within the scope of the present invention in light of all the teachings disclosed. The full scope of the invention is given by the appended claims and any equivalents thereof.

In the description of this specification, reference to the terms "one embodiment," "some embodiments," "exemplary embodiment," "example," "specific example," or "some examples," or the like, is meant to incorporate the embodiment. A particular feature, structure, material, or characteristic described by an example or example is included in at least one embodiment or example of the present invention. In this specification, schematic representations of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the particular features, structures, materials or characteristics described may be combined in any suitable manner in any one or more embodiments or examples.

Claims (32)

1. A method for determining a probe sequence based on a reference sequence, comprising:

(1) constructing a first candidate probe set based on a plurality of discrete high-frequency SNP sites, wherein the first candidate probe set is composed of a plurality of candidate probes, each of the plurality of candidate probes corresponds to at least one discrete high-frequency SNP site, and the allele frequency of each of the plurality of discrete high-frequency SNP sites is at least 10 percent respectively;

(2) aligning the plurality of candidate probes in the first candidate probe set with a reference sequence to obtain an alignment result;

(3) performing a first screening on the first candidate probe set based on the comparison result to obtain a second candidate probe set consisting of a plurality of candidate probes,

wherein the first screening comprises retaining candidate probes that satisfy at least one of the following conditions:

a candidate probe that is uniquely aligned with the reference sequence;

candidate probes aligned to a plurality of positions of the reference sequence, and at least two of the plurality of positions each having a mismatch ratio of less than 10%;

(4) dividing the reference sequence into a plurality of windows with predetermined lengths respectively, and distributing a plurality of candidate probes in the second candidate probe set to the matched windows respectively so as to determine the position information of the candidate probes respectively;

(5) performing a second screening of the second candidate probe set based on the positional information and the allele frequencies of the discrete high frequency SNP sites to determine the probe sequences,

wherein the probe is determined according to the following steps:

(a) if a plurality of candidate probes are positioned in the same window, determining the candidate probe with the highest allele frequency of the corresponding discrete high-frequency SNP locus;

(b) and if the same window only has one candidate probe with the highest allele frequency of the corresponding discrete high-frequency SNP locus, selecting the candidate probe with the highest allele frequency of the corresponding discrete high-frequency SNP locus as the probe, and if the same window has a plurality of candidate probes with the highest allele frequency of the corresponding discrete high-frequency SNP locus, selecting the candidate probe which is closest to the center of the window from the candidate probes with the highest allele frequency of the corresponding discrete high-frequency SNP locus as the probe.

2. The method of claim 1, wherein the allele frequency of each of the plurality of discrete high frequency SNP sites is no more than 90% respectively.

3. The method of claim 1, wherein any two adjacent discrete high frequency SNP sites in the plurality of discrete high frequency SNP sites are not physically closer to the reference sequence than the length of the candidate probe.

4. The method of claim 1, wherein the candidate probe is 50-250 mers in length.

5. The method of claim 4, wherein the candidate probe is 100 mers in length.

6. The method of claim 1, wherein the candidate probe corresponds to one of the discrete high frequency SNP sites, and wherein the discrete high frequency SNP site corresponds to a mid-section of the candidate probe.

7. The method of claim 6, wherein the discrete high frequency SNP sites correspond to midpoints of the candidate probes.

8. The method of claim 1, wherein the candidate probe is truncated from the reference sequence.

9. The method of claim 1, wherein prior to performing the alignment, the first set of candidate probes is pre-screened in advance based on at least one of GC content and number of single base repeats of the candidate probes;

the prescreening includes retaining candidate probes that satisfy at least one of:

the GC content is 35 to 65 percent; and

the single base gravity was less than 7.

10. The method of claim 1, wherein in step (4), the reference sequence is divided into a plurality of windows each having the same predetermined length.

11. The method of claim 10, wherein the reference sequence is partitioned into a plurality of windows of 10Kb in length.

12. The method of claim 1, wherein after the second candidate probe set is subjected to the second screening, when the distance between two adjacent candidate probes in the second candidate probe set respectively falling into two adjacent windows on the reference genome is greater than the length of either of the two adjacent windows, the short tandem repeat sequence or a part of the short tandem repeat sequence located between the two adjacent candidate probes on the reference genome is further added to the second candidate probe set subjected to the second screening to form the probe sequence together.

13. The method of claim 1, wherein the reference sequence is a reference genome or a portion thereof.

14. A method of detecting a genomic structural variation comprising at least one of a chromosomal aneuploidy, a copy number variation, and an indel, for a non-diagnostic purpose, the method comprising,

(1) sequencing the genomic nucleic acid of the target sample to obtain a genomic sequencing result, the genomic sequencing result being comprised of a plurality of reads, wherein the sequencing comprises screening with a probe, wherein the probe is obtained by the method of any one of claims 1 to 13;

(2) dividing the reference genome into m regions, calculating the depth of coverage TD of region i using the number of reads falling into region i_iM and i are natural numbers, i represents the number of the region, i is more than or equal to 1 and less than or equal to m, 10<m；

(3) And determining whether the region i has structural variation or not based on the difference degree between the coverage depth of the region i and the coverage depth of the regions i of k reference samples, wherein k is a natural number and is more than or equal to 2.

15. The method of claim 14, wherein the depth of coverage of the region i is determined using the following equation:

or

Where i represents the number of the region.

16. The method of claim 14, wherein the test for the degree of difference between the depth of coverage of the genomic region i of the target sample and the depth of coverage of the region i of the k reference samples is performed by a t-test.

17. The method according to claim 14, wherein the comparison of the degree of difference between the depth of coverage of region i and the depth of coverage of region i of the k reference samples is performed by comparing the depth of coverage coefficients of genomic region i of the target sample and the reference sample, wherein the depth of coverage coefficient R of region i is_iThe determination of (a) comprises the steps of,

(a) for TD_iPerforming a first correction to obtain a first corrected coverage depth TD_aiThe first correction is implemented by performing linear regression on the depth-covered values of 2n consecutive areas including an area i, where n is a natural number, 10<n≤m/2；

(b) For TD_aiIs homogenized to obtain

Thereby obtaining

18. The method of claim 17, wherein in step (a), the first correction coverage depth TD is determined based on the following formula_ai：

Wherein, TD_jAnd j is a natural number, and j is more than or equal to 1 and less than or equal to n.

19. The method of claim 18, wherein in step (b), the TD is identified based on the following formula_aiIs homogenized to obtain

20. The method of any one of claims 16 to 19, wherein R is measured at the time of obtaining the target sample_iThen further comprises the step of reacting with R_iPerforming a second correction to obtain r_i，

Wherein R is_aiIs the average of the depth of coverage coefficients for k reference sample genomic regions i,

y is a natural number representing the reference sample number, R_i，yThe coverage depth coefficient for genomic region i of the reference sample y is shown.

21. The method of any one of claims 16 to 19, wherein R is measured at the time of obtaining the target sample_iThen further comprises the step of reacting with R_iPerforming a second correction to obtain r_i，

Wherein R is_aiThe mean of the depth of coverage coefficients for genomic region i for k reference samples and one target sample,

y is a natural number representing the reference sample number, R_i，yThe coverage depth coefficient for genomic region i of the reference sample y is shown.

22. The method of claim 20, wherein the t-test is performed such that the t-statistic for the genomic region i of the target sample is calculated as

Wherein,

r representing k reference samples_i,yAverage value of r_i,yFor reference to the second corrected depth of coverage coefficient of genomic region i of sample ygenome,

s is the standard deviation of k reference samples,

23. the method of claim 22, wherein t is based on the genomic region i of the target sample_iValue, obtaining a significance level P_iWhen P is_i<0.05, judging that the structural variation exists in the region i; otherwise, judging that the structural variation does not exist in the region i.

24. The method of claim 22, wherein t is based on the genomic region i of the target sample_iValue and predetermined significance level P_i0Obtaining t_iTheoretical value t_i0When t is_i≥t_i0Judging that the region i has structural variation, otherwise, judging that the region i has no structural variation; the predetermined P_i0≤0.05。

25. The method according to any one of claims 14 to 19, wherein after performing step (3), W regions that are co-directional and continuous are merged to obtain a primary merged region, when the two primary merged regions are co-directional and span no more than L regions, the two primary merged regions are merged to obtain a secondary merged region, and structural variation of the secondary merged region is detected based on the degree of difference between the coverage depth of the secondary merged region of the target sample genome and the coverage depth of the corresponding regions on the plurality of reference sample genomes; wherein, the equidirectional region refers to a region in which the t statistics of the region are both greater than 0 or both less than 0, W and L are both natural numbers, W is greater than or equal to 2, and L-W is less than or equal to 1.

26. A method for detecting loss of heterozygosity for non-diagnostic purposes comprising,

(1) sequencing the genomic nucleic acid of the target sample to obtain a genomic sequencing result, the genomic sequencing result being comprised of a plurality of reads, wherein the sequencing comprises screening with a probe, wherein the probe is obtained by the method of any one of claims 1 to 13;

(2) dividing a reference genome into m' regions, obtaining SNP sites shared by a target sample genome region i and a population region i to form a shared SNP set based on read information falling in the region i and data of the population region i in the genome sequencing result, respectively calculating the heterozygosity of fragments of the SNP sites in the target sample genome region i and the population shared SNP set, and obtaining a heterozygosity set U of the target sample genome region i_iAnd heterozygosity set U of population region i_0iComparing the target samples U_iAnd group U_0iTo determine whether there is loss of heterozygosity in the target sample region i; wherein, the segment where the SNP locus is located is a boundary point of two upstream SNPs and downstream SNPs adjacent to the SNP, m 'and i are natural numbers, m' is not less than 1 and not less than 6.

27. The method of claim 26, wherein each SNP in the common set of SNPs has an allele frequency greater than 0.1.

28. The method according to claim 26, wherein the heterozygosity of the fragment containing the SNP site is represented by a frequency coefficient R of a sub-allele of the SNP site_hetMAF/(1-MAF), which is the sub-allele frequency of the SNP.

29. The method of claim 28, wherein the comparison target sample U is_iAnd group U_0iTo determine whether loss of heterozygosity in the target sample region i has occurred, comprises determining U using an F-test_iVariance of (2)

And U_0iVariance of (2)

Whether there is a significant difference, if U_iAnd U_0iIf the variance difference is significant, it is determined that the target sample region i has loss of heterozygosity, otherwise, it is determined that the target sample region i has no loss of heterozygosity.

30. The method of claim 29, wherein the F-test comprises separately calculating U_iAnd U_i0Using the obtained target sample U_iVariance of (2)

And group U_i0Variance of (2)

Calculating to obtain two statistics F reciprocal to each other_upperAnd F_underObtaining significance level p using said reciprocal statistics_FComparison of p_FWith a predetermined significance level p_F0The size of (a), including the calculation formula,

p_F＝p_upper+(1-p_under) Wherein v is the number of SNPs in the high frequency SNP set shared by the target sample genomic region i and the population region i, q is the number of SNPs in the high frequency SNP set shared by the target sample genomic region i and the population region i, and R_het,i,vThe sub-allele frequency coefficient of the v-th SNP in the common high-frequency SNP set of the target sample genome region i,

is the average value of the sub-allele frequency coefficients, R, of q SNPs in the common high-frequency SNP set of the target sample genomic region i_het,i0,vThe sub-allele frequency coefficient of the v-th SNP in the shared high-frequency SNP set of the genome region i of the population sample,

is the average of the sub-allelic gene frequency coefficients of q SNPs in a common high-frequency SNP set of a population sample genomic region i, p_upperAnd p_underAre respectively according to F_upperAnd F_underObtaining of p_F0≤0.05。

31. The method according to any one of claims 26 to 30, wherein after step (2), W ' regions with loss of heterozygosity and continuity are merged to obtain a three-level merged region, when the span between the two three-level merged regions does not exceed L ' regions, the two three-level merged regions are merged to obtain a four-level merged region, a heterozygosity set of the target sample four-level merged region and a heterozygosity set of the same region of the population are respectively obtained, and the two heterozygosity sets are compared to determine whether loss of heterozygosity occurs in the target sample four-level merged region, wherein W ' and L ' are both natural numbers, W ' is not less than 2, and W '/2 is not less than L '.

32. A method for detecting an monadic diploid, said method being used for non-diagnostic purposes, characterized in that when the loss of heterozygosity is detected in a genomic region of a target sample, the copy number of the genomic region is calculated, and when the copy number of the genomic region is the same as that of the genomic region of a normal genome of the same species, the genomic region of the target sample is determined to be the monadic diploid; the determination of loss of heterozygosity in a genomic region of a target sample is carried out by the method of any one of claims 26 to 30.

Images