inst/doc/mixnorm_Vignette.R

### 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-142
###################################################
#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, qc.sd.outliers=2)


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


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


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


###################################################
### code chunk number 12: mixnorm_Vignette.Rnw:167-168
###################################################
head(euMetabNorm$conv)


###################################################
### code chunk number 13: fig2
###################################################
plot.list<-lapply(ynames, metabplot, batch="batch", raw.obs.data=euMetabData, raw.cont.data=euMetabCData, 
			norm.obs.data=euMetabNorm$obsNorm, norm.cont.data=euMetabNorm$ctlNorm,
			 color.var="pheno", cont.outlier.sd.thresh=2, norm.outlier.sd.thresh=4)

#just show plot for one of the metabolites 
plot.list[[2]]


###################################################
### code chunk number 14: fig2
###################################################
plot.list<-lapply(ynames, metabplot, batch="batch", raw.obs.data=euMetabData, raw.cont.data=euMetabCData, 
			norm.obs.data=euMetabNorm$obsNorm, norm.cont.data=euMetabNorm$ctlNorm,
			 color.var="pheno", cont.outlier.sd.thresh=2, norm.outlier.sd.thresh=4)

#just show plot for one of the metabolites 
plot.list[[2]]


###################################################
### code chunk number 15: mixnorm_Vignette.Rnw:198-201
###################################################
euMetabNormRC <- mixnorm(ynames, 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 16: mixnorm_Vignette.Rnw:206-207
###################################################
euMetabNormRC$normParamsZ[rownames(euMetabNormRC$normParamsZ)=="pyruvic_acid", ]


###################################################
### code chunk number 17: fig4
###################################################
metabplot("pyruvic_acid",  batch="batch", raw.obs.data=euMetabData, raw.cont.data=euMetabCData, 
		norm.obs.data=euMetabNormRC$obsNorm, norm.cont.data=euMetabNormRC$ctlNorm, 
		color.var="pheno", cont.outlier.sd.thresh=2, norm.outlier.sd.thresh=4)


###################################################
### code chunk number 18: fig4
###################################################
metabplot("pyruvic_acid",  batch="batch", raw.obs.data=euMetabData, raw.cont.data=euMetabCData, 
		norm.obs.data=euMetabNormRC$obsNorm, norm.cont.data=euMetabNormRC$ctlNorm, 
		color.var="pheno", cont.outlier.sd.thresh=2, norm.outlier.sd.thresh=4)


###################################################
### code chunk number 19: mixnorm_Vignette.Rnw:232-235
###################################################
norm.with.outliers <- 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, qc.sd.outliers=Inf)


###################################################
### code chunk number 20: fig5
###################################################
metabplot("malonic_acid",  batch="batch", raw.obs.data=euMetabData, raw.cont.data=euMetabCData, 
		norm.obs.data=norm.with.outliers$obsNorm, norm.cont.data=norm.with.outliers$ctlNorm, 
		color.var="pheno", cont.outlier.sd.thresh=Inf, norm.outlier.sd.thresh=4)


###################################################
### code chunk number 21: fig5
###################################################
metabplot("malonic_acid",  batch="batch", raw.obs.data=euMetabData, raw.cont.data=euMetabCData, 
		norm.obs.data=norm.with.outliers$obsNorm, norm.cont.data=norm.with.outliers$ctlNorm, 
		color.var="pheno", cont.outlier.sd.thresh=Inf, norm.outlier.sd.thresh=4)


###################################################
### code chunk number 22: mixnorm_Vignette.Rnw:257-262
###################################################
cData.missing.batch<-euMetabCData
cData.missing.batch[cData.missing.batch$batch==2, "malonic_acid"]<-NA				
norm.with.missing.batch <- mixnorm(ynames, batch="batch", mxtrModel=~pheno+batch|pheno+batch, 
				batchTvals=c(10.76,11.51,11.36,10.31,11.90), cData=cData.missing.batch,
				data=euMetabData)	


###################################################
### code chunk number 23: mixnorm_Vignette.Rnw:267-268
###################################################
norm.with.missing.batch$normParamsZ[rownames(norm.with.missing.batch$normParamsZ)=="malonic_acid", ]


###################################################
### code chunk number 24: mixnorm_Vignette.Rnw:273-274
###################################################
norm.with.missing.batch$conv


###################################################
### code chunk number 25: fig6
###################################################
metabplot("malonic_acid",  batch="batch", raw.obs.data=euMetabData, raw.cont.data=cData.missing.batch,
		 norm.obs.data=norm.with.missing.batch$obsNorm, 
		 norm.cont.data=norm.with.missing.batch$ctlNorm,  color.var="pheno")


###################################################
### code chunk number 26: fig6
###################################################
metabplot("malonic_acid",  batch="batch", raw.obs.data=euMetabData, raw.cont.data=cData.missing.batch,
		 norm.obs.data=norm.with.missing.batch$obsNorm, 
		 norm.cont.data=norm.with.missing.batch$ctlNorm,  color.var="pheno")


###################################################
### code chunk number 27: mixnorm_Vignette.Rnw:299-308
###################################################
cData.missing.pheno<-euMetabCData
cData.missing.pheno[cData.missing.pheno$pheno=="BABY", "malonic_acid"]<-NA				
norm.missing.pheno <- mixnorm(ynames, batch="batch", mxtrModel=~pheno+batch|pheno+batch, 
				batchTvals=c(10.76,11.51,11.36,10.31,11.90), cData=cData.missing.pheno,
				data=euMetabData)	

norm.missing.pheno$normParamsZ[rownames(norm.missing.pheno$normParamsZ)=="malonic_acid", ]
norm.missing.pheno$conv



###################################################
### code chunk number 28: fig7
###################################################
metabplot("malonic_acid",  batch="batch", raw.obs.data=euMetabData, raw.cont.data=cData.missing.pheno, 
		norm.obs.data=norm.missing.pheno$obsNorm, norm.cont.data=norm.missing.pheno$ctlNorm,
		 color.var="pheno")


###################################################
### code chunk number 29: fig7
###################################################
metabplot("malonic_acid",  batch="batch", raw.obs.data=euMetabData, raw.cont.data=cData.missing.pheno, 
		norm.obs.data=norm.missing.pheno$obsNorm, norm.cont.data=norm.missing.pheno$ctlNorm,
		 color.var="pheno")


###################################################
### code chunk number 30: mixnorm_Vignette.Rnw:329-330
###################################################
toLatex(sessionInfo())

Try the metabomxtr package in your browser

Any scripts or data that you put into this service are public.

metabomxtr documentation built on Nov. 8, 2020, 6:50 p.m.