I am performing a genomic study where I am interested in determining the heterogeneity of epigenetic alleles (DNA methylation patterns) in specific genomic locations. Previous work on this subject has lead to several different statistical measures of heterogeneity/entropy being developed, but the underlying assumptions for these statistical measures is not well defined in any source that I can locate, and I am not confident in my statistical knowledge to understand what they might be.

The type of statistic I am speaking of is easiest explained through example: Given a sequence of length 4 with each position in the sequence capable of taking on the value of either 1 or 0, there are $2^4$ distinct possible patterns that this sequence could take (1111, 1110 .... 1000, 0000, etc.). We can sample this sequence a number of times, and say that $N$ is the total number of observations, while $n_i$ is the number of times each pattern is observed. $p_i$ is the frequency that pattern $i$ is observed, calculated as $n_i / N$.

One statistic that has been used as a measure of heterogeneity of the patterns observed is called epipolymorphism which was to my knowledge first used in this paper, and is calculated as below:

$$ epipolymorphism = 1-\sum_{i=1}^n p_i^2 $$

$p_i$ = frequency that pattern $i$ is observed ($n_i/N$)

$n = $ the number of different patterns observed

For reasons that are not immediately obvious, the researchers who devised this statistic limited their analysis to sequences that they had $N > 40$ observations . I am curious if there are statistical reasons why using an $N > 40 $ would be important? I understand that to observe all possible $2^4 = 16$ patterns you would need to have a minimum of 16 observations, but I am not sure I understand why it is important to have 2.5 times this many observations, and how one would decide where to draw the line for minimum $N$ for this calculation?

Please let me know if I can clarify the question at all. I had trouble creating an accurate title, so if you have suggestions for a more descriptive title or tags please let me know.

  • 1
    $\begingroup$ One ought to refer to the original definition for the reasons. Do you have a reference? One thing that is apparent, though, is your formula expresses the opposite of heterogeneity: it attains its maximum when the $p_i$ are as non-heterogeneous as possible (they are all equal) and is minimal when the range of the $p_i$ is as large as possible. It might be worth noting that $1-n\operatorname{Var}(p)=1-\sum p_i^2+O(1/n^2),$ thereby relating this statistic to the variance of the $p_i$ for large $n$. $\endgroup$ – whuber Sep 28 '16 at 20:32
  • 1
    $\begingroup$ It often happens in statistics that a concept has many, many names. Your research might be aided by the observation that your expression for epipolymorphism is exactly the definition of the gini impurity. IAlso, notice that your definition is distinct from entropy which is $-\sum_i p_i\ln(p_i)$ in the resources with which I am familiar. $\endgroup$ – Sycorax says Reinstate Monica Sep 28 '16 at 20:34
  • $\begingroup$ @whuber I've added a link to the original paper in which I came across this statistic. I may have explained the goal of the epipolymorphism incorrectly. I believe the goal is to determine which regions have the least bias in regards to a specific pattern, in which case if the frequencies of each pattern observed are equal, that would show that there was no bias, and hence the highest epipolymorphism. It is not measuring whether the frequencies are heterogeneous, but rather the patterns themselves. Does that make sense? $\endgroup$ – Reilstein Sep 28 '16 at 20:56
  • $\begingroup$ @Sycorax Thanks a lot for linking to the gini impurity, that might help me better understand the goal of this statistic. I am familiar with entropy, as it is another statistic that has been used to capture heterogeneity of epialleles. I have used Shannon's Information Entropy which looks like the equation you listed, as well as a modified version which multiplies the Shannon's Entropy times a scaling factor $1/b$ where b is the length of the sequence. epipolymorphism has been used most frequently in the studies of "heterogeneity" that I have seen, but I'm exploring entropy as well, thanks $\endgroup$ – Reilstein Sep 28 '16 at 21:02

The statistical issue might be easiest to think about in terms of finding the true frequency of heads for a coin that might be biased (binomial distribution). If you only flip a fair coin 3 times per experiment, then in one out of eight experiments you will get all heads, and in one out of eight experiments you will get all tails. That is, in one-quarter of the experiments with a fair coin you will find an empirical frequency that is as far away as possible from the true frequency of 1/2 for heads, if you only sample 3 times in your experiment.

You might think analogously about the case of 16 distinct possibilities (a multinomial distribution) sampled only 40 times. With not even 3 times as many samples as there are possible outcomes, you will not have much assurance that the observed frequencies for each possibility are close to the true underlying frequencies of the population. For example, S. K. Thompson examined sampling from multinomial distributions (The American Statistician, Vol. 41, No. 1, 1987, pp. 42-46; available through JSTOR if you have that access). Based on his table, with only 40 samples, half the time in a worst-case scenario at least one of the observed probabilities $p_i$ will be more than 0.1 unit away from the true probability. (If each of the 16 possibilities is equally likely, each true probability would be 0.0625.)

