阅读说明：本技术 一种基于滤波降噪的基因组结构变异检测方法 (Genome structure variation detection method based on filtering noise reduction ) 是由 刘志岩 刘珍 王海宁 姜玥 于 2021-07-07 设计创作，主要内容包括：本发明提出了在一种基于滤波降噪的基因组拷贝数变异检测方法,该方法考虑了染色体序列本身的特性以及序列中碱基GC含量对读深的影响,在更好的对数据进行预处理的基础上,将读深数据与高斯核函数进行卷积运算得到尺度空间图像函数,并对尺度空间图像进行边缘检测及基准检测,确定候选拷贝数区域,及检测拷贝数变异类型和位置,提高了拷贝数变异检测的精度。(The invention provides a filtering noise reduction-based genome copy number variation detection method, which considers the characteristics of a chromosome sequence and the influence of base GC content in the sequence on reading depth, performs convolution operation on the reading depth data and a Gaussian kernel function to obtain a scale space image function on the basis of better preprocessing the data, performs edge detection and reference detection on the scale space image, determines a candidate copy number region, detects the type and the position of copy number variation, and improves the accuracy of copy number variation detection.)

1. A genome copy number variation detection method based on filtering noise reduction is characterized by comprising the following steps:

s1, preprocessing data;

extracting a read depth signal from the bam file by using a SAMtools tool, wherein the read depth signal consists of the following two signals: r_m＝r_m+E_m，R_mRepresenting the actual value of the read depth signal observed, r_mRepresenting the expected read depth signal in the chromosomal sequence, E_mRepresents a noise signal;

using haar functionsAnd (3) removing noise:

removing the influence of the base GC content on a read depth signal by GC correction;

s2, obtaining a scale space image;

read the depth data r [ i-j]Performing convolution operation with Gaussian kernel function K (I, j) to obtain scale space image function I_SS[i，l]：

Wherein the content of the first and second substances,

σ_lrepresents the scale parameter of the l-th layer, and m represents the window value size of the Gaussian kernel function K (i, j);

s3, detecting the edges of the scale space image;

scale space image function I_SS[i，l]The components are obtained under different scales x and y, and the components areAndat each scale, the module value MI of each pixel point is calculated_SS[i，l]And phase angle AI_SS[i，l]：

MI_SS[i，l]At phase angle AI_SS[i，l]Points with maximum values corresponding to discontinuities in the scale space image, from MI_SS[i，l]And AI_SS[i，l]The extreme point can be obtained, so that the edge detection is carried out on the scale space image;

s4, detecting the scale space graph reference;

setting three reference standards m_t(l)、m_t(l)+λδ_t(l) And m_t(l)-λδ_t(l) Wherein m is_t(l) And delta_t(l) Is a function of the image in scale space I_SS[i，l]At a zero-crossing point where there are two non-zerosFunction Z_SS[i，l]At the mean and standard deviation of layer i, λ is the reference check coefficient, the normal range of the scale space function values is m (k) ± 2 δ (k), scale space function values outside the normal range will be filtered out;

s5, determining a candidate copy number area;

if Z is_SS[s_m，l，l]·Z_SS[e_m，l，l]< 0, l-th layer middle zone { i | s_m，l≤i≤e_m，lAll points in satisfy Z_SS[i，l]0, and interval { i | s_m，l≤i≤e_m，lUpper scale space image function I_SS[i，l]Mean value ofAt m_t(l)+λδ_t(l) And m_t(l)-λδ_t(l) In between, then [ s_m，l，e_m，l]Is a region of candidate copy number variation; where i is at zero crossing point function Z_SS[i，l]Within the corresponding position interval on the l layer;

s6, detecting copy number variation;

mean of scale space image functionAt m_t(l)+λδ_t(l) Above, copy number variation increases; mean of scale space image functionAt m_t(l)-λδ_t(l) In contrast, the copy number variation is absent.

2. The method for detecting genomic copy number variation based on filtering noise reduction according to claim 1, wherein: in step S1, noise signal E_mThe mathematical expectation of (A) is E (X)_K) μ, (k ═ 1,2, … …), variance D (X)_K)＝σ²Not equal to 0, (k ═ 1,2, … …); where μ is the expected value of the random variable, σ²Is variance of。

