Method and system for detecting gene copy number variation based on deep information of population samples

Abstract

The present invention belongs to the field of biotechnology, and more specifically, the present invention discloses a method and system for detecting gene copy number variation based on deep information of population samples. The method comprises: (1) obtaining sequencing data of multiple samples; (2) calculating the copy ratio of the region to be detected relative to the whole genome; (3) sorting the copy ratios of the multiple samples from small to large and grouping them; (4) determining the nearest copy number reference ratio based on the distance between its copy ratio and the copy number reference ratio, and determining the copy number for each group. The method of the present invention accurately detects different types of CNV based on the deep information of population samples and calculates the copy number value; can accurately detect CNVs with common frequencies, and normally calculate the copy number values ​​including DUP; can detect different types of CNVs of population samples and correctly type them.

Description

Technical Field

[0001] The present invention relates to the field of biotechnology, and more specifically to a method and system for detecting gene copy number variation based on deep information of population samples. Background Art

[0002] Gene copy number variation (CNV) detection is generally based on several existing sequencing information: 1. Based on sequence assembly information, the assembled sequence is compared with the reference genome sequence; 2. Based on sequencing depth information, the sequencing depth of the target area is compared with the surrounding area and the area of ​​the control population to see if there is a difference; 3. Based on sequencing read length information, the length of the inserted fragment based on the paired-end sequencing (PE) read length comparison result and the situation where the read length can be cut and compared to different positions. The prior art article "PSCC: Sensitive and Reliable Population-Scale Copy Number Variation Detection Method Based on Low Coverage Sequencing" proposes a CNV detection method (referred to as PSCC method), which calculates CNV based on sequencing depth.

[0003] However, the methods of the prior art can only detect rare CNVs, but cannot detect CNVs with common frequencies. The reason is that the algorithms used to detect CNVs by sequencing depth methods will treat population samples as normal controls, that is, as if there is no variation, in order to ensure accuracy. Although this improves accuracy, it sacrifices sensitivity, especially for CNVs with common frequencies; or it can detect CNVs with common frequencies but with poor accuracy; the methods of the prior art cannot accurately type population genetic CNVs; moreover, the methods of the prior art cannot refine the specific copy number for copy number amplification (DUP). Summary of the invention

[0004] In view of the problems existing in the prior art, the present invention proposes a method and system for detecting gene copy number variation, which effectively solves the problems of common frequency CNV detection, population genetic CNV typing and refinement of DUP copy number.

[0005] Therefore, in one aspect, the present invention provides a method for detecting gene copy number variation based on population sample depth information, the method comprising:

[0006] (1) Obtain sequencing data for multiple samples;

[0007] (2) For each sample, for the region to be detected, the average sequencing depth of the region to be detected is calculated, and the average sequencing depth of the whole genome is calculated, and the average depth of the region to be detected is divided by the average depth of the whole genome to obtain the copy ratio of the region relative to the whole genome;

[0008] (3) sorting the copy ratios of the multiple samples from small to large, and comparing each value with the next value in turn starting from the minimum value, and if the difference between the two values ​​is less than a first threshold, merging the samples into a group, and grouping all the samples;

[0009] (4) For each group, the nearest copy number reference ratio is determined based on the distance between its copy ratio and the copy number reference ratio, and the copy number of each group is determined as the copy number.

[0010] In one embodiment, in (1), the plurality of samples is greater than 50 samples.

[0011] In one embodiment, in (2), the average sequencing depth of the region to be detected is corrected for GC.

[0012] In one embodiment, in (3), the first threshold is less than 0.15, more preferably less than 0.12, and most preferably less than 0.1.

[0013] In one embodiment, in (3), groups whose sample numbers are smaller than the second threshold value and groups whose sample numbers are larger than the second threshold value but whose distribution does not conform to the normal distribution are removed.

[0014] In one embodiment, determining whether a grouping conforms to a normal distribution includes calculating the mean, maximum, minimum, and variance of the copy ratio of each grouping.

