[BioC] Differential gene expression: EdgeR / DESeq and identifying noise/outliers

[Gordon K Smyth]

Hi Michael,

Here's a link to a previous reply that I made to a similar question last 
month:

  https://www.stat.math.ethz.ch/pipermail/bioconductor/2012-May/045483.html

The short answer is to set the argument prior.n of the 
estimateTagwiseDisp() function in edgeR to a smaller value.  The default 
in edgeR is to use a largish value for the prior df, which means that 
edgeR squeezes the tagwise dispersions strongly towards the global value, 
meaning that it isn't able to adapt sufficiently to individual genes with 
outliers such as yours.

Your data has 14 libraries and two groups, so there are 12 residual 
degrees of freedom for each gene.  The prior degrees of freedom are set to 
20 by default, so the prior number of observations defaults to prior.n = 
20/12 = 1.67 for your data.  Try instead a smaller value like prior.n = 
6/12.  The smaller prior.n is, the more edgeR will de-prioritize those 
hyper variable genes.

It would preferable if edgeR did this for you, adapting to the 
characteristics of your data automatically.  A new version of edgeR should 
be able to do that in a few months.

Best wishes
Gordon

> Date: Thu, 28 Jun 2012 17:43:09 -0500
> From: "Michael Salbaum" <Michael.Salbaum at pbrc.edu>
> To: <bioconductor at r-project.org>
> Subject: [BioC] Differential gene expression: EdgeR / DESeq and
> 	identifying	noise/outliers
>
> 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

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:5}}