### R code from vignette source 'mixnorm_Vignette.Rnw'

###################################################
### code chunk number 1: mixnorm_Vignette.Rnw:34-35
###################################################
  library(metabomxtr)


###################################################
### code chunk number 2: mixnorm_Vignette.Rnw:40-44
###################################################
 data(euMetabData)
 class(euMetabData)
 dim(euMetabData)
 head(euMetabData)


###################################################
### code chunk number 3: mixnorm_Vignette.Rnw:51-55
###################################################
data(euMetabCData)
class(euMetabCData)
dim(euMetabCData)
head(euMetabCData)


###################################################
### code chunk number 4: mixnorm_Vignette.Rnw:60-61
###################################################
ynames<-c('pyruvic_acid','malonic_acid')


###################################################
### code chunk number 5: mixnorm_Vignette.Rnw:66-91
###################################################
#get data in suitable format for plotting 
library(plyr)
control.data.copy<-euMetabCData[,c("batch","pheno","pyruvic_acid","malonic_acid")]
control.data.copy$type<-paste("Run Order=",substr(rownames(control.data.copy),6,6),sep="")
revalue(control.data.copy$pheno,c("MOM"="Mother Control","BABY"="Baby Control"))->control.data.copy$pheno
library(reshape2)
control.plot.data<-melt(control.data.copy,measure.vars=c("pyruvic_acid","malonic_acid"),variable.name="Metabolite",value.name="Abundance")
revalue(control.plot.data$Metabolite,c("malonic_acid"="Malonic Acid","pyruvic_acid"="Pyruvic Acid"))->control.plot.data$Metabolite

#get annotation of mean values for mother and baby samples 
control.mean.annotation<-ddply(control.plot.data,c("Metabolite","pheno","batch"),summarize,Abundance=mean(Abundance,na.rm=T))
control.mean.annotation$type<-"Batch Specific Mean"

#merge on to control data 
control.plot.data<-rbind(control.plot.data,control.mean.annotation[,c("batch","pheno","type","Metabolite","Abundance")])

#check for missing values in metabolite abundance
sum(is.na(control.plot.data$Abundance))
control.plot.data[which(is.na(control.plot.data$Abundance)),]

#now for plotting purposes, fill in missing metabolite values with something below the minimum detectable threshold 
#Based on experimental evidence not available here, the minimum detectable threshold for batch 1 was 10.76
min.threshold<-10.76
impute.val<-min.threshold-0.5
control.plot.data[which(is.na(control.plot.data$Abundance)),"Abundance"]<-impute.val


###################################################
### code chunk number 6: fig1
###################################################
#plot abundance by batch and run order
library(ggplot2)

#beginning of plot 
control.plot<-ggplot(control.plot.data,aes(x=batch,y=Abundance,color=pheno,shape=type,group=pheno))+
						geom_point(size=7)+geom_line(data=control.plot.data[which(control.plot.data$type=="Batch Specific Mean"),])+
						scale_color_discrete(name="Sample Type")+scale_shape_manual(name="Point Type",values=c(8,15,16,17))+
						xlab("Analytic Batch")+ylab("Metabolite Abundance")+facet_wrap(~Metabolite)
						
#now add in line segments indicating batch specific thresholds of detectability for the first 4 batches
batch.thresholds1.to.4<-c(10.76,11.51,11.36,10.31)
for (x in 1:4){
	threshold<-batch.thresholds1.to.4[x]
	batch.lower<-x-0.3
	batch.upper<-x+0.3
	control.plot<-control.plot+geom_segment(x=batch.lower,y=threshold,xend=batch.upper,yend=threshold,color="purple",linetype=2)
}

#add in threshold of detectability for the last batch and also add a legend to describe the threshold lines
batch.thresholds5<-11.90
control.plot<-control.plot+geom_segment(aes(x=4.7,y=batch.thresholds5,xend=5.3,yend=batch.thresholds5,linetype="Batch Specific\nThreshold of\nDetectability"),color="purple")+
										scale_linetype_manual(name="",values=2)+guides(color=guide_legend(order=1),shape=guide_legend(order = 2),linetype=guide_legend(order = 3))

