4
$\begingroup$

I begin with a matrix of features (~30000 rows) vs cases (columns). I have binary - TRUE or FALSE - data in each cell indicating whether there is an event spanning the feature for that case. Next, I sum each row to get the total number of cases from the sample set that have an event for each feature. I will call this c for each feature. Next, I retrieve a list of features. How can I test whether c (the total number of cases with an event for a given feature) is significantly greater than c among the feature population in general? Furthermore, is there some way of calculating some uncertainty measure of my result based on the fact that only known features associated with a particular process are in the list.

NB. My attempt (probably wrong and overly complicated):

  1. Bootstrap with replacement 100000 times to get 100000 samples each of size 1023 from the 30000 features. The 1023 features (for example) are excluded from the population from which the features for each sample are drawn.
  2. Use var.test across c for each sample vs the 1023 feature list. If the test does not produce a significant value, that sample can be used in step 3. Samples that produce a significant value are discarded. This is done to satisfy the assumption of the Wilcoxon test that the variances between the two test samples are similar.
  3. Do a Kolmogorov–Smirnov test (ks.test) for each sample brought forward vs the 1023 feature list. If the test does not produce a significant value, that sample can be used in step 4. Samples that produce a significant value are discarded. This is done to satisfy the assumption of the Wilcoxon test that the distributions between the two test samples are similar.
  4. Do a Wilcoxon test (wilcox.test) for each sample brought forward vs the 1023 feature list. If the mean (or median?) p-value for these tests is significant, I can say that c is significantly greater for this set of features than it is for features that do not belong to this set.
$\endgroup$
4
  • $\begingroup$ You say that the sample is not random but the statistical procedures you are using including the bootstrap implicitly assume random samples. So any statistical testing procedrue may be invalid. I don't know how to address the specifics of your problem but I can recommend Brad Efron's Large Scale Inference book as a source for techniques to handle these types of problems. $\endgroup$ Commented May 3, 2012 at 14:58
  • $\begingroup$ The bootstraps are meant to be randomly sampled. The only "sample" that is not random is the gene set to which all the random samples are compared. $\endgroup$ Commented May 3, 2012 at 15:03
  • $\begingroup$ Bootstrapping is random sampling with replacement from the original data. It is intended to get approximations for the sampling distribution so that you can make inferences without requiring parametric assumptions to give you the sampling distribution. It is not used to generate data. The way you describe the use of bootstrap make me wonder how you intend to apply the bootstrap. Could you explain? $\endgroup$ Commented May 3, 2012 at 20:09
  • $\begingroup$ @MichaelChernick I took the liberty of moving your post + comment where I believe they should have been posted (and clean up nice comments you got on them from Macro and cardinal). As they kindly pointed out, this site shall be viewed as a wiki (with full support for edits, history, $\LaTeX$, etc.) where we intend to collect long-lasting responses to specific questions about statistical analysis, with minimal clutter in the response thread. Should you want to ask for clarifications or comment on the present question, you can always make use of the comment functionality. Welcome to this site! $\endgroup$
    – chl
    Commented May 3, 2012 at 21:03

1 Answer 1

1
$\begingroup$

Your idea of prescreening samples to fit the assumptions of the Wilcoxon test, while throwing away anything that does not fit is very wrong, with every step messing up the inference in an additional way. And it is also unnecessarily complicated. The goal of looking at randomly selected gene sets is to figure out how c would look like for a typical set of genes that has nothing to do with copy number aberrations. It would be much simpler if you could define a summary function of the 1023 c values that would capture what you are interested in. Perhaps the mean, or the median, etc. I'll assume that the median c value would be meaningful.

  1. Generate $B$ random gene sets of size 1023.
  2. For each of the sets compute the median c value, denoted $mc_i, i=1,\ldots,B$
  3. Compute the median c value for the actual list of interest, $mc_0$
  4. The p-value would be the proportion of the $mc_i$ values that exceed $mc_0$: $$p=\frac{\# (mc_i \geq mc_0)+1}{B+1}$$ (the "$+1$" terms are often used to avoid 0 p-values, but are not really important).

But you have to be careful with the interpretation of the results. If the test is significant, it does mean that the selected gene set has higher median c value than a random gene set, but it might be due to a variety of reasons. For example, if the stability genes are more likely to be on chromosomes with more CNA's, or if they are more likely to be located close to each other (assuming CNA's cover many neighboring genes simultaneously), etc. You could address some of these alternative explanations by being careful in step 1 above. For example, you could select random gene sets that have the same distribution over the chromosomes as the stability set, or you could impose some restrictions on the spacing of the genes, etc.

$\endgroup$
5
  • $\begingroup$ Thanks,I'll get going trying to implement this. I believe the median is indeed more apt. Sorry for my ignorance, but is there a name for this test? Also, should I drop the 1023 stability genes from the population from which the $B$ random gene sets are drawn? $\endgroup$ Commented May 3, 2012 at 19:27
  • $\begingroup$ I am not sure what the name would be, perhaps something about a Monte-Carlo estimation of the null distribution. And do not drop the 1023 stability genes - you want the actual set of interest be a possible random set. $\endgroup$
    – Aniko
    Commented May 3, 2012 at 20:24
  • $\begingroup$ Ok thanks, may I ask why I shouldn't drop the actual set of interest? $\endgroup$ Commented May 3, 2012 at 20:32
  • $\begingroup$ Also, should replacement be allowed. I assume it should be. $\endgroup$ Commented May 3, 2012 at 20:52
  • $\begingroup$ It should not matter a lot whether you do or don't allow replacement. However you certainly should not drop the actual set of interest, because you want to have the actual set as part of your sampling space: you should be able to "randomly" select this set. $\endgroup$
    – Aniko
    Commented May 4, 2012 at 12:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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