BACKGROUND
dna copy number variation  has long been known as a source of genetic variation, but its importance has only been recognized recently  <cit> .
in a landmark study in  <dig>  redon and colleagues found that  <dig>  cnv regions cover at least 12% of the human genome, with no large stretches exempt from cnv  <cit> .
the cnv regions cover more nucleotide content per genome than single nucleotide polymorphisms , suggesting the importance of cnv in genetic diversity  <cit> .
a common way to detect cnv is to utilize microarray-based methods  <cit> .
the most commonly used method, array comparative genomic hybridization  was first used to detect cnv a decade ago  <cit> .
microarray-based methods have revolutionized the way of how large-scale genome studies are carried out.
today, the next-generation sequencing technologies are transforming biology research  <cit> .
the rapid development of new sequencing technologies is continuously increasing the speed of sequencing and decreasing the cost.
the next-generation sequencing, such as  <dig>  <cit> , solexa  <cit>  and solid  <cit>  have already showed advantages over microarrays in several aspects.
apart from being rapid and cheap, data produced by sequencing can be re-used for varied purposes as opposed to data from microarray-based methods that can usually solely be used by one specific study.
in addition, reproducibility has been one of the major challenges for microarray technology  <cit> .
the once revolutionizing microarray-based chip-chip technology is being replaced by chip-seq, in which the dna fragments are sequenced instead of being hybridized to an array  <cit> .
sequencing-based methods are also used to produce genome-wide dna methylation profiles, detect snp, study chromosome translocations and rna transcriptome profiling  <cit> .
variation in sequencing coverage in genome assemblies has been used as an indicator for potential cnv between an assembled genome and shotgun data from another genome  <cit> .
this is analogous to a comparison of copy number between microarray probes and a single set of dna fragments.
there are two major problems with this kind of approach.
given a certain hybridization condition, hybridization efficiency varies among microarray probes.
likewise, given a certain alignment threshold, sequencing errors in combination with differences between genomes may result in erroneous distribution of the reads.
secondly, the number of probes on a microarray does not represent the real copy number of probe sequences in a genome.
likewise, the copy number of dna segments in an assembled genome may not represent the true one.
notably, the regions containing multiple copies are the most difficult to assemble correctly and is still the key unsolved problem in shotgun assembly  <cit> .
assembly errors like these cause false variation in the sequencing coverage and thus yield erroneous indication of cnv.
in this paper we describe an efficient solution based on a robust model that combines the advantages of acgh and high-throughput sequencing.
we also assessed cnv between two individuals .
an implementation of our method is freely available at .
RESULTS
the model
we have developed a method to detect cnv by shotgun sequencing, cnv-seq.
the method is based on a robust statistical model that allows confidence assessment of observed copy number ratios and is conceptually derived from acgh .
the microarray-based procedure, acgh involves a whole genome microarray where two sets of labeled genomic fragments are hybridized.
instead of a microarray, cnv-seq uses a sequence as a template and two sets of shotgun reads, one set from each target individual, x and y .
the two sets of shotgun reads are mapped by sequence alignment on a template genome.
we use a sliding window approach to analyze the mapped regions and cnvs are detected by computing the number of reads for each individual in each of the windows, yielding ratios.
these observed ratios are assessed by the computation of a probability of a random occurrence, given no copy number variation.
the random sampling in shotgun sequencing results in uneven coverage that may lead to observed coverage ratios that falsely imply cnv.
therefore, a statistical model is essential for the assessment of the probability of false positive ratios.
the average number of reads in a region of a genome is dependent on the total number of reads sampled, the length of the genome and the length of the region.
we use this relationship to compute a minimum sliding window size to achieve a desired minimum confidence level of the observations.
the mean number of reads for x and y genomes in a sliding window determines the distribution of the ratios.
the number of reads in a window is approximately distributed according to poisson, po, where the mean number of reads per window is λ, given by
λ=nwg
where n is the total number of sequenced reads, g is the size of the genome and w is the size of the sliding window, and w < <g.
we use the gaussian distribution to approximate the poisson distribution with mean and variance λ = μ = σ <dig>  this approximation is good when the mean number of reads per window is greater than  <dig> with continuity correction.
the predicted copy number ratio, r in each sliding window can be computed by
r=z×nynx
where z is the ratio of read counts in the window and nx and ny are the total number of reads in the genomes x and y respectively.
assuming an independent distribution of the read counts, we can obtain a probability, p of a copy number ratio being r or divergent from 1: <dig> ratio by a random chance.
for this purpose, we need the distribution of the read count ratio z. this distribution is given by gaussian ratio distribution  <cit> .
the computation of this distribution is cumbersome, but it can be transformed to another variable, t, by geary-hinkley transformation  <cit> :
t=μyz−μxσy2z2+σx <dig>
where μx, σx <dig>  μy and σy <dig> are the means and the variances for x and y respectively.
the new variable t approximately have a standard gaussian distribution when the mean number of reads per window is greater than  <dig> in y and less than  <dig>  in x. the p-value can be computed by
p={2×)if r≥ <dig> ×Φif r< <dig>
where Φ  is the cumulative standard gaussian distribution function.
the probability p decreases with increasing sliding window size  and we would like p to be as low as possible.
conversely, increased sliding window size leads to a decreased resolution of cnv regions.
therefore it is advantageous to compute a window size, which yields the best possible resolution according to a preset threshold of p for r. based on the above equations, we can calculate the best possible resolution, or the theoretical minimum window size, w by
w=gt22nxny
and
t={Φ−1)for det⁡ecting cnv with ratio≥r′,Φ−1for det⁡ecting cnv with ratio≤r′.
where p' is the desired significance level, and r' is the cnv detection threshold ratio.
Φ- <dig> is the inverse function of Φ. the number of reads sampled will affect the minimum window size.
for example, if one wants to detect cnv with ratio ≥  <dig> :  <dig> at significance level  <dig> , a genome size of  <dig> g bases and  <dig> m reads in both genomes will yield the minimum window size of  <dig>  bases, while  <dig> m reads will yield the window size of  <dig>  bases.
the use larger number of reads allows detection of ten times shorter cnv.
an alternative approach is to calculate the range of copy number ratios that can be detected at a certain significance level p', with a predefined window size w':
r≥1+1−abb, and, r≤1−1−abb
where
a=1−)2gnxw′
and
b=1−)2gnyw'
validation
in order to assess the performance of cnv-seq, we used simulated and real human data.
for the simulation of shotgun data, in total of  <dig> genomes were constructed, containing varied number, sizes and locations of cnv regions, snp and short insertions/deletions .
we simulated three sequencing methods, solexa,  <dig> and sanger  <cit>  for each constructed genome on  <dig> × to 8× coverage.
this resulted in the total of  <dig>  simulations.
the figure  <dig> shows the results of the simulations on varied coverage and varied p' for constant log <dig> =  <dig> .
each dot represents an average of  <dig> simulations and the sizes of the dots reflect the sizes of the lengths of the sliding windows that are the theoretical minimum lengths, given by equation .
the overall specificity for our method is between  <dig>  –  <dig> %, the sensitivity between  <dig>  –  <dig> % with the median of  <dig> % and  <dig> % respectively.
the mean sequence length is dependent on the technology simulated.
thus, in order to reach the same coverage, a larger number of fragments need to be sequenced when sequencing is performed with solexa, which produces short reads compared to the sanger and  <dig> methods.
according to our model, the largest number of sequenced reads yields the shortest length of the sliding window and thus the best resolution.
the range of window sizes in our simulations varies from  <dig>  bases to  <dig> , <dig> bases, decreasing with increasing average sequencing coverage.
the results show that our model performs well in the presence of errors.
despite of increased resolution due to shortening of the sliding window size, the sensitivity is increased together with increased sequencing coverage.
slight drop in specificity with increasing sequencing coverage can be observed .
this is likely to be due to snps, short indels, and read mapping errors, that are not considered in our statistical model and have a more profound effect on small windows.
the specificity does not drop in error free data.
the effect of errors may be reduced by using a window size that is larger than the theoretical minimum.
for example, the theoretical minimum window for 8× solexa sequencing at p =  <dig>  is  <dig> bases.
this window size gives a specificity of  <dig> %, while a  <dig> times larger window yields specificity of  <dig> % .
analysis of human data
the genomes of two individuals, dr. craig j. venter and dr. j. watson were recently sequenced on  <dig> × and  <dig> × coverage respectively  <cit> .
the genome of dr. craig j. venter is sequenced using sanger method and dr. j. watson's genome using  <dig> technology.
we compared the two genomes using cnv-seq .
the thresholds p' = 10- <dig> and log <dig> =  <dig>  yield sliding window size of  <dig>  bases for autosomal chromosomes.
the sex chromosomes have a lower sequencing coverage than autosomal chromosomes, therefore larger window sizes are used:  <dig>  bases for chromosome x and  <dig>  bases for chromosome y. we identified  <dig> contiguous regions covered by four or more consecutive windows.
the sizes of these regions range from  <dig>  bases to  <dig>  bases.
the comparison of the  <dig> cnv calls with those in the database of genomic variants   <cit>  revealed  <dig> of the cnv calls to overlap more than 50% with previously reported cnv regions.
in order to asses the significance of cnv calls, we performed  <dig>  permutation tests, using  <dig> randomly distributed cnv regions of the same sizes as in the original experiment.
in average, only  <dig> and maximum  <dig> of  <dig> regions overlap more than 50% with cnv in dgv   <dig>  random sets.
the real cnv calls have significantly larger overlap with dgv .
we also intersected the cnv calls with the cnvs identified by acgh in the two genomes.
there are  <dig> and  <dig> cnv regions reported in watson's and venter's genome respectively  <cit> .
we found  <dig> of our cnv calls overlap with  <dig> of previously reported watson's cnv regions, and only  <dig> of our cnv calls overlap with  <dig> of venter's.
the low overlap with venter's cnv calls made by acgh is not surprising, for the reason that the majority of the cnv regions were detected by only one of three microarray platforms  <cit> .
there are  <dig> cnv calls that made by cnv-seq but not acgh and overlap with dgv data, suggesting that cnv-seq can detect cnv regions that were missed by acgh.
one of these regions is shown in figure  <dig> , a  <dig> kb region  containing two genes  and one mirna .
we have used stringent thresholds in our analysis, thus by lowering thresholds, such as p-value and the number of consecutive windows, will increase the number of reported cnv calls.
a major assumption in cnv-seq is that shotgun sampling of dna fragments is random, and therefore the cnv calls made by cnv-seq are not due to different sequencing bias between the two sets of data compared.
when the two sets of data are prepared in the same way, this assumption is valid.
however, when the shotgun sequences are generated using two different sequencing methods, the assumption may not hold any more.
solexa sequencing reads are recently reported to be gc-biased dependent on a library preparation procedure  <cit> .
venter's genome was sequenced using  <dig> and watson's genome was sequenced using the sanger method.
we compared the distribution of gc frequencies in the shotgun reads in both genomes.
there are no significant differences between the two distributions .
CONCLUSIONS
we have developed a method to detect cnv using shotgun data.
our approach not only combines the advantages of microarray methods and high-throughput sequencing, but is also based on a robust statistical model allowing confidence assessment.
we tested the approach on both simulated and real data and the results show that the method can be applied to relatively low sequencing coverage with good specificity and sensitivity.
we have also developed a model to compute the theoretical limit of resolution for given data at a desired confidence level.
we expect the continued rapid development of sequencing technologies to further lower the cost and increase the speed of sequencing.
thus, sequencing-based approaches are likely to gain increased advantage over microarrays.
next-generation sequencing methods mostly produce a large number of short reads and our results show that the number of reads sequenced – not the length of the reads, is the most important factor that determines the resolution, i.e.
larger number of sequenced fragments results in increased resolution.
alternatively, given a constant resolution an increase in the number of sequenced reads will result in increased sensitivity and specificity.
therefore, a large number of short reads is an advantage as opposed to a small number of long reads.

 dna copy number variation has long been known as a source of genetic variation, but its importance has only been recognized recently. in a landmark study in redon and colleagues found that cnv regions cover at least 12% of the human genome, with no large stretches exempt from cnv. the cnv regions cover more nucleotide content per genome than single nucleotide polymorphisms, suggesting the importance of cnv in genetic diversity. a common way to detect cnv is to utilize microarray-based methods. the most commonly used method, array comparative genomic hybridization was first used to detect cnv a decade ago. microarray-based methods have revolutionized the way of how large-scale genome studies are carried out. today, the next-generation sequencing technologies are transforming biology research. the rapid development of new sequencing technologies is continuously increasing the speed of sequencing and decreasing the cost. the next-generation sequencing, such as, solexa and solid have already showed advantages over microarrays in several aspects. apart from being rapid and cheap, data produced by sequencing can be re-used for varied purposes as opposed to data from microarray-based methods that can usually solely be used by one specific study. in addition, reproducibility has been one of the major challenges for microarray technology. the once revolutionizing microarray-based chip-chip technology is being replaced by chip-seq, in which the dna fragments are sequenced instead of being hybridized to an array. sequencing-based methods are also used to produce genome-wide dna methylation profiles, detect snp, study chromosome translocations and rna transcriptome profiling. variation in sequencing coverage in genome assemblies has been used as an indicator for potential cnv between an assembled genome and shotgun data from another genome. this is analogous to a comparison of copy number between microarray probes and a single set of dna fragments. there are two major problems with this kind of approach. given a certain hybridization condition, hybridization efficiency varies among microarray probes. likewise, given a certain alignment threshold, sequencing errors in combination with differences between genomes may result in erroneous distribution of the reads. secondly, the number of probes on a microarray does not represent the real copy number of probe sequences in a genome. likewise, the copy number of dna segments in an assembled genome may not represent the true one. notably, the regions containing multiple copies are the most difficult to assemble correctly and is still the key unsolved problem in shotgun assembly. assembly errors like these cause false variation in the sequencing coverage and thus yield erroneous indication of cnv. in this paper we describe an efficient solution based on a robust model that combines the advantages of acgh and high-throughput sequencing. we also assessed cnv between two individuals. an implementation of our method is freely available at.  the model we have developed a method to detect cnv by shotgun sequencing, cnv-seq. the method is based on a robust statistical model that allows confidence assessment of observed copy number ratios and is conceptually derived from acgh. the microarray-based procedure, acgh involves a whole genome microarray where two sets of labeled genomic fragments are hybridized. instead of a microarray, cnv-seq uses a sequence as a template and two sets of shotgun reads, one set from each target individual, x and y. the two sets of shotgun reads are mapped by sequence alignment on a template genome. we use a sliding window approach to analyze the mapped regions and cnvs are detected by computing the number of reads for each individual in each of the windows, yielding ratios. these observed ratios are assessed by the computation of a probability of a random occurrence, given no copy number variation. the random sampling in shotgun sequencing results in uneven coverage that may lead to observed coverage ratios that falsely imply cnv. therefore, a statistical model is essential for the assessment of the probability of false positive ratios. the average number of reads in a region of a genome is dependent on the total number of reads sampled, the length of the genome and the length of the region. we use this relationship to compute a minimum sliding window size to achieve a desired minimum confidence level of the observations. the mean number of reads for x and y genomes in a sliding window determines the distribution of the ratios. the number of reads in a window is approximately distributed according to poisson, po, where the mean number of reads per window is λ, given by λ=nwg where n is the total number of sequenced reads, g is the size of the genome and w is the size of the sliding window, and w < <g. we use the gaussian distribution to approximate the poisson distribution with mean and variance λ = μ = σ this approximation is good when the mean number of reads per window is greater than with continuity correction. the predicted copy number ratio, r in each sliding window can be computed by r=z×nynx where z is the ratio of read counts in the window and nx and ny are the total number of reads in the genomes x and y respectively. assuming an independent distribution of the read counts, we can obtain a probability, p of a copy number ratio being r or divergent from 1: ratio by a random chance. for this purpose, we need the distribution of the read count ratio z. this distribution is given by gaussian ratio distribution. the computation of this distribution is cumbersome, but it can be transformed to another variable, t, by geary-hinkley transformation: t=μyz−μxσy2z2+σx where μx, σx μy and σy are the means and the variances for x and y respectively. the new variable t approximately have a standard gaussian distribution when the mean number of reads per window is greater than in y and less than in x. the p-value can be computed by p={2×)if r≥ ×Φif r< where Φ is the cumulative standard gaussian distribution function. the probability p decreases with increasing sliding window size and we would like p to be as low as possible. conversely, increased sliding window size leads to a decreased resolution of cnv regions. therefore it is advantageous to compute a window size, which yields the best possible resolution according to a preset threshold of p for r. based on the above equations, we can calculate the best possible resolution, or the theoretical minimum window size, w by w=gt22nxny and t={Φ−1)for det⁡ecting cnv with ratio≥r′,Φ−1for det⁡ecting cnv with ratio≤r′. where p' is the desired significance level, and r' is the cnv detection threshold ratio. Φ- is the inverse function of Φ. the number of reads sampled will affect the minimum window size. for example, if one wants to detect cnv with ratio ≥: at significance level, a genome size of g bases and m reads in both genomes will yield the minimum window size of bases, while m reads will yield the window size of bases. the use larger number of reads allows detection of ten times shorter cnv. an alternative approach is to calculate the range of copy number ratios that can be detected at a certain significance level p', with a predefined window size w': r≥1+1−abb, and, r≤1−1−abb where a=1−)2gnxw′ and b=1−)2gnyw' validation in order to assess the performance of cnv-seq, we used simulated and real human data. for the simulation of shotgun data, in total of genomes were constructed, containing varied number, sizes and locations of cnv regions, snp and short insertions/deletions. we simulated three sequencing methods, solexa, and sanger for each constructed genome on × to 8× coverage. this resulted in the total of simulations. the figure shows the results of the simulations on varied coverage and varied p' for constant log =. each dot represents an average of simulations and the sizes of the dots reflect the sizes of the lengths of the sliding windows that are the theoretical minimum lengths, given by equation. the overall specificity for our method is between – %, the sensitivity between – % with the median of % and % respectively. the mean sequence length is dependent on the technology simulated. thus, in order to reach the same coverage, a larger number of fragments need to be sequenced when sequencing is performed with solexa, which produces short reads compared to the sanger and methods. according to our model, the largest number of sequenced reads yields the shortest length of the sliding window and thus the best resolution. the range of window sizes in our simulations varies from bases to, bases, decreasing with increasing average sequencing coverage. the results show that our model performs well in the presence of errors. despite of increased resolution due to shortening of the sliding window size, the sensitivity is increased together with increased sequencing coverage. slight drop in specificity with increasing sequencing coverage can be observed. this is likely to be due to snps, short indels, and read mapping errors, that are not considered in our statistical model and have a more profound effect on small windows. the specificity does not drop in error free data. the effect of errors may be reduced by using a window size that is larger than the theoretical minimum. for example, the theoretical minimum window for 8× solexa sequencing at p = is bases. this window size gives a specificity of %, while a times larger window yields specificity of %. analysis of human data the genomes of two individuals, dr. craig j. venter and dr. j. watson were recently sequenced on × and × coverage respectively. the genome of dr. craig j. venter is sequenced using sanger method and dr. j. watson's genome using technology. we compared the two genomes using cnv-seq. the thresholds p' = 10- and log = yield sliding window size of bases for autosomal chromosomes. the sex chromosomes have a lower sequencing coverage than autosomal chromosomes, therefore larger window sizes are used: bases for chromosome x and bases for chromosome y. we identified contiguous regions covered by four or more consecutive windows. the sizes of these regions range from bases to bases. the comparison of the cnv calls with those in the database of genomic variants revealed of the cnv calls to overlap more than 50% with previously reported cnv regions. in order to asses the significance of cnv calls, we performed permutation tests, using randomly distributed cnv regions of the same sizes as in the original experiment. in average, only and maximum of regions overlap more than 50% with cnv in dgv random sets. the real cnv calls have significantly larger overlap with dgv. we also intersected the cnv calls with the cnvs identified by acgh in the two genomes. there are and cnv regions reported in watson's and venter's genome respectively. we found of our cnv calls overlap with of previously reported watson's cnv regions, and only of our cnv calls overlap with of venter's. the low overlap with venter's cnv calls made by acgh is not surprising, for the reason that the majority of the cnv regions were detected by only one of three microarray platforms. there are cnv calls that made by cnv-seq but not acgh and overlap with dgv data, suggesting that cnv-seq can detect cnv regions that were missed by acgh. one of these regions is shown in figure, a kb region containing two genes and one mirna. we have used stringent thresholds in our analysis, thus by lowering thresholds, such as p-value and the number of consecutive windows, will increase the number of reported cnv calls. a major assumption in cnv-seq is that shotgun sampling of dna fragments is random, and therefore the cnv calls made by cnv-seq are not due to different sequencing bias between the two sets of data compared. when the two sets of data are prepared in the same way, this assumption is valid. however, when the shotgun sequences are generated using two different sequencing methods, the assumption may not hold any more. solexa sequencing reads are recently reported to be gc-biased dependent on a library preparation procedure. venter's genome was sequenced using and watson's genome was sequenced using the sanger method. we compared the distribution of gc frequencies in the shotgun reads in both genomes. there are no significant differences between the two distributions.  we have developed a method to detect cnv using shotgun data. our approach not only combines the advantages of microarray methods and high-throughput sequencing, but is also based on a robust statistical model allowing confidence assessment. we tested the approach on both simulated and real data and the results show that the method can be applied to relatively low sequencing coverage with good specificity and sensitivity. we have also developed a model to compute the theoretical limit of resolution for given data at a desired confidence level. we expect the continued rapid development of sequencing technologies to further lower the cost and increase the speed of sequencing. thus, sequencing-based approaches are likely to gain increased advantage over microarrays. next-generation sequencing methods mostly produce a large number of short reads and our results show that the number of reads sequenced – not the length of the reads, is the most important factor that determines the resolution, i.e. larger number of sequenced fragments results in increased resolution. alternatively, given a constant resolution an increase in the number of sequenced reads will result in increased sensitivity and specificity. therefore, a large number of short reads is an advantage as opposed to a small number of long reads. 
array comparative genomic hybridization is commonly used for cnv detection, but the microarray platform has a number of inherent limitations. our results show that the number of reads, not the length of the reads is the key factor determining the resolution of detection. we also show the results for assessment of cnv between two individual human genomes. the method is based on a robust statistical model that describes the complete analysis procedure and allows the computation of essential confidence values for detection of cnv. dna copy number variation has been recognized as an important source of genetic variation.
 dna copy number variation has been recognized as an important source of genetic variation. array comparative genomic hybridization is commonly used for cnv detection, but the microarray platform has a number of inherent limitations.  here, we describe a method to detect copy number variation using shotgun sequencing, cnv-seq. the method is based on a robust statistical model that describes the complete analysis procedure and allows the computation of essential confidence values for detection of cnv. our results show that the number of reads, not the length of the reads is the key factor determining the resolution of detection. this favors the next-generation sequencing methods that rapidly produce large amount of short reads.  simulation of various sequencing methods with coverage between × to 8× show overall specificity between – %, and sensitivity between – %. we also show the results for assessment of cnv between two individual human genomes. 