Dear Mark, Probably, there is not much difference, maybe I will check it. I will let you know but it will take some time.
Best regards Christian On Mar 17, 11:46 am, Mark Robinson <[email protected]> wrote: > 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/... > > 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 <[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 > >> ------------------------------ > > ------------------------------ > 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 -~----------~----~----~----~------~----~------~--~---
