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
-~----------~----~----~----~------~----~------~--~---
