1
$\begingroup$

I have 3 experimental conditions, each generate a series of values for some genes, like the following:

    c1  c2  c3

g1 n11 n12 n13
g2 n21 n22 n23 . . .

c1, c2, c3 are conditions; g1, g2, ... are gene names; n11, ... are the expression levels. c1 & c2 belong to one group, and c3 is the other group by itself.

The biological question is: which genes are differentially expressed?

Thank you!

$\endgroup$
2
  • $\begingroup$ I understand the title of the question. But, I have no idea about the detail of your question. If you can find a way to ask the question which does not require the reader to understand the biology and your notation, more people are going to be able to answer it. $\endgroup$
    – Tim
    Commented Feb 21, 2017 at 21:43
  • $\begingroup$ OK, let me try to rephrase. c1 & c2 are two samples which fall into the same category (control), c3 is another sample which falls into the 2nd category (treatment). nij is the gene expression of gene i in sample cj. What I want to do is to find out whether ni3 is statistically different from ni1 & ni2 for gene i. $\endgroup$ Commented Feb 22, 2017 at 13:54

3 Answers 3

1
$\begingroup$

This is a classical application of statistical biology where you have small sample size. Your example is even worse, you don't have any technical replicate!! Without any biological replicate, you can't model the volatility for each gene. This is a bit like doing t-test between two samples, each with a single respondent.

All modern statistical techniques rely on borrowing information from similar abundant genes. Although we don't have lots of samples, we have lots of genes.

Unfortunately, those packages also rely on replicates to estimate the negative binomial dispersion correct. Without replicates, you may:

  • Calculate the fold-ratio for each gene and use it to compare with your chosen threshold. One disadvantage is that you don't get a p-value.

  • Apply Fischer's exact test to each gene against all the remaining genes. This implicitly assumes most genes are not statistically differentiated.

enter image description here

The paper is here.

  • Follow the limma package described by @Gordon_Smyth.

EDITED

I apology I thought your experiment had no replicate. I misread the question.

EDITED

Just to prove Fischer's exact test on 2x2 table is a possibility. Statistical Methods for Detecting Differentially Abundant Features in Clinical Metagenomic Samples has:

enter image description here

$\endgroup$
6
  • $\begingroup$ Not sure about the "don't get a p-value" part of fold-ratio -- you can use a distribution on the fold-ratio to get the p-value or even just use rank of fold-ratio? The Fischer's exact test sounds interesting but as you said, it assumes most genes are not statistically differentiated, and it probably also assumes the genes follow IID? This could be a big assumption. Mr. Gordon Smith mentioned above about limma. This seems to be the most powerful test. Thank you so mcuh for your answer! $\endgroup$ Commented Feb 23, 2017 at 14:49
  • $\begingroup$ @MoushengXu using fold change is simplest but there is no distribution assumed and thus no p value. $\endgroup$
    – SmallChess
    Commented Feb 23, 2017 at 14:52
  • $\begingroup$ The assumption of most genes are not differentiated is a common assumption and work very well in practice. $\endgroup$
    – SmallChess
    Commented Feb 23, 2017 at 14:54
  • 1
    $\begingroup$ Yes, Fisher's test is a possibility, and it is implemented in the sage.test() function of the statmod package for this sort of data. But Fisher's test ignores biological variation between the replicates and so gives p-values that are much too small. People haven't done it like that for quite a few years. If the data are counts, then using a modern method like edgeR would be far better. $\endgroup$ Commented Feb 23, 2017 at 23:30
  • $\begingroup$ +1 @GordonSmyth Agreed. But what if the experiment has no biological replicate? The question doesn't say there're replicates. $\endgroup$
    – SmallChess
    Commented Feb 23, 2017 at 23:33
2
$\begingroup$

The experiment you have outlined is a minimal example of what is called the differential expression problem in genomic research. RNA samples are collected from two or more cell types or treatment conditions and gene expression profiling is conducted on each RNA sample. This provides an abundance (or expression) measure for each gene in each sample. The aim is to do a statistical test for each gene to find which genes show a change in mean expression level between the conditions. The challenge of course is that the sample sizes are tiny while the number of tests to be done is very large. Genewise t-tests would be on 1 df! In typical genomic applications, the number of rows (genes or genomic ids) is usually over 10,000 but can be much larger, even into the millions.

