Hi Christian. > > 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?
Short answer is I haven't investigated this too thoroughly. But, my guess is that it wouldn't be "too" different. That prediction is based on the fact that the chip effects are in the same ballpark, as you can see from the Aroma_vs_Affy (Aroma=R_rlm_rma, Affy=medpol) plot in the following thread: http://groups.google.com/group/aroma-affymetrix/browse_thread/thread/1b0ab11fad9b4df3/f745ed0860546313 But, I'd be interested to hear more details if you do look into it more. Cheers, Mark > 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 >> ------------------------------ > > ------------------------------ 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 -~----------~----~----~----~------~----~------~--~---