[0015] In one embodiment, the second threshold value is greater than 25, preferably greater than 30.

[0016] In one embodiment, in (4), the copy number reference ratio is a set of arithmetic progressions starting from 0, the arithmetic progression is 1 / N, N is the species ploidy value, preferably, the sample is from a diploid species, and the discrete number starts from 0 and advances by 0.5. For example, the copy number reference ratio includes 0, 0.5, 1, 1.5, and 2.

[0017] In one embodiment, in (4), the copy number to which each group belongs is determined based on the distance between the average value of the copy ratio of the group and the copy number reference ratio. Preferably, the formula for determining the copy number to which each group belongs is: copy number = N*the nearest copy number reference ratio, where N is the species ploidy value.

[0018] In another aspect, the present invention provides a system for detecting gene copy number variation based on population sample deep information, the system comprising:

[0019] A sequencing data acquisition module, used to acquire sequencing data of samples;

[0020] A copy ratio calculation module, used to calculate the average sequencing depth of the region to be detected based on the sequencing data of the sample, calculate the average sequencing depth of the whole genome, and divide the average depth of the region to be detected by the average depth of the whole genome to obtain the copy ratio of the region relative to the whole genome;

[0021] A sample copy ratio grouping module, used to sort the copy ratios of multiple samples from small to large, and compare each value with the next value in turn starting from the minimum value, and if the difference between the two values ​​is less than a first threshold, the samples are merged into one group, and all samples are grouped;

[0022] The copy number determination module is used to determine the nearest copy number reference ratio according to the distance between the copy ratio of the group and the copy number reference ratio, and determine it as the copy number of the group. Preferably, the formula for determining the copy number of each group is: copy number = N*the nearest copy number reference ratio, where N is the species ploidy value.

[0023] In one embodiment, the system further comprises a group verification module for performing a normal distribution test on groups whose sample numbers are greater than a second threshold.

[0024] In one embodiment, the sequencing data acquisition module includes a sequencing instrument.

[0025] In one embodiment, the copy ratio calculation module is also used to perform GC correction on the average sequencing depth of the region to be detected.

[0026] In one embodiment, the first threshold is less than 0.15, more preferably less than 0.12, most preferably less than 0.1.

[0027] In one embodiment, the group verification module is used to calculate the average value, maximum value, minimum value and variance of the copy ratio of the group.

[0028] In one embodiment, the second threshold value is greater than 25, preferably greater than 30, more preferably greater than 50.

[0029] In one embodiment, the copy number reference ratio is a set of arithmetic progressions starting from 0, the arithmetic progression is 1 / N, N is the species ploidy value, preferably, the sample is from a diploid species, the discrete number starts from 0 and advances by 0.5. For example, the copy number reference ratio includes 0, 0.5, 1, 1.5, 2.

[0030] In one embodiment, the copy number determination module determines the copy number of each grouping based on the distance between the average value of the copy ratio of the grouping and the copy number reference ratio. Preferably, the formula for determining the copy number of each grouping is: copy number = N*the nearest copy number reference ratio, where N is the species ploidy value.

[0031] The method of the present invention can be closer to the actual situation, accurately detect different types of CNV based on the deep information of the population sample, and calculate the copy number value. The method of the present invention can accurately detect CNVs with common frequencies and normally calculate the specific values ​​of the copy number including DUP. The method of the present invention can detect different types of CNVs in population samples and correctly type them. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] The present invention is explained by means of the following drawings.

[0033] Figure 1 The distribution of copy ratios in this region for 790 samples is shown as an example.

[0034] Figure 2 The depth and copy ratio analysis of the chr8:39226335-39388919 region of sample SZCH0056 is shown as an example. DETAILED DESCRIPTION

[0035] Without wishing to be bound by any theory, the inventors assume that the target region may have different types of copy numbers in the test samples and the population samples, such as CN0, CN1, CN2, etc. For the same type such as CN2, the sequencing depth ratio calculated for the sample of this type in the target region should obey the normal distribution, and there are obvious differences in the distribution between samples with different copy numbers.

