[BioC] Limma - RNA-Seq DE genes - Quantile normalized log transformed RPKM data

[Gordon K Smyth]

Dear Avinash,

Briefly, you have used

  design <- model.matrix(~0+Group)

when you needed

  design <- model.matrix(~Group)

However, I do not recommend a limma analysis of RPKM values, because it 
does not respect the mean-variance relationship inherent in the count 
data.  Nicole Cloonan's excellent Nature Methods paper was written four 
years ago, and our understanding of the statistical analysis RNA-Seq data 
has moved on considerably since then.  Please see the last case study in 
the limma User's Guide, which analyses the Nigerian HapMap data, for how I 
recommend limma be used to analyse RNA-Seq data.

Best wishes
Gordon


> Date: Fri, 2 Dec 2011 10:31:26 -0600
> From: Avinash S <avins.s at googlemail.com>
> To: bioconductor at r-project.org
> Subject: [BioC] Limma - RNA-Seq DE genes - Quantile normalized log
> 	transformed RPKM data
>
> Hello All,
>
> I'm trying to understand the method of differential expression analysis
> described in :
> *
> *
> *Stem cell transcriptome profiling via massive-scale mRNA sequencing*
> *Nicole Cloonan et al*
> *NATURE METhODS | VOL.5 NO.7 | JULY 2008 | 613*
> http://www.nature.com/nmeth/journal/v5/n7/abs/nmeth.1223.html
> *Section: Statistical Analysis *
> To calculate differential expression of SQRL tag data we analyzed the
> normalized gene signals (tags per Refseq transcript, length-normalized) for
> each library using the Limma package in R. After Quantile normalization, we
> used Limma to fit a linear model to the log2-transformed data using an
> empirical Bayes method32 to moderate standard errors. Differentially
> expressed genes were defined as those with a B statistic > zero.
>
> I have quantile normalized log2-transformed RPKM data and I wanted to find
> DE genes based on B statistic and log2foldchange.
> *Sample Rawdata:*
>
> GeneModel CON1 CON2 TR1 TR2
> 1s00200.1 2.723945276 3.775939811 3.623211245 3.717795434
> 1s00210.1 4.999354495 3.738129253 3.268468778 3.822220668
> 1s00220.1 1.450861588 1.265013193 0.942706046 0.744551693
>
> I'm using the following R-script:
>
> library(limma)
> raw.data <-
> read.delim("INPUT-QuantileNormalizedLog2Transformed-RPKM-Data.txt")
> attach(raw.data)
> names(raw.data)
> d <- raw.data[, 2:5]
> rownames(d) <- raw.data[, 1]
> Group <- factor(c("CON","CON","MYC","MYC"), levels=c("CON","MYC"))
> design <- model.matrix(~0+Group)
> colnames(design) <- c("CON","MYC")
> fit <- lmFit(d, design)
> fit <- eBayes(fit)
> *Warning message:*
> *Zero sample variances detected, have been offset*
> topTable(fit,coef=2,number=1000,genelist=fit$genes)
> - THIS LISTS THE GENES ALL POSITIVE LOGFC VALUES - NON -NEGATIVE.
>
> R version 2.14.0 (2011-10-31)
> Platform: i386-apple-darwin9.8.0/i386 (32-bit)
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] limma_3.10.0
>
>
> Can someone please let me know what I'm doing wrong here. Is it in the
> array design or input data?
>
> Thank you,
> Avinash

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}