Hi everyone, 

I am working on a differential gene expression paradigm; n=7, 3’-expression tag sequencing, AB SOLiD5500XL. We look for expression of ~26,000 RefSeq genes; as input, I’m using a bottom-shaved counts table (~16,500 genes), i.e. rows with a total count sum of 12 or less have been trimmed off. Using edgeR, I find 2569 genes (1224 up / 1345 down) to be differentially expressed at FDR better than 0.05.
I’ve noticed that in this list are many genes where the differential expression call appears to be driven by one ‘outlier’ on one side of the paradigm, for instance:
Nodal	WT: 7 7 5 4 19 6 1 Het: 320 2 16 8 6 1 13 logFC: 2.65 FDR: 0.00814198

Of course, this is not desirable for follow-up studies, and I’m wondering whether there’s a way to filter out such situations.

Right now, I’ve resorted to looking at the coefficient of variation calculated from normalized data as a potential discriminator.

I’ve also used DESeq on the same count data set, which identifies 1569 genes (719 up / 850 down) at padj<0.05. Of these 1569, 1544 are also called by edgeR, so, excellent agreement. Relaxing the cutoff to padj<0.1 gives another 761 genes, of which 585 make the FDR<0.05 cutoff in edgeR.

Plotting –log10(FDR) against the coefficient of variation shows that ~90% the genes called by both DESeq and edgeR (at padj<0.05) have a CV of 0.3 or better. The same plot for the genes called only by edgeR (at padj<0.05) seems to show that this group contains many genes at higher CV – I suspect these are the outlier-driven ones mentioned above.

I’m perfectly happy with the outcome (…and with the procedure of using two programs and continue with the intersect between the two…), but those single outlier-driven genes were irksome – particularly so because they were of a nature that got us excited, in biological terms ?. So, to avoid the letdown, I’d appreciate advice how not to get those in the first place.

And I apologize for the long post, but as a newbie, I figured I better include what info I have.

I’m not sure this is proper etiquette here, but here is a link to the plots:
http://inlinethumb04.webshots.com/51523/2373514640050256648S600x600Q85.jpg

R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows"
‘edgeR’ version 2.6.7

library(edgeR)
Loading required package: limma
x <- read.delim("/Users/jms/Desktop/SAGE3/SAGE3F.txt",row.names="Gene")
> head(x)
              WT_1 WT_2 WT_3 WT_4 WT_5 WT_6 WT_7 HET_1 HET_2 HET_3 HET_4 HET_5 HET_6 HET_7
0610005C13Rik    5    1    6    0    5    2   18    34     7    13    18    28     1    13
> group <- factor(c(1,1,1,1,1,1,1,2,2,2,2,2,2,2))
> y <- DGEList(counts=x,group=group)
Calculating library sizes from column totals.
> y <- estimateCommonDisp(y, verbose=TRUE)
Disp = 0.06254 , BCV = 0.2501 
> y <- estimateTagwiseDisp(y)
> et <- exactTest(y)
> topTags(et)
> summary(de <- decideTestsDGE(et, p=0.05, adjust="BH"))
   [,1] 
-1  1345
0  14155
1   1224



‘DESeq’ version 1.8.3
library(DESeq)
Loading required package: Biobase
Loading required package: BiocGenerics
Attaching package: ‘BiocGenerics’
The following object(s) are masked from ‘package:stats’:
    xtabs	
The following object(s) are masked from ‘package:base’:
    anyDuplicated, cbind, colnames, duplicated, eval, Filter, Find, get, intersect, lapply, Map, mapply,
    mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rbind, Reduce, rep.int, rownames, sapply,
    setdiff, table, tapply, union, unique
Welcome to Bioconductor
    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.
Loading required package: locfit
locfit 1.5-8 	 2012-04-25
> f = "/Users/jms/Desktop/SAGE3/SAGE3F.txt"
> countsTable <- read.table( f, header=TRUE, row.names=1, stringsAsFactors=TRUE )
> head(countsTable)
              WT_1 WT_2 WT_3 WT_4 WT_5 WT_6 WT_7 HET_1 HET_2 HET_3 HET_4 HET_5 HET_6 HET_7
0610005C13Rik    5    1    6    0    5    2   18    34     7    13    18    28     1    13
> conds <- c( "WT", "WT", "WT", "WT", "WT", "WT", "WT", "HET", "HET", "HET", "HET", "HET", "HET", "HET" )
> cds <- newCountDataSet( countsTable, conds )
> cds <- estimateSizeFactors( cds )
> sizeFactors( cds )
     WT_1      WT_2      WT_3      WT_4      WT_5      WT_6      WT_7     HET_1     HET_2     HET_3     HET_4     HET_5 
1.3598468 0.9759626 1.0788248 1.0080744 1.0097703 0.5502473 0.8747237 1.0877753 0.5544828 0.9594141 1.6036028 1.4868868 
    HET_6     HET_7 
0.8422686 1.5328149 
> cds <- estimateDispersions( cds )
> res <- nbinomTest( cds, "WT", "HET" )
> write.table (res, sep = "\t", file = "/Users/jms/Desktop/SAGE3/SAGE3F_DE.txt", col.names = NA)
> write.table (counts( cds, normalized=TRUE ), sep = "\t", file = "/Users/jms/Desktop/SAGE3/SAGE3F_NORM.txt", col.names = NA)
	
Cheers, michael


J. Michael Salbaum, Ph.D.
Associate Professor
Pennington Biomedical Research Center
Louisiana State University System
6400 Perkins Road
Baton Rouge, LA 70808

(225) 763-2782

	[[alternative HTML version deleted]]

