# Load one or the other. You'll need to adjust the path to fit where you saved the data
# The first part splits the data based on arbor type, cell type, and/or brain region.
# The second part does the analysis and produces the resulting analysis csv file.
# The third part has code for doing correlations
# The fourth part can be used to generate histograms of whatever data you're interested in
caulDatFull = read.csv("C:/caulescence/analysis/allCaulescence-sum.csv")

maxCaulDatFull = read.csv("C:/caulescence/analysis/allCaulescence-max.csv")

# First part
caulDat = subset(caulDatFull,Bifurcations > 4)
apicCaul = subset(caulDat,ArborType=="ApicalDendrite")
dendCaul = subset(caulDat,ArborType=="Dendrite")
axonCaul = subset(caulDat,ArborType=="Axon")

maxcaulDat = subset(maxCaulDatFull,Bifurcations > 4)
maxapicCaul = subset(maxcaulDat,ArborType=="ApicalDendrite")
maxdendCaul = subset(maxcaulDat,ArborType=="Dendrite")
maxaxonCaul = subset(maxcaulDat,ArborType=="Axon")

#summary(dendCaul$CellType) # gives counts for each cell type

pyrBasal = dendCaul[dendCaul$CellType=="Pyramidal cell",]
pyrCortBasal = pyrBasal[pyrBasal$BrainRegion=="Cerebral Cortex",]
pyrHipBasal = pyrBasal[pyrBasal$BrainRegion=="Hippocampus",]
maxpyrBasal = maxdendCaul[maxdendCaul$CellType=="Pyramidal cell",]
maxpyrCortBasal = maxpyrBasal[maxpyrBasal$BrainRegion=="Cerebral Cortex",]
maxpyrHipBasal = maxpyrBasal[maxpyrBasal$BrainRegion=="Hippocampus",]

pyrApical = apicCaul[apicCaul$CellType=="Pyramidal cell",]
pyrCortApical = pyrApical[pyrApical$BrainRegion=="Cerebral Cortex",]
pyrHipApical = pyrApical[pyrApical$BrainRegion=="Hippocampus",]
maxpyrApical = maxapicCaul[maxapicCaul$CellType=="Pyramidal cell",]
maxpyrCortApical = maxpyrApical[maxpyrApical$BrainRegion=="Cerebral Cortex",]
maxpyrHipApical = maxpyrApical[maxpyrApical$BrainRegion=="Hippocampus",]

pyrRodentApical = pyrApical[pyrApical$Species%in%c("Rat","Mouse","Rabbit","Guinea-pig"),]

pyrCortRodentApical = pyrCortApical[pyrCortApical$Species%in%c("Rat","Mouse","Rabbit","Guinea-pig"),]
pyrCortPrimateApical = pyrCortApical[pyrCortApical$Species%in%c("Human","Monkey"),]
pyrHipRodentApical = pyrHipApical[pyrHipApical$Species%in%c("Rat","Mouse","Rabbit","Guinea-pig"),]
pyrHipPrimateApical = pyrHipApical[pyrHipApical$Species%in%c("Human","Monkey"),]
plot(pyrCortRodentApical$L_Caulescence,pyrCortRodentApical$L_GlobalAsym,col=1,xlim=c(.1,.9),ylim=c(.1,.9))
points(pyrHipRodentApical$L_Caulescence,pyrHipRodentApical$L_GlobalAsym,col=4)

points(pyrCortPrimateApical$L_Caulescence,pyrCortPrimateApical$L_GlobalAsym,col=5,pch=3)

