Hi Wolfgang,
Thanks for the apologies.

In terms of the data quality, I am pretty confident in  the sample
collection.As a background,  this pair of samples (FC v CC)  was extracted
from a much bigger timecourse done by microarrays. These are mice, the
tissue is brain and I am measuring response to learning, there's a
timecourse of controls and a timecourse of learning animals. The PCA of the
big matrix of 72 samples does indeed extract sources of variability in gene
expression in the brain that are not my treatment. Most of those can be
accounted for by reasons that I somewhat expected and thus were randomized
beforehand.

The pair of samples I used for RNAseq were deliberately selected because it
is the time point when learning has its bigger effect and it is the main
source of global variability in gene expression. These samples did not have
big contributions in any of the PC that were not explained by the
treatment.

I have now done some further data exploration on the RNASeq data following
your suggestions. Using MDS I get a good separation between treatment and
control  data in the second dimension while the samples align in the first
dimension, (one sample being a huge outlier, which I removed from further
analysis) the separation gets better when you filter out low counts. The MA
plot looks normal as far as I can tell, as well as the distribution of
p-values.

The DE gene list I now get from DESeq is a subset of the microarray list
almost completely contained in it with a very good rate of true positives,
the main problem is the false negatives and the fact that there seems to be
a problem detecting down-regulation. So the data now looks like this (at
FDR <0.1) :


Up                                                  Down

RNAseq (DESeq)

90 (69)

2

Microarray

104

79

The numbers for RNASeq are Refseq gene models with the number of unique
gene names in parenthesis, while it is probesets for the microarray. The
main difference in the analysis is the multiple testing correction.

The reason we settled on using locfdr instead of BH in our microarray
analysis was the fact that  >90% of the genes don't change with our
treatment, thus estimating an empirical null makes a lot of sense and
should in theory be more accurate, that may not be true in other
situations. Judging by independent qPCR, the locfdr procedure returns more
true positives and the GO enrichment analysis makes more sense.

I am going to feed the vector of uncorrected p-values from DESeq to locfdr
and see what happens.
I  am also going to try generating my count matrix using HTSeq and run the
analysis again, I want to use HTSeq to generate counts for DEXSeq as well.

In any case, I have not been able to find a study in which microarrays and
RNAseq are compared head to head in multiple biological replicates of the
same samples. I am not that familiar with the RNASeq literature but,
can it be possible that when dealing with biological (not technical ) noise
at the gene level it is still better to use microarrays?

Lucia








On Tue, May 21, 2013 at 6:49 PM, Wolfgang Huber <whuber@embl.de> wrote:

> Dear Lucia
>
> I just googled one of the words I had used and realised it has a much more
> sinister meaning than I ever intended. My apologies for that very
> unfortunate choice. A spectacularly incompetent attempt at using
> informal/slang language in a place where it does not belong.
>
> You reported a large difference between the results of the 'locfdr' and
> 'BH' methods for multiple testing correction on the microarray data. This
> could happen (among a few other reasons) if there are problems with data
> quality, such as batch effects, that affect the distribution of the test
> statistic for non differentially expressed genes. How does the histogram of
> unadjusted p-values look like? For the BH-method, it should look like a
> mixture of a uniform between 0 and 1 (for the true nulls) and a peak at the
> left. See also, e.g., Fig 3c in
> http://openi.nlm.nih.gov/detailedresult.php?img=2865860_btq118f3&req=4
>
> The bias towards up-regulation in the RNA-Seq data is curious. It does not
> seem to be caused by unequal library sizes. One way to get to the bottom of
> this could be the MA-plots - how do they look like?
>
> Best wishes
>         Wolfgang
>
> On 21 May 2013, at 23:19, Lucia Peixoto <luciap@iscb.org> wrote:
>
> > Dear Simon,
> >
> > My apologies for not being concise in explaining the problem. I am pretty
> > new to the list and I can see that it can be frustrating to try to help
> > someone that is not giving all the details you need to do so. I
> appreciate
> > your apology, since some of the initial responses to my question had a
> very
> > unfortunate choice of words. Going through the suggestions of the thread
> > has been very helpful. So here is a more thorough description of what I
> > did, as well as some observations made after taking into account
> > suggestions:
>



-- 
Lucia Peixoto PhD
Postdoctoral Research Fellow
Laboratory of Dr. Ted Abel
Department of Biology
School of Arts and Sciences
University of Pennsylvania

"Think boldly, don't be afraid of making mistakes, don't miss small
details, keep your eyes open, and be modest in everything except your
aims."
Albert Szent-Gyorgyi

	[[alternative HTML version deleted]]