[0036] First, the sum of the sequencing depth of the target region is calculated, and the sum of the sequencing depth of the target region is divided by the length of the target region to obtain the average sequencing depth of the region, and preferably the average sequencing depth of the region is corrected by GC. The target region is at least 500 bp long, and can be tens of thousands of bp long. It is preferably a region known to have copy number variation, a region predicted to have copy number variation, or a candidate region with copy number variation. For a long target region, the target region can be segmented, for example, each segment is 500 bp or 1k bp, and the method or system of the present invention is applied to the segmented target region respectively, and adjacent segments with the same copy number variation are merged in the calculation result to obtain the full length of the segment with copy number variation. The method of GC correction can refer to the correction method of PSCC, which is briefly described as follows: divide the window to count the depth value and GC base ratio of each window, and for windows with the same GC ratio, take the median value of the depth value of these windows as the depth value under the GC ratio. Lasso regression is performed for all GC ratios and corresponding depth values ​​to obtain the GC ratio and the corresponding regression depth value, and the actual sequencing depth is standardized based on the regression depth of the GC ratio of the region. At the same time, the average sequencing depth of the whole genome of each sample is calculated. The average sequencing depth of the whole genome can be all the sequencing sequences of a sample, and it is not required to cover the whole genome of the species. In fact, there are a lot of repeated sequences in the genome, which are not suitable for sequencing. Then, the average sequencing depth of the region to be tested is divided by the average sequencing depth of the whole genome to obtain the copy ratio of the target region relative to the whole genome.

[0037] Then, sort the copy ratios of all samples from small to large, starting from the sample copy ratio with the smallest value. If the difference between the second sample copy ratio and the first sample copy ratio does not exceed 0.1, the two sample copy ratios are merged into one group; then compare the third sample copy ratio with the second sample copy ratio. If the size does not exceed 0.1, they are merged into one group, and so on, merging one by one; when the difference between the latter sample copy ratio and the previous sample copy ratio is greater than 0.1, the current grouping is terminated and the next grouping is restarted until all sample copy ratios are grouped.

[0038] Afterwards, normal distribution test quality control was performed for all groups with more than 30 samples; the mean, maximum, minimum and variance of the copy ratio of each group of samples were calculated. The distance between the mean values ​​of the copy ratios of each group of samples and the distance between the mean value of each group and the reference point (for example, 0, 0.5, 1, 1.5 and 2 for diploid species) were also calculated. For Nploid species, the advancement of discrete numbers can be different from that of diploid species, for example, in the form of 0 / N, 1 / N, 2 / N, 3 / N... N can be 1. For sperm or eggs, their genomes can be understood as haploid, in which case N can be regarded as 1, and their reference points can include 0, 1, 2, 3, 4...

[0039] Finally, the copy number of the target region of the sample in each group is determined according to the value of the reference point that is closest to the average value of each group; for N-ploid species, the copy number can be determined according to the multiple of its chromosomes and the attribution of the average value. Preferably, the formula for determining the copy number to which each group belongs is: copy number = N*the nearest copy number reference ratio, where N is the species ploidy value. For example, for diploid species, if the average value of a group is close to 0, the copy number is 0, if it is close to 0.5, the copy number is 1, and so on.

[0040] In another aspect, the present invention can be presented in a system manner. For example, the present invention relates to a system for detecting gene copy number variation based on population sample depth information, the system comprising a sequencing data acquisition module, a copy ratio calculation module, a sample copy ratio grouping module, a copy number determination module and a preferred grouping verification module.

[0041] In the present invention, the sequencing data acquisition module is used to obtain the sequencing data of the sample. The sequencing data acquisition module includes a sequencing instrument, and the sequencing instrument can sequence the sample to obtain sequencing data. The sequencing data acquisition module can also obtain the sequencing data of the sample from other places, such as sequencing data stored in a local computer or a remote server, or sequencing data stored in a medium such as a CD, a floppy disk, or a hard disk.

