1
$\begingroup$

I have an RNA_Seq dataset (FPKM values) obtained from the ENCODE Caltech datasets for K562 cells; I also have 6 different lists of genes and I want to show that the expression values of my panel of genes is not random. Therefore, I decided to use R and perform bootstrapping.

Bootstrapping, as far as I know, requires normally distributed data. My RNA-Seq data are not normally distributed and they are strongly skewed on the left.

I used BoxCox transformation:

box.cox.powers(RNA_Seq[,2])

to normalize them.

I obtained a negative value (-0.03701075 to be precise) and raised my FPKM values to that. I then used:

scale()

to obtain the z-scores and see if the z-scores of my genes are bigger/smaller than +/-1.96. that's because I am interested to see if my panel of genes is more highly expressed -or it has a lower expression- than what you would expect by chance.

However, the values with the highest FPKM became those with the lowest z-score and vice versa since I elevated to a negative number. I understand why this happens mathematically; however, I don't really know how to handle the bootstrap now because the data show the exact opposite of what I would reasonably expect!

So, how would you guys handle this situation to test my genes of interest?
I looked up online but I could not find anything specific for RNA-Seq, or I didn't know how to apply it to RNA-Seq data.

I asked this question on Stack Overflow but it was marked as off-topic and I should have posted it here. So here we go :)

$\endgroup$
3
  • 3
    $\begingroup$ Why do you think bootstrapping requires normal data? Also, what exactly are you trying to show? "Not random" is not a testable hypothesis, without specifying what distribution you want to show they are not random on. $\endgroup$
    – Peter Flom
    Commented Aug 31, 2013 at 20:17
  • $\begingroup$ Well, I want to show that the way we selected the group of genes led to a cohort of genes whose expression signature is significantly different that what you would expect by chance...I thought that bootstrap would be the correct approach to show that, in that way we can calculate an empirical p-value associated with it...which way would you suggest? $\endgroup$
    – TriRook
    Commented Sep 1, 2013 at 3:18
  • $\begingroup$ I think it would be best to start with what you want to do (in terms of what questions you want to ask of your data) free of statistical terminology and your own guesses about what might be suitable things to do. $\endgroup$
    – Glen_b
    Commented Oct 1, 2013 at 5:53

1 Answer 1

2
$\begingroup$

ENCODE Caltech dataset contains two replicates for K562 cell line. Is that the experiment you are using? And how did you select the gene lists that you currently have at hand? Further, do you want to test all the genes in the list as a "set", or do you want to have a single p-value for each gene?

It is very common to select the gene lists on the basis of some statistical test. For example, R add-on packages edgeR, DESeq and limma (from the Bioconductor project) offer suitable methods for rna-seq data, and especially for small sample sizes. There are also other ways to do this, such as, simple filtering based on FPKM values or their standard deviation. Using some statistical test to arrive to list of differentially expressed genes will also simultaneously give you a p-value for each gene, which might be what you are after. Now that you have a single condition only (just one cell line, K562), the statistical test is equivalent to testing whether the gene's expression is different from 1, or 0 if the FPKM values are also log-transformed.

Each gene's statistical significance can also be tested by using a permutation test, where usually the sample labels are shuffled a large number of times. This approach is not really applicable here, since there are only two samples. For more information on a possible implementation, see this page.

In addition, it is also possible to test the whole genelist as a set, and possibly compare it to all the other genes in the experiment. This can be accomplished with, e.g., the package globaltest (from the Bioconductor project) in R.

$\endgroup$
2
  • $\begingroup$ JTT, thanks for the answer! yes, that's the experiment I'm using and the selected list of genes was selected based on epigenetic marks on the TSS of those genes. The approach that I ended up using was: 1) starting from FPKM values, sample 100k times (with replacement)the whole RNASeq dataset (let's call it A)and take the median of the 100k vectors 2)sample my genes dataset (let's call it B) and take the median 3)measure the 95th percentile and the mean of both 4)see if the average of B falls outside the 95th percentile of A 5)use Wilcoxon Rank Test to test the difference..makes sense? $\endgroup$
    – TriRook
    Commented Oct 1, 2013 at 12:26
  • $\begingroup$ Unfortunately that does seem to be quite right. I would strongly suggest you to take a loog at the globaltest package from the Bioconductor project, since it addresses this type of situations, and includes, among other, permutation testing. The original paper is bioinformatics.oxfordjournals.org/content/20/1/93.long, and the package documentation is bioconductor.org/packages/2.13/bioc/vignettes/globaltest/inst/…. See especially Gene Set Testing chapter of the manual. Also, use the argument method="poisson" in the function gt(). $\endgroup$
    – anon
    Commented Oct 2, 2013 at 13:55

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.