Question: How many tests does it take to warrant a microarray/RNA-seq approach to differential gene expression analysis?
gravatar for lucas.james.wager
4.0 years ago by
lucas.james.wager10 wrote:

Background Information:

I have 53 genes of interest. I have two conditions Tumour and Normal samples at 24, 48 and 72 hours from three mice (i.e. 3xTumour and 3x Normal for each timepoint). The Tumours and Normal skin are sampled repeatedly (i.e. each mouse has Tumor and Normal samples taken at each time point). I have measured the expression of all my genes of interest (53) in each sample using qRT-PCR.

If these data were RNA-seq or microarray I would turn to something like DEseq2 to perform a differential expression analysis. If I only had up to 10 genes of interest I would have no concern to use a repeated measures ANOVA or linear models to investigate differential expression between my two conditions.

However I feel like I am somewhere in between each method, as correcting for multiple testing for 53 ANOVAs is diminishing my Pvalues.

Here are 20 of the genes of interest over the 24, 48, and 72 hours (for now ignore the green 0H data points). I am interested in the difference between conditions.

I am aware that N=3 is a very small sample, but this is the nature of biological sciences (especially during a PhD). With this in mind humour me in your answer to my questions.


After how many measurements does one decide that they are performing enough tests to warrant a software package such as DEseq2 (i.e. a microarray/RNA-seq approach) to investigate differential expression and how many tests are not enough?

Are there tests or packages that better deal with low to moderate numbers of statistical tests?

In the case that I am not raising a valid concern (quite possible), what would be the most powerful approach to discerning any difference in gene expression that exists between tumour and normal tissue, at any point in time, in my study?

Thanks in advance


R gene • 1.9k views
ADD COMMENTlink modified 4.0 years ago by Devon Ryan98k • written 4.0 years ago by lucas.james.wager10

As a remark: the assumption in tools such as DESeq2 and edgeR is that most genes in the dataset are not differentially expressed. I don't think that is the case for your selection? Did you quantify housekeeping genes using qRT-PCR?

ADD REPLYlink written 4.0 years ago by WouterDeCoster45k

I did quantify housekeeping genes, however these are not shown in my selection, but I use them to create these graphs. I have plotted the DELTA Ct of my genes of interest calculated by subtracting the average Ct of three housekeeping genes at each time point for each condition, converting to fold change and then normalising to the 0H point. To comment on your remark regarding the presence of many changes in my genes of interest, this difference may not be as prominent as it seems upon first glance. I have allowed a free y scale to exemplify the patterns in the genes individually (patterns I would like to see). However these patterns diminish when the Y scale is kept constant. I will try to test for differential expression with DEseq2. As it will consider all the genes at once perhaps the changes that seem apparent here will not have a detrimental effect on the assumption of consistency for the majority of genes. I have attached the same genes graphed with a consistent Y scale.

Same genes with a consistent Y scale

ADD REPLYlink written 4.0 years ago by lucas.james.wager10

Did you try something like bioconductor HTqPCR? I never used it myself, but if I had data like yours I would probably try to find something in bioconductor first.

ADD REPLYlink written 4.0 years ago by Benn8.1k

No I haven't yet. Now it's installed and I'll give it a go. Thanks for the suggestion.

ADD REPLYlink written 4.0 years ago by lucas.james.wager10

So I tried HTqPCR, but found that I was unsure of what was happening behind the scenes, so I resorted to limma - which I understand was implemented in HTqPCR anyway.

I have come to a point of confusion. There is an example of applying limma to differentiate between the expression of 3 replicates forming two groups (dummy data) given by Gordon Smyth in the limma package vignette (bottom of page). In this example there is a comment "Would like to consider original two estimates plus difference between first 3 and last 3 arrays" that I do not understand. This comment and following method of constructing a contrast matrix seems at odds with another BiostarUser for investigating the gene DE between groups that do not "consider original two estimates" and examples I have followed in the limma user guide (section 9.2). I'm sure this comes down to my lack of understanding when it comes to contrast matrices. However I have been unable to resolve this myself.

Help with understanding the difference between the two methods is greatly appreciated.

Smyth method