[0042] In the present invention, the copy ratio calculation module is used to calculate the average sequencing depth of the region to be detected based on the sequencing data of the sample, calculate the average sequencing depth of the whole genome, and divide the average depth of the region to be detected by the average depth of the whole genome to obtain the copy ratio of the region relative to the whole genome. The sequencing depth can be understood as the ratio of the total amount of bases (bp) obtained by sequencing to the size of the genome (Genome), which is one of the indicators for evaluating the amount of sequencing. The sequencing depth can be calculated for a sequence on the genome, that is, the ratio of the amount of bases (bp) obtained by sequencing a sequence to the size of the sequence.

[0043] In the present invention, the sample copy ratio grouping module is used to sort the copy ratios of multiple samples from small to large, and each value is compared with the next value in turn starting from the minimum value. If the difference between the two values ​​is less than the first threshold, the samples are merged into one group and all samples are grouped. Grouping the sample copy ratios is actually to separate samples with different copy numbers in the target area so that samples with the same copy number are in the same group. The choice of threshold is very important. If the threshold is too small, it will be divided into too many groups, and samples originally in one group will be divided into different groups; if the threshold is too large, samples originally in different groups will be divided into one group. In either case, the detected gene copy number variation will be biased. In one embodiment of the present invention, the first threshold is less than 0.15, more preferably less than 0.12, and most preferably less than 0.1; preferably, the above threshold is applicable to diploid species. For sperm or egg, its genome can be understood as haploid, and the first threshold can be less than 0.2. For Nploid species, for example, N is greater than 2, the first threshold is different from diploid species, and can be less than 0.05, or even less than 0.03.

[0044] In the present invention, the group verification module is used to perform a normal distribution test for groups whose sample numbers are greater than the second threshold value. In the present invention, the inventors found that the copy ratio of the population sample calculated in depth is significantly different under the normal distribution of the same genotype and different genotypes, based on which the copy number variation can be detected or judged. The method of the present invention can be calculated based on the copy ratio or based on the copy number. The method of the present invention can be used to detect CNV or to identify the true or false CNV. The normal distribution requires a large sample size, and it is difficult to observe the normal distribution of the grouping with a sample size in the single digit. In one embodiment, the second threshold value is greater than 25, preferably greater than 30, and more preferably greater than 50. The index of the normal distribution test is used to assist in the accuracy of the judgment result. If the grouping is less than 25, the grouping is not tested, and the auxiliary index output is empty.

[0045] In the present invention, the copy number determination module is used to determine the nearest copy number reference ratio according to the distance between the copy ratio of the group and the copy number reference ratio, and determine it as the copy number of the group. The mean value, maximum value, minimum value and variance of the copy ratio of each group of samples are calculated. The distance between the mean values ​​of the copy ratios of each group of samples and the distance between the mean value of each group and the reference point (for example, for diploid species, 0, 0.5, 1, 1.5 and 2, etc.) are also calculated. Finally, the copy number of the target area of ​​the sample in each group is judged according to the value of the reference point that is closest to the mean value of each group; for Nploid species, the copy number can be determined according to the multiple of its chromosomes and the attribution of the mean value. For example, for diploid species, if the mean value of a group is close to 0, the copy number is 0, if it is close to 0.5, the copy number is 1, and so on.

[0046] The copy ratio calculation module, the sample copy ratio grouping module, the grouping verification module and the copy number determination module can be implemented by a computer program, for example, by programming voice, writing a computer executable program according to the calculation formula of the sequencing depth, and then implementing it with the help of computer hardware.

[0047] Those skilled in the art should understand that the division and order of the steps in the method for detecting gene copy number variation based on deep information of population samples of the present invention are merely illustrative and not restrictive, and those skilled in the art may make deletions, additions, substitutions, modifications and changes without departing from the spirit and scope of the present invention as set forth in the attached claims and their equivalent technical solutions.

