Dear List,

 

I'm using genefilter and Limma to analyze the human exon array.  I did
my preprocessing with aroma.affymetrix and then loaded it into limma and
used genefilter to isolate the significant genes.  At least that is what
I think I am doing.  Here is my script starting at the function I have
the question about

 

 

##Filter for differentialy expressed genes

f1<-kOverA(2,7)

ffun<-filterfun(f1)

which<-genefilter(exprs(RMASet), ffun)

rma_no<-summary(which)

rmafiltset<-RMASet[which,]

 

 

 

grps <- paste(pData(rmafiltset)$hours)

lev <- c("time0", "time3", "time6", "time12", "time24", "time48")

f <- factor(grps, levels = lev)

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

colnames(design) <- lev

rmafit <- lmFit(rmafiltset, design)

 

 

 

contrast.matrix <- makeContrasts("time3-time0", "time6-time3",
"time12-time6", "time24-time12", "time48-time24", "time6-time0",
"time12-time0", "time24-time0", "time48-time0", levels=design)

rmafit2 <- contrasts.fit(rmafit, contrast.matrix)

rmafit3 <- eBayes(rmafit2)

 

 

pval<-0.05

l2foldch<-  0.58496250072

 

rmaresults<-decideTests(rmafit3, p.value=pval, lfc=l2foldch)

 

 

rmalist1 <- rmaresults[,1]!=0

rmalist2 <- rmaresults[,2]!=0

rmalist3 <- rmaresults[,3]!=0

rmalist4 <- rmaresults[,4]!=0

rmalist5 <- rmaresults[,5]!=0

rmalist6 <- rmaresults[,6]!=0

rmalist7 <- rmaresults[,7]!=0

rmalist8 <- rmaresults[,8]!=0

rmalist9 <- rmaresults[,9]!=0

 

rmagroup1 <- rmafiltset[rmalist1,]

rmagroup2 <- rmafiltset[rmalist2,]

rmagroup3 <- rmafiltset[rmalist3,]

rmagroup4 <- rmafiltset[rmalist4,]

rmagroup5 <- rmafiltset[rmalist5,]

rmagroup6 <- rmafiltset[rmalist6,]

rmagroup7 <- rmafiltset[rmalist7,]

rmagroup8 <- rmafiltset[rmalist8,]

rmagroup9 <- rmafiltset[rmalist9,]

 

rmaglist1 <- featureNames(rmagroup1)

rmaglist2 <- featureNames(rmagroup2)

rmaglist3 <- featureNames(rmagroup3)

rmaglist4 <- featureNames(rmagroup4)

rmaglist5 <- featureNames(rmagroup5)

rmaglist6 <- featureNames(rmagroup6)

rmaglist7 <- featureNames(rmagroup7)

rmaglist8 <- featureNames(rmagroup8)

rmaglist9 <- featureNames(rmagroup9)

 

rmaunion <- union(rmaglist1, union(rmaglist2, union(rmaglist3,
union(rmaglist4, union(rmaglist5, union(rmaglist6, union(rmaglist7,
union(rmaglist8, rmaglist9))))))))

 

rmaunionset <- rmafiltset[rmaunion,]

 

plotdata2 <- exprs(rmaunionset)

colnames(plotdata2) <- paste(pData(RMASet)$label)

 

 

I go on to create a heatmap and html output of significant genes.   My
problem is that when I change kOverA(2, 7) to, say, kOverA(2, 6), I get
fewer genes, I thought this should increase the number of genes I'm
seeing in the final output.  If I increase from 7 to 8 I still get a
shorter list.  It seems like 7 is the value that generates the maximum
output.  If anyone could explain why this is happening I would really
appreciate it.   If I should post more information please just let me
know.

 

 

Ray Roberts

 

Lovelace Respiratory Research Institute

 


	[[alternative HTML version deleted]]

