Ryan -- that's extremely helpful, thanks! I'll try PCA and look at my dispersion estimates and see what light that sheds. Much appreciated. Jessica Jessica P. Hekman, DVM, MS PhD student, University of Illinois, Urbana-Champaign Animal Sciences / Genetics, Genomics, and Bioinformatics On 07/15/2014 05:23 PM, Ryan wrote: > Hi Jessica, > > Have you done some exploratory analysis on your dataset? A good place to > start would be to generate PCA plots using plotMDS for edgeR and plotPCA > for DESeq2. Do your samples cluster into two groups as expected on the > PCA plots? Secondly, you should have a look at the dispersion estimates > from both methods and compare them to typically observed values. You can > do this with plotBCV for edgeR and plotDispEsts for DESeq2 (but remember > that BCV is the square root of dispersion, so pay attention to the > y-axis label). If your dispersions are too high, this indicates that the > variation within groups is large, which means that detecting significant > differences between groups is difficult and you will get fewer genes. > The edgeR User's Guide says in section 2.10 that typical BCV values are > 0.1 for genetically identical animals (e.g. lab mice) and 0.4 for human > samples. The latter value is probably closer to what you should expect. > If your BCV is a lot higher than that, your experiment may not be > well-controlled, or there may be some other problem in the methods or > the data that you need to track down. > > Home this helps, > > -Ryan Thompson. > > On Tue Jul 15 14:54:47 2014, Jessica Perry Hekman wrote: >> I'm getting only a few dozen differentially expressed genes when I >> analyze my RNA-Seq data with DESeq2 (79) and EdgeR (34) (even fewer >> when I use EBSeq). I had expected many more -- hundreds or even a >> thousand. If this is the real answer, I'm fine with it, but I'm >> concerned that I'm doing something wrong. What are the ranges of >> numbers of differentially expressed genes that one would expect from >> DESeq2 or EdgeR? >> >> More information: >> >> I'm in the midst of my first RNA-seq project (as many of you have >> probably surmised from my frequent postings to a variety of lists). My >> initial goal is to get a list of differentially expressed (DE) genes. >> >> I have 24 samples, 12 from each of 2 treatment groups. >> >> My species is fox (Vulpes vulpes), which aligns very nicely to dog >> (Canis familiaris). >> >> My current approach is to use the dog reference genome (to which my >> fox reads align at about 83%) + GTF with location of exons. >> >> Can I feel confident about DESeq2 and EdgeR's calls? >> >> Thanks very much for any insights, >> >> Jessica >>

On 07/15/2014 05:32 PM, Steve Lianoglou wrote: > We can only help you to answer the last point, but to do so we will > need to see the code you used and an explanation of your dataset (ie. > which samples are which condition). > > Also: why do you expect to get upwards of a thousand differentially > expressed genes? Have you (or others) done this experiment before and > verified those results, or? Steve -- thanks for your feedback. This is my first time doing RNA-Seq and I had been under the impression that there would be a large number of DE genes just by chance in most experiments. I think to a large extent I just need to adjust my expectations. I just don't have a feel for what normal results look like! However, as you asked whether this had been done before, I sheepishly recalled that we do have some analysis from a different lab using different tissues from the same animals and a slightly different approach to dealing with the data (de novo transcriptome instead of aligning to dog). It turns out that they found 350-650 DE genes, depending on the tissue. So they did get more DE genes than I did. > What does the enrichment of reads that align to exons vs. intergenic > reads look like? I haven't done this analysis, but I'll take a look and see what it looks like. I'm guessing that I should expect 80-90% of reads to align to exons, since this is RNA? You suggested that seeing my code might be helpful, so here it is, in case you're curious: For EdgeR: x <- read.delim("cf3ncbi.genes.matrix", row.names="GeneName") group <- factor(c(1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2)) y <- DGEList(counts=x,group=group) y <- calcNormFactors(y) y <- estimateCommonDisp(y) y <- estimateTagwiseDisp(y) et <- exactTest(y) tt <- topTags(et, n=100) write.table(tt, file = "edgeR-cf3ncbi-topTags100.txt") ...where cf3ncbi.genes.matrix looks like: GeneName 1_481 2_124 3_23 4_124 5_125 6_126 7_127 8_128 9_129 10_130 11_131 12_132 49_169 50_170 51_171 52_172 53_173 54_174 55_535 56_176 57_177 58_178 59_179 60_180 A1BG 27 39 55 41 50 42 33 66 63 31 50 34 52 53 34 44 35 20 45 43 75 47 43 45 A1CF 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 and so forth (as suggested in the code, the first 12 samples are from one treatment group and the second 12 samples are from another). For DESeq2: sample.table <- data.frame(read.table("DESeq2SampleTable.txt", header=TRUE)) dds <- DESeqDataSetFromHTSeqCount(sample.table, directory = "cf3ncbi/", formula(~ type)) dds2 <- DESeq(dds) res <- results(dds2) write.table (res, "cf3ncbi.deseq2.results.csv") sum( res$padj < 0.05, na.rm=TRUE) ...where DESeq2SampleTable.txt looks like: sampleName sampleFile type 1_481 cf3ncbi-1_481.genes.counts tame 2_124 cf3ncbi-2_124.genes.counts tame ... Thanks, Jessica