[0048] The present invention may be implemented as a computer-readable storage medium having a computer program stored thereon, which causes the steps of the method of the present invention to be performed when executed by a processor. In one embodiment, the computer program is distributed on a plurality of computer devices or processors coupled to a network so that the computer program is stored, accessed and executed in a distributed manner by one or more computer devices or processors. A single method step / operation, or two or more method steps / operations, may be performed by a single computer device or processor or by two or more computer devices or processors. One or more method steps / operations may be performed by one or more computer devices or processors, and one or more other method steps / operations may be performed by one or more other computer devices or processors. One or more computer devices or processors may perform a single method step / operation, or perform two or more method steps / operations.

[0049] Example 1

[0050] The present invention is described exemplarily by taking human chr20 chromosome 1561287-1594229 as the target region.

[0051] First, for the sequencing data of 790 samples, for each sample, the sum of the sequencing depths of the target region was calculated, and the sum of the sequencing depths of the target region was divided by the length of the target region to obtain the average sequencing depth of the region, and the average sequencing depth of the region was corrected for GC. At the same time, the average sequencing depth of the whole genome of each sample was calculated, and the average sequencing depth of the target region to be detected was divided by the average sequencing depth of the whole genome to obtain the copy ratio of the target region relative to the whole genome. The copy ratio distribution of the target region in the 790 samples is as follows: Figure 1 As shown;

[0052] Then, the copy ratios of all samples are sorted from small to large, starting from the sample copy ratio with the smallest value. If the difference between the second sample copy ratio and the first sample copy ratio does not exceed 0.1, the two sample copy ratios are merged into one group; then the third sample copy ratio is compared with the second sample copy ratio. If the size does not exceed 0.1, they are merged into one group, and so on, and merged one by one; when the difference between the latter sample copy ratio and the previous sample copy ratio is greater than 0.1, the current grouping is terminated and the next grouping is restarted until all sample copy ratios are grouped;

[0053] Then, normal distribution test quality control was performed for all groups with more than 30 samples; the mean, maximum, minimum and variance of the copy ratio of each group of samples were calculated. The distance between the mean values ​​of the copy ratio of each group of samples and the distance between the mean value of each group and the reference point (0, 0.5, 1, 1.5 and 2, etc.) were also calculated;

[0054] Finally, the copy ratios of the target area of ​​790 samples were divided into three groups; the first group had a mean value of 0.067, a standard deviation of 0.013, a distance of 0.561 from the second group, and was close to 0 in the reference point, so all samples in this group had 0 copies of the target area; the second group had a mean value of 0.627, a standard deviation of 0.031, a distance of 0.561 from the first group, and a distance of 0.560 from the third group, which was close to 0.5 in the reference point, so all samples in this group had 1 copy in this area; the third group had a mean value of 1.187, a standard deviation of 0.569, a distance of 0.560 from the second group, and was close to 1 in the reference point, so all samples in this group had 2 copies in this area.

[0055] Example 2

[0056] The sequencing depth and copy ratio of the chr8:39226335-39388919 region of sample SZCH0056 were analyzed, and the distribution of the samples in the population was as follows Figure 2 The copy ratio of this sample is 0.022, and the distribution in the population is around the 0 reference system, so the copy number is 0. The present invention uses the PSCC algorithm to detect this area but does not detect it, which is verified by the CMA chip of Affymetrix CytoScan750K.

[0057] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0058] Although the present invention has been described with reference to exemplary embodiments, it should be understood that the present invention is not limited to the configuration and method of the above-described embodiments. On the contrary, the present invention is intended to cover various modifications and equivalent configurations. In addition, although various elements and method steps disclosed in the present invention are shown in various exemplary combinations and configurations, other combinations including more or less elements or methods also fall within the scope of the present invention.

Claims

