Hi everybody,
I would like to ask some general questions about the limma analysis.
Suppose you have an experiment with two groups (n_1 and n_2 arrays for
each group) and you have selected a specific gene g with n_g
probesets. If you would like to test a probeset of that gene for
differential expression/alternative splicing you could just do a t-
test for that probeset, right?
Now you want to test the whole gene for alternative splicing. Since
you have n_g probesets you have a multiple comparison and therefor you
should use a multicomparison test, right?
My impression was, the limma-package can do this multicomparison test
by executing:
#fs is the 'standard' FirmaSet-object
fsDF <- extractDataFrame(fs, addNames=TRUE)
fsDF <- fsDF[(1:1000),] #only for test purpose to speed things up
fsDF[,-c(1:5)] <- log2(fsDF[,-c(1:5)])
design <- cbind(Grp1=1,Grp2=c(rep(0,n_1),rep(1,n_2)))
fit<-lmFit(fsDF[,-c(1:5)],design)
fit<-eBayes(fit)
fit$genes<-fsDF[,1]
(x<-topTable(fit,coef="Grp2",adjust="BH",number=10,genelist=fsDF[,
1]))
...
ID logFC t P.Value
adj.P.Val B
964 ENSG00000004455 -5.397440e+00 -1.075052e+01 5.735220e-07
0.0005654927 6.0790934
961 ENSG00000004455 3.835360e+00 8.000724e+00 9.142608e-06
0.0045073059 3.8392302
963 ENSG00000004455 2.538539e+00 6.312477e+00 7.306077e-05
0.0240126401 1.9954011
962 ENSG00000004455 2.896122e+00 5.319056e+00 2.944042e-04
0.0706149086 0.7042825
...
Obviously this provides test-statistics per probeset and not per gene.
Then what is the difference to the normal t-test per exon mentioned
above (moderated t-statistic?)?
The User's Guide, Chapter 10 gives some explanations for the topTable-
statistics. For example the B-statistic = log-odds that the >GENE< is
differentially expressed. But with the above input I get a B-statistic
for every >PROBESET<. (Or maybe I got the meaning of "gene" in the
User's Guide wrong?)
So my problem is: shouldn't I somehow group my probesets into genes or
'tell' limma which probesets relate to which gene, if I want to draw
any conclusions on a gene-level and not on a probeset-level. It seems
that limma (if executed like above) does not, what I expected. And so
I got big problems interpreting the limma results. Could someone help
me out and tell me, how to interpret the given output of topTable to
me and especially direct me how I can come from probeset-based
conclusions to gene-based conclusions.
Furthermore I got no idea what the difference is between
design <- cbind(Grp1=1,Grp2=c(rep(0,n_1),rep(1,n_2)))
design <- cbind(Grp1=c(rep(1,n_1),rep(0,n_2)),Grp2=c(rep(0,n_1),rep
(1,n_2)))
Both designs are given in the User's Guide and were already discussed
here in another thread, but I can't tell the difference. Could it be
explained in more practical words? I guess that you also have to
adjust the coef-variable in topTable according to the used design-
variable.
I know these are very general questions, but maybe someone can help me
out. Thanks a lot,
Frank
--~--~---------~--~----~------------~-------~--~----~
When reporting problems on aroma.affymetrix, make sure 1) to run the latest
version of the package, 2) to report the output of sessionInfo() and
traceback(), and 3) to post a complete code example.
You received this message because you are subscribed to the Google Groups
"aroma.affymetrix" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/aroma-affymetrix?hl=en
-~----------~----~----~----~------~----~------~--~---
