Search
Question: Limma- Filtering After Statistical Analysis
0
gravatar for Sally
9.8 years ago by
Sally250
Sally250 wrote:
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=TRU E,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="no ne",p.value=1) write.table(s24vss48,file="s24vss48nofdr.txt",sep="\t") s48vss96<-topTable(fit2,coef="s48vss96",number=17328,adjust.method="no ne",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="no ne",p.value=1) write.table(c24vsc48,file="c24vsc48nofdr.txt",sep="\t") c48vsc96<-topTable(fit2,coef="c48vsc96",number=17328,adjust.method="no ne",p.value=1) write.table(c48vsc96,file="c48vsc96nofdr.txt",sep="\t") [[alternative HTML version deleted]]
ADD COMMENTlink modified 9.8 years ago by Gordon Smyth35k • written 9.8 years ago by Sally250
0
gravatar for Gordon Smyth
9.8 years ago by
Gordon Smyth35k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth35k wrote:
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=T RUE,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="no ne",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="no ne",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")
ADD COMMENTlink written 9.8 years ago by Gordon Smyth35k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 351 users visited in the last hour