An enormous amount of attention has been devoted to this problem in the statistical bioinformatics literature with (I guess) probably more than 1000 methodological papers. The consensus is that univariate statistical tests perform very badly if they are simply applied genewise and one must instead use more advanced statistical methods that borrow information between genes while allowing for the possibility that genes may have different variances. One of the most popular methods is that of the limma package, which uses empirical Bayes moderated t-tests. (I'm the limma author, so I may be biased.)

Here is an example run using R and the limma package. I will assume that your expression levels are log2-expression values. If the entries are actually RNA-seq counts or something else, then a little pre-processing is required. First install the package:

source("http://www.bioconductor.org/biocLite.R")
biocLite("limma")
library(limma)

Generate some example data with 1000 genes. In this simulation, the first 5 genes are differentially expressed:

y <- matrix(rnorm(3000),1000,3)
rownames(y) <- paste0("Gene",1:1000)
y[1:5,3] <- y[1:5,3]+10

Now analyse with limma:

Group <- factor(c(1,1,2))
design <- model.matrix(~Group)
fit <- lmFit(y, design)
fit <- eBayes(fit)

Show the top 10 DE genes by p-value:

> topTable(fit)
        logFC AveExpr     t  P.Value adj.P.Val      B
Gene4    9.54   2.408  7.75 6.51e-08  5.91e-05  7.208
Gene2    9.22   2.870  7.46 1.24e-07  5.91e-05  6.715
Gene1    8.98   3.287  7.30 1.77e-07  5.91e-05  6.436
Gene3    8.70   3.823  7.05 3.16e-07  7.89e-05  5.986
Gene5    8.81   3.200  6.68 7.50e-07  1.50e-04  5.301
Gene414 -4.65  -0.988 -3.80 8.95e-04  1.32e-01 -0.586
Gene507  4.64  -0.241  3.79 9.21e-04  1.32e-01 -0.611
Gene566  4.05   0.702  3.18 4.11e-03  5.14e-01 -1.882
Gene628 -3.49  -0.133 -2.86 8.85e-03  9.84e-01 -2.530
Gene131 -3.26  -0.337 -2.67 1.36e-02  9.94e-01 -2.888

The t column is the moderated t-statistic. The adj.P.value column here gives an upper bound for the expected false discovery rate if that number of genes are chosen as DE. The B column is the log posterior odds of being DE. In this example, all the truly DE genes have positive log posterior odds.

$\endgroup$
5
  • $\begingroup$ Thank you so much! Very helpful. Then I assume I could use edgeR? It uses limma underneath, I think. $\endgroup$ Commented Feb 23, 2017 at 14:31
  • 1
    $\begingroup$ edgeR would be the ideal solution if your data is RNA-seq (although you haven't told us that). You should ask questions about limma and edgeR on the Bioconductor support site support.bioconductor.org . Cross Validated is not an appropriate forum to ask for help on these specialized tools. Indeed the whole topic raised by your question is better dealt with on that forum. $\endgroup$ Commented Feb 23, 2017 at 22:42
  • $\begingroup$ Without biological replicates (as stated in the question), edgeR could be doggy (read section 2.11 in the edgeR user guide). $\endgroup$
    – SmallChess
    Commented Feb 23, 2017 at 22:54
  • 1
    $\begingroup$ @Student T. The question states that there are two replicates in the first group. edgeR works with that quite happily. Even if there were no replicates, edgeR would be still more more reliable than Fisher's exact test. Fisher's test is equivalent to the limiting case in edgeR when the dispersion is set to zero. So Fisher's test is always known to be overly liberal. I wrote section 2.11 of the edgeR User's Guide, so I know what it says. $\endgroup$ Commented Feb 23, 2017 at 23:48
  • $\begingroup$ @GordonSmyth You asked the right question -- I should not have said it's gene expression level -- it's ChIP-seq signal. Although it's not RNA-seq, can I still assume the counts fall into a negative binomial distribution? $\endgroup$ Commented Feb 24, 2017 at 15:29
0
$\begingroup$

Two options: 1. Pool together all of the data from c1 and c2, and perform an independent samples t-test. However, this only tells you that c3 is different from the joint c1 and c2. 2. Use Dunnett's test, but, what you are calling the 'treatment' will be the 'control' in the terminology of this test. This will test whether c1 is different from c3 and whether c2 is different from c3.

$\endgroup$
1
  • $\begingroup$ Sounds an appropriate test. How do you think about "borrowing" information from other genes to increase the power as mentioned in Gordon Smith's & Student T's replies? $\endgroup$ Commented Feb 23, 2017 at 14:38

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.