#add in annotation for the point below the threshold of detectability 
annotate.df<-data.frame(Metabolite=factor("Malonic Acid",levels=c("Pyruvic Acid","Malonic Acid")),batch=factor(1,levels=1:5),Abundance=impute.val-0.5)
control.plot.final<-control.plot+geom_text(aes(x=batch,y=Abundance,label="(Unknown\nAbundance)",shape=NULL,group=NULL),color="Red",annotate.df,size=7)+
					theme(axis.title=element_text(size=35), axis.text=element_text(size=25),  strip.text=element_text(size=30), legend.title=element_text(size=25), legend.text=element_text(size=25))+
					guides(color=guide_legend(override.aes=list(size=2),order=1, keywidth=0.7, keyheight=0.5,default.unit="inch"), shape=guide_legend(override.aes=list(size=5), order=2,  keywidth=0.7, keyheight=0.5,default.unit="inch"), linetype=guide_legend(order=3),  keywidth=0.7, keyheight=0.5,default.unit="inch")
control.plot.final


###################################################
### code chunk number 7: fig1
###################################################
#plot abundance by batch and run order
library(ggplot2)

#beginning of plot 
control.plot<-ggplot(control.plot.data,aes(x=batch,y=Abundance,color=pheno,shape=type,group=pheno))+
						geom_point(size=7)+geom_line(data=control.plot.data[which(control.plot.data$type=="Batch Specific Mean"),])+
						scale_color_discrete(name="Sample Type")+scale_shape_manual(name="Point Type",values=c(8,15,16,17))+
						xlab("Analytic Batch")+ylab("Metabolite Abundance")+facet_wrap(~Metabolite)
						
#now add in line segments indicating batch specific thresholds of detectability for the first 4 batches
batch.thresholds1.to.4<-c(10.76,11.51,11.36,10.31)
for (x in 1:4){
	threshold<-batch.thresholds1.to.4[x]
	batch.lower<-x-0.3
	batch.upper<-x+0.3
	control.plot<-control.plot+geom_segment(x=batch.lower,y=threshold,xend=batch.upper,yend=threshold,color="purple",linetype=2)
}

#add in threshold of detectability for the last batch and also add a legend to describe the threshold lines
batch.thresholds5<-11.90
control.plot<-control.plot+geom_segment(aes(x=4.7,y=batch.thresholds5,xend=5.3,yend=batch.thresholds5,linetype="Batch Specific\nThreshold of\nDetectability"),color="purple")+
										scale_linetype_manual(name="",values=2)+guides(color=guide_legend(order=1),shape=guide_legend(order = 2),linetype=guide_legend(order = 3))

#add in annotation for the point below the threshold of detectability 
annotate.df<-data.frame(Metabolite=factor("Malonic Acid",levels=c("Pyruvic Acid","Malonic Acid")),batch=factor(1,levels=1:5),Abundance=impute.val-0.5)
control.plot.final<-control.plot+geom_text(aes(x=batch,y=Abundance,label="(Unknown\nAbundance)",shape=NULL,group=NULL),color="Red",annotate.df,size=7)+
					theme(axis.title=element_text(size=35), axis.text=element_text(size=25),  strip.text=element_text(size=30), legend.title=element_text(size=25), legend.text=element_text(size=25))+
					guides(color=guide_legend(override.aes=list(size=2),order=1, keywidth=0.7, keyheight=0.5,default.unit="inch"), shape=guide_legend(override.aes=list(size=5), order=2,  keywidth=0.7, keyheight=0.5,default.unit="inch"), linetype=guide_legend(order=3),  keywidth=0.7, keyheight=0.5,default.unit="inch")
control.plot.final


###################################################
### code chunk number 8: mixnorm_Vignette.Rnw:138-143
###################################################
#execute normalization
euMetabNorm <- mixnorm(ynames, batch="batch", mxtrModel=~pheno+batch|pheno+batch, 
				batchTvals=c(10.76,11.51,11.36,10.31,11.90), cData=euMetabCData,
				data=euMetabData)
