[BioC] limma: topTable

Gordon K Smyth smyth at wehi.EDU.AU
Fri Mar 19 13:03:42 MET 2004


> Sorry, here comes my code:
>
> library(affy)
> library(limma)
> library(vsn)
> data  	<- ReadAffy()
> normalize.AffyBatch.methods <- c(normalize.AffyBatch.methods, "vsn")
> Set		<- expresso(data, normalize.method="vsn",
> bgcorrect.method="none",pmcorrect.method="pmonly",
> summary.method="medianpolish")
>
> # R and T stand for different treatment, R1 and R2 as well as T1 and T2
> are  biological reps.
> design      <- model.matrix(~ -1+factor(c(1,1,2,2,3,3,4,4)))
> colnames(design) <- c("R1","R2","T1","T2")
> fit              <- lmFit(Set, design)
> contrast.matrix  <- makeContrasts(((R1+R2)-(T1
> +T2))/2,R1-T1,R2-T2,levels=design);
> fit2             <- contrasts.fit(fit, contrast.matrix)
> fit2             <- eBayes(fit2)
> data          <- read.table("ATH1-121501_annot.csv", header=TRUE,
> sep=",") res              <- topTable(fit2, number=10,coef=1,
> adjust="fdr",genelist=data,sort.by="P",resort.by="M")

Mmmm, it appears that the genes are a different order in your data object
'Set' compared to your annotation file read into 'data'.  Just use

  genelist=geneNames(Set)

in the topTable argument list to ensure you're using the right names.

>From limma 1.5.8, topTable will extract the IDs directly from the exprSet
object so you won't have to input the names yourself.

Gordon

> Thanks a lot again,
> Julia
>
>> >Hi Bioconductor folks,
>> >
>> >I have really been enjoying using the limma package, but I have just
>> come across a problem.
>> >When I use the topTable-command, the slot "result$Probe.Set.ID" does
>> not
>> > seem to match the other entries I am interested in, namely M, P, t.
>> >I am using R 1.8.0, limma 1.5.5 on Affymetrix ATH1-121501 chips.
>> >
>> >An example would be:
>> >Output of the topTable-result:
>> >
>> >                 Probe.Set.ID            M               t
>> > P.Value         B
>> >5598    259302_at       4.593339 38.9793 1.126260e-05 11.88238
>> >
>> >If I check on gene number 5598 it says
>> >geneNames(myExprSet)[5598]
>> >[1] "250498_at"
>> >
>> >Am I misinterpreting 5598 as the index of my ExpressionSet?
>>
>> I think you've done something non-standard, anyway you can't have
>> simply used lmFit() and eBayes() and topTable() on our exprSet object.
>> Please show us enough of your code so that we can see how topTable is
>> getting the probe set ID's.
>>
>> Gordon
>>
>> >Thanks a lot for any suggestions!
>> >Julia
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor



More information about the Bioconductor mailing list