Differential Gene Expression

Question

I have a matrix two gene expression datasets, one from 10 normal individuals and ten from diseased individuals, which are normalised log fold change values--, i.e., already processed.

My question is: how do I obtain a list of differentially expressed genes? That is the up and down expressed genes. I was thinking of getting means for each group and then sorting the average fold change values would work, but I am not sure...

Thanks.

Answer 1

It hard to figure out what sort of data you have. qPCR? Let me give you an example with tools involved on how I get differential genes for RNAseq data, that is whole transcriptom and for mouse that transcriptom that's about 23000 genes.

• Step 1. Map raw reads to the reference genome.
• Step 2. Count number of reads within the feature. Feature can be several things, but lets just take protein coding feature (feature == gene).
• Step 3. Do differential expressions analysis

From step 2 you get a matrix of counts.

      sampleA_rep1 sampleA_rep2 sampleB_rep1 sampleB_rep2
gene1      2             5        20             23 
gene2      10            12       50             50
gene3      20            23       4              6
..
geneN 

The packages that I use are limma and/or edgeR. You can find them on bioconductor Both packages take raw counts and give you topTable of DEG, default is top 100 genes.

I recommend reading edgeR or limma manuals, there come with great explanations as well as examples. Also A scaling normalization method for differential expression analysis of RNA-seq data is original paper that explains TMM normalisation (mention a bit more in depth below)

Two things; one you need to filter genes with low counts, convert to cpm() counts per million before filtering as this takes library size (depth of sequencing) into consideration. two calcNormFactors() will do TMM

The hardest thing (I understand) in RNAseq is to get total RNA output for each sample. Because a treatment will effect the number of transcripts expressed in that sample, increasing or decreasing amount of transcript, which give you unbalanced counts. You won't know just by looking at raw counts if geneA has less counts because it's been down regulated OR because that sample after treatment has more transcripts present in the cell and geneA simply has less coverage, due to other transcripts also getting reads now. I can't really explain that too technically because I don't understand this myself, but TMM normalisation takes care of this and gives you an estimate of degree of variability.

I know that I'm not giving that much info on actual differential testing, which was your question. edgeR has exactTest() which is Exact Tests for Differences between Two Groups of Negative-Binomial Counts , which might help you, but you might have normalise your counts differently as mentioned above. The general idea (I understand) is you need to fit a model (lmFit()) and then do the test eBayes()

Other

• Degust this is very nice tool that uses limma and edgeR in the backend, but simplifies a lot of things as well as gives you nice interface to work with your data. You upload raw counts! in the format shown above
• DESeq one of the other rather popular R packages that can also do differential gene testing

Answer 2

There is software for this, but if you want a very primitive way to analyze it, do a t-test (which is basically looking at both the difference in means and the STDEV of each group)