3. The method for detecting genomic copy number variation based on filtering noise reduction according to claim 1, wherein: in the step S1, in the step S,the function after shifting the unit k to the right or left becomesWherein k is an integer.

4. The method for detecting genomic copy number variation based on filtering noise reduction according to claim 1, wherein: in step S3, setting a threshold TH for each scale of edge image, and obtaining edge points; reserving points which are larger than or equal to the threshold TH as edge points, and setting the points which are smaller than the threshold TH to be zero; and solving the edge of the image in the scale space, and linking edge points under different scales to obtain the edge of the image in the scale space under different scales.

5. The method for detecting genomic copy number variation based on filtering noise reduction according to claim 1, wherein: in step S6, when the current time is in the m-th interval [ S ] of the l-th layer_m，l，e_m，1]When the copy number variation is found, searching from the l-1 layer to the m-th interval [ s ]_m，l-1，e_m，l-1]And circulating the steps until the original read-depth signal is iterated.

6. The method for detecting genomic copy number variation based on filtering noise reduction according to claim 1, wherein: in step S2, [ sigma ]_lThe scale parameter representing the ith layer, the scale space image is smoother and smoother as the sigma increases, and the contour of the scale space image is kept unchanged in the process of smoothing.

7. The method for genomic copy number variation detection based on filtering noise reduction according to claim 1,the method is characterized in that: in step S2, m represents the window value size of the gaussian kernel K (i, j), and m is 3 σ_l。

Technical Field

The invention relates to the field of bioinformatics, in particular to a genome copy number variation detection method based on filtering and noise reduction.

Background

The development of New Generation Sequencing (NGS) technology is more and more mature, each sequencing platform is endless, the sequencing cost of gene sequences is greatly reduced, and the sequencing speed is higher and higher, so that the DNA sequence data generated by sequencing is very huge, and the accuracy of data processing becomes urgent.

Wavelet transform is an advantageous instrument for removing data noise in signal processing, and scale space filtering is an important means for visually describing signals. Sequencing data in the field of bioinformatics often need to consider the influence of noise on a model when modeling data due to high repeatability of biological sequences and inevitable errors in a sequencing process. In the detection of copy number variation, because the existence of repeated sequence fragments of human chromosome sequences can cause inevitable errors in sequencing, in order to describe the copy number variation by using a proper scale, multi-scale spatial filtering is adopted to find out a zero crossing point in a read-depth signal, and finally a copy number variation region is determined.

In addition, with the implementation and development of the human genome project and 1000genome project, the scale of sequence data of protein, DNA and RNA is increasing, and it has not been possible to meet the practical needs to study biological genetic variation and disease generation only by means of biological experiments, so that it is necessary to study and clarify biological problems from mass data by means of theoretical and conceptual methods of computer, mathematics and other subjects. Copy number variation detection is one of the effective methods for studying gene structure changes in bioinformatics.

The current techniques applied to copy number variation detection mainly include:

1. comparative Genomic Hybridization (CGH): this technology has been developed to date, and has been combined with chip technology (Microarray) to derive chip comparative genomic hybridization technology (Array-CGH). The technology can detect the copy number of DNA sequences among different genomes on the level of all chromosomes or chromosome sub-bands, thereby finding copy number variation. However, this technique has resolution at the Mb level, and copy number fragments of smaller fragments are not easily detected. Meanwhile, the technology is complex to operate, low in flux, long in time consumption, high in cost, and not beneficial to large-scale popularization, and a large amount of template DNA is needed.

