Nothing
### 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())
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.