I've been wondering a similar thing: suppose I use transcripts per million (TPM) as a coherent estimate of abundance and feed it to limma/voom in order to fit a multi-factorial blocked model in a poorly annotated organism (draft genome and txome, but orthologs poorly characterized thus far). It appears to work properly in human and mouse samples (i.e. the top hits are reasonable); is it sensible to generalize to non-model or poorly- annotated organisms in this fashion? Probably something that could be puzzled out given time, but since the authors have likely addressed this before, a reference would be great. Statistics is the grammar of science. Karl Pearson <http: en.wikipedia.org="" wiki="" the_grammar_of_science=""> On Mon, Aug 25, 2014 at 2:50 AM, assaf www <assafwww at="" gmail.com=""> wrote: > Dear Edger developers and users, > > I would like to compare transcription levels of orthologous genes belonging > to different species, in order to find significant species-dependent > changes in transcription levels. I though of using Edger for such > analysis. > Specifically, I have the read-counts data for several RNA-Seq samples, for > 2 different species (e.g., read counts produced by Htseq-count, and Rsem). > > I would like to ask: > 1) because Edger uses CPM values, which are not normalized by gene- length, > and because the length of orthologous genes differ, it would lead to a > serious length-dependet bias, and I would ask how to normalize for that. > 2) if the above length-bias can be eliminated, and the compared genes are > true orthologs, are you aware of any other major problems that should be > considered in the above case ? > > > Thanks in advance, > Assaf > > [[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 > [[alternative HTML version deleted]]