# Comparing ratio of GlobalAsym to Caul in higher caul pyrApical vs lower caul pyrApical (Figure 4, top-left inset)
pyrApicalHi = pyrApicalT[pyrApicalT$L_Caulescence > .5,]
pyrApicalLow = pyrApicalT[pyrApicalT$L_Caulescence <= .5,]
hiRatio = pyrApicalHi$L_GlobalAsym/pyrApicalHi$L_Caulescence
lowRatio = pyrApicalLow$L_GlobalAsym/pyrApicalLow$L_Caulescence
hiSE = sd(hiRatio)/sqrt(length(hiRatio))
hiSE = sd(lowRatio)/sqrt(length(lowRatio))
xdat = c(mean(lowRatio),mean(hiRatio))
xErr = c(lowSE,hiSE)
ymax = (max(xdat) + max(xErr))*1.1
#ymin = (min(xdat) - min(xErr))*0.7
ymin = 0
cols=c("grey40","grey70")
mp = barplot(xdat,names.arg=c("Caulescence <= .5","Caulescence > .5"),ylab="Mean Ratio of Global Asymmetry to Caulescence",ylim=c(ymin,ymax),col=cols)
segments(mp, xdat, mp, xdat + xErr, col = 1, lwd = 1.5)
segments(mp, xdat, mp, xdat - xErr, col = 1, lwd = 1.5)

pyrCortRodentBasal = pyrCortBasal[pyrCortBasal$Species%in%c("Rat","Mouse","Rabbit","Guinea-pig"),]
pyrCortPrimateBasal = pyrCortBasal[pyrCortBasal$Species%in%c("Human","Monkey"),]
pyrHipRodentBasal = pyrHipBasal[pyrHipBasal$Species%in%c("Rat","Mouse","Rabbit","Guinea-pig"),]
pyrHipPrimateBasal = pyrHipBasal[pyrHipBasal$Species%in%c("Human","Monkey"),]

pyrCortApicalHi = pyrCortApical[pyrCortApical$L_Caulescence > .5,]
pyrCortApicalLow = pyrCortApical[pyrCortApical$L_Caulescence <= .5,]
pyrHipApicalHi = pyrHipApical[pyrHipApical$L_Caulescence > .5,]
pyrHipApicalLow = pyrHipApical[pyrHipApical$L_Caulescence <= .5,]

pyrCortPrimateApicalHi = pyrCortPrimateApical[pyrCortPrimateApical$L_Caulescence > .5,]
pyrCortPrimateApicalLow = pyrCortPrimateApical[pyrCortPrimateApical$L_Caulescence <= .5,]
pyrRodentApicalHi = pyrRodentApical[pyrRodentApical$L_Caulescence > .5,]
pyrRodentApicalLow = pyrRodentApical[pyrRodentApical$L_Caulescence <= .5,]
pyrCortRodentApicalHi = pyrCortRodentApical[pyrCortRodentApical$L_Caulescence > .5,]
pyrCortRodentApicalLow = pyrCortRodentApical[pyrCortRodentApical$L_Caulescence <= .5,]
pyrHipRodentApicalHi = pyrHipRodentApical[pyrHipRodentApical$L_Caulescence > .5,]
pyrHipRodentApicalLow = pyrHipRodentApical[pyrHipRodentApical$L_Caulescence <= .5,]

plot(pyrHipRodentApical$L_Caulescence,pyrHipRodentApical$L_GlobalAsym,col=3,xlim=c(.1,.9),ylim=c(.1,.9))
points(pyrCortRodentApical$L_Caulescence,pyrCortRodentApical$L_GlobalAsym,col=4,pch=4,cex=.6)
points(pyrCortPrimateApical$L_Caulescence,pyrCortPrimateApical$L_GlobalAsym,col=5,pch=5,cex=.6)
lines(c(0,1),c(0,1))
legend(.45,.9,c("Hippocampal","Cortical"),col=c(3,4),pch=c(1,4))

pyrApicalHiT = pyrCortPrimateApicalHi
pyrApicalLowT = pyrCortPrimateApicalLow
pyrApicalHiT = pyrCortRodentApicalHi
pyrApicalLowT = pyrCortRodentApicalLow
pyrApicalHiT = pyrHipRodentApicalHi
pyrApicalLowT = pyrHipRodentApicalLow
hiRatio = pyrApicalHiT$L_GlobalAsym/pyrApicalHiT$L_Caulescence
lowRatio = pyrApicalLowT$L_GlobalAsym/pyrApicalLowT$L_Caulescence
hiSE = sd(hiRatio)/sqrt(length(hiRatio))
hiSE = sd(lowRatio)/sqrt(length(lowRatio))
xdat = c(mean(lowRatio),mean(hiRatio))
xErr = c(lowSE,hiSE)
ymax = (max(xdat) + max(xErr))*1.1
#ymin = (min(xdat) - min(xErr))*0.7
ymin = 0
cols=c("grey60","grey60")
mp = barplot(xdat,names.arg=c("Caulescence <= .5","Caulescence > .5"),ylab="Mean Ratio of Global Asymmetry to Caulescence",ylim=c(ymin,ymax),col=cols)
segments(mp, xdat, mp, xdat + xErr, col = 1, lwd = 1.5)
segments(mp, xdat, mp, xdat - xErr, col = 1, lwd = 1.5)

