# Author Sucheta

# Has 2 different versions of background correction using GCRMA

# The first one returns the probeset values and the second one gets the

# probe values

#SCRIPT -1

#----------

#-----------

library(gcrma)

celpath <- "YOUR PATH"

eset <- justGCRMA(celfile.path=celpath, normalize=FALSE)

write.exprs(eset, "out.txt");

#-------- This would return the background corrected probeset values ---------------------------

#SCRIPT-2

#--------

#--------

library(gcrma);

# Get your affinity data else the program can automatically compute affinity

# Time consuming

affinityData <- "YOUR AFFINITY DATA";

load(affinityData);

# Compress CEL files

compress <- getOption("BioC")$affy$compress.cel

fast <- TRUE;

cdfName <- "Soybean";

celPath <- "YOUR PATH";

l <- AllButCelsForReadAffy(celfile.path = celPath);

fileNames <- l$filenames;

# Get file name length

n <- length(fileNames);

# Stop if file name length is zero

if(n == 0)

stop("File name Length is zero EXITING ...");

# Calling the function

pms <- mem.bkg(filenames=fileNames,

pm.affinities = pm.affinities,

mm.affinities = mm.affinities,

index.affinities = index.affinities,

type = "fullmodel",

optical.correct = TRUE,

verbose = TRUE,

minimum = 1,

k = 6 * fast+0.5*(1-fast),

rho = 0.7,

correction = 1,

stretch = 1.15*fast+1*(1-fast),

fast = TRUE,

cdfname = cdfName)

output <- paste(celPath,"_bg.txt", sep="");

write.table(pms, file=output, row.names=FALSE, col.names=TRUE, sep="\t");

# ---- FUNCTION DEFINITION --------------#

## A function to do background correction using less RAM

mem.bkg <- function(filenames, pm.affinities, mm.affinities,

index.affinities, type, minimum, optical.correct,

verbose, k, rho, correction, stretch, fast, cdfname){

pms <- read.probematrix(filenames=filenames, which="pm", cdfname = cdfname)$pm

## tmps used to carry optical correct value to bkg correction loop

if(optical.correct){

if(verbose) cat("Adjusting for optical effect.")

tmps <- NULL

for (i in 1:ncol(pms)){

if(verbose) cat(".")

mm <- read.probematrix(filenames=filenames[i], which="mm")$mm[,1]

tmp <- min(c(pms[,i], mm), na.rm=TRUE)

pms[,i] <- pms[,i]- tmp + minimum

tmps <- c(tmps, tmp)

}

}

if(verbose) cat("Done.\n")

if(type=="fullmodel" | type=="affinities"){

set.seed(1)

Subset <- sample(1:length(pms[index.affinities,]),25000)

y <- log2(pms)[index.affinities,][Subset]

Subset <- (Subset-1)%%nrow(pms[index.affinities,])+1

x <- pm.affinities[Subset]

fit1 <- lm(y~x)

}

if(verbose) cat("Adjusting for non-specific binding")

for(i in 1:ncol(pms)){

if(verbose) cat(".")

mm <- read.probematrix(filenames=filenames[i], which="mm")$mm[,1]

if(optical.correct)

mm <- mm - tmps[i] + minimum

if(type=="fullmodel"){

pms[,i] <- bg.adjust.fullmodel(pms[,i],mm,ncs=NULL,

pm.affinities,mm.affinities,anc=NULL,

index.affinities,k=k,

rho=rho,fast=fast)

pms[index.affinities,i] <- 2^(log2(pms[index.affinities,i])-

fit1$coef[2]*pm.affinities+mean(fit1$coef[2]*pm.affinities))

}

if(type=="affinities"){

pms[,i] <- bg.adjust.affinities(pms[,i],ncs=mm,pm.affinities,mm.affinities,

index.affinities, k=k,

fast=fast)

pms[index.affinities,i] <- 2^(log2(pms[index.affinities,i])-

fit1$coef[2]*pm.affinities +

mean(fit1$coef[2]*pm.affinities))

}

if(type=="mm") pms[,i] <- bg.adjust.mm(pms[,i],correction*mm,k=k,fast=fast)

if(type=="constant"){

pms[,i] <- bg.adjust.constant(pms[,i],

k=k,Q=correction*mean(pms[,i]< mm ),

fast=fast)

}

if(stretch!=1){

mu <- mean(log(pms[,i]))

pms[,i] <- exp(mu + stretch*(log(pms[,i])- mu))

}

}

if(verbose) cat("Done.\n")

return(pms)

}

## 2 comments:

Hello tsucheta,

How can I write you an e-mail? Thank you,

IĆ±igo

I don't think you need to write to me, just post here and I will reply back.

Thanks

Sucheta

Post a Comment