Thanks Gordon.  Would your advice about X still hold if the smallest group
size is 5 out of 500 total samples (i.e. an experimental design profiling
100 treatments with 5 replicates each)?  5 out of 5000 samples?  It seems
to me that the available df should somehow impact your choice of X, and
that for very large df, you'd wish to raise X somewhat.

What I'm trying to get at is that it seems, at least to my intuition, that
the desire to keep "at least X nonzero" stems from a possibility of a
transcript being entirely "off" (zero counts) in all but the single,
smallest group where it is found "on" (non-zero counts) -- but the models
we use don't really test/allow for such "binary" switching behavior (i.e.
infinite fold changes); so instead we're really just trying to weed out
those transcripts that tend to have zero counts more often than others
(which may not be due to expression at all, but other artifacts of the
experiment, including vagaries of read mapping), that negatively impact the
statistics of what we can model (non-infinite fold changes).

Thus, we have a desire to come up with a metric or decision criterion for
which the "rule of thumb" is a good proxy in small/modest-sized experiments
but which can also be extended in larger many-factorial experiments ...
i.e. something akin to getPriorN().  I imagine some tuning strategy (much
like lasso/ridge) where the decrease in error associated with more
aggressive filtering is offset by a loss in power, until a balance/crossing
point is found; in an edgeR (not voom) scenario, this might be some aspect
of the BCV distribution/loess fitt vs. some multivariate F-statistic of the
experimental design?

Thanks to all for brainstorming this,

-Aaron

On Mon, Sep 10, 2012 at 7:59 PM, Gordon K Smyth <smyth@wehi.edu.au> wrote:

> Dear Mark,
>
> I think that voom() should be pretty tolerant of the amount of filtering
> that is done, so you can feel free to be more inclusive.
>
> Note that our recommended filtering is
>
>   keep <- rowSums(cpm(dge) > k) >= X
>
> where X is the sample size of the smallest group size.  Since X is usually
> smaller than half the number of arrays, our recommended filtering is
> usually more inclusive than the filter you give.
>
> You are also free to vary k, depending on your sequencing depth.  The idea
> is to filter low counts.
>
> Best wishes
> Gordon
>
> -------------- original message -------------
> [BioC] A metric to determine best filtration in the limma package
> Aaron Mackey amackey at virginia.edu
> Mon Sep 10 16:27:21 CEST 2012
>
> Hello Bioconductor Gurus!
>
> (I apologize if this goes through more than once)
>
> We are currently using limma (through the voom() function) to analyze
> RNA-seq data, represented as RSEM counts. We currently have 246 samples
> (including replicates) and our design matrix has 65 columns.
>
> My question is in regard to how much we should be filtering our data
> before running it through the analysis pipeline. Our current approach is to
> look for a CPM of greater than 2 in at least half of the samples. The code
> is:
>
> keep <- rowSums(cpm(dge) > 2) >= round(ncol(dge)/2)
>
> This brings down our transcript count from 73,761 to less than 20,000.
> While we do see groupings and batch effects we expect to see in the MDS
> plots, we are afraid we might be filtering too severely.
>
> So finally my question: What is a good metric for determining how well we
> have filtered the data?
>
> Thank you,
> Mark Lawson, PhD
>
> ______________________________**______________________________**__________
> The information in this email is confidential and inte...{{dropped:10}}