This is super helpful. Just to be clear, the most robust solution is to use edgeR and offset for putative gene length, TMM & library size while using raw counts (not effective counts) estimated by e.g. RSEM, eXpress, or the like? Also re: cross-species comparisons, while my experience is that it is indeed a can of worms, Mark Gerstein's group recently published a method that might interest others working on non-model or incompletely annotated organisms: http://genomebiology.com/2014/15/8/R100 Any thoughts on applicability of the method for kooky experiments such as comparing Drosophila hemocytes, zebrafish vascular endothelial progenitors and the same in mice? Or for that matter, alligator differentiation. I never realized how hard RNAseq in non-model organisms was until I tried it. Statistics is the grammar of science. Karl Pearson <http: en.wikipedia.org="" wiki="" the_grammar_of_science=""> On Thu, Aug 28, 2014 at 8:37 AM, assaf www <assafwww at="" gmail.com=""> wrote: > I checked, it's true, the results look the same. > As for FPKM, of course you are right. > > Thanks a lot > Assaf > > > > On Thu, Aug 28, 2014 at 2:47 AM, Gordon K Smyth <smyth at="" wehi.edu.au=""> wrote: > > > The code should have been: > > > > offset <- expandAsMatrix(getOffset(y),dim(y)) > > offset <- offset + gl > > > > This should give same result as your code. > > > > rpkm() corrects for gene length as well as library size -- that's the > > whole purpose of RPKMs: > > > > rpkm(y, gene.length=geneLength) > > > > should work for you without any modification. > > > > Gordon > > > > > > On Wed, 27 Aug 2014, assaf www wrote: > > > > This is very helpful for me, thanks. > >> > >> A slight change that I made in the code you sent, to avoid some R erros, > >> is > >> > >> # to replace: > >> offset = offset + gl > >> # with: > >> offset = sweep(gl,2,offset,"+") > >> > >> In addition to differential expression tests, I wanted also to extract > >> FPKMs values (and/or normalized CPM values), that would take into > account > >> all components of the offset (which if I'm not mistaken are: > library_size > >> + > >> TMM + gene_size). > >> I assume rpkm()/cpm() should correct only for library_size + TMM. > >> Is there a possibly "decent" solution for that ? > >> > >> all the best, and thanks, > >> Assaf > >> > >> > >> > >> On Wed, Aug 27, 2014 at 4:45 AM, Gordon K Smyth <smyth at="" wehi.edu.au=""> > >> wrote: > >> > >> It works something like this: > >>> > >>> library(edgeR) > >>> y <- DGEList(counts=counts) > >>> y <- calcNormFactors(y) > >>> > >>> # Column correct log gene lengths > >>> # Columns of gl should add to zero > >>> > >>> gl <- log(geneLength) > >>> gl <- t(t(gl)-colMeans(gl)) > >>> > >>> # Combine library sizes, norm factors and gene lengths: > >>> > >>> offset <- expandAsMatrix(getOffset(y)) > >>> offset <- offset + gl > >>> > >>> Then > >>> > >>> y$offset <- offset > >>> y <- estimateGLMCommonDisp(y,design) > >>> > >>> etc. > >>> > >>> Note that I have not tried this myself. It should work in principle > from > >>> a differential expression point of view. > >>> > >>> On the other hand, there may be side effects regarding dispersion trend > >>> estimation -- I do not have time to explore this. > >>> > >>> Gordon > >>> > >>> --------------------------------------------- > >>> Professor Gordon K Smyth, > >>> Bioinformatics Division, > >>> Walter and Eliza Hall Institute of Medical Research, > >>> 1G Royal Parade, Parkville, Vic 3052, Australia. > >>> http://www.statsci.org/smyth > >>> > >>> On Wed, 27 Aug 2014, assaf www wrote: > >>> > >>> Probably wrong, but the reason I thought of using quantile > normalization > >>> > >>>> is > >>>> the need to correct both for the species-length, and library size. > >>>> > >>>> > >>>> On Wed, Aug 27, 2014 at 2:40 AM, Gordon K Smyth <smyth at="" wehi.edu.au=""> > >>>> wrote: > >>>> > >>>> That doesn't look helpful to me. I suggested that you incorporate > gene > >>>> > >>>>> lengths into the offsets, not do quantile normalization of cpms. > >>>>> > >>>>> Sorry, I just don't have time to develop a code example for you. I > >>>>> hope > >>>>> someone else will help. > >>>>> > >>>>> The whole topic of interspecies differential expression is a can of > >>>>> worms. > >>>>> Even if you adjust for gene length, there will still be differences > in > >>>>> gc > >>>>> content and mappability between the species. > >>>>> > >>>>> Gordon > >>>>> > >>>>> > >>>>> On Wed, 27 Aug 2014, assaf www wrote: > >>>>> > >>>>> Dear Gordon thanks, > >>>>> > >>>>> > >>>>>> Suppose I start with the following matrices: > >>>>>> > >>>>>> # 'counts' is the Rsem filtered counts > >>>>>> > >>>>>> counts[1:4,] > >>>>>> > >>>>>>> > >>>>>>> h0 h1 h2 n0 n1 n2 > >>>>>>> > >>>>>> ENSRNOG00000000021 36 17 20 10 25 38 > >>>>>> ENSRNOG00000000024 1283 615 731 644 807 991 > >>>>>> ENSRNOG00000000028 26 12 11 18 23 28 > >>>>>> ENSRNOG00000000029 22 13 12 16 17 15 > >>>>>> > >>>>>> # 'geneLength' is the species-specific gene lengths, for species 'h' > >>>>>> and > >>>>>> 'n': > >>>>>> > >>>>>> geneLength[1:3,] > >>>>>> > >>>>>>> > >>>>>>> h0.length h1.length h2.length n0.length > n1.length > >>>>>>> > >>>>>> n2.length > >>>>>> ENSRNOG00000000021 1200 1200 1200 1303 1303 > >>>>>> 1303 > >>>>>> ENSRNOG00000000024 1050 1050 1050 3080 3080 > >>>>>> 3080 > >>>>>> ENSRNOG00000000028 1047 1047 1047 1121 1121 > >>>>>> 1121 > >>>>>> > >>>>>> > >>>>>> does the following code look correct, and may allow allows across > >>>>>> species > >>>>>> analysis ?: > >>>>>> (technically it works, I checked) > >>>>>> > >>>>>> # quantile normalization: (as in here: > >>>>>> http://davetang.org/wiki/tiki-index.php?page=Analysing+ > >>>>>> digital+gene+expression > >>>>>> ) > >>>>>> > >>>>>> x1 = counts*1000/geneLength > >>>>>> library(limma) > >>>>>> x2 = normalizeBetweenArrays(data.matrix(x1),method="quantile") > >>>>>> offset = log(counts+0.1)-log(x2+0.1) > >>>>>> > >>>>>> ... > >>>>>> > >>>>>> y <- estimateGLMCommonDisp(y,design,offset=offset) > >>>>>> y <- estimateGLMTrendedDisp(y,design,offset=offset) > >>>>>> y <- estimateGLMTagwiseDisp(y,design,offset=offset) > >>>>>> fit <- glmFit(y,design,offset=offset) > >>>>>> > >>>>>> ... > >>>>>> > >>>>>> > >>>>>> Thanks in advance for any help.., > >>>>>> Assaf > >>>>>> > >>>>>> > >>>>> ____________________________________________________________ > >>> __________ > >>> The information in this email is confidential and intended solely for > the > >>> addressee. > >>> You must not disclose, forward, print or use it without the permission > of > >>> the sender. > >>> ______________________________________________________________________ > >>> > >>> > >> > > ______________________________________________________________________ > > The information in this email is confidential and inte...{{dropped:10}} > > _______________________________________________ > 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 > [[alternative HTML version deleted]]