BasketAxonCb = axonCaul[axonCaul$CellType=="Basket cell" & axonCaul$BrainRegion=="Cerebral Cortex",]
pyrAxon = axonCaul[axonCaul$CellType=="Pyramidal cell",]
upnAxon = axonCaul[axonCaul$CellType=="Uniglomerular projection neuron",]

maxBasketAxonCb = maxaxonCaul[maxaxonCaul$CellType=="Basket cell" & axonCaul$BrainRegion=="Cerebral Cortex",]
maxpyrAxon = maxaxonCaul[maxaxonCaul$CellType=="Pyramidal cell",]
maxupnAxon = maxaxonCaul[maxaxonCaul$CellType=="Uniglomerular projection neuron",]

analysis = getAllAnalysis(axonCaul,mean)
rownames(analysis) = "Axon Means"
analysis["Axon Pyr Means",] = getAllAnalysis(pyrAxon,mean)
analysis["Axon UPN Means",] = getAllAnalysis(upnAxon,mean)
analysis["Axon Cortical Basket Means",] = getAllAnalysis(BasketAxonCb,mean)

maxanalysis = getAllAnalysis(maxaxonCaul,mean)
rownames(maxanalysis) = "Axon Means"
maxanalysis["Axon Pyr Means",] = getAllAnalysis(maxpyrAxon,mean)
maxanalysis["Axon UPN Means",] = getAllAnalysis(maxupnAxon,mean)
maxanalysis["Axon Cortical Basket Means",] = getAllAnalysis(maxBasketAxonCb,mean)

sdAnalysis = getAllAnalysis(axonCaul,sd)
rownames(sdAnalysis) = "Axon SDs"
sdAnalysis["Axon Pyr SDs",] = getAllAnalysis(pyrAxon,sd)
sdAnalysis["Axon UPN SDs",] = getAllAnalysis(upnAxon,sd)
sdAnalysis["Axon Cortical Basket SDs",] = getAllAnalysis(BasketAxonCb,sd)

BasketDend = dendCaul[dendCaul$CellType=="Basket cell",]
GangDend = dendCaul[dendCaul$CellType=="Ganglion cell",]
MotoDend = dendCaul[dendCaul$CellType=="Motoneuron",]
maxBasketDend = maxdendCaul[maxdendCaul$CellType=="Basket cell",]
maxGangDend = maxdendCaul[maxdendCaul$CellType=="Ganglion cell",]
maxMotoDend = maxdendCaul[maxdendCaul$CellType=="Motoneuron",]

analysis = getAllAnalysis(dendCaul,mean)
rownames(analysis) = "Dendrite Means"
analysis["Dend Pyr Means",] = getAllAnalysis(pyrBasal,mean)
analysis["Dend Hip Pyr Means",] = getAllAnalysis(pyrHipBasal,mean)
analysis["Dend Cort Pyr Means",] = getAllAnalysis(pyrCortBasal,mean)
analysis["Dend Ganglion Means",] = getAllAnalysis(GangDend,mean)
analysis["Dend Basket Means",] = getAllAnalysis(BasketDend,mean)
analysis["Dend Motoneuron Means",] = getAllAnalysis(MotoDend,mean)

maxanalysis = getAllAnalysis(maxdendCaul,mean)
rownames(maxanalysis) = "Dendrite Means"
maxanalysis["Dend Pyr Means",] = getAllAnalysis(maxpyrBasal,mean)
maxanalysis["Dend Ganglion Means",] = getAllAnalysis(maxGangDend,mean)
maxanalysis["Dend Basket Means",] = getAllAnalysis(maxBasketDend,mean)
maxanalysis["Dend Motoneuron Means",] = getAllAnalysis(maxMotoDend,mean)