#this produces warnings about NaNs produced, but this is the expected behavior of the function 


###################################################
### code chunk number 9: mixnorm_Vignette.Rnw:148-149
###################################################
euMetabNorm$normParamsZ


###################################################
### code chunk number 10: mixnorm_Vignette.Rnw:154-155
###################################################
head(euMetabNorm$ctlNorm)


###################################################
### code chunk number 11: mixnorm_Vignette.Rnw:161-162
###################################################
head(euMetabNorm$obsNorm)


###################################################
### code chunk number 12: mixnorm_Vignette.Rnw:169-184
###################################################
#get normalized control data in right format for plotting 
normalized.control.data<-control.data.copy
normalized.control.data[,ynames]<-euMetabNorm$ctlNorm[,ynames]
normalized.control.plot.data<-melt(normalized.control.data,measure.vars=ynames,variable.name="Metabolite",value.name="Abundance")
revalue(normalized.control.plot.data$Metabolite,c("malonic_acid"="Malonic Acid","pyruvic_acid"="Pyruvic Acid"))->normalized.control.plot.data$Metabolite

#get annotation of normalized mean values for mother and baby samples 
normalized.control.mean.annotation<-ddply(normalized.control.plot.data,c("Metabolite","pheno","batch"),summarize,Abundance=mean(Abundance,na.rm=T))
normalized.control.mean.annotation$type<-"Batch Specific Mean"

#merge on to control data 
normalized.control.plot.data<-rbind(normalized.control.plot.data,normalized.control.mean.annotation[,c("batch","pheno","type","Metabolite","Abundance")])

#impute value below the minimum detectable threshold for the 1 missing value (as above)
normalized.control.plot.data[which(is.na(normalized.control.plot.data$Abundance)),"Abundance"]<-impute.val


###################################################
### code chunk number 13: fig2
###################################################
#beginning of plot 
normalized.control.plot<-ggplot(normalized.control.plot.data,aes(x=batch,y=Abundance,color=pheno,shape=type,group=pheno))+
						geom_point(size=7)+geom_line(data=normalized.control.plot.data[which(normalized.control.plot.data$type=="Batch Specific Mean"),])+
						scale_color_discrete(name="Sample Type")+scale_shape_manual(name="Point Type",values=c(8,15,16,17))+
						xlab("Analytic Batch")+ylab("Metabolite Abundance")+facet_wrap(~Metabolite)
						
#now add in line segments indicating batch specific thresholds of detectability for the first 4 batches
batch.thresholds1.to.4<-c(10.76,11.51,11.36,10.31)
for (x in 1:4){
	threshold<-batch.thresholds1.to.4[x]
	batch.lower<-x-0.3
	batch.upper<-x+0.3
	normalized.control.plot<-normalized.control.plot+geom_segment(x=batch.lower,y=threshold,xend=batch.upper,yend=threshold,color="purple",linetype=2)
}

#add in threshold of detectability for the last batch and also add a legend to describe the threshold lines
batch.thresholds5<-11.90
normalized.control.plot<-normalized.control.plot+geom_segment(aes(x=4.7,y=batch.thresholds5,xend=5.3,yend=batch.thresholds5,linetype="Batch Specific\nThreshold of\nDetectability"),color="purple")+
										scale_linetype_manual(name="",values=2)+guides(color=guide_legend(order=1),shape=guide_legend(order = 2),linetype=guide_legend(order = 3))

#add in annotation for the point below the threshold of detectability 
normalized.control.plot.final<-normalized.control.plot+geom_text(aes(x=batch,y=Abundance,label="(Unknown\nAbundance)",shape=NULL,group=NULL),color="Red",annotate.df,size=7)+
					theme(axis.title=element_text(size=35), axis.text=element_text(size=25),  strip.text=element_text(size=30), legend.title=element_text(size=25), legend.text=element_text(size=25))+
					guides(color=guide_legend(override.aes=list(size=2),order=1, keywidth=0.7, keyheight=0.5,default.unit="inch"), shape=guide_legend(override.aes=list(size=5), order=2,  keywidth=0.7, keyheight=0.5,default.unit="inch"), linetype=guide_legend(order=3),  keywidth=0.7, keyheight=0.5,default.unit="inch")