In principle, to "draw a line" at $N=40$ might have been based on a judgement by the authors about the error that such limited sampling would make in their calculation of "epipolymorphism." In this case, however, the requirement for $N>40$ might instead represent an attempt to omit technically questionable results rather than to achieve some strictly statistical purpose.

Data analyzed for "epipolymorphism" in the paper were reduced representation bisulfite sequencing (RBBS) results. This method determines whether certain cytosine residues next to guanines (CpG) in DNA have been methylated, one type of "epigenetic" modification of different parts of the genome. As next-generation sequencing (NGS) was used for RBSS, several consecutive CpG sites could be analyzed on each single piece of DNA. The authors examined thousands of sets of 4 consecutive CpGs (4-mers) across the human genome and determined how the variability of methylation patterns for the 4-mers changed as a function of their average CpG methylation, among many biological samples.

If all 4 CpGs in a 4-mer are completely methylated or completely unmethylated then there is no variability (0 epipolymorphism) but at intermediate average methylation there are distinct, biologically interesting, possibilities. For example, 25% average CpG methylation of a 4-mer could come from one out of every 4 DNA pieces having all 4 of its CpGs methylated, or all DNA pieces having one identically placed CpG methylated out of the 4, or something more random providing 25% average methylation. The "epipolymorphism" (noted in a comment by @Sycorax to be the same as Gini impurity) was evidently intended to capture the variability in this epigenetic modification among individual pieces of DNA covering each analyzed 4-mer in the genome, a type of epigenetic heterogeneity among cells.

The supplementary data for the paper suggest that 300 to 600 individual pieces of DNA ("reads") were typically analyzed for each 4-mer for each biological sample. Having no more than 40 reads available for a particular 4-mer suggests substantial under-representation of that 4-mer (only about 1/10 of the average representation) in that DNA sample. RBBS involves several processing steps including the polymerase chain reaction (PCR), which can have different amplification ability on different regions of the genome.

Omitting under-represented genomic regions in NGS can be a way to try to avoid regions of the genome that might have been affected by technical limitations. The hard cutoff at $N>40$ might have been Procrustean but it was a simple way to proceed; it still left tens of thousands of 4-mers to examine in almost all samples. Clearly the 4-mers with low representation would also have been less reliable in terms of identifying their methylation patterns, but I would guess that the authors simply looked at the distribution of 4-mer coverage among all genomic sites and found that those with $N\le40$ appeared to be technical outliers.

The authors of the paper could presumably clarify how they came up with that cutoff.

  • $\begingroup$ Thanks for your insight! Are you referring to figure 11A when you mention that their typical coverage was 300-600 reads per loci? I think that figure is showing the typical read depth for amplicon sequencing, rather than RRBS. Is there another figure I am missing? I see the logic if your argument in observing the coverage distribution and setting the cutoff at a level that avoids outliers, however for the # of reads they got from RRBS samples I would expect their mean coverage to be 5-10x, so the cutoff does not really control for outliers, it enriches for high cov. outliers right? $\endgroup$ – Reilstein Oct 6 '16 at 18:35
  • $\begingroup$ And you're right, the authors can probably help answer this question, I was just hoping for some insight on the statistical backing of different coverage thresholds, because I am having trouble conceptualizing how the statistic may be biased by higher or lower coverage cutoffs. Thanks again! $\endgroup$ – Reilstein Oct 6 '16 at 18:37
  • $\begingroup$ @Reilstein I used supplementary table 1, which contains for each sample a column for the total mapped reads and the number of 4-mers with at least 40X RRBS coverage that were seen in at least 5 samples. Just took the ratio as a rough estimate of coverage per 4-mer, so I suppose that could be off. I'll edit my answer in a day or two to get more directly at the statistical problem, which basically is that with few observations you can't have much confidence that the observed frequencies will be close to the true frequencies in the population. $\endgroup$ – EdM Oct 6 '16 at 19:41
  • $\begingroup$ I see how you calculated now. I compared their reads/library to a different paper where they provide mean coverage per RRBS library, and that paper had a mean coverage of 10-15x, but 3x as many reads per library, so I estimate the mean coverage from the present paper to be 5-10x. Thanks for taking the time to think about and comment on this issue. I think I'm beginning to understand, but I would love to hear the rest of your extended explanation in an edit. Thanks! $\endgroup$ – Reilstein Oct 6 '16 at 19:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.