R/eval.target.mz.R

In kuppal2/xMSanalyzer: xMSanalyzer: R package for data analysis, comparison, and annotation of metabolomics data.

eval.target.mz <-
function(dataA,refMZ,feature.eval.result,mzthresh=10,timethresh=NA,outloc,folderheader=NA,alignment.tool=NA){

feature09<-dataA
if(is.na(refMZ[1,1])==FALSE){
stddata<-refMZ
}else{
	Name<-paste("mz",seq(1,dim(feature09)[1]),sep="")
	stddata<-cbind(feature09[,c(1)],Name)
	
	}
qcresults<-feature.eval.result

rm(dataA)
outloc1<-paste(outloc,"/",folderheader,"targetedeval",mzthresh,"ppm/",sep="")
dir.create(outloc1,showWarnings=FALSE)
setwd(outloc1)

col.names.dataA<-colnames(feature09)

if(is.na(alignment.tool)==FALSE){
	 if(alignment.tool=="apLCMS")
        {
              sample.col.start=5
        }
        else
        {
              if(alignment.tool=="XCMS")
              {
              	
                    sample.col.start=9
                    col.names.dataA[1]="mz"
                    col.names.dataA[2]="time"
		    
                    colnames(data_a)=col.names.dataA
                   
              }
	      
	}}else{
                    #stop(paste("Invalid value for alignment.tool. Please use either \"apLCMS\" or \"XCMS\"", sep=""))
		   
		    col.names.dataA[1]="mz"
		   
		   
		     col.names.dataA[2]="time"
		      sample.col.start=3
		     		      print("Using the 1st columns as \"mz\" and 2nd columsn as \"retention time\"")
		     
		    colnames(feature09)=col.names.dataA
                   
		}
		
par(mfrow=c(2,2))


overlapres5ppm<-getVenn(dataA=feature09, name_a=paste("Expdata",sep=""), dataB=stddata, name_b="Target", mz.thresh = mzthresh, time.thresh=timethresh,alignment.tool="apLCMS",
xMSanalyzer.outloc=outloc1,plotvenn=FALSE)

save(overlapres5ppm,file="overlapRes.Rda")
match5ppmdata<-feature09[overlapres5ppm$common$index.A,]

name_mz<-{}
if(dim(match5ppmdata)[1]<1){
print("No matches found for targeted metabolites.")
return(name_mz);
}


qcresults5ppmdata<-qcresults[overlapres5ppm$common$index.A,]


#m1<-cbind(match5ppmdata[,c(1:(sample.col.start-1))],qcresults5ppmdata[,c(3,7,8)],match5ppmdata[,-c(1:(sample.col.start-1))])

fnames<-paste("../../Stage3b/Targeted_feature_table_",mzthresh,"ppm_filtered.txt",sep="")
#write.table(m1,file=fnames,sep="\t",quote=F,col.name=T,row.names=F)

#write.table(match5ppmdata,file="Matching_data.txt",sep="\t",quote=F,col.name=T,row.names=F)

fnames<-paste("Boxplot_sampleintensity_usingtargetmatches",mzthresh,"ppm.tiff",sep="")


int_data<-log10(match5ppmdata[,-c(1:(sample.col.start-1))]+1)
#int_data<-match5ppmdata[,-c(1:(sample.col.start-1))]
main_lab<-paste("Intensity distribution (log10; all samples) \n of each m/z matching targets (+/- ",mzthresh,"ppm)",sep="")


vennDiagram(circle.col="red",overlapres5ppm$vennCounts,counts.col="blue")

#(fnames,width=2000,height=2000,res=300)
#	boxplot(int_data,cex.names=0.35,cex.axis=1,main=main_lab, ylab="Intensity",xlab="Sample",col="orange",cex.main=0.7)
#	dev.off()

	
	

cv5ppm<-apply(match5ppmdata[,-c(1:(sample.col.start-1))],1,function(x){
	x<-replace(x,which(x==0),NA)
	cvres<-100*sd(x,na.rm=TRUE)/mean(x,na.rm=TRUE)
	return(cvres)
	})
names(cv5ppm)<-round(match5ppmdata[,1],5)

fnames<-paste("TotalCVallsamples_refmz",mzthresh,"ppm.tiff",sep="")

#pdf("Targeted_mz_QC.pdf")
main_name<-paste("Total CV (across all samples) \n of each m/z matching targets (+/- ",mzthresh,"ppm)",sep="")

if(length(cv5ppm)>0){
#tiff(fnames,width=2000,height=2000,res=300)
barplot(cv5ppm,cex.names=0.4,cex.axis=1,main=main_name,col="orange",cex.main=0.7)
}
#dev.off()
fnames<-paste("BarplotTIC_targetmatches",mzthresh,"ppm.tiff",sep="")

main_name<-paste("TIC per sample for matching targets (+/- ",mzthresh,"ppm)",sep="")
#tiff(fnames,width=2000,height=2000,res=300)

tic<-apply(match5ppmdata[,-c(1:(sample.col.start-1))],2,function(x){
	x<-replace(x,which(x==0),NA)
	res<-sum(x,na.rm=TRUE)
	return(res)
})

if(length(tic)>0){
barplot(tic,cex.names=0.4,cex.axis=1,main=main_name,col="orange",cex.main=0.7)
}

fnames<-paste("Pairwiseplot_overall",mzthresh,"ppm.tiff",sep="")
#tiff(fnames,width=2000,height=2000,res=300)
if(length(overlapres5ppm$common$index.A)>0){
plot(cbind(feature09[overlapres5ppm$common$index.A,c(1:2)],qcresults[overlapres5ppm$common$index.A,c(3,7,8)]),main="Pairwise plots of m/z, time, CV, Qscore",cex.main=0.7)
}
#dev.off()



#dev.off()



mean_tic<-mean(tic,na.rm=TRUE)

#cv_tic<-100*sd(tic,na.rm=TRUE)/mean(tic,na.rm=TRUE)

#tic_res<-cbind(mean_tic,cv_tic)

#colnames(tic_res)<-c("Average_TIC","CV_TIC")

#fnames<-paste("TICperfeat_targetmatches",mzthresh,"ppm.txt",sep="")
#write.table(tic_res,file=fnames,sep="\t",quote=F,col.name=T,row.names=F)

tic<-as.data.frame(tic)
names(tic)<-c("sample_TIC")

fnames<-paste("TICpersamp_targetmatches",mzthresh,"ppm.txt",sep="")
write.table(tic,file=fnames,sep="\t",quote=F,col.name=T,row.names=TRUE)


#refMZ<-cbind(refMZ[overlapres5ppm$common$index.B,],match5ppmdata)

#colnames(refMZ)<-c("mz","Name")

name_mz<-cbind(match5ppmdata[,c(1:2)],stddata[overlapres5ppm$common$index.B,],match5ppmdata[,-c(1:2)])

name_mz<-as.data.frame(name_mz)

#name_mz<-merge(refMZ,match5ppmdata,by="mz")

if(folderheader=="raw"){
fnames<-paste("../../Stage3b/",folderheader,"_targetfeaturetable",mzthresh,"ppm_final.txt",sep="")
}else{
	
	if(folderheader=="ComBat"){
	fnames<-paste("../../Stage4b/",folderheader,"_targetfeaturetable",mzthresh,"ppm_final.txt",sep="")
	}
}
#fnames<-paste("Targetmz_matching_data",mzthresh,"ppm.txt",sep="")
write.table(name_mz,file=fnames,sep="\t",quote=F,col.name=T,row.names=F)


#outloc1<-paste(outloc1,"target_barplots/",sep="")
#dir.create(outloc1,showWarnings=FALSE)
#setwd(outloc1)

match5ppmdata<-name_mz
match5ppmdata<-as.data.frame(match5ppmdata)
match5ppmdata_int<-match5ppmdata[,-c(1:(2+dim(refMZ)[2]))]
match5ppmdata_int<-apply(match5ppmdata_int,2,as.numeric)

print(head(match5ppmdata))
print(head(match5ppmdata_int))
print(dim(match5ppmdata))
print(dim(match5ppmdata_int))
if(length(overlapres5ppm$common$index.B)>0){
for(i in 1:dim(match5ppmdata_int)[1]){
	fname<-paste("mz",round(match5ppmdata[i,1],5),".tiff",sep="")
#	tiff(fname,width=2000,height=2000,res=300)
delta_ppm<-10^6*(abs(match5ppmdata[i,1]-match5ppmdata[i,3]))/match5ppmdata[i,3]
delta_ppm<-round(delta_ppm,2)

#int_val<-as.numeric(int_val)
if(length(overlapres5ppm$common$index.B)>1){
cur_vec<-t(match5ppmdata_int[i,])
}else{
cur_vec<-t(match5ppmdata_int[i])
}
	plot(x=t(cur_vec),main=paste("mz: ",round(match5ppmdata[i,1],5),"; time:",round(match5ppmdata[i,2],0), "(s) ;\n",match5ppmdata$Name[i], " (delta=",delta_ppm," ppm)",sep=""),ylab="Intensity",xlab="Sample index",cex=0.6,cex.main=0.7,col="brown")
	#dev.off()
}

}
#dev.off()

return(name_mz)
}