normalized.control.plot.final


###################################################
### code chunk number 14: fig2
###################################################
#beginning of plot 
normalized.control.plot<-ggplot(normalized.control.plot.data,aes(x=batch,y=Abundance,color=pheno,shape=type,group=pheno))+
						geom_point(size=7)+geom_line(data=normalized.control.plot.data[which(normalized.control.plot.data$type=="Batch Specific Mean"),])+
						scale_color_discrete(name="Sample Type")+scale_shape_manual(name="Point Type",values=c(8,15,16,17))+
						xlab("Analytic Batch")+ylab("Metabolite Abundance")+facet_wrap(~Metabolite)
						
#now add in line segments indicating batch specific thresholds of detectability for the first 4 batches
batch.thresholds1.to.4<-c(10.76,11.51,11.36,10.31)
for (x in 1:4){
	threshold<-batch.thresholds1.to.4[x]
	batch.lower<-x-0.3
	batch.upper<-x+0.3
	normalized.control.plot<-normalized.control.plot+geom_segment(x=batch.lower,y=threshold,xend=batch.upper,yend=threshold,color="purple",linetype=2)
}

#add in threshold of detectability for the last batch and also add a legend to describe the threshold lines
batch.thresholds5<-11.90
normalized.control.plot<-normalized.control.plot+geom_segment(aes(x=4.7,y=batch.thresholds5,xend=5.3,yend=batch.thresholds5,linetype="Batch Specific\nThreshold of\nDetectability"),color="purple")+
										scale_linetype_manual(name="",values=2)+guides(color=guide_legend(order=1),shape=guide_legend(order = 2),linetype=guide_legend(order = 3))

#add in annotation for the point below the threshold of detectability 
normalized.control.plot.final<-normalized.control.plot+geom_text(aes(x=batch,y=Abundance,label="(Unknown\nAbundance)",shape=NULL,group=NULL),color="Red",annotate.df,size=7)+
					theme(axis.title=element_text(size=35), axis.text=element_text(size=25),  strip.text=element_text(size=30), legend.title=element_text(size=25), legend.text=element_text(size=25))+
					guides(color=guide_legend(override.aes=list(size=2),order=1, keywidth=0.7, keyheight=0.5,default.unit="inch"), shape=guide_legend(override.aes=list(size=5), order=2,  keywidth=0.7, keyheight=0.5,default.unit="inch"), linetype=guide_legend(order=3),  keywidth=0.7, keyheight=0.5,default.unit="inch")
normalized.control.plot.final


###################################################
### code chunk number 15: mixnorm_Vignette.Rnw:225-251
###################################################
#get experimental data into format suitable for plotting
exp.data.copy<-euMetabData[,c("batch","pheno","pyruvic_acid","malonic_acid")]
exp.data.copy$norm_pyruvic_acid<-euMetabNorm$obsNorm[,"pyruvic_acid"]
exp.data.copy$norm_malonic_acid<-euMetabNorm$obsNorm[,"malonic_acid"]
revalue(exp.data.copy$pheno,c("MOM"="Mother Experimental","BABY"="Baby Experimental"))->exp.data.copy$pheno
exp.plot.data<-melt(exp.data.copy,measure.vars=c("pyruvic_acid","malonic_acid","norm_pyruvic_acid","norm_malonic_acid"),variable.name="raw.metab",value.name="Abundance")
exp.plot.data$Type<-factor(ifelse(exp.plot.data$raw.metab %in% c("pyruvic_acid","malonic_acid"),"Before Normalization","After Normalization"),levels=c("Before Normalization","After Normalization"))
exp.plot.data$Metabolite<-ifelse(exp.plot.data$raw.metab %in% c("pyruvic_acid","norm_pyruvic_acid"),"Pyruvic Acid","Malonic Acid")


#check for missing values in metabolite abundance
sum(is.na(exp.plot.data$Abundance))
exp.plot.data[which(is.na(exp.plot.data$Abundance)),]