#  Simulate gene expression data: 6 microarrays and 100 genes
#  with one gene differentially expressed in first 3 arrays
M <- matrix(rnorm(100*6,sd=0.3),100,6)
M[1,1:3] <- M[1,1:3] + 2
#  Design matrix corresponds to oneway layout, columns are orthogonal
design <- cbind(First3Arrays=c(1,1,1,0,0,0),Last3Arrays=c(0,0,0,1,1,1))
fit <- lmFit(M,design=design)
#  Would like to consider original two estimates plus difference between first 3 and last 3 arrays
contrast.matrix <- cbind(First3=c(1,0),Last3=c(0,1),"Last3-First3"=c(-1,1))
fit2 <-,contrast.matrix)
fit2 <- eBayes(fit2)
#  Large values of eb$t indicate differential expression
results <- classifyTestsF(fit2)


Biostar user method

celFiles <- list.celfiles(path_to_cel_files,full.names=TRUE)
affyExpressionFS <- read.celfiles(celFiles)
eset <- rma(affyExpressionFS, background = TRUE, normalize = TRUE)
stress <- c(1,1,1,0,0,0) # first three replicates (case)
control <- c(0,0,0,1,1,1) # last three replicates (control)
design <- cbind(stress, control) # first term is numerator, second denominator of the ratio
fit <- lmFit(eset, design)
cont.matrix <- makeContrasts(STvsCO=stress-control, levels=design)
fit2 <-, cont.matrix)
fit2 <- eBayes(fit2)
tab <- topTable(fit2, adjust = "BH", confin =T , number = 1000,"logFC")
ADD REPLYlink modified 4.0 years ago • written 4.0 years ago by lucas.james.wager10

I don't exactly understand your question, but I think you're better off asking Gordon Smyth himself for help. He's active on and usually answers questions about limma within a day. Be sure to add limma as a tag.

ADD REPLYlink written 4.0 years ago by Benn8.1k
gravatar for Devon Ryan
4.0 years ago by
Devon Ryan98k
Freiburg, Germany
Devon Ryan98k wrote:
  1. You start gaining power by using limma/DESeq2/et al. once you have more than one gene.
  2. Yes, limma/DESeq2/etc. (for RNAseq, for qPCR there's a different bioconductor package that can probably do most of the analysis for you).
  3. What WouterDeCoster said.
ADD COMMENTlink modified 4.0 years ago • written 4.0 years ago by Devon Ryan98k

Can we really use DESeq2 with continuous data (vs integer count) ? It's a genuine question, I really don't know.

ADD REPLYlink written 4.0 years ago by Carlo Yague5.5k

It's not that it's been designed for qPCR, but tximport for kallisto/salmon pseudocounts isn't exclusively for integers.

ADD REPLYlink written 4.0 years ago by WouterDeCoster45k

It appears that DESeq2 still only accepts integer input (in fact, the DESeqDatasetFromTxImport function rounds them internally ( ).

ADD REPLYlink written 4.0 years ago by Devon Ryan98k

Thank you for the response. The idea of rounding these values is of concern, as I am dealing with such small changes. After internal rounding I am likely to end up with a great deal of '1's.

Would a consistent multiplication of the inputs before analysing them with DEseq be a bad idea? For example, if raw inputs to DEseq, 1.48 and 0.84, were tested without multiplication they would both be rounded to 1 and show no difference, but by multiplying by 10 they would round to 15 and 8, making the tests more sensitive(?).

I guess the question would then be one of the relevance of any changes found. I'm interested in your opinion on this.

ADD REPLYlink written 4.0 years ago by lucas.james.wager10

I'm not sure how nicely that'd play with a negative binomial fit. In any case, you could always use edgeR or limma/voom instead.

ADD REPLYlink written 4.0 years ago by Devon Ryan98k

Using edgeR, DESeq or DESeq (or any NB-based software) for PCR data is not at all appropriate, and for reasons that run much deeper than whether the values are integer or not. Use limma or HTqPCR (a limma wrapper) instead (and ordinary limma, not voom).

ADD REPLYlink modified 2.3 years ago • written 2.3 years ago by Gordon Smyth2.1k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2017 users visited in the last hour