[Bioc] RNAseq less sensitive than microarrays? Is it a statistical issue?
2
0
Entering edit mode
Lucia Peixoto ▴ 330
@lucia-peixoto-4203
Last seen 10.2 years ago
Dear All, I have a dataset for which I have two conditions. I have 9 replicates per group for microarrays, 5 per group for RNAseq (which are a subset of the RNA samples used in the microarrays, couldn't sequence all 9), and 8 per group for qPCR (which is an independent set of experiments). Each n is an independent mouse, in and independent day from and independent experiment, so that one experiment with yield n=1 for each of the groups. The correlation between control and treatments within the same day is not better than across days, however. Theoretically they all measure the same biological phenomenon, which is gene expression changes, so I have been doing some comparisons between them to try to get at the truth of what is really being differentially expressed. In particular I have focused in the 5 samples in each of the three groups in which the only difference is whether the RNA was hybridized by microarray or sequenced. To my surprise the gene lists obtained from analyzing differential expression using RNASeq (using either edgeR or DESeq) is considerably smaller than the one obtained from microarray analysis (using locfdr on pairwise t-statistics) at the same FDR. The RNASeq list is included in the microarray list, but there are several differences I have validated by qPCR that the RNASeq analysis is not able to detect at a reasonable FDR. Moreover, there seems to be an unusual bias towards not being able to detect down-regulated genes. I am a little bit puzzled by this, since one of the reasons we are sequencing is that it is supposed to have a better dynamic range. These are the same RNA samples so this apparent lack of sensitivity has to be related to either library prep or statistical analysis. So these are my questions: - can the inability to distinguish down-regulated genes be related to filtering low count reads? (in order to get good separation between groups in an MDS plot I need to filter cpm >0.1) - Is it possible that I need more coverage to improve sensitivity? I am currently sequencing at 50X pair end, that seemed enough. Is there any published study looking at RNASeq sequencing depth and sensitivity in human or mouse genomes? - Are the multiple testing corrections applied in EdgeR and DESeq too stringent thus rendering the overall analysis less sensitive? For the record my count matrices are of counts of transcripts, averaging counts over all exons from the same gene model for all RefSeq genes. I did this because the microarray data is per transcript. In log scale I have on average 0.7 R2 correlation between microarray intensity and RPKM from the same sample. Thanks for the insight! Lucia -- 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]]
Sequencing RNASeq Microarray qPCR Coverage edgeR DESeq Sequencing RNASeq Microarray qPCR • 3.0k views
ADD COMMENT •
0
Entering edit mode
@steve-lianoglou-2771
Last seen 21 months ago
United States
Hi, Some quick comments: > - can the inability to distinguish down-regulated genes be related to > filtering low count reads? (in order to get good separation between groups > in an MDS plot I need to filter cpm >0.1) It's not clear how you are filtering. If the gene never has a cpm > 0.1 in any sample, then you'll have no hope of detecting differential expression anywhere. If it is 0.1 in some and 5 cpm in others, then that's more like it. Be sure that your filtering strategy keeps the latter. You should be filtering on the number of samples that detect the gene at a minimum threshold and you should not be using the "label" of the samples (which sample/treatment) for counting. > - Is it possible that I need more coverage to improve sensitivity? I am > currently sequencing at 50X pair end, that seemed enough. Is there any > published study looking at RNASeq sequencing depth and sensitivity in human > or mouse genomes? The important number you need to look at is the number of reads you are getting per sample. > - Are the multiple testing corrections applied in EdgeR and DESeq too > stringent thus rendering the overall analysis less sensitive? Unlikely. > For the record my count matrices are of counts of transcripts, averaging > counts over all exons from the same gene model for all RefSeq genes. I did > this because the microarray data is per transcript. In log scale I have on > average 0.7 R2 correlation between microarray intensity and RPKM from the > same sample. It's not clear to me how you are doing the counting, actually, but "averaging" anything suggests to me that you are doing it incorrectly. Look at the vignettes from a variety of packages (edgeR, DESeq2, easRNAseq, QuasR) for suggested counting methods. -steve -- Steve Lianoglou Computational Biologist Department of Bioinformatics and Computational Biology Genentech
ADD COMMENT •
0
Entering edit mode
> Look at the vignettes from a variety of packages (edgeR, DESeq2, > easRNAseq, QuasR) for suggested counting methods. Or HTSeqGenie, even ... Point is that there are many existing ways to counts reads per gene (or transcript) and there is likely not a need for you to invent another one. I believe HTSeqGenie gives you the results for many different ways (including at the transcript level) (see page 7 of their vignette), so it's especially useful. -steve -- Steve Lianoglou Computational Biologist Department of Bioinformatics and Computational Biology Genentech
ADD REPLY •
0
Entering edit mode
@wolfgang-huber-3550
Last seen 3 months ago
EMBL European Molecular Biology Laborat…
Dear Lucia there are many reasons why you might see less differentially expressed genes from your RNA-Seq data. Firstly - the post doesn't say what the number of fragments per sample is ("50X pair end" is not a useful statement for RNA), and it might indeed just be too low. - the data quality of the microarray libraries might happen to be better than the RNA-Seq library in your experiment (e.g. do you see many PCR duplicates in the latter?) - the asymmetry in power to detect up vs down regulated genes might be a problem of unequal library sizes. > For the record my count matrices are of counts of transcripts, averaging > counts over all exons from the same gene model for all RefSeq genes. I'm afraid this is not correct. The documentation of the DESeq and edgeR packages explicitly states [1] that you need to use the sum of the counts mapping to the gene (or transcript), not some average of per-exon values. It could well be that this is the major reason for your getting retarded results. As far as I can see, differences such as you report have nothing to do with independent filtering or multiple testing adjustment. The former, if done correctly, will increase, not decrease the number of hits; the latter is independent of technology (i.e. the same as with microarrays). Hope this helps - Wolfgang [1] For the record: http://bioconductor.org/packages/release/bioc/vign ettes/DESeq/inst/doc/DESeq.pdf --> Section 1.1 on page 2. On May 15, 2013, at 10:01 pm, Lucia Peixoto <luciap at="" iscb.org=""> wrote: > Dear All, > > I have a dataset for which I have two conditions. I have 9 replicates per > group for microarrays, 5 per group for RNAseq (which are a subset of the > RNA samples used in the microarrays, couldn't sequence all 9), and 8 per > group for qPCR (which is an independent set of experiments). > Each n is an independent mouse, in and independent day from and independent > experiment, so that one experiment with yield n=1 for each of the groups. > The correlation between control and treatments within the same day is not > better than across days, however. > > Theoretically they all measure the same biological phenomenon, which is > gene expression changes, so I have been doing some comparisons between them > to try to get at the truth of what is really being differentially > expressed. In particular I have focused in the 5 samples in each of the > three groups in which the only difference is whether the RNA was hybridized > by microarray or sequenced. > > To my surprise the gene lists obtained from analyzing differential > expression using RNASeq (using either edgeR or DESeq) is considerably > smaller than the one obtained from microarray analysis (using locfdr on > pairwise t-statistics) at the same FDR. The RNASeq list is included in the > microarray list, but there are several differences I have validated by qPCR > that the RNASeq analysis is not able to detect at a reasonable FDR. > Moreover, there seems to be an unusual bias towards not being able to > detect down-regulated genes. I am a little bit puzzled by this, since one > of the reasons we are sequencing is that it is supposed to have a better > dynamic range. > > These are the same RNA samples so this apparent lack of sensitivity has to > be related to either library prep or statistical analysis. So these are my > questions: > > - can the inability to distinguish down-regulated genes be related to > filtering low count reads? (in order to get good separation between groups > in an MDS plot I need to filter cpm >0.1) > - Is it possible that I need more coverage to improve sensitivity? I am > currently sequencing at 50X pair end, that seemed enough. Is there any > published study looking at RNASeq sequencing depth and sensitivity in human > or mouse genomes? > - Are the multiple testing corrections applied in EdgeR and DESeq too > stringent thus rendering the overall analysis less sensitive? > > For the record my count matrices are of counts of transcripts, averaging > counts over all exons from the same gene model for all RefSeq genes. I did > this because the microarray data is per transcript. In log scale I have on > average 0.7 R2 correlation between microarray intensity and RPKM from the > same sample. > > Thanks for the insight! > > Lucia > > > > > > > > > -- > 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]] > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
ADD COMMENT •
0
Entering edit mode
Thanks everyone for the helpful input My counts are in fact total counts per transcript and not averages, my mistake in the original post I have an average of 70 million read pairs of 100bp, with very even number of reads for all my samples I am curious regarding how many clonal reads is too many though. I always see quite a bit of it on my samples. However, I have heard mixed advice regarding filtering duplicates in RNAseq, since with high coverage you expect some of the duplication to actually reflect abundance. In my samples filtering duplicates reduces the dynamic range dramatically and for highly expressed transcripts a high proportion of the reads that map to them have 2 or more copies Thanks very much for all the help and advice Lucia Sent from my iPhone On May 15, 2013, at 4:33 PM, Wolfgang Huber <whuber at="" embl.de=""> wrote: > Dear Lucia > > there are many reasons why you might see less differentially expressed genes from your RNA-Seq data. Firstly > - the post doesn't say what the number of fragments per sample is ("50X pair end" is not a useful statement for RNA), and it might indeed just be too low. > - the data quality of the microarray libraries might happen to be better than the RNA-Seq library in your experiment (e.g. do you see many PCR duplicates in the latter?) > - the asymmetry in power to detect up vs down regulated genes might be a problem of unequal library sizes. > >> For the record my count matrices are of counts of transcripts, averaging >> counts over all exons from the same gene model for all RefSeq genes. > > I'm afraid this is not correct. The documentation of the DESeq and edgeR packages explicitly states [1] that you need to use the sum of the counts mapping to the gene (or transcript), not some average of per-exon values. It could well be that this is the major reason for your getting retarded results. > > As far as I can see, differences such as you report have nothing to do with independent filtering or multiple testing adjustment. The former, if done correctly, will increase, not decrease the number of hits; the latter is independent of technology (i.e. the same as with microarrays). > > Hope this helps - > Wolfgang > > [1] For the record: http://bioconductor.org/packages/release/bioc/vi gnettes/DESeq/inst/doc/DESeq.pdf --> Section 1.1 on page 2. > > > On May 15, 2013, at 10:01 pm, Lucia Peixoto <luciap at="" iscb.org=""> wrote: > >> Dear All, >> >> I have a dataset for which I have two conditions. I have 9 replicates per >> group for microarrays, 5 per group for RNAseq (which are a subset of the >> RNA samples used in the microarrays, couldn't sequence all 9), and 8 per >> group for qPCR (which is an independent set of experiments). >> Each n is an independent mouse, in and independent day from and independent >> experiment, so that one experiment with yield n=1 for each of the groups. >> The correlation between control and treatments within the same day is not >> better than across days, however. >> >> Theoretically they all measure the same biological phenomenon, which is >> gene expression changes, so I have been doing some comparisons between them >> to try to get at the truth of what is really being differentially >> expressed. In particular I have focused in the 5 samples in each of the >> three groups in which the only difference is whether the RNA was hybridized >> by microarray or sequenced. >> >> To my surprise the gene lists obtained from analyzing differential >> expression using RNASeq (using either edgeR or DESeq) is considerably >> smaller than the one obtained from microarray analysis (using locfdr on >> pairwise t-statistics) at the same FDR. The RNASeq list is included in the >> microarray list, but there are several differences I have validated by qPCR >> that the RNASeq analysis is not able to detect at a reasonable FDR. >> Moreover, there seems to be an unusual bias towards not being able to >> detect down-regulated genes. I am a little bit puzzled by this, since one >> of the reasons we are sequencing is that it is supposed to have a better >> dynamic range. >> >> These are the same RNA samples so this apparent lack of sensitivity has to >> be related to either library prep or statistical analysis. So these are my >> questions: >> >> - can the inability to distinguish down-regulated genes be related to >> filtering low count reads? (in order to get good separation between groups >> in an MDS plot I need to filter cpm >0.1) >> - Is it possible that I need more coverage to improve sensitivity? I am >> currently sequencing at 50X pair end, that seemed enough. Is there any >> published study looking at RNASeq sequencing depth and sensitivity in human >> or mouse genomes? >> - Are the multiple testing corrections applied in EdgeR and DESeq too >> stringent thus rendering the overall analysis less sensitive? >> >> For the record my count matrices are of counts of transcripts, averaging >> counts over all exons from the same gene model for all RefSeq genes. I did >> this because the microarray data is per transcript. In log scale I have on >> average 0.7 R2 correlation between microarray intensity and RPKM from the >> same sample. >> >> Thanks for the insight! >> >> Lucia >> >> >> >> >> >> >> >> >> -- >> 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]] >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at r-project.org >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor >
ADD REPLY •
0
Entering edit mode
Hi Lucia On 16/05/13 19:25, Lucia wrote: > My counts are in fact total counts per transcript and not averages, > my mistake in the original post I hope you mean counts per _gene_, not per transcript. Otherwise, it is easy to guess what is going wrong. Maybe explain in more detail how you obtained the counts. There are surprisingly many ways of getting this wrong. You are right that filtering for PCR duplicates in RNA-Seq is rarely a good idea due to the dramatic impact this has on the dynamic range. Simon
ADD REPLY •
0
Entering edit mode
Dear Simon, To obtain my count matrix I use each Ucsc gene model as one " transcript", I limit my comparisons with the microarray data to those Ucsc gene models that have a unique RefSeq match with an Affy probe set Regarding multiple testing, I guess there is no reason why I can't use the uncorrected pvalues produced by DESeq and run locfdr right? Thanks for the help Lucia Sent from my iPad On May 16, 2013, at 3:02 PM, Simon Anders <anders at="" embl.de=""> wrote: > Hi Lucia > > On 16/05/13 19:25, Lucia wrote: >> My counts are in fact total counts per transcript and not averages, >> my mistake in the original post > > I hope you mean counts per _gene_, not per transcript. Otherwise, it is easy to guess what is going wrong. > > Maybe explain in more detail how you obtained the counts. There are surprisingly many ways of getting this wrong. > > You are right that filtering for PCR duplicates in RNA-Seq is rarely a good idea due to the dramatic impact this has on the dynamic range. > > Simon > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
ADD REPLY •
0
Entering edit mode
Dear Lucia On 18/05/13 14:01, Lucia wrote: > To obtain my count matrix I use each Ucsc gene model as one " > transcript", I limit my comparisons with the microarray data to those > Ucsc gene models that have a unique RefSeq match with an Affy probe > set You have to understand that with such piecewise information, this thread is not very helpful for anybody. So, for the record and for others who might read this thread: Lucia's observation that RNA-Seq shows less sensitivity than microarrays is unusual and not what one usually sees. The reason is likely either extraordinarily low read counts in the RNA-Seq experiment or a mistake in the generation of the count table. Lucia: If you want the help of the list to get to the bottom of it, you have to explain what you did. This means giving a precise and complete description of your workflow, including all commands, scripts, parameters, outputs etc. I don't understand why you don't want to tell us _how_ you created your count table (which command-line tool, R package, or self-made script have you used and how?), but without that, and some useful statistics about read counts (number of total and mapped read counts, column sum of count matrix, size factors reported by DESeq), this discussion won't lead anywhere. > Regarding multiple testing, I guess there is no reason why I can't > use the uncorrected pvalues produced by DESeq and run locfdr right? Sure, this is fine, even though the use of _local_ FDRs is a it unusual in transcriptomics, where people prefer to use tail-based FDRs such as those reported by the Benjamini-Hochberg method. Simon
ADD REPLY •
0
Entering edit mode
Dear Lucia and list On second reading, I noticed that my previous post sounded quite aggressive, which was not my intention. Sorry. I shouldn't write e-mails that late at night. Anyway: We had a lot of discussion on this list and others recently about how to correctly obtain a count table for differential expression analysis from aligned RNA-Seq reads. From these discussions, it has become clear that this is a task with many more pitfalls than one might expect at first. In microarray analysis, there is no need to do this, and so read counting is a likely culprit when such discrepancies are noted. This is why exact details on the procedure are so important. Simon
ADD REPLY •
0
Entering edit mode
So you still think the read counting problem is the first issue on this topic shan [[alternative HTML version deleted]]
ADD REPLY •
0
Entering edit mode
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: My experiment has 2 conditions, with 9 replicates of each being hybridized to independent micro arrays, 5 of which were also sequenced by RNA- seq. The comparisons I am making are between those RNA samples that were both hybridized and sequenced. I have also an independent experiment in which I have looked at the same samples by qPCR in several genes, so I have a decent set of positive controls. RNA-seq was obtained after polyA selection, with an average of 70 million pairs of 100bp reads per sample, each sample having very similar number of reads. I did my alignments using RUM ( http://www.ncbi.nlm.nih.gov/pubmed/21775302?dopt=Abstract). I have ~65 million unique pairs, which cover about 10% of the genome, and I am using for DE analysis using both edgeR and DESeq. In terms of my count table, the counts are generated by RUM as described in the paper, I just parse them out into a table format using custom scripts. I do not have much experience with this, but in principle I do not see anything wrong with the way the counts are generated, neither do the developers which are down the hall. It produces very high correlation to Ct values by qPCR after logRPKM transformation based on the data I have seen from them. Here are the DESeq size factors CC3 1.1395037 CC5 1.2704539 CC6 0.7973817 CC7 1.1593902 CC8 0.9147930 FC3 0.9944328 FC5 0.8701476 FC6 0.9405719 FC7 0.8326854 FC8 1.1706489 Now, using BH instead of locfdr shrinks the DE list produced by the microarray considerably and then they are more similar to those produced by DESeq, so the choice for multiple testing correction does have a big effect The locfdr procedure is better able to detect true positives. It can be something particular about my data set, or maybe it is a good idea to use locfdr always, not sure. even then, my RNASeq data has an obvious bias towards up-regulation, as it detects almost no down-regulated genes. I did filter low counts as follows: >rs <- rowSums ( counts ( cds)) >use <- (rs > 10) >cdsFilt <- cds[ use, ] but my filter is not too stringent and shouldn't really have a big effect thanks very much for all your help and suggestions. Lucia On Mon, May 20, 2013 at 7:15 PM, Simon Anders <anders@embl.de> wrote: > Dear Lucia and list > > On second reading, I noticed that my previous post sounded quite > aggressive, which was not my intention. Sorry. I shouldn't write e-mails > that late at night. > > Anyway: We had a lot of discussion on this list and others recently about > how to correctly obtain a count table for differential expression analysis > from aligned RNA-Seq reads. From these discussions, it has become clear > that this is a task with many more pitfalls than one might expect at first. > In microarray analysis, there is no need to do this, and so read counting > is a likely culprit when such discrepancies are noted. This is why exact > details on the procedure are so important. > > Simon > -- 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]]
ADD REPLY •
0
Entering edit mode
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 at="" 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:
ADD REPLY •
0
Entering edit mode
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]]
ADD REPLY •
0
Entering edit mode
Hi Lucia, On Wed, May 22, 2013 at 12:57 PM, Lucia Peixoto <luciap at="" iscb.org=""> wrote: [snip] > 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? There are likely better ones, but here's the first one that fit your criteria from a google scholar search: Evaluating Gene Expression in C57BL/6J and DBA/2J Mouse Striatum Using RNA-Seq and Microarrays http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.001 7820 The methods section for the RNA-seq data begins like so: """ Total RNA from 21 male mice (10 B6 and 11 D2) was provided to the Oregon Health & Science University Massively Parallel Sequencing Shared Resource facility [25] for transcriptome sequencing (NCBI SRA accession number: SRA026846.1). """ HTH, -steve -- Steve Lianoglou Computational Biologist Department of Bioinformatics and Computational Biology Genentech
ADD REPLY •
0
Entering edit mode
On Wed, May 22, 2013 at 3:57 PM, Lucia Peixoto <luciap@iscb.org> wrote: > 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? > > Perhaps you will enjoy Hansen, K. D., Wu, Z., Irizarry, R. A. & Leek, J. T. Sequencing technology does not eliminate biological variability. Nat Biotechnol 29, 572–573 (2011). For your question, you should also look at the supplementary figures. Kasper > 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]] > > _______________________________________________ > Bioconductor mailing list > Bioconductor@r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > [[alternative HTML version deleted]]
ADD REPLY •
0
Entering edit mode
Hi Lucia, just one comment to one of your questions. On 5/22/13 9:57 PM, Lucia Peixoto wrote: > 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? you can find matching RNA-seq and microarray biological replicates in the following Bioconductor experiment data package: library(tweeDEseqCountData) data(commonPickrell1Huang) in the help page of those data, i.e., do help(commonPickrell1Huang), you will find all details about them. i believe they could help you to explore this question by applying the same pipeline you use on your own data. cheers, robert. > Lucia > > > > > > > > > On Tue, May 21, 2013 at 6:49 PM, Wolfgang Huber <whuber at="" 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 at="" 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: > >
ADD REPLY •
0
Entering edit mode
Hi Lucia On 21/05/13 23:19, Lucia Peixoto wrote: > In terms of my count table, the counts are generated by RUM as described > in the paper, I just parse them out into a table format using custom > scripts. I do not have much experience with this, but in principle I do > not see anything wrong with the way the counts are generated, neither do > the developers which are down the hall. It produces very high > correlation to Ct values by qPCR after logRPKM transformation based on > the data I have seen from them. I did not know that RUM produces a count table, but if this works well, it sounds like a useful feature. However, as I pointed out in my reply to Thomas, counting for expression strength estimation and counting for differential testing are slightly different tasks, and it would be important to know what the RUM developers aimed for when they decided on a strategy to deal with ambiguously mapped reads. If you ask them, let us know what they say, please. To get back to the original issue: The sensitivity of RNA-Seq obviously depends on read depth. If you have few reads, Poisson noise becomes strong, and such an RNA-Seq experiment will be less precise than a microarray assay, especially if overall expression strength of genes is all you care about. The break-even point where RNA-Seq becomes better than microarrays is probably somewhere at 1 to 10 million reads. (This is a rough guess, the papers mentioned in the thread may contain proper figures.) You seem to have more than that. But you also mentioned that you have an unusually high number of PCR duplicates, which means that the information content of your reads is reduced. Also, check the column sums of your count table: Which fraction of your reads actually gets assigned to a gene? It's rather common for things to go wrong here. Simon
ADD REPLY •
0
Entering edit mode
Hi Lucia, I know this post is a little bit redundant to Simons Post but I would also think that it's worth checking logFC and mean Expression Values of the Microarray DE transcripts that have not been flagged as DE by RNA-Seq and maybe as Simon stated, technical noise is just too large to overcome the effect size of your treatment. There is a nice Paper http://bioinformatics.oxfordjournals.org/content/27/13/i383.full that investigates the relationship of transcript coverage (or number of reads per transcript) and measurement precision. As far as I remember the authors are rather pessimistic regarding measurement precision of low expressed transcripts using nowadays common sequencing depth. For instance they could only reliably (with 20% or less error) detect 41% of transcripts in mammary epithelial cells using a whole SOLID Flowcell (331 Mreads). Also there is non neglectable noise occuring during library preparation and finally low amounts of RNA used for library prep could introduce a large bias. Hope that helped a little bit Moritz 2013/5/23 Simon Anders <anders@embl.de> > Hi Lucia > > > On 21/05/13 23:19, Lucia Peixoto wrote: > >> In terms of my count table, the counts are generated by RUM as described >> in the paper, I just parse them out into a table format using custom >> scripts. I do not have much experience with this, but in principle I do >> not see anything wrong with the way the counts are generated, neither do >> the developers which are down the hall. It produces very high >> correlation to Ct values by qPCR after logRPKM transformation based on >> the data I have seen from them. >> > > I did not know that RUM produces a count table, but if this works well, it > sounds like a useful feature. However, as I pointed out in my reply to > Thomas, counting for expression strength estimation and counting for > differential testing are slightly different tasks, and it would be > important to know what the RUM developers aimed for when they decided on a > strategy to deal with ambiguously mapped reads. If you ask them, let us > know what they say, please. > > To get back to the original issue: The sensitivity of RNA-Seq obviously > depends on read depth. If you have few reads, Poisson noise becomes strong, > and such an RNA-Seq experiment will be less precise than a microarray > assay, especially if overall expression strength of genes is all you care > about. > > The break-even point where RNA-Seq becomes better than microarrays is > probably somewhere at 1 to 10 million reads. (This is a rough guess, the > papers mentioned in the thread may contain proper figures.) > > You seem to have more than that. But you also mentioned that you have an > unusually high number of PCR duplicates, which means that the information > content of your reads is reduced. Also, check the column sums of your count > table: Which fraction of your reads actually gets assigned to a gene? It's > rather common for things to go wrong here. > > Simon > > > ______________________________**_________________ > Bioconductor mailing list > Bioconductor@r-project.org > https://stat.ethz.ch/mailman/**listinfo/bioconductor<https: stat.et="" hz.ch="" mailman="" listinfo="" bioconductor=""> > Search the archives: http://news.gmane.org/gmane.** > science.biology.informatics.**conductor<http: news.gmane.org="" gmane.="" science.biology.informatics.conductor=""> > -- *Moritz Heß PhD Candidate * *Research associate Forest Research Institute of Baden Württemberg (FVA) Wonnhalde 4 79100 Freiburg (Germany) phone +49 761 4018 301* [[alternative HTML version deleted]]
ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 548 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6