6
$\begingroup$

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.

enter image description here

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.

$\endgroup$
3
$\begingroup$

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

enter image description here

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 enter image description here

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.

$\endgroup$
  • $\begingroup$ Hmm. I don't understand. Aren't two samples supposed to be (positively) correlated whenever we run a paired t-test? I mean, that's the whole point of running a paired t-test, compared to an unpaired t-test: namely, paired test has higher power precisely when there is correlation due to pairing. Are you saying that one should expect to see such skewed distributions of p-values under the null for a paired t-test? That would be weird. $\endgroup$ – amoeba says Reinstate Monica Jul 8 '17 at 22:13
  • $\begingroup$ I don't have Matlab installed on the laptop I am currently working, but I ran your code in octave-online.net - I only get a matching histogram if I replace ttest with ttest2. Are you sure your plot matches the code? $\endgroup$ – amoeba says Reinstate Monica Jul 8 '17 at 22:13
  • $\begingroup$ @amoeba You're right! The original code block should say test2, i.e. the first figure is from unpaired t-tests. That was a combination of errors and late-night editing. Things are corrected now. I re-ran with paired t-tests and different beta, and that's included now too. $\endgroup$ – R Greg Stacey Jul 9 '17 at 1:12
  • $\begingroup$ Thanks for fixing the code/images. However, I believe that your interpretation of what is happening here is pretty much entirely wrong. E.g. this "This distribution is possible if your two samples, alpha and beta, are correlated." must be wrong as I tried to explain in the comment above. I am thinking that what you see is probably because the assumptions of the t-test are grossly violated in the samples you generate: they are not normal but coming from the uniform distribution, and with very low $n$ t-test behaves weird. $\endgroup$ – amoeba says Reinstate Monica Jul 9 '17 at 7:33

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.