Interpretation of p-value histogram for differential methylation analysis: what can explain prevalence of large p-values?

Question

I am conducting an analysis of methylation data for ~20k+ genes.

For n=3, I am doing a t-test for every gene to see if it has been differentially methylated after treatment. Methylation values range from [0,1].

So, in my data I have 6 rows, with a before and after row for each patient and # columns = # of genes (~20k).

Here is code for generating p-values:

for (i in 1:ncol(df))
{
  alpha = c(df[c(1,3,5),i])
  beta = c(df[c(2,4,6),i])
  df.p[i] = t.test(alpha,beta,paired=TRUE)$p.value
}

hist(df.p)

This churns out the p-values just fine, but when I make a histogram of the resulting p-values, it is skewed strongly to the left, which is strange since you'd at least expect a uniform distribution if you didn't have any significant methylation differences. Below is a screenshot of the distribution.

Side note: I've also used the limma package in BioConductor and got the same results.

Am I conducting the t-test wrong? How am I to interpret these results? Any advice is appreciated, as I am a novice biostatistician.

Answer 1

Edit: As Amoeba pointed out there was an error in my code, and the first plot is from an unpaired t-test. I re-ran with paired t-tests and different alpha and beta, and the results are at the bottom.

I simulated it in Matlab with the code below. I generated 1000 random alphas (1x3 vectors); corresponding betas are made by adding a smaller-magnitude, random 1x3 vector to each alpha. The result is that alpha and beta are correlated. The resulting p-value distribution is given here (unpaired t-test). As you can see it's skewed to near-1 values.

Unpaired t-test

pp = zeros(1000,1);
for iter = 1:1000
  alpha = rand(1,3);
  beta = alpha + rand(1,3)*.1;
  [~,pp(iter)] = ttest2(alpha,beta);
end

Edit: Paired t-test

pp = zeros(10000,1);
for iter = 1:10000
  alpha = randn(1,3);
  beta = alpha + [rand 0.1 -1*rand];
  [~,pp(iter)] = ttest(alpha,beta);
end

To achieve the right-skewed distribution with a paired t-test I had to make different assumptions about beta: beta = alpha + [r1 0.1 r2], where $r1$ and $r2$ are uniform random numbers, and $r1$ is in the range [0 1] and $r2$ is in the range [-1 0]. That is $r1$ is always positive, $r2$ is always negative, and their combined effect on alpha-beta is, on average, 0. The effect is to fairly reproduce your original histogram.

If this actually describes your data, it means that if e.g. gene methylation is high in one replicate, it's low in another replicate. Which is weird. Although it could happen if your data is a ratio and you've swapped numerator and denominator in one replicate. To test that, can you make scatter plots / calculate correlation between each pair of replicates, i.e. between df[1,] and df[3,], df[1,] and df[5,], etc. If your replicates are good they should all have high positive correlation - if any are negatively correlated then that's evidence for swapping the numerator and denominator.