###### FROM BEFORE, FIND THIS SECTION

Index = which(pns %in% keep_pns)
newP = p[Index,]
newchr=chr[Index];newpos=pos[Index]
newpns = clusterMaker(newchr,newpos)
cat("\n")

###############
# NEW BELOW - ADD

coi = pd$Age # what I'm trying, not sure what youre looking at
# make groups into factors
if(is.character(coi) | is.factor(coi)) {
	groups = factor(coi)
	gNames= levels(groups)
}
gIndexes=split(1:length(groups),groups)

coiName = "Aging"
plotGaps=1e6

outPath = getwd()



#############################
# covariate info
cat("Averaging Methylation in Probe Groups.")

Indexes = splitit(newpns)
expMat = matrix(nr = length(Indexes), nc = length(coi))
colnames(expMat) = colNames(newP)

for(i in seq(along=Indexes)) {
	if(i %% 10000 == 0) cat(".")
	Index = Indexes[[i]]
	if(length(Index) > 1) {
		expMat[i,] = colMeans(newP[Index,])
	} else {
		expMat[i,] = as.numeric(newP[Index,])
	}
}

cat("\n")
# make stuff to loop over
chrIndexes = sort(unique(newchr))
# chrIndexes = chrIndexes[-grep("random",chrIndexes)]
# chrIndexes = chrIndexes[-grep("chrM",chrIndexes)]

##########
# make some plots
cat("plotting","\n")
pdf(paste(outPath,"/",coiName,"_blocks.pdf",sep=""),
	height = 6, width = 16)
mypar(1,1)

for(i in 1:length(chrIndexes)) {
	Index = which(newMat$chr %in% chrIndexes[i])
	tmp = newMat[Index,]
	ppos = rowMeans(tmp[,2:3])
		tmpp = expMat[Index,]

	# do the centomeres? or just break it up
	cc = clusterMaker(tmp$chr,ppos,maxGap = plotGaps)
	Indexes = splitit(cc)
	for(j in seq(along=Indexes)) {
		Index2 = Indexes[[j]]
		if(length(Index2) > 20) {
			xx = ppos[Index2]
			pp = tmpp[Index2,]
			matplot(x = xx, y = pp,type = "p",
				col = as.numeric(factor(coi)), pch = 19, cex = 0.7,
				xlab = "position", ylab = "methylation",
				main = paste(chrIndexes[i], "- chunk",
					names(Indexes)[j]), ylim = c(0,1))
			for(k in seq(along = gIndexes)) {
				  if(length(gIndexes[[k]]) == 1) yy=pp[,gIndexes[[k]]]
				  if(length(gIndexes[[k]]) > 1)	  yy=rowMeans(pp[,gIndexes[[k]]])
				fit1=loess(yy~xx,degree=1,span=0.2,	family="symmetric")
				lines(xx,fit1$fitted,col=k,lwd=2)
			}
			legend("topleft",legend=gNames,col=1:length(gNames),
			  lty=1, lwd = 4,cex=1)
		}
	}
}
dev.off()