pcr <- read.csv("/Users/Guest/Public/Shreff_lab/JAX/pcr_data/pcr.csv")
genef <- factor(pcr$GENE)
exp <- tapply(pcr$Ct, genef, fivenum)
quartz(,8,4); boxplot(exp)

# explore using product of housekeeping genes to filter out poor quality data
quartz(,12,4); layout(matrix(c(1:3),1,3))
plot(subset(pcr, GENE=="RPL41")$Ct*subset(pcr, GENE=="RPS9")$Ct, subset(pcr, GENE=="Tbet")$Ct, xlab="RPL X RPS", ylab="Tbet Ct")
abline(v=600,lty=3)
plot(subset(pcr, GENE=="RPL41")$Ct*subset(pcr, GENE=="RPS9")$Ct, subset(pcr, GENE=="CD25")$Ct, xlab="RPL X RPS", ylab="CD25 Ct")
abline(v=600,lty=3)
plot(subset(pcr, GENE=="RPL41")$Ct*subset(pcr, GENE=="RPS9")$Ct, subset(pcr, GENE=="FOXP3")$Ct, xlab="RPL X RPS", ylab="FOXP3 Ct")
abline(v=600,lty=3)

# read in document containing reference back to sample ID
rna_id <- read.csv("/Users/Guest/Public/Shreff_lab/JAX/pcr_data/rna_isolation.csv")
rpl <- subset(pcr, GENE == "RPL41")
rps <- subset(pcr, GENE == "RPS9")
index1 <- match(rpl$SAMPLEID, rna_id$STORED.IN.TUBE)

quartz(,8,4); layout(matrix(c(1,2),1,2))
plot(rps$Ct, rna_id$PURITY[index1],xlim=c(15,35), xlab="Ct")
points(rpl$Ct, rna_id$PURITY[index1],col=3)
legend("topright", c("rps", "rpl"), pch=19,col=c(1,3))
plot(rps$Ct, rna_id$CONCENTRATION[index1],xlim=c(15,35), ylim=c(0,50), xlab="Ct")
points(rpl$Ct, rna_id$CONCENTRATION[index1],col=3)


# explore using only RPL41 to filter out poor quality data
quartz(,12,4); layout(matrix(c(1:3),1,3))
plot(subset(pcr, GENE=="RPL41")$Ct, subset(pcr, GENE=="Tbet")$Ct, xlab="RPL", ylab="Tbet Ct")
abline(v=25,lty=3)
plot(subset(pcr, GENE=="RPL41")$Ct, subset(pcr, GENE=="CD25")$Ct, xlab="RPL", ylab="CD25 Ct")
abline(v=25,lty=3)
plot(subset(pcr, GENE=="RPL41")$Ct, subset(pcr, GENE=="FOXP3")$Ct, xlab="RPL", ylab="FOXP3 Ct")
abline(v=25,lty=3)

# use this to filter

rpl <- subset(rpl, Ct <30)

# ONLY THE TUBE ID'S FOR WHICH THE RPL41 CT WAS <30 WILL BE INCLUDED IN FURTHER ANALYSES

index2 <- pcr$SAMPLEID %in% rpl$SAMPLEID

pcr <- data.frame(SAMPLE = NA, TUBEID = pcr$SAMPLEID[index2], WELLID = pcr$WELL[index2], GENE = pcr$GENE[index2], Ct = pcr$Ct[index2], Exp = pcr$Exp[index2], stim = NA, dCT = NA)

genef <- factor(pcr$GENE)
exp <- tapply(pcr$Ct, genef, fivenum)
quartz(,8,4); boxplot(exp)

# now need to work on analyses by stim, patient, visit

i <- grep("MC", rna_id$STIMULANT)
MC <- pcr$TUBEID %in% rna_id$STORED.IN.TUBE[i]
pcr$stim[MC] <- "MED_CTRL"

i <- grep("MM", rna_id$STIMULANT)
MM <- pcr$TUBEID %in% rna_id$STORED.IN.TUBE[i]
pcr$stim[MM] <- "MUS_M1"

i <- grep("TT", rna_id$STIMULANT)
TT <- pcr$TUBEID %in% rna_id$STORED.IN.TUBE[i]
pcr$stim[TT] <- "TT"


# need to calculate dCt values by tube and grab sample id from rna_isolation file


for (step in 1:length(levels(pcr$TUBEID))) {
	w <- subset(pcr, TUBEID == levels(pcr$TUBEID)[step])
	if (dim(w)[1] != 0) {
		i <- grep(levels(pcr$TUBEID)[step], pcr$TUBEID)
		pcr$SAMPLE[i] <- as.character(subset(rna_id, STORED.IN.TUBE. == levels(pcr$TUBEID)[step])$SAMPLE.ID)
		# calculate dCt using RPL41
		pcr[i,]$dCT = pcr[i,]$Ct - subset(pcr[i,], GENE == "RPL41")$Ct
		}
	}

sum <- by(pcr$dCT, list(pcr$stim,pcr$GENE), fivenum)
quartz(,10,4)

boxplot(sum)
abline(v=c(3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),lty=3)
mtext(dimnames(sum)[[2]],at=c(2,5,8,11,14,17,20,23,26,29))