#for plotting purposes, fill in missing metabolite values with something below the minimum detectable threshold 
thresholds<-c(10.76,11.51,11.36,10.31,11.90)
missing.rows<-row.names(exp.plot.data[which(is.na(exp.plot.data$Abundance)),])
exp.plot.imputed.data<-exp.plot.data
exp.plot.imputed.data[missing.rows,"Abundance"]<-thresholds[exp.plot.imputed.data[missing.rows,"batch"]]-0.5

#add in annotation of mean values for mother and baby samples 
mean.annotation<-ddply(exp.plot.data,c("Type","Metabolite","pheno","batch"),summarize,Abundance=mean(Abundance,na.rm=T))
mean.annotation$point<-"Batch Specific Mean"
exp.plot.imputed.data$point<-"Observed Data"
exp.plot.imp.mean.data<-rbind(exp.plot.imputed.data[,c("batch","pheno","Abundance","Type","Metabolite","point")],
								mean.annotation[,c("batch","pheno","Abundance","Type","Metabolite","point")])


###################################################
### code chunk number 16: fig3
###################################################
#beginning of plot, showing observed data and batch specific means connected by line
exp.plot<-ggplot(exp.plot.imp.mean.data,aes(x=batch,y=Abundance,color=pheno,group=pheno,shape=point,size=point))+
						geom_point()+scale_shape_manual(name="Point Type",values=c(8,16))+
						xlab("Analytic Batch")+ylab("Metabolite Abundance")+facet_grid(Metabolite~Type)+geom_line(aes(x=batch,y=Abundance),size=0.5,data=mean.annotation)+
						scale_size_manual(values=c(7.5,4.5),guide=FALSE)+scale_color_manual(name="Sample Type",values=c( "darkgreen","orange"))					
						
#now add in line segments indicating batch specific thresholds of detectability for the first 4 batches
for (x in 1:4){
	threshold<-batch.thresholds1.to.4[x]
	batch.lower<-x-0.3
	batch.upper<-x+0.3
	exp.plot<-exp.plot+geom_segment(x=batch.lower,y=threshold,xend=batch.upper,yend=threshold,color="purple",linetype=2,size=0.5)
}

#add in threshold of detectability for the last batch and also add a legend to describe the threshold lines
final.exp.plot<-exp.plot+geom_segment(aes(x=4.7,y=batch.thresholds5,xend=5.3,yend=batch.thresholds5,linetype="Batch Specific\nThreshold of\nDetectability"),color="purple",size=0.5)+
				scale_linetype_manual(name="",values=2)+theme(axis.title=element_text(size=35), axis.text=element_text(size=25),  strip.text=element_text(size=30), legend.title=element_text(size=25), legend.text=element_text(size=25))+
				guides(color=guide_legend(override.aes=list(size=2),order=1, keywidth=0.7, keyheight=0.5,default.unit="inch"), shape=guide_legend(override.aes=list(size=5), order=2,  keywidth=0.7, keyheight=0.5,default.unit="inch"), linetype=guide_legend(order=3),  keywidth=0.7, keyheight=0.5,default.unit="inch")
				
				
final.exp.plot


###################################################
### code chunk number 17: fig3
###################################################
#beginning of plot, showing observed data and batch specific means connected by line
exp.plot<-ggplot(exp.plot.imp.mean.data,aes(x=batch,y=Abundance,color=pheno,group=pheno,shape=point,size=point))+
						geom_point()+scale_shape_manual(name="Point Type",values=c(8,16))+
						xlab("Analytic Batch")+ylab("Metabolite Abundance")+facet_grid(Metabolite~Type)+geom_line(aes(x=batch,y=Abundance),size=0.5,data=mean.annotation)+
						scale_size_manual(values=c(7.5,4.5),guide=FALSE)+scale_color_manual(name="Sample Type",values=c( "darkgreen","orange"))					
						