MLPA: the method is a copy number detection method developed in 2002 and is fully called a multiple ligation probe amplification technology. At present, corresponding kits are available for detecting diseases such as SMA, Down syndrome and the like. The technology has a relatively accurate relative quantification function. However, the method has the disadvantages of complex probe preparation, complicated operation steps and long time consumption. And capillary electrophoresis is adopted as an analysis means, so that the flux is low, the cost is high, the method belongs to open operation, and the pollution of PCR products is easily caused.

Disclosure of Invention

The method considers the characteristics of a chromosome sequence and the influence of the GC content of a base in the sequence on the read depth, performs convolution operation on the read depth data and a Gaussian kernel function to obtain a scale space image function on the basis of better preprocessing the data, performs edge detection and reference detection on the scale space image, determines a candidate copy number area, and detects the type and the position of the copy number variation, thereby improving the accuracy of the copy number variation detection.

The method specifically comprises the following steps:

s1, preprocessing data;

extracting a read depth signal from the bam file by using SAMtools, wherein the read depth signal consists of two signals of R_m＝r_m+E_m，R_mRepresenting the actual value of the read depth signal observed, r_mRepresenting the expected read depth signal in the chromosomal sequence, E_mRepresents a noise signal;

using haar functionsAnd (3) removing noise:

removing the influence of the base GC content on a read depth signal by GC correction;

s2, obtaining a scale space image;

read the depth data r [ i-j]Performing convolution operation with Gaussian kernel function K (x, y) to obtain scale space image function I_SS[i，l]:

Wherein the content of the first and second substances,

σ_lrepresents the scale parameter of the l-th layer, and m represents the window value size of the Gaussian kernel function K (i, j);

s3, detecting the edges of the scale space image;

scale space image function I_SS[i，l]The components are obtained under different scales x and y, and the components areAndat each scale, the module value MI of each pixel point is calculated_SS[i，l]And phase angle AI_SS[i，l]：

MI_SS[i，l]At phase angle AI_SS[i，l]Points with maximum values corresponding to discontinuities in the scale space image, from MI_SS[ i, l and AI_SS[i，l]The extreme point can be obtained, so that the edge detection is carried out on the scale space image;

s4, detecting the scale space image reference;

setting three reference standards m_t(l)、m_t(l)+λδ_t(l) And m_t(l)-λδ_t(l) Wherein m is_t(l) And delta_t(l) Is a scale space image function y_SS[i，l]At a zero-crossing point function Z having two non-zeros_SS[i，l]At the mean and standard deviation of layer i, λ is the reference check coefficient, the normal range of scale-space values is m (k) ± 2 δ (k), scale-space function values outside the normal range will be filtered out;

s5, determining a candidate copy number area;

if Z is_SS[s_m，l,l]·Z_SS[e_m，l,l]<0, l-th layer middle area { is_m，l≤i≤e_m，lAll points in satisfy Z_SS[i，l]0, and interval { i | s_m，l≤i≤e_m，lUpper scale space image function I_SS[i，l]Mean value ofAt m_t(l)+λδ_t(l) And m_t(l)-λδ_t(l) In between, then [ s_m，l,e_m，l]Is a region of candidate copy number variation; where i is at zero crossing point function Z_SS[i，l]Within the corresponding position interval on the l layer;

s6, detecting the copy number variation type and position;

mean of scale space image functionAt m_t(l)+λδ_t(l) Above, copy number variation increases; mean of scale space image functionAt m_t(l)-λδ_t(l) Next, the copy number variation is absent.

Further, in step S1, the noise signal E_mThe mathematical expectation of (A) is E (X)_K) μ, (k ═ 1,2, … …), variance D (X)_K)＝σ²Not equal to 0, (k ═ 1,2, … …); where μ is the expected value of the random variable, σ²Is the variance.

Further, in step S1,the function after shifting the unit k to the right or left becomesWherein k is an integer.

Further, in step S3, a threshold TH is set for the edge image of each scale, and an edge point is obtained; reserving points which are larger than or equal to the threshold TH as edge points, and setting the points which are smaller than the threshold TH to be zero; and solving the edge of the image in the scale space, and linking edge points under different scales to obtain the edge of the image in the scale space under different scales.

Further, in step S6, when the current time is in the m-th interval [ S ] of the l-th layer_m，l,e_m，l]When the copy number variation is found, searching from the l-1 layer to the m-th interval [ s ]_m，l-1,e_m，l-1]And circulating the steps until the original read-depth signal is iterated.

