[BioC] Ideal Mismatch value
James W. MacDonald
jmacdon at med.umich.edu
Wed Feb 28 20:47:29 CET 2007
I think I over-thought this problem.
You probably wouldn't have to do much work to get the values you want.
abatch <- ReadAffy()
abatch <- bg.correct.mas(abatch)
pmcorrect.mas.ims(pm(abatch, "someprobeID"), mm(abatch, "someprobeID"))
should do it.
pmcorrect.mas.ims <- function (pms, mms, contrast.tau = 0.03,
scale.tau = 10, delta = 2^(-20)){
diff <- log2(pms) - log2(mms)
delta <- rep(delta, nrow(diff))
out <- matrix(NA, ncol = dim(diff)[2], nrow = dim(diff)[1])
for (i in 1:ncol(diff)) {
sb <- tukey.biweight(diff[, i])
pps.pm <- pms[, i]
pps.mm <- mms[, i]
pps.im <- pps.mm
j <- (pps.mm >= pps.pm) & (sb > contrast.tau)
pps.im[j] <- pps.pm[j]/2^sb
j <- (pps.mm >= pps.pm) & (sb <= contrast.tau)
pps.im[j] <- pps.pm[j]/2^(contrast.tau/(1 + (contrast.tau -
sb)/scale.tau))
out[,i] <- pps.im
}
out
}
> library(affy)
> data(affybatch.example)
> affybatch.example <- bg.correct.mas(affybatch.example)
> pmcorrect.mas.ims(pm(affybatch.example, "A28102_at"),
mm(affybatch.example, "A28102_at"))
[,1] [,2] [,3]
[1,] 146.5810 116.02368 122.06893
[2,] 141.1703 122.70979 114.68573
[3,] 129.8570 109.14092 103.36482
[4,] 127.8895 96.35865 113.20909
[5,] 113.9200 112.38565 108.77917
[6,] 128.8732 140.60497 120.59229
[7,] 129.6603 121.92319 120.10008
[8,] 137.7271 126.83944 105.62900
[9,] 120.0194 169.11926 101.88818
[10,] 131.1359 175.01877 105.62900
[11,] 133.7921 143.35807 97.26137
[12,] 121.0031 135.68871 99.42711
[13,] 130.3489 128.80595 108.58228
[14,] 127.8895 111.89402 104.84146
[15,] 122.9707 117.00693 102.18351
[16,] 131.3327 161.25325 108.28695
> pm(affybatch.example, "A28102_at")
20A 20B 10A
A28102_at1 149.0 118.0 124.0
A28102_at2 143.5 124.8 116.5
A28102_at3 132.0 111.0 105.0
A28102_at4 130.0 98.0 115.0
A28102_at5 115.8 114.3 110.5
A28102_at6 131.0 143.0 122.5
A28102_at7 131.8 124.0 122.0
A28102_at8 140.0 129.0 107.3
A28102_at9 122.0 172.0 103.5
A28102_at10 133.3 178.0 107.3
A28102_at11 136.0 145.8 98.8
A28102_at12 123.0 138.0 101.0
A28102_at13 132.5 131.0 110.3
A28102_at14 130.0 113.8 106.5
A28102_at15 125.0 119.0 103.8
A28102_at16 133.5 164.0 110.0
> mm(affybatch.example, "A28102_at")
20A 20B 10A
A28102_at1 847.0 694.0 999.0
A28102_at2 860.3 667.3 1084.8
A28102_at3 815.3 650.0 1057.0
A28102_at4 855.0 664.0 1038.0
A28102_at5 729.0 594.5 1096.0
A28102_at6 711.0 619.3 933.0
A28102_at7 798.0 642.0 803.0
A28102_at8 800.0 631.5 734.0
A28102_at9 862.0 663.3 892.0
A28102_at10 835.5 605.3 1070.0
A28102_at11 886.0 614.0 1040.3
A28102_at12 900.0 589.5 981.0
A28102_at13 941.8 629.3 1062.0
A28102_at14 899.0 608.3 1061.5
A28102_at15 846.0 594.8 961.0
A28102_at16 860.0 538.0 927.0
Best,
Jim
James W. MacDonald wrote:
> Wang, Yonghong (NIH/NCI) [C] wrote:
>
>>Hi, All:
>>As part of the suggestions from a reviewer of my submitted paper, I am
>>asked to provide the IM (Ideal Mismatch) values generated from MAS5 for
>>some probes of some probeset. I know BioConductor can provide both PM
>>and MM intensities for each individual probes, would someone please tell
>>me whether I can also get IM values (when applicable) from cel file by
>>using BioConductor? If we can, how?
>
>
> The short answer is no.
>
> The long answer is yes, but it would take some work to do so. The work
> is done by the function pmcorrect.mas(), which is called on each
> probeset individually. You would have to background correct your
> AffyBatch using bgcorrect.mas(), then for whichever probesets you are
> interested in, you could either create your own function based on
> pmcorrect.mas() that returns the IM values, and just run that function
> on the probeset in question. Or you could just debug() the existing
> pmcorrect.mas(), and when the IM values you are interested in appear,
> you could dump them into the .GlobalEnv using the <<- operator.
>
> Note that pmcorrect.mas() takes a ProbeSet object as input. What is done
> in expresso() is to make an empty ProbeSet object, then get the PM and
> MM values for a particular probeset (using a combination of get() and
> intensity(), although you could probably just use pm(abatch,
> "probesetID")) and then stick these values in the pm and mm slots of the
> ProbeSet. This is done in computeExprSet().
>
> If it sounds like a lot of work, it is. If you really need these data,
> you should look at the functions pmcorrect.mas(), and computeExprSet()
> to see how it is done and to figure out how to do it yourself.
>
> To see the functions you can either download the sources from BioC and
> grep them out of the R directory (preferable), or you can just load the
> affy package and type
>
> pmcorrect.mas
>
> and
>
> showMethods(computeExprSet, class = "AffyBatch", includeDefs = TRUE)
>
> Best,
>
> Jim
>
>
>
>>Thanks a lot for your help.
>>
>>Best regards
>>
>>YH Wang
>>
>>_______________________________________________
>>Bioconductor mailing list
>>Bioconductor at stat.math.ethz.ch
>>https://stat.ethz.ch/mailman/listinfo/bioconductor
>>Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>
>
>
--
James W. MacDonald, M.S.
Biostatistician
Affymetrix and cDNA Microarray Core
University of Michigan Cancer Center
1500 E. Medical Center Drive
7410 CCGC
Ann Arbor MI 48109
734-647-5623
**********************************************************
Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues.
More information about the Bioconductor
mailing list