Hi Fraser, the bump at 0.5 stems from probe sets with only one probe. For those is modeling not required ;-) I would recommend not to use them, as they their expression values are probably not reliable. This way you get rip off them: Cheers, Okko SNR <- apply(eset at assayData$se.exprs,1,min) tmp <- which(SNR==0.5) eset2 <- eset[-tmp] inis <- INIcalls(eset) summary(inis) plot(inis) -- Djork-Arn? Clevert, PhD Phone: +49 30 6883 5306 Fax: +49 30 6883 5307 Email: okko at clevert.de Am 15.11.2013 um 16:15 schrieb Fraser Sim <fjsim at="" buffalo.edu="">: > Hi Okko, > > Here?s the plot(INIs) using the latest code you sent and including the ?normalize = TRUE? in the call to farms. > > I now have 4919 features labeled as informative. > > Do you think this is ok? There appears to be strange albeit small bump at 0.5, which I?ve never seen in this type of plot before with Affy&Farms on U133 or similar 3?arrays. > > Cheers, > Fraser > > <image001.png> > > > Fraser Sim, PhD > Assistant Professor of Pharmacology and Toxicology > University at Buffalo > 3435 Main Street > 119 Farber Hall > Buffalo, NY 14214 > > T: (716) 829-2151 > fjsim at buffalo.edu > > Confidentiality Notice: This message, including any attachments, is for the sole use of the intended recipient(s) and may contain confidential and privileged information. Any unauthorized review, use, disclosure, or distribution is strictly prohibited. If you are not the intended recipient, please contact the sender by reply e-mail and destroy all copies of the original message. > > > -----Original Message----- > From: Djork-Arn? Clevert [mailto:okko at clevert.de] > Sent: Friday, November 15, 2013 9:46 AM > To: Fraser Sim > Cc: bioconductor at r-project.org > Subject: Re: [BioC] Farms/Call.INI using oligo package? > > Hi Fraser, > > I just realized that my previous email was everything else than ?self-explanatory? and on top of that - the workaround contained a typo? > > Cheers, > Okko > > Please find below the revised code including an example at the end: > > library(oligo) > library(farms) > > farms <- function(object, background=FALSE, normalize=TRUE, subset=NULL, target="core", weight=.5, > mu=0, weighted.mean=FALSE, laplacian=FALSE, robust=TRUE, > correction=0, centering=c("median","mean"), spuriousCorrelation=0, verbose=TRUE ) > { > target <- match.arg(target, c("core", "probeset")) > featureInfo <- oligo:::stArrayPmInfo(object, target=target) > theClass <- class(exprs(object)) > pmi <- featureInfo[["fid"]] > pnVec <- as.character(featureInfo[["fsetid"]]) > if ("matrix" %in% theClass){ > pms <- exprs(object)[pmi,, drop=FALSE] > dimnames(pms) <- NULL > colnames(pms) <- sampleNames(object) > theExprs <- basicFARMS(pms, pnVec, normalize, background, weight=weight, mu=mu, weighted.mean=weighted.mean, > laplacian=laplacian, robust=robust, correction=correction, centering=centering, spuriousCorrelation=spuriousCorrelation) > rm(pms) > } > else{ > stop("farms is currently not implemented for '", theClass, "' objects.") > } > > > colnames(theExprs$exprs) <- colnames(theExprs$se.exprs) <- sampleNames(object) > rownames(theExprs$exprs) <- rownames(theExprs$se.exprs) <- unique(pnVec) > > out <- new("ExpressionSet", > phenoData=phenoData(object), > experimentData=experimentData(object), > exprs=theExprs$exprs, > se.exprs=theExprs$se.exprs, > annotation=annotation(object), > protocolData=protocolData(object)) > > if (validObject(out)){ > return(out) > }else{ > stop("Resulting object is invalid.") > } > } > > > > > basicFARMS <- function(pmMat, pnVec, normalize, background, weight, mu, weighted.mean, > laplacian, robust, correction, centering, spuriousCorrelation){ > pns <- unique(pnVec) > nPn <- length(unique(pnVec)) > pnVec <- split(1:(length(pnVec)), pnVec) > narray <- ncol(pmMat) > if(background){ > pmMat <- oligo:::backgroundCorrect(pmMat) > } > if(normalize){ > pmMat <- oligo:::normalize(pmMat) > } > ini.mat <- matrix(0.5, nPn, narray ) > exprs.mat <- matrix(0, nPn, narray ) > message("Summarizing... ", appendLF=FALSE) > for(i in 1:length(pnVec)){ > if(length(pnVec[[i]])<2){ > exprs.mat[i,] <- log2(pmMat[pnVec[[i]],]) > } > else{ > res <- generateExprVal.method.farms(probes=pmMat[pnVec[[i]],], weight=weight, mu=mu, cyc=30, tol=0.00001, > weighted.mean=weighted.mean, robust=robust, minNoise=minNoise, correction=correction, > laplacian=laplacian, centering=centering, spuriousCorrelation=spuriousCorrelation) > exprs.mat[i,] <- res$exprs > ini.mat[i,] <- res$se.exprs > } > } > message("OK") > return(list(exprs=exprs.mat, se.exprs=ini.mat)) } > > > > ### example to run > # library(oligo) > # library(farms) > # celFiles <- list.celfiles("/path/to/celfiles/", full.name=TRUE) # raw <- read.celfiles(celFiles) # eset <- farms(raw, normalize=TRUE) #### > > > > > > -- > Djork-Arn? Clevert, PhD > > Phone: +49 30 6883 5306 > Fax: +49 30 6883 5307 > Email: okko at clevert.de > > Am 13.11.2013 um 14:32 schrieb Fraser Sim <fjsim at="" buffalo.edu="">: > > > Hi, > > > > > > > > I am analyzing Affymetrix hugene2.0st arrays using the oligo package and had > > hoped to be able to use informative/non-informative probe calling (FARMS, > > Tolloen et 2007) as a means on filtering data. I had used this on older 3' > > biased arrays analyzed with affy but the farms package only works with an > > affyBatch object. Is this possible with another package or should I be > > considering a different approach to filter array data on the newer Affy > > arrays? > > > > > > > > Cheers, > > > > Fraser > > > > > > [[alternative HTML version deleted]] > > > > _______________________________________________ > > Bioconductor mailing list > > Bioconductor at r-project.org > > https://stat.ethz.ch/mailman/listinfo/bioconductor > > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor