#for analysis of MIT MIT data
rm(list=ls())
PMIT_BATCH <-function(directory, exportname="", export=FALSE) {
dir <- readLines(pipe(paste("ls", directory)))
paths <- paste(directory, dir, sep="/")
FILENUM <- 0
#non_data <- 0
#cv <- function(x) 100*((sqrt(var(x)/length(x)))/mean(x)) #define CV
non_baso <- 0

cat("processing . . .\n")

for (step in 1:length(paths)) {
	subdir <- readLines(pipe(paste("ls", paths[step])))
	for (instep in 1:length(subdir)) {
		filename <- paste(paths[step], subdir[instep], sep="/")
		shortname <- subdir[instep]; corename <- substring(shortname, 1, 6)
		visit_num <- substring(shortname, 8, 8)
		baso_test <- substring(filename, nchar(filename)-4, nchar(filename))
		boo <- any(as.logical(grep(baso_test, "-CD63")))
		if (boo) {
			if (FILENUM == 0) FILENUM<-1 else FILENUM <- 2 # used below for define of data_array
			con <- file(filename, encoding="macroman")
			TABLE <- read.table(con, skip=1, sep="\t", check.names=FALSE, fill=TRUE, strip.white=TRUE)
			con <- file(filename, encoding="macroman")
			ROWS <- length(TABLE[[1]])		#GET TABLE LENGTH TO LEAVE OFF LAST TWO ROWS
			TABLE <- read.table(con, skip=1, sep="\t", check.names=FALSE, fill=TRUE, strip.white=TRUE, nrows=ROWS-2)
			cat(paste(filename, "\n"))
			# remove additional text from stim field
			TABLE[,1] <- substring(TABLE[,1],4,15)
			TABLE[,1] <- substring(TABLE[,1],1,length(TABLE[,1])-1)
			TABLE[,1] <- gsub(" ","",TABLE[,1])

			data_sum <- data.frame(experiment=corename, visit=visit_num, stim=TABLE[,1], count=TABLE[,2], baso_count=TABLE[,3], baso_CD69_mfi=TABLE[,4], baso_CD69_CV=TABLE[,5], baso_CD203_mfi=TABLE[,6], baso_CD203_CV=TABLE[,7], percent63low=TABLE[,8], count63low=TABLE[,9], CD63neg_63MFI=TABLE[,10], CD63neg_63CV=TABLE[,11], CD63neg_203MFI=TABLE[,12], CD63neg_203CV=TABLE[,13], CD63neg_123MFI=TABLE[,14], CD63neg_123CV=TABLE[,15], percent63hi=TABLE[,16], count63hi=TABLE[,17], CD63pos_63MFI=TABLE[,18], CD63pos_63CV=TABLE[,19], CD63pos_203MFI=TABLE[,20], CD63pos_203CV=TABLE[,21], CD63pos_123MFI=TABLE[,22], CD63pos_123CV=TABLE[,23], pDC_percent63hi=TABLE[,24], pDC_count=TABLE[,25], pDC_123MFI=TABLE[,26], pDC_123CV=TABLE[,27], pDC_CD69_MFI=TABLE[,28], pDC_CD69_CV=TABLE[,29],CD123posDRhi=TABLE[,32], CD123posDRlo=TABLE[,33])
				if (FILENUM==1) {data_array <<- data_sum}
				else {data_array <<- rbind(data_array, data_sum)}			}
		else non_baso <- non_baso + 1	
	}

}
cat(paste("total number of folders:", paste(step, "\n")))

#remove funny character from flowjo files
for (x in 4:31) {
	data_array[,x] <- as.numeric(gsub("•","NA",data_array[,x]))
	}


if (export) {
	exportpath <- paste(readLines(pipe("pwd")),"/",sep="")
	con <- file(paste(exportpath, exportname, sep = ""), encoding="macroman")
	write.csv(data_array, con)
	cat(paste("writing file", paste(exportpath, paste(exportname, "\n"), sep="")))
	}
}

ARK_BATCH <-function(directory, exportname="", export=FALSE) {
dir <- readLines(pipe(paste("ls", directory)))
paths <- paste(directory, dir, sep="/")
FILENUM <- 0
#non_data <- 0
#cv <- function(x) 100*((sqrt(var(x)/length(x)))/mean(x)) #define CV
non_baso <- 0

cat("processing . . .\n")

for (step in 1:length(paths)) {
	subdir <- readLines(pipe(paste("ls", paths[step])))
	for (instep in 1:length(subdir)) {
		filename <- paste(paths[step], subdir[instep], sep="/")
		shortname <- subdir[instep]; corename <- substring(shortname, 1, 4)
		visit_num <- substring(shortname, 6, 6)
		baso_test <- substring(filename, nchar(filename)-4, nchar(filename))
		boo <- any(as.logical(grep(baso_test, "-CD63")))
		if (boo) {
			if (FILENUM == 0) FILENUM<-1 else FILENUM <- 2 # used below for define of data_array
			con <- file(filename, encoding="macroman")
			TABLE <- read.table(con, skip=1, sep="\t", check.names=FALSE, fill=TRUE, strip.white=TRUE)
			con <- file(filename, encoding="macroman")
			ROWS <- length(TABLE[[1]])		#GET TABLE LENGTH TO LEAVE OFF LAST TWO ROWS
			TABLE <- read.table(con, skip=1, sep="\t", check.names=FALSE, fill=TRUE, strip.white=TRUE, nrows=ROWS-2)
			cat(paste(filename, "\n"))
			# remove additional text from stim field
			TABLE[,1] <- substring(TABLE[,1],4,15)
			TABLE[,1] <- substring(TABLE[,1],1,length(TABLE[,1])-1)
			TABLE[,1] <- gsub(" ","",TABLE[,1])

			data_sum <- data.frame(experiment=corename, visit=visit_num, stim=TABLE[,1], count=TABLE[,2], baso_count=TABLE[,3], baso_CD69_mfi=TABLE[,4], baso_CD69_CV=TABLE[,5], baso_CD203_mfi=TABLE[,6], baso_CD203_CV=TABLE[,7], percent63low=TABLE[,8], count63low=TABLE[,9], CD63neg_63MFI=TABLE[,10], CD63neg_63CV=TABLE[,11], CD63neg_203MFI=TABLE[,12], CD63neg_203CV=TABLE[,13], CD63neg_123MFI=TABLE[,14], CD63neg_123CV=TABLE[,15], percent63hi=TABLE[,16], count63hi=TABLE[,17], CD63pos_63MFI=TABLE[,18], CD63pos_63CV=TABLE[,19], CD63pos_203MFI=TABLE[,20], CD63pos_203CV=TABLE[,21], CD63pos_123MFI=TABLE[,22], CD63pos_123CV=TABLE[,23], pDC_percent63hi=TABLE[,24], pDC_count=TABLE[,25], pDC_123MFI=TABLE[,26], pDC_123CV=TABLE[,27], pDC_CD69_MFI=TABLE[,28], pDC_CD69_CV=TABLE[,29],CD123posDRhi=TABLE[,32], CD123posDRlo=TABLE[,33])
				if (FILENUM==1) {data_array <<- data_sum}
				else {data_array <<- rbind(data_array, data_sum)}			}
		else non_baso <- non_baso + 1	
	}

}
cat(paste("total number of folders:", paste(step, "\n")))

#remove funny character from flowjo files
for (x in 4:31) {
	data_array[,x] <- as.numeric(gsub("•","NA",data_array[,x]))
	}


if (export) {
	exportpath <- paste(readLines(pipe("pwd")),"/",sep="")
	con <- file(paste(exportpath, exportname, sep = ""), encoding="macroman")
	write.csv(data_array, con)
	cat(paste("writing file", paste(exportpath, paste(exportname, "\n"), sep="")))
	}
}

ARK_BATCH("/Users/Guest/Public/Shreff_lab/Duke_peanut/Out_box_analysis2/Arkansas_analysis2", "Arkansas.csv", TRUE)
PMIT_BATCH("/Users/Guest/Public/Shreff_lab/Duke_peanut/Out_box_analysis2/PMIT_analysis2", "Duke.csv", TRUE)
# after running batch having saved each data_array as indicated below:
rm(list=ls())
mit <- read.csv("Arkansas.csv", as.is=TRUE) # note that as.is prevents read.table from creating factors
mit2 <- read.csv("Duke.csv", as.is=TRUE)
MIT <- rbind(mit, mit2)

# extensive cleanup needed

i <- grep("Peanut", MIT$stim, TRUE)
MIT$stim[i] <- substring(MIT$stim[i],1,7)
MIT$stim[i] <- gsub("Peanut", "PN", MIT$stim[i])

i <- grep("ige", MIT$stim, TRUE)
MIT$stim[i] <- "aIgE"

i <- grep("fmlp", MIT$stim, TRUE)
MIT$stim[i] <- "fMLP"

i <- grep("eggwhite", MIT$stim, TRUE)
MIT$stim[i] <- substring(gsub("EggWhite", "EGG", MIT$stim[i]), 1, 4)
i <- grep("egg", MIT$stim, TRUE)
MIT$stim[i] <- toupper(substring(MIT$stim[i],1,4))
MIT$stim <- gsub("EW4", "EGG1", MIT$stim)
MIT$stim <- gsub("EW3", "EGG2", MIT$stim)
MIT$stim <- gsub("EW2", "EGG3", MIT$stim)
MIT$stim <- gsub("EW1", "EGG4", MIT$stim)


i <- grep("IL-3", MIT$stim, TRUE)
MIT$stim[i] <- "IL3"

i <- grep("PMI", MIT$stim, TRUE)
MIT$stim[i] <- "RPMI"

MIT$stim <- gsub("PN0", "PN", MIT$stim)
MIT$stim <- gsub("PNE", "PN", MIT$stim)

i <- grep("PN", MIT$stim)
MIT$stim[i] <- substring(MIT$stim[i], 1, 3)

MIT$stim <- gsub("E2", "PN1", MIT$stim)
MIT$stim <- gsub("E3", "PN2", MIT$stim)
MIT$stim <- gsub("E4", "PN3", MIT$stim)
MIT$stim <- gsub("E5", "PN4", MIT$stim)
MIT$stim <- gsub("E6", "PN5", MIT$stim)
MIT$stim <- gsub("E7", "PN6", MIT$stim)

MIT$stim <- gsub("P1", "PN1", MIT$stim)
MIT$stim <- gsub("P2", "PN2", MIT$stim)
MIT$stim <- gsub("P3", "PN3", MIT$stim)
MIT$stim <- gsub("P4", "PN4", MIT$stim)
MIT$stim <- gsub("P5", "PN5", MIT$stim)
MIT$stim <- gsub("P6", "PN6", MIT$stim)

# convert to factors
MIT$stim <- as.factor(MIT$stim)
MIT$experiment <- as.factor(MIT$experiment)
MIT$visit <- as.factor(MIT$visit)
	
# NEED TO MASH UP WITH CLINICAL DATA TO ANALYZE BY TREATMENT GROUP WITH TIME AND DOSE

clinical <- read.csv("MIT_clinical.csv")
# add variables that I need to MIT
MIT <- data.frame(MIT, treatment = NA, day = NA, dose = NA)

# loop through MIT_clinical by PatientID and VISIT to extract

for (step in 1:length(clinical$PatientID)) {
	i <- MIT$experiment %in% clinical$PatientID[step] & MIT$visit %in% clinical$VISIT[step]
	MIT$treatment[i] <- as.character(clinical$TREATMENT[step])
	MIT$day[i] <- clinical$DAYS[step]
	MIT$dose[i] <- clinical$DOSE[step]
	}

MIT$treatment <- as.factor(MIT$treatment)


# NOW WE CAN PRODUCE SOME PLOTS!!

#-------Start with Peanut-------------
quartz(,8,6); lo<-layout(matrix(c(1:8),2,4, byrow=TRUE)); par(cex.lab=1.4, cex.axis=2)

treated <- subset(MIT, treatment != "PLACEBO")

#PN1-------------------------------------
plot(subset(treated, stim == "PN1")$day, subset(treated, stim == "PN1")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("PN1", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "PN1")$experiment))) {
	lines(subset(subset(treated, stim == "PN1"), experiment == levels(subset(treated, stim == "PN1")$experiment)[step])$day, subset(subset(treated, stim == "PN1"), experiment == levels(subset(treated, stim == "PN1")$experiment)[step])$percent63hi, col="green")
	}

#PN2--------------------------------------
plot(subset(treated, stim == "PN2")$day, subset(treated, stim == "PN2")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("PN2", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "PN2")$experiment))) {
	lines(subset(subset(treated, stim == "PN2"), experiment == levels(subset(treated, stim == "PN2")$experiment)[step])$day, subset(subset(treated, stim == "PN2"), experiment == levels(subset(treated, stim == "PN2")$experiment)[step])$percent63hi, col="green")
	}

#PN3--------------------------------------
plot(subset(treated, stim == "PN3")$day, subset(treated, stim == "PN3")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("PN3", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "PN3")$experiment))) {
	lines(subset(subset(treated, stim == "PN3"), experiment == levels(subset(treated, stim == "PN3")$experiment)[step])$day, subset(subset(treated, stim == "PN3"), experiment == levels(subset(treated, stim == "PN3")$experiment)[step])$percent63hi, col="green")
	}

#PN4--------------------------------------
plot(subset(treated, stim == "PN4")$day, subset(treated, stim == "PN4")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("PN4", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "PN4")$experiment))) {
	lines(subset(subset(treated, stim == "PN4"), experiment == levels(subset(treated, stim == "PN4")$experiment)[step])$day, subset(subset(treated, stim == "PN4"), experiment == levels(subset(treated, stim == "PN4")$experiment)[step])$percent63hi, col="green")
	}

#PLACEBO-----------------------------------


placebo <- subset(MIT, treatment != "PEANUT")

#PN1-------------------------------------
plot(subset(placebo, stim == "PN1")$day, subset(placebo, stim == "PN1")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("PN1", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "PN1")$experiment))) {
	lines(subset(subset(placebo, stim == "PN1"), experiment == levels(subset(placebo, stim == "PN1")$experiment)[step])$day, subset(subset(placebo, stim == "PN1"), experiment == levels(subset(placebo, stim == "PN1")$experiment)[step])$percent63hi, col="red")
	}

#PN2--------------------------------------
plot(subset(placebo, stim == "PN2")$day, subset(placebo, stim == "PN2")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("PN2", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "PN2")$experiment))) {
	lines(subset(subset(placebo, stim == "PN2"), experiment == levels(subset(placebo, stim == "PN2")$experiment)[step])$day, subset(subset(placebo, stim == "PN2"), experiment == levels(subset(placebo, stim == "PN2")$experiment)[step])$percent63hi, col="red")
	}

#PN3--------------------------------------
plot(subset(placebo, stim == "PN3")$day, subset(placebo, stim == "PN3")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("PN3", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "PN3")$experiment))) {
	lines(subset(subset(placebo, stim == "PN3"), experiment == levels(subset(placebo, stim == "PN3")$experiment)[step])$day, subset(subset(placebo, stim == "PN3"), experiment == levels(subset(placebo, stim == "PN3")$experiment)[step])$percent63hi, col="red")
	}

#PN4--------------------------------------
plot(subset(placebo, stim == "PN4")$day, subset(placebo, stim == "PN4")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("PN4", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "PN4")$experiment))) {
	lines(subset(subset(placebo, stim == "PN4"), experiment == levels(subset(placebo, stim == "PN4")$experiment)[step])$day, subset(subset(placebo, stim == "PN4"), experiment == levels(subset(placebo, stim == "PN4")$experiment)[step])$percent63hi, col="red")
	}

# Controls
#-----------------------------------
quartz(,8,6); lo<-layout(matrix(c(1:8),2,4, byrow=TRUE)); par(cex.lab=1.4, cex.axis=2)

treated <- subset(MIT, treatment != "PLACEBO")

#RPMI-------------------------------------
plot(subset(treated, stim == "RPMI")$day, subset(treated, stim == "RPMI")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("RPMI", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "RPMI")$experiment))) {
	lines(subset(subset(treated, stim == "RPMI"), experiment == levels(subset(treated, stim == "RPMI")$experiment)[step])$day, subset(subset(treated, stim == "RPMI"), experiment == levels(subset(treated, stim == "RPMI")$experiment)[step])$percent63hi, col="green")
	}

#IL3--------------------------------------
plot(subset(treated, stim == "IL3")$day, subset(treated, stim == "IL3")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("IL-3", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "IL3")$experiment))) {
	lines(subset(subset(treated, stim == "IL3"), experiment == levels(subset(treated, stim == "IL3")$experiment)[step])$day, subset(subset(treated, stim == "IL3"), experiment == levels(subset(treated, stim == "IL3")$experiment)[step])$percent63hi, col="green")
	}

#ANTI-IGE--------------------------------------
plot(subset(treated, stim == "aIgE")$day, subset(treated, stim == "aIgE")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("α-IgE", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "aIgE")$experiment))) {
	lines(subset(subset(treated, stim == "aIgE"), experiment == levels(subset(treated, stim == "aIgE")$experiment)[step])$day, subset(subset(treated, stim == "aIgE"), experiment == levels(subset(treated, stim == "aIgE")$experiment)[step])$percent63hi, col="green")
	}

#fMLP--------------------------------------
plot(subset(treated, stim == "fMLP")$day, subset(treated, stim == "fMLP")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("fMLP", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "fMLP")$experiment))) {
	lines(subset(subset(treated, stim == "fMLP"), experiment == levels(subset(treated, stim == "fMLP")$experiment)[step])$day, subset(subset(treated, stim == "fMLP"), experiment == levels(subset(treated, stim == "fMLP")$experiment)[step])$percent63hi, col="green")
	}

#RPMI-------------------------------------
plot(subset(placebo, stim == "RPMI")$day, subset(placebo, stim == "RPMI")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("RPMI", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "RPMI")$experiment))) {
	lines(subset(subset(placebo, stim == "RPMI"), experiment == levels(subset(placebo, stim == "RPMI")$experiment)[step])$day, subset(subset(placebo, stim == "RPMI"), experiment == levels(subset(placebo, stim == "RPMI")$experiment)[step])$percent63hi, col="red")
	}

#IL3--------------------------------------
plot(subset(placebo, stim == "IL3")$day, subset(placebo, stim == "IL3")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("IL-3", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "IL3")$experiment))) {
	lines(subset(subset(placebo, stim == "IL3"), experiment == levels(subset(placebo, stim == "IL3")$experiment)[step])$day, subset(subset(placebo, stim == "IL3"), experiment == levels(subset(placebo, stim == "IL3")$experiment)[step])$percent63hi, col="red")
	}

#ANTI-IGE--------------------------------------
plot(subset(placebo, stim == "aIgE")$day, subset(placebo, stim == "aIgE")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("α-IgE", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "aIgE")$experiment))) {
	lines(subset(subset(placebo, stim == "aIgE"), experiment == levels(subset(placebo, stim == "aIgE")$experiment)[step])$day, subset(subset(placebo, stim == "aIgE"), experiment == levels(subset(placebo, stim == "aIgE")$experiment)[step])$percent63hi, col="red")
	}

#fMLP--------------------------------------
plot(subset(placebo, stim == "fMLP")$day, subset(placebo, stim == "fMLP")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("fMLP", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "fMLP")$experiment))) {
	lines(subset(subset(placebo, stim == "fMLP"), experiment == levels(subset(placebo, stim == "fMLP")$experiment)[step])$day, subset(subset(placebo, stim == "fMLP"), experiment == levels(subset(placebo, stim == "fMLP")$experiment)[step])$percent63hi, col="red")
	}

#------Now Egg

quartz(,8,6); lo<-layout(matrix(c(1:8),2,4, byrow=TRUE)); par(cex.lab=1.4, cex.axis=2)

treated <- subset(MIT, treatment != "PLACEBO")

#EGG1-------------------------------------
plot(subset(treated, stim == "EGG1")$day, subset(treated, stim == "EGG1")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("EGG1", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "EGG1")$experiment))) {
	lines(subset(subset(treated, stim == "EGG1"), experiment == levels(subset(treated, stim == "EGG1")$experiment)[step])$day, subset(subset(treated, stim == "EGG1"), experiment == levels(subset(treated, stim == "EGG1")$experiment)[step])$percent63hi, col="green")
	}

#EGG2--------------------------------------
plot(subset(treated, stim == "EGG2")$day, subset(treated, stim == "EGG2")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("EGG2", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "EGG2")$experiment))) {
	lines(subset(subset(treated, stim == "EGG2"), experiment == levels(subset(treated, stim == "EGG2")$experiment)[step])$day, subset(subset(treated, stim == "EGG2"), experiment == levels(subset(treated, stim == "EGG2")$experiment)[step])$percent63hi, col="green")
	}

#EGG3--------------------------------------
plot(subset(treated, stim == "EGG3")$day, subset(treated, stim == "EGG3")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("EGG3", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "EGG3")$experiment))) {
	lines(subset(subset(treated, stim == "EGG3"), experiment == levels(subset(treated, stim == "EGG3")$experiment)[step])$day, subset(subset(treated, stim == "EGG3"), experiment == levels(subset(treated, stim == "EGG3")$experiment)[step])$percent63hi, col="green")
	}

#EGG4--------------------------------------
plot(subset(treated, stim == "EGG4")$day, subset(treated, stim == "EGG4")$percent63hi, main="PEANUT", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("EGG4", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(treated, stim == "EGG4")$experiment))) {
	lines(subset(subset(treated, stim == "EGG4"), experiment == levels(subset(treated, stim == "EGG4")$experiment)[step])$day, subset(subset(treated, stim == "EGG4"), experiment == levels(subset(treated, stim == "EGG4")$experiment)[step])$percent63hi, col="green")
	}

#PLACEBO-----------------------------------


placebo <- subset(MIT, treatment != "PEANUT")

#EGG1-------------------------------------
plot(subset(placebo, stim == "EGG1")$day, subset(placebo, stim == "EGG1")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("EGG1", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "EGG1")$experiment))) {
	lines(subset(subset(placebo, stim == "EGG1"), experiment == levels(subset(placebo, stim == "EGG1")$experiment)[step])$day, subset(subset(placebo, stim == "EGG1"), experiment == levels(subset(placebo, stim == "EGG1")$experiment)[step])$percent63hi, col="red")
	}

#EGG2--------------------------------------
plot(subset(placebo, stim == "EGG2")$day, subset(placebo, stim == "EGG2")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("EGG2", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "EGG2")$experiment))) {
	lines(subset(subset(placebo, stim == "EGG2"), experiment == levels(subset(placebo, stim == "EGG2")$experiment)[step])$day, subset(subset(placebo, stim == "EGG2"), experiment == levels(subset(placebo, stim == "EGG2")$experiment)[step])$percent63hi, col="red")
	}

#EGG3--------------------------------------
plot(subset(placebo, stim == "EGG3")$day, subset(placebo, stim == "EGG3")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("EGG3", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "EGG3")$experiment))) {
	lines(subset(subset(placebo, stim == "EGG3"), experiment == levels(subset(placebo, stim == "EGG3")$experiment)[step])$day, subset(subset(placebo, stim == "EGG3"), experiment == levels(subset(placebo, stim == "EGG3")$experiment)[step])$percent63hi, col="red")
	}

#EGG4--------------------------------------
plot(subset(placebo, stim == "EGG4")$day, subset(placebo, stim == "EGG4")$percent63hi, main="PLACEBO", xlab="time (days)", ylab="%CD63 bright", ylim=c(0,100)); mtext("EGG4", 1, line=6, cex=1.5, at=0)

for (step in 1:length(levels(subset(placebo, stim == "EGG4")$experiment))) {
	lines(subset(subset(placebo, stim == "EGG4"), experiment == levels(subset(placebo, stim == "EGG4")$experiment)[step])$day, subset(subset(placebo, stim == "EGG4"), experiment == levels(subset(placebo, stim == "EGG4")$experiment)[step])$percent63hi, col="red")
	}


#-----------------------------------------
evivo_pl <- subset(MIT, stim == "RPMI" & treatment != "PEANUT")
evivo_pn <- subset(MIT, stim == "RPMI" & treatment != "PLACEBO")

quartz()
plot(evivo_pn$day, evivo_pn$baso_CD203_mfi)
for (step in 1:length(levels(evivo_pn$experiment))) {
	lines(subset(evivo_pn, experiment == levels(evivo_pn$experiment)[step])$day,  subset(evivo_pn, experiment == levels(evivo_pn$experiment)[step])$baso_CD203_mfi, col="blue")
	}

quartz()
plot(evivo_pl$day, evivo_pl$baso_CD203_mfi)
for (step in 1:length(levels(evivo_pl$experiment))) {
	lines(subset(evivo_pl, experiment == levels(evivo_pl$experiment)[step])$day,  subset(evivo_pl, experiment == levels(evivo_pl$experiment)[step])$baso_CD203_mfi, col="blue")
	}

quartz()
plot(evivo_pn$day, evivo_pn$baso_CD69_mfi)
for (step in 1:length(levels(evivo_pn$experiment))) {
	lines(subset(evivo_pn, experiment == levels(evivo_pn$experiment)[step])$day,  subset(evivo_pn, experiment == levels(evivo_pn$experiment)[step])$baso_CD69_mfi, col="blue")
	}
	
quartz()
plot(evivo_pl$day, evivo_pl$baso_CD69_mfi)
for (step in 1:length(levels(evivo_pl$experiment))) {
	lines(subset(evivo_pl, experiment == levels(evivo_pl$experiment)[step])$day,  subset(evivo_pl, experiment == levels(evivo_pl$experiment)[step])$baso_CD69_mfi, col="blue")
	}



quartz(); plot(MIT$baso_CD203_mfi, MIT$percent63hi, type = "n", xlim = c(300, 10000), ylim = c(0,100))
for (step in 1:length(levels(MIT$experiment))) {
	sub <- subset(MIT, experiment == levels(MIT$experiment)[step])
	points(subset(sub, stim == "IL3")$baso_CD203_mfi, subset(sub, stim == "aIgE")$percent63hi, xlim = c(300, 10000), ylim = c(0,100), col = step)
	}
	
	
	
quartz(); plot(MIT$baso_CD69_mfi, MIT$percent63hi, type = "n", xlim = c(0, 1000), ylim = c(0,100))
for (step in 1:length(levels(MIT$experiment))) {
	sub <- subset(MIT, experiment == levels(MIT$experiment)[step])
	points(subset(sub, stim == "IL3")$baso_CD69_mfi, subset(sub, stim == "PN3")$percent63hi, xlim = c(300, 10000), ylim = c(0,100), col = step)
	}