sdAnalysis = getAllAnalysis(dendCaul,sd)
rownames(sdAnalysis) = "Dendrite SDs"
sdAnalysis["Dend Pyr SDs",] = getAllAnalysis(pyrBasal,sd)
sdAnalysis["Dend Ganglion SDs",] = getAllAnalysis(GangDend,sd)
sdAnalysis["Dend Basket SDs",] = getAllAnalysis(BasketDend,sd)
sdAnalysis["Dend Motoneuron SDs",] = getAllAnalysis(MotoDend,sd)

analysis = getAllAnalysis(apicCaul,mean)
rownames(analysis) = "Apical Means"
analysis["Apical Cort Means",] = getAllAnalysis(pyrCortApical,mean)
analysis["Apical Hip Means",] = getAllAnalysis(pyrHipApical,mean)

maxanalysis = getAllAnalysis(maxapicCaul,mean)
rownames(maxanalysis) = "Apical Means"
maxanalysis["Apical Cort Means",] = getAllAnalysis(maxpyrCortApical,mean)
maxanalysis["Apical Hip Means",] = getAllAnalysis(maxpyrHipApical,mean)

# Second part - produces means and standard deviations for each group
getAnalysis = function(dat,metric,analysis=mean){
	valNames = c("TotalSize","UniformAsym","MinDegAsym","GlobalAsym","ExcessAsym","Caulescence")
	valNames = paste(metric,valNames,sep="")
	data.frame(t(analysis(dat[,names(dat)%in%valNames])))
}
getAllAnalysis = function(dat,analysis=mean){
	merge(
		merge(getAnalysis(dat,"N_",analysis),getAnalysis(dat,"L_",analysis)),
		merge(getAnalysis(dat,"SA_",analysis),getAnalysis(dat,"V_",analysis)))
}
analysis = getAllAnalysis(pyrBasal,mean)
rownames(analysis) = "Basal Pyr Means"
analysis["Basal Pyr StdDevs",] = getAllAnalysis(pyrBasal,sd)
analysis["Basal Cort Pyr Means",] = getAllAnalysis(pyrCortBasal,mean)
analysis["Basal Cort Pyr StdDevs",] = getAllAnalysis(pyrCortBasal,sd)
analysis["Basal Hip Pyr Means",] = getAllAnalysis(pyrHipBasal,mean)
analysis["Basal Hip Pyr StdDevs",] = getAllAnalysis(pyrHipBasal,sd)
analysis["Apical Pyr Means",] = getAllAnalysis(pyrApical,mean)
analysis["Apical Pyr StdDevs",] = getAllAnalysis(pyrApical,sd)
analysis["Apical Cort Pyr Means",] = getAllAnalysis(pyrCortApical,mean)
analysis["Apical Cort Pyr StdDevs",] = getAllAnalysis(pyrCortApical,sd)
analysis["Apical Hip Pyr Means",] = getAllAnalysis(pyrHipApical,mean)
analysis["Apical Hip Pyr StdDevs",] = getAllAnalysis(pyrHipApical,sd)
analysis["Axon Cortical Basket Means",] = getAllAnalysis(BasketAxonCb,mean)
analysis["Axon Cortical Basket StdDevs",] = getAllAnalysis(BasketAxonCb,sd)
analysis["Dend Moto Means",] = getAllAnalysis(MotoDend,mean)
analysis["Dend Moto StdDevs",] = getAllAnalysis(MotoDend,sd)

write.csv(analysis,'C:/caulescence/sumCaulAnalysis.csv')
write.csv(analysis,'C:/caulescence/maxCaulAnalysis.csv')

#2.5 significance
wilcox.test(BasketAxonCb$L_Caulescence,BasketAxonCb$L_GlobalAsym)
wilcox.test(upnAxon$L_Caulescence,upnAxon$L_GlobalAsym)
wilcox.test(apicCaul$L_Caulescence,pyrBasal$L_Caulescence)
wilcox.test(apicCaul$L_Caulescence,apicCaul$L_GlobalAsym)

