I tried bioconductor packages for differential gene expression analysis such as EdgeR, Deseq2, Limma and obtained expressed genes by these methods. I want to compare my results with Two Sample t-test with bootstrapping, but I can not understand this method very well. For example in my table, there are 5 control, 5 treatment sample (columns) with 1000 genes (rows). Can we find differentially expressed genes by applying bootstrap two sample t-test for each gene? I scanned literature, but I could not find good solutions for gene expression analysis. I tried below codes as t-test process with "boot" package:

boot.tee <- function(data, i){
 data <- as.matrix(data)
  for (i in 1:data) {
    t.test(sample(data[i,1:5], 5, replace=T ),sample(data[i,6:10], 5, replace=T), paired = FALSE)$p.value
boot.out <- boot(data=LogT_matrix, statistic=boot.tee , R=10)

then I recieved an error :In 1:data : numerical expression has 10000 elements: only the first used.

In this page http://ww2.coastal.edu/kingw/statistics/R-tutorials/resample.html, there are some examples, but I want to obtain p values for all genes in my table such as toptable, toptags tables in EdgeR, limma packages. I can obtain standart t-test for my data, but I could not use it for bootstraping t-test. Can Bootstrap Statistics be applied to each gene? Thank you.


1 Answer 1


The error is in coding, which is technically off-topic here, but the question raises important statistical issues.

Repeated t-tests like this, even if done by bootstrapping, aren't very useful for gene expression data. Restricting analysis to one gene at a time makes the results very susceptible to the vagaries of the particular experiment. Just by random sampling, some genes might have very low standard errors among samples while others might have very high standard errors. In the first case you might find some "statistically significant" differential expression results of an insignificant magnitude; in the second, you might miss some truly important differences hidden in the excess noise. Also, you still have to take into account the testing on 1000 genes with some type of multiple-comparison correction.

The Bioconductor packages that you cite combine information across genes, treatments and samples to get pooled estimates of variance in expression values, with multiple-comparison correction built in. That approach typically gives much more useful estimates of differential expression. You can start to think about it similarly to the way that classical analysis of variance uses all the observations to get a pooled error estimate for statistical tests, although the pooling with those packages is more sophisticated and designed specifically for gene-expression analysis.

If you do want nevertheless to do bootstrapped-based multiple tests on individual genes, you must be careful with syntax. In for (i in 1:data), R interprets data to be all 10000 entries of your data matrix (1000 rows and 10 columns). You probably wanted to use 1:nrow(data). Also, I assume that you only asked for 10 bootstrap samples (R = 10) because you were still testing your code; many more would be highly advised, and particularly important if you want to get reliable estimates of things like 95% confidence intervals.


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.