Mercurial > repos > mmonsoor > camera_combinexsannos

comparison lib.r @ 4:87570e9b71f5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 9d47e3b467dbbe0af0d90a937c5e9f2c4b958c4b
author lecorguille
date Mon, 25 Apr 2016 11:07:23 -0400
parents
children 198b035d4848
comparison
equal deleted inserted replaced
3:4ca5c7bbc6cf 4:87570e9b71f5
1 # lib.r version="2.2.1"
2
3 #The function create a pdf from the different png generated by diffreport
4 diffreport_png2pdf <- function(filebase, new_file_path) {
5
6 pdfEicOutput = paste(new_file_path,filebase,"-eic_visible_pdf",sep="")
7 pdfBoxOutput = paste(new_file_path,filebase,"-box_visible_pdf",sep="")
8
9 system(paste("gm convert ",filebase,"_eic/*.png ",filebase,"_eic.pdf",sep=""))
10 system(paste("gm convert ",filebase,"_box/*.png ",filebase,"_box.pdf",sep=""))
11
12 file.copy(paste(filebase,"_eic.pdf",sep=""), pdfEicOutput)
13 file.copy(paste(filebase,"_box.pdf",sep=""), pdfBoxOutput)
14 }
15
16 #The function annotateDiffreport without the corr function which bugs
17 annotatediff <- function(xset=xset, listArguments=listArguments, variableMetadataOutput="variableMetadata.tsv", dataMatrixOutput="dataMatrix.tsv",new_file_path=NULL) {
18 # Resolve the bug with x11, with the function png
19 options(bitmapType='cairo')
20
21 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
22 res=try(is.null(xset@filled))
23
24 # ------ annot -------
25 listArguments[["calcCiS"]]=as.logical(listArguments[["calcCiS"]])
26 listArguments[["calcIso"]]=as.logical(listArguments[["calcIso"]])
27 listArguments[["calcCaS"]]=as.logical(listArguments[["calcCaS"]])
28
29 #graphMethod parameter bugs where this parameter is not defined in quick=true
30 if(listArguments[["quick"]]==TRUE) {
31 xa= annotate(object=xset,nSlaves=1,sigma=listArguments[["sigma"]],perfwhm=listArguments[["perfwhm"]],maxcharge=listArguments[["maxcharge"]],maxiso=listArguments[["maxiso"]],minfrac=listArguments[["minfrac"]],ppm=listArguments[["ppm"]],mzabs=listArguments[["mzabs"]],quick=listArguments[["quick"]],polarity=listArguments[["polarity"]],max_peaks=listArguments[["max_peaks"]],intval=listArguments[["intval"]])
32 }
33 else {
34 xa= annotate(object=xset,nSlaves=1,sigma=listArguments[["sigma"]],perfwhm=listArguments[["perfwhm"]],graphMethod=listArguments[["graphMethod"]],cor_eic_th=listArguments[["cor_eic_th"]],pval=listArguments[["pval"]],calcCiS=listArguments[["calcCiS"]],calcIso=listArguments[["calcIso"]],calcCaS=listArguments[["calcCaS"]],multiplier=listArguments[["multiplier"]],maxcharge=listArguments[["maxcharge"]],maxiso=listArguments[["maxiso"]],minfrac=listArguments[["minfrac"]],ppm=listArguments[["ppm"]],mzabs=listArguments[["mzabs"]],quick=listArguments[["quick"]],polarity=listArguments[["polarity"]],max_peaks=listArguments[["max_peaks"]],intval=listArguments[["intval"]])
35
36 }
37 peakList=getPeaklist(xa,intval=listArguments[["intval"]])
38 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
39
40
41 # --- Multi condition : diffreport ---
42 diffrep=NULL
43 if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))>=2) {
44 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
45 res=try(is.null(xset@filled))
46 classes=levels(sampclass(xset))
47 x=1:(length(classes)-1)
48 for (i in seq(along=x) ) {
49 y=1:(length(classes))
50 for (n in seq(along=y)){
51 if(i+n <= length(classes)){
52 filebase=paste(classes[i],class2=classes[i+n],sep="-vs-")
53
54 diffrep=diffreport(object=xset,class1=classes[i],class2=classes[i+n],filebase=filebase,eicmax=listArguments[["eicmax"]],eicwidth=listArguments[["eicwidth"]],sortpval=TRUE,value=listArguments[["value"]],h=listArguments[["h"]],w=listArguments[["w"]],mzdec=listArguments[["mzdec"]])
55 #combines results
56 diffreportTSV=merge(peakList, diffrep[,c("name","fold","tstat","pvalue")], by.x="name", by.y="name", sort=F)
57 diffreportTSV=cbind(diffreportTSV[,!(colnames(diffreportTSV) %in% c(sampnames(xa@xcmsSet)))],diffreportTSV[,(colnames(diffreportTSV) %in% c(sampnames(xa@xcmsSet)))])
58
59 if(listArguments[["sortpval"]]){
60 diffreportTSV=diffreportTSV[order(diffreportTSV$pvalue), ]
61 }
62
63 if (listArguments[["convert_param"]]){
64 #converting the retention times (seconds) into minutes
65 diffreportTSV$rt=diffreportTSV$rt/60;diffreportTSV$rtmin=diffreportTSV$rtmin/60; diffreportTSV$rtmax=diffreportTSV$rtmax/60;
66 }
67 write.table(diffreportTSV, sep="\t", quote=FALSE, row.names=FALSE, file=paste(new_file_path,filebase,"-tabular_visible_tabular",sep=""))
68
69 if (listArguments[["eicmax"]] != 0) {
70 diffreport_png2pdf(filebase, new_file_path)
71 }
72 }
73 }
74 }
75 }
76
77
78
79
80 # --- variableMetadata ---
81 variableMetadata=peakList[,!(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
82 # if we have 2 conditions, we keep stat of diffrep
83 if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))==2) {
84 variableMetadata = merge(variableMetadata, diffrep[,c("name","fold","tstat","pvalue")],by.x="name", by.y="name", sort=F)
85 if(exists("listArguments[[\"sortpval\"]]")){
86 variableMetadata=variableMetadata[order(variableMetadata$pvalue), ]
87 }
88 }
89
90 variableMetadataOri=variableMetadata
91 if (listArguments[["convert_param"]]){
92 #converting the retention times (seconds) into minutes
93 print("converting the retention times into minutes in the variableMetadata")
94 variableMetadata$rt=variableMetadata$rt/60;variableMetadata$rtmin=variableMetadata$rtmin/60; variableMetadata$rtmax=variableMetadata$rtmax/60;
95 }
96 #Transform metabolites name
97 variableMetadata$name= paste("M",round(variableMetadata$mz,digits=listArguments[["num_digits"]]),"T",round(variableMetadata$rt),sep="")
98 write.table(variableMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=variableMetadataOutput)
99
100 # --- dataMatrix ---
101 dataMatrix = peakList[,(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
102 dataMatrix=cbind(peakList$name,dataMatrix); colnames(dataMatrix)[1] = c("name");
103
104 if (listArguments[["convert_param"]]){
105 #converting the retention times (seconds) into minutes
106 print("converting the retention times into minutes in the dataMatrix ids")
107 peakList$rt=peakList$rt/60
108 }
109 dataMatrix$name= paste("M",round(peakList$mz,digits=listArguments[["num_digits"]]),"T",round(peakList$rt),sep="")
110 write.table(dataMatrix, sep="\t", quote=FALSE, row.names=FALSE, file=dataMatrixOutput)
111
112 return(list("xa"=xa,"diffrep"=diffrep,"variableMetadata"=variableMetadataOri));
113
114 }
115
116
117 combinexsAnnos_function <- function(xaP, xaN, listOFlistArgumentsP,listOFlistArgumentsN, diffrepP=NULL,diffrepN=NULL,convert_param=FALSE,pos=TRUE,tol=2,ruleset=NULL,keep_meta=TRUE, variableMetadataOutput="variableMetadata.tsv"){
118
119 #Load the two Rdata to extract the xset objects from positive and negative mode
120 cat("\tObject xset from positive mode\n")
121 print(xaP)
122 cat("\n")
123
124 cat("\tObject xset from negative mode\n")
125 print(xaN)
126 cat("\n")
127
128 cat("\n")
129 cat("\tCombining...\n")
130 #Convert the string to numeric for creating matrix
131 row=as.numeric(strsplit(ruleset,",")[[1]][1])
132 column=as.numeric(strsplit(ruleset,",")[[1]][2])
133 ruleset=cbind(row,column)
134 #Test if the file comes from an older version tool
135 if ((!is.null(xaP)) & (!is.null(xaN))) {
136 #Launch the combinexsannos function from CAMERA
137 cAnnot=combinexsAnnos(xaP, xaN,pos=pos,tol=tol,ruleset=ruleset)
138 } else {
139 stop("You must relauch the CAMERA.annotate step with the lastest version.")
140 }
141
142
143
144 if(pos){
145 xa=xaP
146 listOFlistArgumentsP=listOFlistArguments
147 mode="neg. Mode"
148 } else {
149 xa=xaN
150 listOFlistArgumentsN=listOFlistArguments
151 mode="pos. Mode"
152 }
153 intval = "into"; for (steps in names(listOFlistArguments)) { if (!is.null(listOFlistArguments[[steps]]$intval)) intval = listOFlistArguments[[steps]]$intval }
154 peakList=getPeaklist(xa,intval=intval)
155 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
156 variableMetadata=cbind(peakList, cAnnot[, c("isotopes", "adduct", "pcgroup",mode)]);
157 variableMetadata=variableMetadata[,!(colnames(variableMetadata) %in% c(sampnames(xa@xcmsSet)))]
158
159 #Test if there are more than two classes (conditions)
160 if ( nlevels(sampclass(xaP@xcmsSet))==2 & (!is.null(diffrepN)) & (!is.null(diffrepP))) {
161 diffrepP = diffrepP[,c("name","fold","tstat","pvalue")]; colnames(diffrepP) = paste("P.",colnames(diffrepP),sep="")
162 diffrepN = diffrepN[,c("name","fold","tstat","pvalue")]; colnames(diffrepN) = paste("N.",colnames(diffrepN),sep="")
163
164 variableMetadata = merge(variableMetadata, diffrepP, by.x="name", by.y="P.name")
165 variableMetadata = merge(variableMetadata, diffrepN, by.x="name", by.y="N.name")
166 }
167
168 rownames(variableMetadata) = NULL
169 #TODO: checker
170 #colnames(variableMetadata)[1:2] = c("name","mz/rt");
171
172 #If the user want to convert the retention times (seconds) into minutes.
173 if (listArguments[["convert_param"]]){
174 #converting the retention times (seconds) into minutes
175 cat("\tConverting the retention times into minutes\n")
176 variableMetadata$rtmed=cAnnot$rt/60; variableMetadata$rtmin=cAnnot$rtmin/60; variableMetadata$rtmax=cAnnot$rtmax/60;
177 }
178
179 #If the user want to keep only the metabolites which match a difference
180 if(keep_meta){
181 variableMetadata=variableMetadata[variableMetadata[,c(mode)]!="",]
182 }
183
184 #Write the output into a tsv file
185 write.table(variableMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=variableMetadataOutput)
186 return(variableMetadata);
187
188 }