wilcox.test(dendCaul$L_Caulescence,dendCaul$L_GlobalAsym)
wilcox.test(GangDend$L_Caulescence,GangDend$L_GlobalAsym)
wilcox.test(MotoDend$L_Caulescence,MotoDend$L_GlobalAsym)
wilcox.test(BasketDend$L_Caulescence,BasketDend$L_GlobalAsym)

wilcox.test(pyrApical$L_Caulescence,pyrApical$L_GlobalAsym)
wilcox.test(pyrBasal$L_Caulescence,pyrBasal$L_GlobalAsym)

wilcox.test(BasketAxonCb$L_Caulescence,BasketAxonCb$V_Caulescence)
wilcox.test(upnAxon$L_Caulescence,upnAxon$V_Caulescence)
wilcox.test(BasketAxonCb$N_Caulescence,BasketAxonCb$V_Caulescence)
wilcox.test(upnAxon$N_Caulescence,upnAxon$V_Caulescence)
wilcox.test(BasketAxonCb$L_Caulescence,BasketAxonCb$N_Caulescence)
wilcox.test(upnAxon$L_Caulescence,upnAxon$N_Caulescence)

# Third part, Correlations: Use names(pyrBasal) to see what the field names are
names(pyrBasal)
# They'll be the same for all of the data frames
cor(pyrBasal$N_Caulescence,pyrBasal$N_GlobalAsym)
cor(MotoDend$N_Caulescence,MotoDend$N_GlobalAsym)
cor(BasketAxonCb$N_Caulescence,BasketAxonCb$N_GlobalAsym)
cor(pyrApical$N_Caulescence,pyrApical$N_GlobalAsym)

cor(pyrBasal$L_Caulescence,pyrBasal$L_GlobalAsym)
cor(MotoDend$L_Caulescence,MotoDend$L_GlobalAsym)
cor(BasketAxonCb$L_Caulescence,BasketAxonCb$L_GlobalAsym)
cor(pyrApical$L_Caulescence,pyrApical$L_GlobalAsym)

cor(pyrBasal$N_Caulescence,pyrBasal$L_Caulescence)
cor(pyrBasal$N_Caulescence,pyrBasal$V_Caulescence)
cor(pyrBasal$L_Caulescence,pyrBasal$V_Caulescence)

plot(type='n')
plot(BasketAxonCb$N_Caulescence,BasketAxonCb$N_GlobalAsym,col="red")
points(pyrApical$N_Caulescence,pyrApical$N_GlobalAsym,col="blue")
lines(c(0,1),c(0,1))

# Fourth part, Histograms
#brks = 1:20/20
#windows(width=4.2,height=3.6)
#hist(axonCaul$Exhaustive,main="Axon Caulescence Histogram",include.lowest=TRUE,breaks=brks)
#windows(width=4.2,height=3.6)
#hist(dendCaul$Exhaustive,main="Dendrite Caulescence Histogram",include.lowest=TRUE,breaks=brks)
#windows(width=4.2,height=3.6)
#hist(apicCaul$Exhaustive,main="ApicalDendrite Caulescence Histogram",include.lowest=TRUE,breaks=brks)
#windows(width=4.2,height=3.6)
#hist(axonCaul$UniformAsym,main="Axon Asymmetry Histogram",include.lowest=TRUE,breaks=brks)
#windows(width=4.2,height=3.6)
#hist(dendCaul$UniformAsym,main="Dendrite Asymmetry Histogram",include.lowest=TRUE,breaks=brks)
#windows(width=4.2,height=3.6)
#hist(apicCaul$UniformAsym,main="ApicalDendrite Asymmetry Histogram",include.lowest=TRUE,breaks=brks)

joinDataFrames = function(dat1,dat2){
	if (is.null(dat1)){
		return(dat2);
	}
	dat1[(dim(dat1)[1]+1):(dim(dat1)[1]+dim(dat2)[1]),] = dat2
	return(dat1)
}

caData = read.csv("C:/caulescence/CA.csv",stringsAsFactors=F)

