[Bio] / FigKernelScripts / EcoliRMA.R Repository:
ViewVC logotype

View of /FigKernelScripts/EcoliRMA.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (download) (annotate)
Fri Apr 30 20:20:43 2010 UTC (9 years, 6 months ago) by overbeek
Branch: MAIN
CVS Tags: mgrast_dev_08112011, mgrast_dev_08022011, rast_rel_2014_0912, myrast_rel40, mgrast_dev_05262011, mgrast_dev_04082011, rast_rel_2010_0928, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, rast_rel_2014_0729, mgrast_dev_02212011, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2011_0119, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, mgrast_dev_04012011, rast_rel_2010_0827, myrast_33, rast_rel_2011_0928, mgrast_dev_04052011, mgrast_dev_02222011, mgrast_dev_10262011, HEAD
Changes since 1.1: +2 -2 lines
from Matt

args = commandArgs(TRUE);
if (length(args) != 2) { print("usage:  Rscript EcoliRMA.R peg_probe_table experiment_dir"); q() }
peg_probe_tbl = args[1];
experiments_dir = args[2];

library(affy)

# read old Shewanella CDF environment for chip dimensions
library(ecoli2cdf)

# create the new Shewanella CDF environment from Ross's table
newshewcdf = new.env(hash=TRUE);
probesets = readLines(peg_probe_tbl);
probesets = strsplit(probesets, "\t");

for (i in 1:length(probesets))
{
	probeset = probesets[[i]];
	peg = probeset[1]
	x_y = strsplit(probeset[2], "_")
	index = xy2indices(as.numeric(x_y[[1]][1]),
					   as.numeric(x_y[[1]][2]),cdf="ecoli2cdf")
	
	if (! exists(peg, envir= newshewcdf)) {
		# we're not assigning a mismatch index for now
		pm_mm = matrix(c(index,0), nrow=1,ncol=2, dimnames=list(c(), c("pm", "mm")))
		assign(peg,pm_mm,envir=newshewcdf)
	}
	else {
		pm_mm = get(peg,envir=newshewcdf)
		pm_mm = rbind(pm_mm,c(index,0))
		assign(peg,pm_mm,envir=newshewcdf)
	}
}

# Read the Data and perform RMA
Data = ReadAffy(celfile.path=experiments_dir, cdfname="newshewcdf", compress=TRUE, verbose=TRUE)
eset = rma(Data)
geneexp=exprs(eset)
write.table(geneexp, file="raw_data.tab", row.names=TRUE, col.names=TRUE, sep="\t", quote=FALSE)

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3