Hi, Dr. Gordon and all:
 
Just follow up on my RNAseq analysis on the TCGA RSEM processed so-called raw counts using limma-voom package as Dr. Gordon suggested due to the issue of so-called RSEM raw counts but not true raw counts that edgeR required. just try to get your opinion on it. Here are the commands I used (only show most relevant ones)
 
> library(edgeR);
> library(limma)
> show(head(targets));
  SampleName      Subject   Type
1  T4626_01A TCGA_38_4626  Tumor
2  N4626_11A TCGA_38_4626 Normal
3  T2668_01A TCGA_44_2668  Tumor
4  N2668_11A TCGA_44_2668 Normal
> Patient<-factor(targets$Subject);
> Type<-factor(targets$Type);
> y<-DGEList(count=rawdata[,4:ncol(rawdata)],genes=rawdata[,1:3]);
> show(dim(y));
[1] 20531    58

##Keep genes with least 1 count-per-million reads (cpm) in at least 20 samples:
> isexpr <- rowSums(cpm(y)>1) >= 20;
> y2 <- y[isexpr,]
> show(dim(y2));
[1] 14735    58
> #Recompute the library sizes:
> y2$samples$lib.size <- colSums(y2$counts)
> #Use Entrez Gene IDs as row names:
> rownames(y2$counts) <- rownames(y2$genes) <- y2$genes$entrezID;
 
> ##Apply TMM normalization
> y2 <- calcNormFactors(y2,method=c("TMM"));
> plotMDS(y2,main="y2");

> show(data.frame(Sample=colnames(y2),Patient,Type));
      Sample      Patient   Type
1  T4626_01A TCGA_38_4626  Tumor
2  N4626_11A TCGA_38_4626 Normal
3  T2668_01A TCGA_44_2668  Tumor
4  N2668_11A TCGA_44_2668 Normal
...
> design<-model.matrix(~Patient+Type);
> rownames(design) <- colnames(y2);
> show(design[1:6,]);
          (Intercept) PatientTCGA_38_4626 PatientTCGA_38_4627
T4626_01A           1                   1                   0
N4626_11A           1                   1                   0
T2668_01A           1                   0                   0
...
> v <- voom(y2,design,plot=TRUE);
> plotMDS(v, top=50,main="common for v of y2",labels=paste(substr(Type,1,1),substr(Patient,9,12)),col=ifelse(Type=="Tumor","red","blue"),gene.selection="common");
> plotMDS(v, top=50,main="pairwise for v of y2",labels=paste(substr(Type,1,1),substr(Patient,9,12)),col=ifelse(Type=="Tumor","red","blue"),gene.selection="pairwise");
> fit <- lmFit(v,design)
> fit <- eBayes(fit)
> options(digits=3)
> Rslt<-topTable(fit,coef="TypeTumor",n=Inf,sort="p");
> show(summary(de <- decideTests(fit)));
   (Intercept) PatientTCGA_38_4626 PatientTCGA_38_4627 PatientTCGA_38_4632
-1         126                   1                   3                   1
0         1805               14732               14725               14720
1        12804                   2                   7                  14
...
   PatientTCGA_91_6836 TypeTumor
-1                   4      3917
0                14718      5426
1                   13      5392
 
 
Here are my questions:
 
1. I used isexpr <- rowSums(cpm(y)>1) >= 20 (y2 filtering) to do the filtering (based on one example in the package guide). I have 29 paired of matched tumor and normal samples (total 58 samples), is this reasonble to use for filtering, I asked this last time as well? I have a bit concern, since after filtering, the number of genes went from 20531 to 14735. Shall I use something like y3<-y[rowSums(y$counts)>=50,] (y3 filtering) to filter (this filtering leaves about 18215 genes, about 4k more genes than the 14735 from the other filtering way)?
 
2. I did plot with plotMDS(y2,main="y2"); v <- voom(y2,design,plot=TRUE); plotMDS(v, top=50,main...). The plotMDS all show very nice separation of tumors vs normals in dim1. However, the v <- voom(y2,design,plot=TRUE) plot slightly different between filtering with  isexpr <- rowSums(cpm(y)>1) >= 20  vs filtering with y3<-y[rowSums(y$counts)>=50,]. The first filtering way seems got similar plot as the one on page 116 of the user guide, but the second way seems got bending toward to the down side around log2(count size +0.5) between 0 to 5, and the data points also shift to 0 side due to more of lower reads data points (the first option only have very few data points between 0 to 5 in x-axis). Does this suggest the first filtering way (y2) is better?  

3. the above command show(summary(de <- decideTests(fit))); showed y2 filtering got  
3917 +5392=9309 DEGs at FDR5% (7661 at FDR 1%), whereas the y3 filtering (similar commands as y2, not shown) got final 4917+6022=10939 DEGs at FDR5%, both filtering seem getting too many DEGs, any issues here?   
 
4. using edgeR and the same data, I got more than 13k DEGs (total 18215 genes in the dataset after filtering) with FDR <0.05 (filtering like y3), now using limma voom, I still got quite a lot of DEGs (10939 DEGs at FDR5%) as shown in my question 3,  although slight less than edgeR. Any potenital issue?
 
Thanks so much in advance for your input and advice!
 
Ming 

 
 
 
 
 
> From: yi02@hotmail.com
> To: smyth@wehi.edu.au
> Date: Thu, 30 Jan 2014 03:07:56 +0000
> CC: bioconductor@r-project.org
> Subject: Re: [BioC] limma-voom
> 
> Hi, Dr. Smyth:
> 
> Thanks for the info. I used limma a lot but not much on voom. I will give a shot on voom. Thanks again and best Ming
> 
> > Date: Thu, 30 Jan 2014 13:05:57 +1100
> > From: smyth@wehi.EDU.AU
> > To: yi02@hotmail.com
> > CC: bioconductor@r-project.org
> > Subject: limma-voom
> > 
> > Hi Ming,
> > 
> > voom is part of the limma package. There is a voom case study in the 
> > limma User's Guide with complete working code. It is the last case study 
> > in the users guide.
> > 
> > voom works fine with either counts, or fractional counts, or scaled 
> > counts.
> > 
> > We are currently finalizing additional voom code to detect differential 
> > splicing. That code isn't in the public package, but will be soon.
> > 
> > Best wishes
> > Gordon
> > 
> > On Thu, 30 Jan 2014, Ming Yi wrote:
> > 
> > > I found the following documents
> > >
> > > http://stuff.mit.edu/afs/athena/software/r/r_v2.15.1/lib/R/library/limma/html/voom.html
> > > http://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf
> > > http://www.statsci.org/smyth/pubs/VoomPreprint.pdf
> > >
> > > But I could not find the real limma-voom package document. Is voom just 
> > > as normalization function prior to using limma linear model? Any good 
> > > examples of case study of working commands available as example?
> > >
> > > Thanks again!
> > > Ming
> > 
> > ______________________________________________________________________
> > The information in this email is confidential and inte...{{dropped:9}}
> 
> _______________________________________________
> 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]]