Thanks, Tom. Yes, you summarized my dilemma well, although I am more concerned with false negatives right now than false positives (as we do intend to do PCR to validate any positives we get, but any false negatives lost are lost forever :). Ryan also suggested PCA and I'll definitely be trying that next. But getting the feedback that others have used these techniques and feel moderately confident in them (with followup work) is very helpful. Thanks, Jessica Jessica P. Hekman, DVM, MS PhD student, University of Illinois, Urbana-Champaign Animal Sciences / Genetics, Genomics, and Bioinformatics On 07/15/2014 06:52 PM, Thomas H. Hampton wrote: > Dear Jessica, > > The answer to this sort of question depends on the question you are actually asking. > > You no doubt know a) that in this sort of test many genes are expected to reach > nominal significance by chance and b) that software packages include algorithms > intended to differentiate between these predicted "false positives" and "true positives". > > So maybe you are asking : do these strategies work? > > The short answer is that they definitely do not work perfectly for every experiment or for every gene in any experiment. > > I look at a lot of data where effects are pretty subtle and the usual sort of multiple hypothesis testing corrections suggest that > there are *no* differentially expressed genes in my data. > > I have yet to run into a situation where follow up measurements and experiments validated the assertion that *no* genes were > in fact differentially expressed. > > If I were in your shoes, as I have been on numerous occasions, I would do at least some of the following. > > 1) Establish that samples from your experimental groups are more like samples from the same group than > they are like samples from other groups using comparisons like PCA, multi dimensional scaling or hierarchical clustering. > > 2) Consider your a priori knowledge of genes that you thought would go some place. Did they go in the directions you thought they > should go? > > 3) Use qPCR or some low throughput method on single genes to reassure yourself that genes of interest (say genes with a fold change > of 2 and a nominal, unadjusted p value of 0.05) are in fact regulated, even though multiple hypothesis testing correction suggests that > they may be the result of chance. > > Good luck, > > Tom > > > On Jul 15, 2014, at 5:54 PM, Jessica Perry Hekman wrote: > >> I'm getting only a few dozen differentially expressed genes when I analyze my RNA-Seq data with DESeq2 (79) and EdgeR (34) (even fewer when I use EBSeq). I had expected many more -- hundreds or even a thousand. If this is the real answer, I'm fine with it, but I'm concerned that I'm doing something wrong. What are the ranges of numbers of differentially expressed genes that one would expect from DESeq2 or EdgeR? >> >> More information: >> >> I'm in the midst of my first RNA-seq project (as many of you have probably surmised from my frequent postings to a variety of lists). My initial goal is to get a list of differentially expressed (DE) genes. >> >> I have 24 samples, 12 from each of 2 treatment groups. >> >> My species is fox (Vulpes vulpes), which aligns very nicely to dog (Canis familiaris). >> >> My current approach is to use the dog reference genome (to which my fox reads align at about 83%) + GTF with location of exons. >> >> Can I feel confident about DESeq2 and EdgeR's calls? >> >> Thanks very much for any insights, >> >> Jessica >> >> -- >> Jessica P. Hekman, DVM, MS >> PhD student, University of Illinois, Urbana-Champaign >> Animal Sciences / Genetics, Genomics, and Bioinformatics >> >> _______________________________________________ >> 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 >