#now add in line segments indicating batch specific thresholds of detectability for the first 4 batches
for (x in 1:4){
	threshold<-batch.thresholds1.to.4[x]
	batch.lower<-x-0.3
	batch.upper<-x+0.3
	exp.plot<-exp.plot+geom_segment(x=batch.lower,y=threshold,xend=batch.upper,yend=threshold,color="purple",linetype=2,size=0.5)
}

#add in threshold of detectability for the last batch and also add a legend to describe the threshold lines
final.exp.plot<-exp.plot+geom_segment(aes(x=4.7,y=batch.thresholds5,xend=5.3,yend=batch.thresholds5,linetype="Batch Specific\nThreshold of\nDetectability"),color="purple",size=0.5)+
				scale_linetype_manual(name="",values=2)+theme(axis.title=element_text(size=35), axis.text=element_text(size=25),  strip.text=element_text(size=30), legend.title=element_text(size=25), legend.text=element_text(size=25))+
				guides(color=guide_legend(override.aes=list(size=2),order=1, keywidth=0.7, keyheight=0.5,default.unit="inch"), shape=guide_legend(override.aes=list(size=5), order=2,  keywidth=0.7, keyheight=0.5,default.unit="inch"), linetype=guide_legend(order=3),  keywidth=0.7, keyheight=0.5,default.unit="inch")
				
				
final.exp.plot


###################################################
### code chunk number 18: mixnorm_Vignette.Rnw:290-293
###################################################
euMetabNormRC <- mixnorm("pyruvic_acid", batch="batch", mxtrModel=~pheno+batch|pheno+batch, 
				batchTvals=c(10.76,11.51,11.36,10.31,11.90), cData=euMetabCData,
				removeCorrection="pheno",data=euMetabData)


###################################################
### code chunk number 19: mixnorm_Vignette.Rnw:298-299
###################################################
head(euMetabNormRC$normParamsZ)


###################################################
### code chunk number 20: mixnorm_Vignette.Rnw:304-326
###################################################
#get experimental data into format suitable for plotting
exp.data.copy2<-euMetabData[,c("batch","pheno","pyruvic_acid")]
exp.data.copy2$norm_pyruvic_acid<-euMetabNormRC$obsNorm[,"pyruvic_acid"]
revalue(exp.data.copy2$pheno,c("MOM"="Mother Experimental","BABY"="Baby Experimental"))->exp.data.copy2$pheno
exp.plot.data2<-melt(exp.data.copy2,measure.vars=c("pyruvic_acid","norm_pyruvic_acid"),variable.name="raw.metab",value.name="Abundance")
exp.plot.data2$Type<-factor(ifelse(exp.plot.data2$raw.metab=="pyruvic_acid","Before Normalization","After Normalization"),levels=c("Before Normalization","After Normalization"))

#check for missing values in metabolite abundance
sum(is.na(exp.plot.data2$Abundance))
exp.plot.data2[which(is.na(exp.plot.data2$Abundance)),]

#for plotting purposes, fill in missing metabolite values with something below the minimum detectable threshold 
missing.rows2<-row.names(exp.plot.data2[which(is.na(exp.plot.data2$Abundance)),])
exp.plot.imputed.data2<-exp.plot.data2
exp.plot.imputed.data2[missing.rows2,"Abundance"]<-thresholds[exp.plot.imputed.data2[missing.rows2,"batch"]]-0.5

#add in annotation of mean values for mother and baby samples 
mean.annotation2<-ddply(exp.plot.data2,c("Type","pheno","batch"),summarize,Abundance=mean(Abundance,na.rm=T))
mean.annotation2$point<-"Batch Specific Mean"
exp.plot.imputed.data2$point<-"Observed Data"
exp.plot.imp.mean.data2<-rbind(exp.plot.imputed.data2[,c("batch","pheno","Abundance","Type","point")],
								mean.annotation2[,c("batch","pheno","Abundance","Type","point")])