Further, in step S2, σ_lThe scale parameter representing the ith layer, the scale space image is smoother and smoother as the sigma increases, and the contour of the scale space image is kept unchanged in the process of smoothing.

Further, in step S2, m represents the window value size of the gaussian kernel function K (i, j), and m is 3 σ_l。

The terms used herein are explained as follows:

copy Number Variation (CNV): refers to a change in copy number of at least a portion of a nucleic acid molecule in a sample to be tested, compared to the corresponding nucleic acid sequence in a normal sample, wherein said portion has a length of more than 1 kb. The circumstances and reasons for copy number variation may include: deletions, such as microdeletions; insertions, such as microinsertions, microreplications, replications; inversion, transposition and complex multi-site variation.

Sequencing: refers to the process of obtaining information on the nucleic acid sequence of a sample. Sequencing can be performed by various methods, including but not limited to dideoxy chain termination; preferably, the high throughput sequencing method includes, but is not limited to, next generation sequencing techniques or single molecule sequencing techniques. The higher the sequencing depth, the higher the detection sensitivity, i.e., the smaller the length of deletion and repeat fragments that can be detected.

Reading a deep signal: refers to a nucleic acid sequence of a certain length (usually longer than 20bp), such as the sequencing result of a sequence generated by a sequencer, which can be aligned with a specific region or position of a reference sequence by a sequence alignment method.

Indexing: refers to a nucleic acid sequence having a specific length and performing a labeling function. When the DNA molecules to be tested are derived from a plurality of samples to be tested, each of the plurality of samples may be added with a different index for distinguishing the plurality of samples during sequencing.

GC content deviation: there is some GC bias between batches or within a batch, which may result in copy number bias appearing in regions of the genome with high or low GC content. CG correction was performed using control set-based sequencing data to obtain the corrected relative number of reads in each window, whereby such bias can be eliminated and the accuracy of detecting copy number variation can be improved.

Drawings

FIG. 1 is a flow chart of the method for detecting genomic copy number variation based on filtering and noise reduction according to the present invention.

FIG. 2 is a diagram illustrating a signal comparison after data preprocessing in the mutation detection method.

FIG. 3 is a diagram of an exemplary zero crossing point and inflection point of the smoothed signal in the mutation detection method.

Detailed Description

The present invention is further illustrated by the following examples, which are provided to illustrate and not to represent all the possibilities of the present invention, and the present invention is not limited to the materials, reaction conditions or parameters mentioned in the examples, and any person skilled in the relevant art can use other similar materials or reaction conditions to perform the detection of the variation in the copy number of the gene described in the present invention according to the principles of the present invention. Without thereby departing from the basic concept described in the present invention.

The words "include," "have," or any other variation thereof, as used herein, are intended to cover a non-exclusive inclusion. For example, a process, method, article, or apparatus that comprises a list of elements is not necessarily limited to those elements, but may include other elements not expressly listed or inherent to such process, method, article, or apparatus.

Referring to fig. 1, a flow chart of the method for detecting genome copy number variation based on filtering noise reduction according to the present invention is shown: the variation detection method specifically comprises the following steps:

s1, preprocessing data

The read depth signal is extracted by the tool SAMtools from the bam file, where the information that the sequencing reads match into the reference sequence is stored. Extracting a read count file from the bam file by using a SAMtools tool, wherein the file comprises read counts values and corresponding position information, and considering a reading depth signal as being composed of the following two signals:

R_m＝r_m+E_m (1)

