[BioC] Limma- Filtering After Statistical Analysis

[Gordon K Smyth]

Dear Sally,

I haven't checked your R script, but why not simply

   topTable(fit2[i,],...)

where 'i' indexes the genes you want to keep?

Best wishes
Gordon


> Date: Sun, 8 Feb 2009 13:29:31 -0800
> From: "Sally" <sagoldes at shaw.ca>
> Subject: Re: [BioC] Limma- Filtering After Statistical Analysis
> To: <bioconductor at stat.math.ethz.ch>
> Message-ID: <9004B797DD0E49ACBF9B692CC2A67BE8 at sghome>
> Content-Type: text/plain
>
> Hi,
>
> I am relatively new to Limma.  I want to filter my gene list after my 
> statistical analysis.  What object do I filter on (fit2?)?  I looked at 
> genefilter and multtest and they both filter on the ExpressionSet.  But 
> this occurs before the statistical analysis.
>
> Sally
>
> Here is my Limma script:
>
>
> #Load libraries
> source("http://bioconductor.org/biocLite.R")
> biocLite()
> library(limma)
> library(Biobase)
> #change directory to folder where files are (c:/limmadegenes)
> #Change to directory with original data files
> #read in expression data file and phenotypic data.
> #Note that row.names=1 means that row names are #in column 1,
> exprdata<-read.table("exprsData.txt", header=TRUE,sep="\t",row.names=1,as.is=TRUE,fill=TRUE,)
> class(exprdata)
> #[1] "data.frame"
> dim(exprdata)
> #[1] 17328    28
> colnames(exprdata)
> head(exprdata)
> #printout too long to paste
> phenotypicdata<-read.table("phenotypicdata.txt",row.names=1,header=TRUE,sep="\t")
> class(phenotypicdata)
> #returns: [1] "data.frame"
> dim(phenotypicdata)
> #returns: [1] 28  2
> colnames(phenotypicdata)
> #returns: [1] "Species" "Time"
> rownames(phenotypicdata)
> #Coerse exprdata into a matrix
> myexprdata<-as.matrix(exprdata)
> write.table(myexprdata,file="myexprdata.txt",sep="\t",col.names=NA)
> class(myexprdata)
> #[1] "matrix"
> rownames(myexprdata)
> colnames(myexprdata)
> #Coerse phenotypicdata into a data frame
> myphenotypicdata<-as.data.frame(phenotypicdata)
> write.table(myphenotypicdata,file="myphenotypicdatacheck.txt",sep="\t",col.names=NA)
> rownames(myphenotypicdata)
> colnames(myphenotypicdata)
> #[1] "species" "time"
> summary(myphenotypicdata)
> all(rownames(myphenotypicdata)==colnames(myexprdata))
> #[1] TRUE
> #Create annotated Data Frame
> adf<-new("AnnotatedDataFrame",data=phenotypicdata)
> #dim means: dimension of an object.
> dim(adf)
> #rowNames columnNames
> #     28           2
> rownames(adf)
> #NULL
> #Create eset object
> eset<-new("ExpressionSet",exprs=myexprdata,phenoData=adf)
> #Read in targets file
> targets <- readTargets("targets.txt")
> targets
> # Set up character list defining your arrays, include replicates
> TS <- paste(targets$Species, targets$Time, sep=".")
> #This script returns the following:
> TS
> # Turn TS into a factor variable which facilitates fitting
> TS <- factor(TS)
> #This script returns the following
> design <- model.matrix(~0+TS)
> #write design object to text file
> write.table(design,file="design.txt",sep="\t",col.names=NA)
> colnames(design) <- levels(TS)
> #for eset put in your M values - see ?lmFit for object types
> fit <- lmFit(eset, design)
> cont.matrix<-makeContrasts(s0vss24=s.0-s.24, s24vss48=s.24-s.48, s48vss96=s.48-s.96, c0vsc24=c.0-c.24, c24vsc48=c.24-c.48, c48vsc96=c.48-c.96, levels=design)
> write.table(cont.matrix,file="cont.matrix.txt",sep="\t",col.names=NA)
> # estimate the contrasts and put in fit2
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)
>
> #print fit2 table
> write.table(fit2,file="fit2.txt",sep="\t")
>
> #print MArrayLM table
> write.fit(fit, file="MArrayLM.txt", adjust="none")
>
> s0vss24<-topTable(fit2,coef="s0vss24",number=17328,adjust.method="none",p.value=1)
>
> write.table(s0vss24,file="s0vss24nofdr.txt",sep="\t")
>
>
> s24vss48<-topTable(fit2,coef="s24vss48",number=17328,adjust.method="none",p.value=1)
>
> write.table(s24vss48,file="s24vss48nofdr.txt",sep="\t")
>
>
> s48vss96<-topTable(fit2,coef="s48vss96",number=17328,adjust.method="none",p.value=1)
>
> write.table(s48vss96,file="s48vss96nofdr.txt",sep="\t")
>
>
> c0vsc24<-topTable(fit2,coef="c0vsc24",number=17328,adjust.method="none",p.value=1)
>
> write.table(c0vsc24,file="c0vsc24nofdr.txt",sep="\t")
>
>
> c24vsc48<-topTable(fit2,coef="c24vsc48",number=17328,adjust.method="none",p.value=1)
>
> write.table(c24vsc48,file="c24vsc48nofdr.txt",sep="\t")
>
>
> c48vsc96<-topTable(fit2,coef="c48vsc96",number=17328,adjust.method="none",p.value=1)
>
> write.table(c48vsc96,file="c48vsc96nofdr.txt",sep="\t")