#Fifth part, going through all combinations to find slopes & correlations
# This is useful to quickly skim large amounts of data
correlations = function(dat,writeCSV=T,plotCorrelations=F,minN=14){
	analysisDF = NULL
	for (arborType in c("ApicalDendrite","Dendrite","Axon")){
		arborTypeSubset = dat[dat["ArborType"]==arborType,]
		tmpDF = datasetAnalysis(arborType,arborTypeSubset,plotCorrelations)
		analysisDF = joinDataFrames(analysisDF,tmpDF)
		for (cellType in levels(arborTypeSubset[["CellType"]])){
			cellTypeSubset = arborTypeSubset[arborTypeSubset[["CellType"]]==cellType,]
			if (dim(cellTypeSubset)[1] > minN){
				subsetName = paste(arborType,cellType)
				print(subsetName)
				tmpDF = datasetAnalysis(subsetName,cellTypeSubset,plotCorrelations)
				analysisDF = joinDataFrames(analysisDF,tmpDF)
				for (brainRegion in levels(cellTypeSubset[["BrainRegion"]])){
					BR_CT_Subset = cellTypeSubset[cellTypeSubset[["BrainRegion"]]==brainRegion,]
					if (dim(BR_CT_Subset)[1] > minN){
						subsetName = paste(arborType,cellType,brainRegion)
						print(subsetName)
						tmpDF = datasetAnalysis(subsetName,BR_CT_Subset,plotCorrelations)
						analysisDF = joinDataFrames(analysisDF,tmpDF)
						# Go into CA1/CA3
						for (subRegion in c("CA1","CA3")){
							subRegionSubset = BR_CT_Subset[BR_CT_Subset[["Neuron"]] %in% caData[[subRegion]],]
							if (dim(subRegionSubset)[1] > minN){
								subsetName = paste(arborType,cellType,brainRegion,subRegion)
								print(subsetName)
								tmpDF = datasetAnalysis(subsetName,subRegionSubset,plotCorrelations)
								analysisDF = joinDataFrames(analysisDF,tmpDF)
							}
						}
					}
				}
			}
		}
	}
	if (writeCSV){
		names(analysisDF) = c("Subset","N","Caulescence Mean","Caulescence SD","MinDegAsym Mean","MinDegAsym SD",
			"GlobalAsym Mean","GlobalAsym SD","MainPathRatio Mean","MainPathRatio SD",
			"Caul-GlobalAsym Cor","Caul-GlobalAsym LmIntercept","Caul-GlobalAsym LmSlope",
			"Caul-MinDegAsym Cor","Caul-MinDegAsym LmIntercept","Caul-MinDegAsym LmSlope",
			"Caul-MainPathRatio Cor","Caul-MainPathRatio LmIntercept","Caul-MainPathRatio LmSlope")
		write.csv(analysisDF,"C:/caulescence/deeperCaulAnalysis.csv",row.names=F)
	}
}
correlations(caulDat)