1. A method for detecting gene copy number variation based on population sample deep information, the method comprising: (1) Obtain sequencing data for multiple samples; (2) For each sample, for the region to be detected, calculate the average sequencing depth of the region to be detected, and calculate the average sequencing depth of the whole genome, and divide the average depth of the region to be detected by the average depth of the whole genome to obtain the copy ratio of the region relative to the whole genome; The average sequencing depth is the sum of the regional sequencing depths divided by the region length; (3) sorting the copy ratios of the multiple samples from small to large, and comparing each value with the next value in turn starting from the minimum value. If the difference between the two values ​​is less than a first threshold, merging the samples into a group and grouping all the samples; (4) For each group, the nearest copy number reference ratio is determined based on the distance between its copy ratio and the copy number reference ratio, and the copy number of each group is calculated: copy number = N * the nearest copy number reference ratio, where N is the species ploidy value; the copy number reference ratio is a set of arithmetic progressions starting from 0, with the arithmetic progression being 1 / N, and N being the species ploidy value.

2. The method according to claim 1, wherein in (3), the first threshold is less than 0.

15.

3. The method according to claim 1, wherein in (3), the first threshold is less than 0.

12.

4. The method according to claim 1, in (3), the first threshold is less than 0.

1.

5. The method according to any one of claims 1 to 4, wherein in (3), groups whose sample numbers are smaller than a second threshold value and groups whose sample numbers are larger than the second threshold value but whose distribution does not conform to a normal distribution are removed. The method according to claim 5 , wherein the second threshold is greater than 25. The method according to claim 5 , wherein the second threshold is greater than 30.

8. The method according to any one of claims 1 to 4, wherein in (2), the average sequencing depth of the region to be detected is corrected for GC.

9. The method according to any one of claims 1 to 4, wherein the sample is from a diploid species.

10. The method according to claim 9, wherein the copy number reference ratio includes 0, 0.5, 1, 1.5, and 2.

11. The method according to any one of claims 1 to 4, wherein in (4), the copy number to which each group belongs is determined based on the distance between the average copy ratio of the group and the copy number reference ratio.

12. A system for detecting gene copy number variation based on population sample deep information, the system comprising: A sequencing data acquisition module, used to acquire sequencing data of samples; A copy ratio calculation module, used to calculate the average sequencing depth of the region to be detected based on the sequencing data of the sample, calculate the average sequencing depth of the whole genome, and divide the average depth of the region to be detected by the average depth of the whole genome to obtain the copy ratio of the region relative to the whole genome; The average sequencing depth is the sum of the regional sequencing depths divided by the region length; A sample copy ratio grouping module, used to sort the copy ratios of multiple samples from small to large, and compare each value with the next value in turn starting from the minimum value, and if the difference between the two values ​​is less than a first threshold, the samples are merged into one group, and all samples are grouped; The copy number determination module is used to determine the nearest copy number reference ratio according to the distance between the copy ratio of the group and the copy number reference ratio, and calculate the copy number of the group: copy number = N * the nearest copy number reference ratio, N is the species ploidy value; the copy number reference ratio is a set of arithmetic progressions starting from 0, the arithmetic progression is 1 / N, and N is the species ploidy value.

13. The system of claim 12, wherein in (3), the first threshold is less than 0.

15.

14. The system of claim 12, wherein in (3), the first threshold is less than 0.

12.

15. The system of claim 12, wherein in (3), the first threshold is less than 0.

1.

16. The system according to any one of claims 12 to 15, further comprising a group verification module for performing a normal distribution test on groups whose sample numbers are greater than a second threshold. The system according to claim 16 , wherein the second threshold is greater than 25. The system according to claim 16 , wherein the second threshold is greater than 30. The system of claim 16 , wherein the second threshold is greater than 50.

20. The system according to any one of claims 12 to 15, wherein the copy ratio calculation module is further used to perform GC correction on the average sequencing depth of the region to be detected.

21. The system of any one of claims 12-15, wherein the sample is from a diploid species.

22. The system of claim 21, wherein the copy number reference ratios include 0, 0.5, 1, 1.5, and 2.

23. The system according to any one of claims 12 to 15, wherein the copy number determination module determines the copy number to which each group belongs according to the distance between the average copy ratio of the group and the copy number reference ratio.

Images