in the formula (1), wherein R_mRepresenting the actual value of the read depth signal observed, r_mRepresenting the expected read depth signal in the chromosomal sequence, E_mThe representative noise signal is generally considered to be a white gaussian noise signal. White noise is noise with a uniform power distribution over an infinitely wide frequency range, and is only an idealized noise model. White gaussian noise follows a gaussian distribution with statistical regularity of its amplitude, and "white" in the definition means that its power spectrum is constant in the entire frequency domain. Noise signal E_mThe mathematical expectation of (A) is E (X)_K) μ, (k ═ 1,2, … …), variance D (X)_K)＝σ²Not equal to 0, (k ═ 1,2, … …). Where μ is the expectation (or mean) of a random variable, σ²Is its variance.

We can use a signal-based approach to find breakpoints in the read-deep signal. Wavelet theory is useful in both denoising and detecting breakpoint information of a signal. Wavelet de-noising is to find the optimal mapping of wavelet function space from noisy signals to filter and analyze the signals, so that high-frequency interference signals can be filtered, the original characteristics of the signals can be successfully reserved, and the obtained characteristic signals and the low-pass filtered signals are combined and reconstructed.

We use haar functionsAnd (3) removing noise:

R_mactual values representing the observed read depth signal:

for the integer number k of the series of integers,is as followsThe result after shifting the unit k to the right or left.

Fig. 2 is a diagram showing a comparison between signals after data preprocessing. Our use of haar's function for noise removal is mainly due to the nature of its structural fit to read depth data, the copy number variation along the chromosome is present in blocks, and the adjacent loci that are marked have the same copy number increase or decrease.

After the wavelet denoising is carried out, due to the influence of the base GC content on the read depth signal, a GC correction is needed, and the GC correction adopts a GC correction method in the prior art, and the detailed description is omitted here.

S2, obtaining a scale space image

Read the depth data r [ i-j]Performing convolution operation with Gaussian kernel function K (I, j) to obtain scale space image function I_SS[i，l]:

σ_lThe scale parameter representing the ith layer, the scale space image will become smoother and smoother as σ increases, and its contour will remain unchanged during the smoothing process. m represents the window value size of the gaussian kernel function K (i, j), and is 3 σ as default_l，σ_lThe range of variation of (a) determines the size of the length of the copy number variation interval that can be detected. The ratio between two adjacent scale parameters determines the temporal complexity and the accuracy of the copy number variation interval that can be detected. If a smaller scale is used, the detection accuracy is traded for higher detection accuracy, whereas if a larger scale is used, the detection accuracy is traded for lower time complexity. Therefore, a proper proportion is selected to obtain high detection precision and small time complexity.

S3, edge detection of image in scale space

The filtering process of the scale space function is to smooth the image signal to be detected under different scales by using a smoothing function and find out the mutation point of the signal according to the first order or the second order derivative of the wavelet transform coefficient model of the smoothed signal. The extreme points of the first derivative correspond to the zero crossing points of the second derivative and the inflection points of the smoothed signal, see fig. 2. Image edges can be detected by local maxima of the wavelet transform model.

Scale space image function I_SS[i，l]The components are obtained under different scales x and y, and the components areAndat each rulerUnder the degree, the module value MI of each pixel point is calculated_SS[i，l]And phase angle AI_SS[i，l]：