###################################################
### code chunk number 21: fig4
###################################################
#beginning of plot, showing observed data and batch specific means connected by line
exp.plot2<-ggplot(exp.plot.imp.mean.data2,aes(x=batch,y=Abundance,color=pheno,group=pheno,shape=point,size=point))+
						geom_point(size=7)+scale_color_discrete(name="Sample Type")+scale_shape_manual(name="Point Type",values=c(8,16))+
						xlab("Analytic Batch")+ylab("Pyruvic Acid Abundance")+facet_grid(~Type)+geom_line(aes(x=batch,y=Abundance),size=0.5,data=mean.annotation2)+
						scale_size_manual(values=c(10,4.5),guide=FALSE)+scale_color_manual(name="Sample Type",values=c( "darkgreen","orange"))
						
						
#now add in line segments indicating batch specific thresholds of detectability for the first 4 batches
for (x in 1:4){
	threshold<-batch.thresholds1.to.4[x]
	batch.lower<-x-0.3
	batch.upper<-x+0.3
	exp.plot2<-exp.plot2+geom_segment(x=batch.lower,y=threshold,xend=batch.upper,yend=threshold,color="purple",linetype=2,size=0.5)
}

#add in threshold of detectability for the last batch and also add a legend to describe the threshold lines
final.exp.plot2<-exp.plot2+geom_segment(aes(x=4.7,y=batch.thresholds5,xend=5.3,yend=batch.thresholds5,linetype="Batch Specific\nThreshold of\nDetectability"),color="purple",size=0.5)+
				scale_linetype_manual(name="",values=2)+theme(axis.title=element_text(size=35), axis.text=element_text(size=25),  strip.text=element_text(size=30), legend.title=element_text(size=25), legend.text=element_text(size=25))+
				guides(color=guide_legend(override.aes=list(size=2),order=1, keywidth=0.7, keyheight=0.5,default.unit="inch"), shape=guide_legend(override.aes=list(size=5), order=2,  keywidth=0.7, keyheight=0.5,default.unit="inch"), linetype=guide_legend(order=3),  keywidth=0.7, keyheight=0.5,default.unit="inch")
final.exp.plot2


###################################################
### code chunk number 22: fig4
###################################################
#beginning of plot, showing observed data and batch specific means connected by line
exp.plot2<-ggplot(exp.plot.imp.mean.data2,aes(x=batch,y=Abundance,color=pheno,group=pheno,shape=point,size=point))+
						geom_point(size=7)+scale_color_discrete(name="Sample Type")+scale_shape_manual(name="Point Type",values=c(8,16))+
						xlab("Analytic Batch")+ylab("Pyruvic Acid Abundance")+facet_grid(~Type)+geom_line(aes(x=batch,y=Abundance),size=0.5,data=mean.annotation2)+
						scale_size_manual(values=c(10,4.5),guide=FALSE)+scale_color_manual(name="Sample Type",values=c( "darkgreen","orange"))
						
						
#now add in line segments indicating batch specific thresholds of detectability for the first 4 batches
for (x in 1:4){
	threshold<-batch.thresholds1.to.4[x]
	batch.lower<-x-0.3
	batch.upper<-x+0.3
	exp.plot2<-exp.plot2+geom_segment(x=batch.lower,y=threshold,xend=batch.upper,yend=threshold,color="purple",linetype=2,size=0.5)
}

#add in threshold of detectability for the last batch and also add a legend to describe the threshold lines
final.exp.plot2<-exp.plot2+geom_segment(aes(x=4.7,y=batch.thresholds5,xend=5.3,yend=batch.thresholds5,linetype="Batch Specific\nThreshold of\nDetectability"),color="purple",size=0.5)+
				scale_linetype_manual(name="",values=2)+theme(axis.title=element_text(size=35), axis.text=element_text(size=25),  strip.text=element_text(size=30), legend.title=element_text(size=25), legend.text=element_text(size=25))+
				guides(color=guide_legend(override.aes=list(size=2),order=1, keywidth=0.7, keyheight=0.5,default.unit="inch"), shape=guide_legend(override.aes=list(size=5), order=2,  keywidth=0.7, keyheight=0.5,default.unit="inch"), linetype=guide_legend(order=3),  keywidth=0.7, keyheight=0.5,default.unit="inch")
final.exp.plot2


###################################################
### code chunk number 23: mixnorm_Vignette.Rnw:364-365
###################################################
toLatex(sessionInfo())