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 <mrobin...@wehi.edu.au> 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: m.robin...@garvan.org.au > e: mrobin...@wehi.edu.au > 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 aroma-affymetrix@googlegroups.com To unsubscribe from this group, send email to aroma-affymetrix-unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/aroma-affymetrix?hl=en -~----------~----~----~----~------~----~------~--~---