datasetAnalysis = function(subsetName,dat,plotCorrelations=F){
	metrics = c("N_","L_","SA_","V_")
	names(metrics) = c("Nodal","Length","Surface","Volume")
	analysisDF = NULL
	
	for (metricIndex in 1:length(metrics)){
		metric = metrics[metricIndex]
		caul = dat[paste(metric,"Caulescence",sep="")]
		minDegAsym = dat[paste(metric,"MinDegAsym",sep="")]
		globalAsym = dat[paste(metric,"GlobalAsym",sep="")]
		mainPathRatio = dat[paste(metric,"MainPathSize",sep="")]/dat[paste(metric,"TotalSize",sep="")]

		vals = c("Caulescence","MinDegAsym","GlobalAsym")
		tmpDat = dat[,paste(metric,vals,sep="")]
		#Add main path ratio
		tmpDat = cbind(tmpDat,mainPathRatio)
		names(tmpDat) = c(vals,"MainPathRatio")

		caulMean = mean(caul)
		caulSd = sd(caul)
		minDegAsymMean = mean(minDegAsym)
		minDegAsymSd = sd(minDegAsym)
		globalAsymMean = mean(globalAsym)
		globalAsymSd = sd(globalAsym)
		mainPathRatioMean = mean(mainPathRatio)
		mainPathRatioSd = sd(mainPathRatio)
		caulGlobalAsymCor = cor(tmpDat$Caulescence,tmpDat$GlobalAsym)
		caulGlobalAsymLm = lm(Caulescence~GlobalAsym,tmpDat)
		caulMinDegAsymCor = cor(tmpDat$Caulescence,tmpDat$MinDegAsym)
		caulMinDegAsymLm = lm(Caulescence~MinDegAsym,tmpDat)
		caulMainPathCor = cor(tmpDat$Caulescence,tmpDat$MainPathRatio)
		caulMainPathLm = lm(Caulescence~MainPathRatio,tmpDat)
		caulMinDegAsymCor = cor(tmpDat$Caulescence,tmpDat$MinDegAsym)
		caulMinDegAsymLm = lm(Caulescence~MinDegAsym,tmpDat)
		
		if (plotCorrelations){
			par(ask=T)
			plot(tmpDat$Caulescence,tmpDat$GlobalAsym,xlab=paste(metric,"Caulescence",sep=""),
				ylab=paste(metric,"GlobalAsym",sep=""),main=paste(subsetName,"Scatter Plot"))
			lines(c(0,1),c(0,1))
			plot(tmpDat$Caulescence,tmpDat$MinDegAsym,xlab=paste(metric,"Caulescence",sep=""),
				ylab=paste(metric,"MinDegAsym",sep=""),main=paste(subsetName,"Scatter Plot"))
			lines(c(0,1),c(0,1))
			plot(tmpDat$Caulescence,tmpDat$MainPathRatio,xlab=paste(metric,"Caulescence",sep=""),
				ylab=paste(metric,"MainPathRatio",sep=""),main=paste(subsetName,"Scatter Plot"))
			lines(c(0,1),c(0,1))
		}

		name = paste(subsetName,names(metrics)[metricIndex])
		tmpDF = data.frame(name,dim(tmpDat)[1],caulMean,caulSd,minDegAsymMean,minDegAsymSd,
			globalAsymMean,globalAsymSd,mainPathRatioMean,mainPathRatioSd,
			caulGlobalAsymCor,caulGlobalAsymLm$coefficients[1],caulGlobalAsymLm$coefficients[2],
			caulMinDegAsymCor,caulMinDegAsymLm$coefficients[1],caulMinDegAsymLm$coefficients[2],
			caulMainPathCor,caulMainPathLm$coefficients[1],caulMainPathLm$coefficients[2],stringsAsFactors=F);
		analysisDF = joinDataFrames(analysisDF,tmpDF)
	}
	return(analysisDF)
}

correlations(caulDat)
correlations(caulDat,writeCSV=F,plotCorrelations=T)

all metrics
	correlation, slope
	means, variances
		caulescence global asym, minDeg asym

#Insert CA1/CA3 (subregion) into datafile
addSubregion = function(datafile,caFile){
	caData = read.csv(caFile,stringsAsFactors=F)
	dat = read.csv(datafile)
	ca1Set = dat[dat[["Neuron"]] %in% caData$CA1,]
	ca1Set = addCol(ca1Set,"SubRegion","CA1")
	ca3Set = dat[dat[["Neuron"]] %in% caData$CA3,]
	ca3Set = addCol(ca3Set,"SubRegion","CA3")
	allOther = dat[!(dat[["Neuron"]] %in% caData$CA3 | dat[["Neuron"]] %in% caData$CA1),]
	allOther = addCol(allOther,"SubRegion","")

	recombined = joinDataFrames(joinDataFrames(allOther,ca1Set),ca3Set)
	write.csv(recombined,datafile)
}

addCol = function(x,columnName,value){
	xSize = dim(x)[1]
	newCol = data.frame(rep(value,xSize),stringsAsFactors=F)
	names(newCol) = columnName
	return(cbind(x,newCol))
}

datafile = "C:/caulescence/allCaulescence-sum.csv"
caFile = "C:/caulescence/CA.csv"

addSubregion(datafile,caFile)