[BioC] Interspecies differential expression of orthologs with Edger

Tim Triche, Jr. tim.triche at gmail.com
Thu Aug 28 18:42:21 CEST 2014


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]]



More information about the Bioconductor mailing list