Dear Mark,
Thank you for your fast reply. You are right:
ad 2, "estimate of standard error", when I do:
> u.mad <- mad(unlist(log2(eps)), center=0)
> u.mad
[1] 0.3797971
I get for the 1.probeset
> y1 <- log2(eps[1:4,])
> F1 <- apply(y1/u.mad, 2, median)
HeartA HeartB HeartC MuscleA MuscleB
MuscleC
-0.37628338 0.13465003 -0.63489610 0.06787988 -0.30099279
0.99667714
and for the 2.probeset
> y2 <- log2(eps[5:8,])
> F2 <- apply(y2/u.mad, 2, median)
HeartA HeartB HeartC MuscleA
MuscleB MuscleC
-9.719282e-02 -7.810612e-01 -3.951703e-01 -4.654384e-01
-2.844947e-16
-1.108233e+0
ad 1, Yes, for the moment I am using median polish as an illustration.
However, in this respect I have also the following question:
How does using "median polish" compare to using
"R_rlm_rma_default_model"?
Are the final scores still of some use if you use medpol?
Best regards
Christian
On Mar 16, 9:13 am, Mark Robinson <[email protected]> wrote:
> Hi Christian.
>
> From what I can tell looking at your code (rather quickly, i must
> admit), there will be 2 differences between aroma.affymetrix and what
> you have:
>
> 1. We use the 'preprocessCore' codebase for the robust fitting of the
> linear model (... but maybe you are just using median polish as an
> illustration). For example, you might try:
>
> library(preprocessCore)
> f <- .Call("R_rlm_rma_default_model", log2(yTr), 0,
> 1.345,PACKAGE="preprocessCore")
> [... and piece together the alpha, beta, etc ...]
>
> 2. The "estimate of standard error" is calculated genewise, over
> residuals from all probes/samples (i.e. u.mad should be a scalar not a
> vector).
>
> Hope that helps.
> Mark
>
> On 16/03/2009, at 6:32 PM, cstratowa wrote:
>
>
>
>
>
> > Dear all,
>
> > After reading the FIRMA paper I would like to understand the
> > implementation, but this is not easy since the source code is hard to
> > read. Nevertheless, I tried and would like to know if this is correct.
>
> > According to the page on exon array analysis you do the following:
>
> > I, fit a summary of the entire transcript
> >> plmTr <- ExonRmaPlm(csN, mergeGroups=TRUE)
> >> fit(plmTr, verbose=verbose)
>
> > II, fit the FIRMA model for each exon
> >> firma <- FirmaModel(plmTr)
> >> fit(firma, verbose=verbose)
>
> > However, I would like to understand the underlying source code.
>
> > For this example let us assume that we have quantile-normalized
> > intensities yTr for a transcript containing two exons:
> >> yTr
> > HeartA HeartB HeartC MuscleA MuscleB MuscleC
> > 1 5.74954 18.0296 2.50436 15.5857 26.1744 31.0075
> > 2 9.59819 23.0093 22.01120 70.1742 32.8408 102.0080
> > 3 114.50800 87.1742 70.34080 312.3410 266.1740 601.3410
> > 4 66.34080 52.0075 67.34080 184.1740 266.1740 147.0080
> > 5 210.17400 142.0080 173.34100 514.5080 659.1740 509.6740
> > 6 104.00800 84.3408 70.34080 333.5080 324.1740 231.0080
> > 7 194.00800 124.5080 234.00800 443.6740 767.5080 716.8410
> > 8 319.34100 282.6740 283.50800 656.0080 807.6740 954.6740
>
> > Here rows 1:4 code for exon 1 and rows 5:8 code for exon 2.
>
> > I, fit a summary of the entire transcript
> > To simplify issues I will fit the data using median polish:
> > # 1. fit median polish
> >> mp <- medpolish(log2(yTr))
>
> > # 2. data set specific estimates (probe affinities)
> >> beta <- mp$overall+mp$col
> >> thetaTr <- 2^beta
>
> > # 3. array-specific estimates
> >> alpha <- mp$row
> >> alpha[length(alpha)] <- -sum(alpha[1:(length(alpha)-1)])
> >> phiTr <- 2^alpha
>
> > II, fit FIRMA model for each exon
> > # 1. calculate residuals
> >> phi <- matrix(phiTr, nrow=nrow(yTr), ncol=ncol(yTr))
> >> theta <- matrix(thetaTr, nrow=nrow(yTr), ncol=ncol(yTr),
> > byrow=TRUE)
> >> yhat <- phi *theta
> >> eps <- yTr/yhat # rma uses y/yhat
>
> > # 2. estimate of standard error
> >> u.mad <- apply(log2(eps), 2, mad, center=0)
>
> > # 3. compute final score statisitc
> > # for 1. exon
> >> y1 <- log2(eps[1:4,])
> >> F1 <- apply(y1/u.mad, 2, median)
> >> F1
> > HeartA HeartB HeartC MuscleA MuscleB
> > MuscleC
> > -0.89938777 -0.03792624 -0.69409936 0.11536565 -0.61385296
> > 1.08709568
>
> > # for 2. exon
> >> y2 <- log2(eps[5:8,])
> >> F2 <- apply(y2/u.mad, 2, median)
> >> F2
> > HeartA HeartB HeartC MuscleA MuscleB
> > MuscleC
> > -0.02899616 -1.64645153 -0.70048533 -0.39996057 0.02666064
> > -1.46657055
>
> > Now my question is:
> > Is this calculation of the final score statistic F1 for exon 1 and F2
> > for exon 2 correct?
> > Did I miss something?
>
> > Best regards
> > Christian
>
> ------------------------------
> Mark Robinson
> Epigenetics Laboratory, Garvan
> Bioinformatics Division, WEHI
> e: [email protected]
> e: [email protected]
> p: +61 (0)3 9345 2628
> f: +61 (0)3 9347 0852
> ------------------------------
--~--~---------~--~----~------------~-------~--~----~
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
-~----------~----~----~----~------~----~------~--~---