MI_SS[i，l]At phase angle AI_SS[i，l]Points with maximum values corresponding to discontinuities in the scale space image, from MI_SS[ i, l and AI_SS[i，l]Extreme points can be found, and therefore edge detection can be conducted on the scale space image.

Preferably, a threshold TH may be set for the edge image of each scale, an edge point is obtained, a point greater than or equal to TH is reserved as an edge point, and a point less than TH is set to zero; and solving the image edge, and linking the edge points under different scales to obtain the image edge under different scales.

S4, detecting scale space graph reference

In order to further remove the influence of outliers on the detection accuracy, it is necessary to perform reference detection on the scale space pattern. Here we set three reference standards, m respectively_t(l)、m_t(l)+λδ_t(l) And m_t(l)-λδ_t(l) In that respect Wherein m is_t(l) And delta_t(l is a scale space image function I_SS[i，l]At a zero-crossing point function Z having two non-zeros_SS[i，l]In the mean and standard deviation of the l-th layer, λ is the reference check coefficient, preferably 3. To filter out outliers, the normal range of scale-space values is m (k) ± 2 δ (k), outside of which scale-space function values are to be filtered out;

s5, determining candidate copy number area

When the zero crossing points of each layer are found, the candidate region for copy number variation can start from the zero crossing point of each layer with two non-zero values, and the interval [ s ]_m，l,e_m，l]In the mth zone of the l-th layer, in the region { i | s_m，l≤i≤e_m，lWhere i is at the zero crossing point function Z_SS[i，l]Within the corresponding location interval on the l-th layer.

If Z is_SS[s_m，l,l]·Z_SS[e_m，l,l]<0, l-th layer middle zone { i | s_m，l≤i≤e_m，lAll points in satisfy Z_SS[i，l]0, and interval { i | s_m，l≤i≤e_m，lUpper scale space image function I_SS[i，l]Mean value ofAt m_t(l)+λδ_t(l) And m_t(l)-λδ_t(l) In between, then [ s_m，l,e_m，l]Is a region of candidate copy number variation; where i is at zero crossing point function Z_SS[i，l]Within the corresponding position interval on the l layer;

s6, copy number variation detection

Copy number variation detection includes detecting the type of copy number variation (addition and deletion) and detecting the precise location of the copy number variation region. The increase in copy number variation is defined as: mean of scale space functionAt m_t(l)+λδ_t(l) Above, the copy number variation deletion (LOSS) is defined as: mean of scale space function At m_t(l)-λδ_t(l) The following is a description.

When in the m-th interval of the l-th layer [ s ]_m，l,e_m，l]When the copy number variation is found, searching from the l-1 layer to the m-th interval [ s ]_m，l-1,e_m，l-1]And circulating the steps until the original read-depth signal is iterated.

In another aspect of the present invention, a computer-readable storage medium may be provided, the computer-readable storage medium having stored thereon machine-executable instructions, which when executed, cause a machine to perform the steps of the method for detecting genomic structural variation based on filtering and noise reduction according to the present invention.

In the present invention, a computer readable storage medium may be a tangible device that can hold and store instructions for use by an instruction execution device. The computer readable storage medium includes, for example, but is not limited to, an electronic memory device, a magnetic memory device, an optical memory device, an electromagnetic memory device, a semiconductor memory device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: a portable computer diskette, a hard disk, a Random Access Memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), a Static Random Access Memory (SRAM), a portable compact disc read-only memory (CD-ROM), a Digital Versatile Disc (DVD), a memory stick, a floppy disk, a mechanical coding device, such as punch cards or in-groove projection structures having instructions stored thereon, and any suitable combination of the foregoing. Computer-readable storage media as used herein is not to be construed as transitory signals per se, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through a waveguide or other transmission medium (e.g., optical pulses through a fiber optic cable), or electrical signals transmitted through electrical wires.

The machine-executable instructions described herein may be downloaded from a computer-readable storage medium to a respective computing/processing device, or to an external computer or external storage device via a network, such as the internet, a local area network, a wide area network, and/or a wireless network. The network may include copper transmission cables, fiber optic transmission, wireless transmission, routers, firewalls, switches, gateway computers and/or edge servers. A network adapter card or network interface in each computing/processing device receives machine-executable instructions from the network and forwards the machine-executable instructions for storage in a computer-readable storage medium in the respective computing/processing device.

The foregoing is a more detailed description of the invention in connection with specific preferred embodiments and it is not intended that the invention be limited to these specific details. For those skilled in the art to which the invention pertains, several simple deductions or substitutions can be made without departing from the spirit of the invention, and all shall be considered as belonging to the protection scope of the invention.