annotate BC/batch_correction_3Lfct.R @ 4:23314e1192d4 draft default tip

Uploaded
author melpetera
date Thu, 14 Jan 2021 09:56:58 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
4
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
1 # Author: jfmartin
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
2 # Modified by : mpetera
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
3 ###############################################################################
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
4 # Correction of analytical effects inter and intra batch on intensities using quality control pooled samples (QC-pools)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
5 # according to the algorithm mentioned by Van der Kloet (J Prot Res 2009).
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
6 # Parameters : a dataframe of Ions intensities and an other of samples? metadata which must contains at least the three following columns :
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
7 # "batch" to identify the batches of analyses ; need at least 3 QC-pools for linear adjustment and 8 for lo(w)ess adjustment
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
8 # "injectionOrder" integer defining the injection order of all samples : QC-pools and analysed samples
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
9 # "sampleType" indicates if defining a sample with "sample" or a QC-pool with "pool"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
10 # NO MISSING DATA are allowed
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
11 # Version 0.91 insertion of ok_norm function to assess correction feasibility
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
12 # Version 0.92 insertion of slope test in ok_norm
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
13 # Version 0.93 name of log file define as a parameter of the correction function
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
14 # Version 0.94 Within a batch, test if all QCpools or samples values = 0. Definition of an error code in ok_norm function (see function for details)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
15 # Version 0.99 include non linear lowess correction.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
16 # Version 1.00 the corrected result matrix is return transposed in Galaxy
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
17 # Version 1.01 standard deviation=0 instead of sum of value=0 is used to assess constant data in ok_norm function. Negative values in corrected matrix are converted to 0.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
18 # Version 1.02 plotsituation create a result file with the error code of non execution of correction set by function ok_norm
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
19 # Version 1.03 fix bug in plot with "reg" option. suppression of ok_norm=4 condition if ok_norm function
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
20 # Version 2.00 Addition of loess function, correction indicator, plots ; modification of returned objects' format, some plots' displays and ok_norm ifelse format
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
21 # Version 2.01 Correction for pools negative values earlier in norm_QCpool
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
22 # Version 2.10 Script refreshing ; vocabulary adjustment ; span in parameters for lo(w)ess regression ; conditionning for third line ACP display ; order in loess display
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
23 # Version 2.11 ok1 and ok2 permutation (ok_norm) ; conditional display of regression (plotsituation) ; grouping of linked lignes + conditioning (normX) ; conditioning for CVplot
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
24 # Version 2.20 acplight function added from previous toolBox.R [# Version 1.01 "NA"-coding possibility added in acplight function]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
25 # Version 2.30 addition of suppressWarnings() for known and controlled warnings ; suppression of one useless "cat" message ; change in Rdata names ; 'batch(es)' in cat
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
26 # Version 2.90 change in handling of generated negative and Inf values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
27 # Version 2.91 Plot improvement
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
28 # Version 3.00 - handling of sample tags' parameters
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
29 # - accepting sample types beyond "pool" and "sample"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
30 # - dealing with NA
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
31 # - changes in the normalisation strategy regarding mean values to adjust for NA or 0 values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
32 # - changes in the normalisation strategy regarding unconsistant values (negative or Inf)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
33
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
34 ok_norm=function(qcp,qci,spl,spi,method,normref=NA,valimp="0") {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
35 # Function used for one ion within one batch to determine whether or not batch correction is possible
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
36 # ok_norm values :
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
37 # 0 : no preliminary-condition problem
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
38 # 1 : standard deviation of QC-pools or samples = 0
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
39 # 2 : insufficient number of QC-pools within a batch (n=3 for linear, n=8 for lowess or loess)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
40 # 2.5 : less than 2 samples within a batch
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
41 # 3 : significant difference between QC-pools' and samples' means
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
42 # 4 : denominator =0 when on 1 pool per batch <> 0
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
43 # 5 : (linear regression only) the slopes ratio ?QC-pools/samples? is lower than -0.2
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
44 # 6 : (linear regression only) none of the pool or sample could be corrected if negative and infinite values are turned into NA
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
45 # Parameters:
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
46 # qcp: intensity of a given ion for pools
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
47 # qci: injection numbers for pools
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
48 # spl: intensity of a given ion for samples
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
49 # spi: injection numbers for samples
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
50 # method: to provide specific checks for "linear"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
51
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
52 ok=0
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
53 if (method=="linear") {minQC=3} else {minQC=8}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
54 if (length(qcp[!is.na(qcp)])<minQC) { ok=2 } else { if (length(spl[!is.na(spl)])<2) { ok=2.5
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
55 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
56 if (sd(qcp,na.rm=TRUE)==0 | sd(spl,na.rm=TRUE)==0) { ok=1
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
57 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
58 cvp= sd(qcp,na.rm=TRUE)/mean(qcp,na.rm=TRUE); cvs=sd(spl,na.rm=TRUE)/mean(spl,na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
59 rttest=t.test(qcp,y=spl)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
60 reslsfit=lsfit(qci, qcp)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
61 reslsfitSample=lsfit(spl, spi)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
62 ordori=reslsfit$coefficients[1]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
63 penteB=reslsfit$coefficients[2]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
64 penteS=reslsfitSample$coefficients[2]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
65 # Significant difference between samples and pools
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
66 if (rttest$p.value < 0.01) { ok=3
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
67 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
68 # to avoid denominator =0 when on 1 pool per batch <> 0
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
69 if (method=="linear" & length(which(((penteB*qci)+ordori)==0))>0 ){ ok=6
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
70 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
71 # different sloop between samples and pools
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
72 if (method=="linear" & penteB/penteS < -0.20) { ok=5
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
73 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
74 #
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
75 if (method=="linear" & !is.na(normref) & valimp=="NA") {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
76 denom = (penteB * c(spi,qci) + ordori)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
77 normval = c(spl,qcp)*normref / denom
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
78 if(length(which((normval==Inf)|(denom<1)))==length(normval)){ok=6}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
79 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
80 }}}}}}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
81 ok_norm=ok
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
82 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
83
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
84 plotsituation <- function (x, nbid,outfic="plot_regression.pdf", outres="PreNormSummary.txt",fact="batch",span="none",
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
85 sm_meta=list(batch="batch", injectionOrder="injectionOrder", sampleType="sampleType",
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
86 sampleTag=list(pool="pool",blank="blank",sample="sample"))) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
87 # Checks for all ions in every batch if linear or lo(w)ess correction is possible.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
88 # Uses ok_norm function and creates a file (PreNormSummary.txt) with the corresponding error codes.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
89 # Also creates a pdf file with plots of linear and lo(w)ess regression lines.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
90 # Parameters:
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
91 # x: dataframe with ions in columns and samples in rows ; x is the result of concatenation of sample metadata file and ions file
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
92 # nbid: number of samples description columns (id and factors) with at least : "batch","injectionOrder","sampleType"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
93 # outfic: name of regression plots pdf file
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
94 # outres: name of summary table file
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
95 # fact: factor to be used as categorical variable for plots and PCA
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
96 # span: span value for lo(w)ess regression; "none" for linear or default values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
97 # sm_meta: list of information about sample metadata coding
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
98 indfact=which(dimnames(x)[[2]]==fact)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
99 indtypsamp=which(dimnames(x)[[2]]==sm_meta$sampleType)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
100 indbatch=which(dimnames(x)[[2]]==sm_meta$batch)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
101 indinject=which(dimnames(x)[[2]]==sm_meta$injectionOrder)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
102 lastIon=dim(x)[2]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
103 nbi=lastIon-nbid # Number of ions = total number of columns - number of identifying columns
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
104 nbb=length(levels(x[[sm_meta$batch]])) # Number of batch = number of levels of "batch" comlumn (factor)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
105 nbs=length(x[[sm_meta$sampleType]][x[[sm_meta$sampleType]] %in% sm_meta$sampleTag$sample])# Number of samples = number of rows with "sample" value in sampleType
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
106 pdf(outfic,width=27,height=7*ceiling((nbb+2)/3))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
107 cat(nbi," ions ",nbb," batch(es) \n")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
108 cv=data.frame(matrix(0,nrow=nbi,ncol=2))# initialisation de la dataset qui contiendra les CV
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
109 pre_bilan=matrix(0,nrow=nbi,ncol=3*nbb) # dataset of ok_norm function results
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
110 for (p in 1:nbi) {# for each ion
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
111 par (mfrow=c(ceiling((nbb+2)/3),3),ask=F,cex=1.2)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
112 labion=dimnames(x)[[2]][p+nbid]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
113 indpool=which(x[[sm_meta$sampleType]] %in% sm_meta$sampleTag$pool) # QCpools subscripts in x
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
114 pools1=x[indpool,p+nbid]; cv[p,1]=sd(pools1,na.rm=TRUE)/mean(pools1,na.rm=TRUE)# CV before correction
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
115 for (b in 1:nbb) {# for each batch...
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
116 xb=data.frame(x[(x[[sm_meta$batch]]==levels(x[[sm_meta$batch]])[b]),c(indtypsamp,indinject,p+nbid)])
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
117 indpb = which(xb[[sm_meta$sampleType]] %in% sm_meta$sampleTag$pool)# QCpools subscripts in the current batch
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
118 indsp = which(xb[[sm_meta$sampleType]] %in% sm_meta$sampleTag$sample)# samples subscripts in the current batch
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
119 normLinearTest=ok_norm(xb[indpb,3],xb[indpb,2], xb[indsp,3],xb[indsp,2],method="linear",normref=mean(xb[c(indpb,indsp),3],na.rm=TRUE),valimp="NA")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
120 normLoessTest=ok_norm(xb[indpb,3],xb[indpb,2], xb[indsp,3],xb[indsp,2],method="loess")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
121 normLowessTest=ok_norm(xb[indpb,3],xb[indpb,2], xb[indsp,3],xb[indsp,2],method="lowess")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
122 pre_bilan[ p,3*b-2]=normLinearTest
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
123 pre_bilan[ p,3*b-1]=normLoessTest
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
124 pre_bilan[ p,3*b]=normLowessTest
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
125 if(length(indpb)>1){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
126 if(span=="none"){span1<-1 ; span2<-2*length(indpool)/nbs}else{span1<-span ; span2<-span}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
127 if(normLoessTest!=2){resloess=loess(xb[indpb,3]~xb[indpb,2],span=span1,degree=2,family="gaussian",iterations=4,surface="direct")}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
128 if(length(which(!(is.na(xb[indsp,3]))))>1){resloessSample=loess(xb[indsp,3]~xb[indsp,2],span=2*length(indpool)/nbs,degree=2,family="gaussian",iterations=4,surface="direct") }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
129 if(normLowessTest!=2){reslowess=lowess(xb[indpb,2],xb[indpb,3],f=span2)}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
130 if(length(which(!(is.na(xb[indsp,3]))))>1){reslowessSample=lowess(xb[indsp,2],xb[indsp,3])}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
131 liminf=min(xb[,3],na.rm=TRUE);limsup=max(xb[,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
132 firstinj=min(xb[,2],na.rm=TRUE);lastinj=max(xb[,2],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
133 plot(xb[indsp,2],xb[indsp,3],pch=16, main=paste(labion,"batch ",b),ylab="intensity",xlab="injection order",ylim=c(liminf,limsup),xlim=c(firstinj,lastinj))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
134 if(nrow(xb)>(length(indpb)+length(indsp))){points(xb[-c(indpb,indsp),2], xb[-c(indpb,indsp),3],pch=18,col="grey")}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
135 points(xb[indpb,2], xb[indpb,3],pch=5)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
136 if(normLoessTest!=2){points(cbind(resloess$x,resloess$fitted)[order(resloess$x),],type="l",col="green3")}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
137 if(length(which(!(is.na(xb[indsp,3]))))>1){points(cbind(resloessSample$x,resloessSample$fitted)[order(resloessSample$x),],type="l",col="green3",lty=2)}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
138 if(normLowessTest!=2){points(reslowess,type="l",col="red")}; if(length(which(!(is.na(xb[indsp,3]))))>1){points(reslowessSample,type="l",col="red",lty=2)}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
139 abline(lsfit(xb[indpb,2],xb[indpb,3]),col="blue")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
140 if(length(which(!(is.na(xb[indsp,3]))))>1){abline(lsfit(xb[indsp,2],xb[indsp,3]),lty=2,col="blue")}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
141 legend("topleft",c("pools","samples"),lty=c(1,2),bty="n")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
142 legend("topright",c("linear","lowess","loess"),lty=1,col=c("blue","red","green3"),bty="n")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
143 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
144 plot.new()
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
145 legend("center","Plot only available when the\nbatch contains at least 2 pools.")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
146 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
147 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
148 # series de plot avant correction
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
149 minval=min(x[p+nbid],na.rm=TRUE);maxval=max(x[p+nbid],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
150 plot( x[[sm_meta$injectionOrder]], x[,p+nbid],col=x[[sm_meta$batch]],ylim=c(minval,maxval),ylab=labion,
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
151 main=paste0("before correction (CV for pools = ",round(cv[p,1],2),")"),xlab="injection order")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
152 suppressWarnings(plot.design( x[c(indtypsamp,indbatch,indfact,p+nbid)],main="factors effect before correction"))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
153 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
154 dev.off()
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
155 pre_bilan=data.frame(pre_bilan)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
156 labion=dimnames(x)[[2]][nbid+1:nbi]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
157 for (i in 1:nbb) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
158 dimnames(pre_bilan)[[2]][3*i-2]=paste("batch",i,"linear")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
159 dimnames(pre_bilan)[[2]][3*i-1]=paste("batch",i,"loess")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
160 dimnames(pre_bilan)[[2]][3*i]=paste("batch",i,"lowess")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
161 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
162 bilan=data.frame(labion,pre_bilan)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
163 write.table(bilan,file=outres,sep="\t",row.names=F,quote=F)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
164 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
165
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
166
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
167 normlowess=function (xb,detail="no",vref=1,b,span=NULL,valneg="none",sm_meta=list(batch="batch", injectionOrder="injectionOrder", sampleType="sampleType",
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
168 sampleTag=list(pool="pool",blank="blank",sample="sample")),min_norm=1){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
169 # Correction function applied to 1 ion in 1 batch.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
170 # Uses a lowess regression computed on QC-pools in order to correct samples intensity values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
171 # xb: dataframe for 1 ion in columns and samples in rows.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
172 # vref: reference value (average of ion)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
173 # b: batch subscript
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
174 # detail: level of detail in the outlog file
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
175 # span: span value for lo(w)ess regression; NULL for default values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
176 # valneg: to determine what to do with generated negative and Inf values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
177 # sm_meta: list of information about sample metadata coding
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
178 # min_norm: minimum value accepted for normalisation term (denominator); should be strictly positive
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
179 indpb = which(xb[[sm_meta$sampleType]] %in% sm_meta$sampleTag$pool) # pools subscripts of current batch
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
180 indsp = which(xb[[sm_meta$sampleType]] %in% sm_meta$sampleTag$sample) # samples of current batch subscripts
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
181 labion=dimnames(xb)[[2]][3]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
182 newval=xb[[3]] # initialisation of corrected values = intial values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
183 ind <- 0 # initialisation of correction indicator
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
184 normTodo=ok_norm(xb[indpb,3],xb[indpb,2], xb[indsp,3],xb[indsp,2],method="lowess")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
185 #cat("batch:",b," dim xb=",dim(xb)," ok=",normTodo,"\n")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
186 if (normTodo==0) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
187 if(length(span)==0){span2<-2*length(indpb)/length(indsp)}else{span2<-span}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
188 reslowess=lowess(xb[indpb,2],xb[indpb,3],f=span2) # lowess regression with QC-pools
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
189 if(length(which(reslowess$y<min_norm))!=0){ # to handle cases where 0<denominator<min_norm or negative
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
190 toajust <- which(reslowess$y<min_norm)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
191 if(valneg=="NA"){ reslowess$y[toajust] <- NA
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
192 } else { if(valneg=="0"){ reslowess$y[toajust] <- -1
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
193 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
194 mindenom <- min(reslowess$y[reslowess$y>=min_norm],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
195 reslowess$y[toajust] <- mindenom
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
196 } } }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
197 for(j in 1:nrow(xb)) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
198 if (j %in% indpb) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
199 newval[j]=(vref*xb[j,3]) / (reslowess$y[which(indpb==j)])
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
200 } else { # for samples other than pools, the correction value "corv" correspond to the nearest QCpools
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
201 corv= reslowess$y[which(abs(reslowess$x-xb[j,2])==min(abs(reslowess$x-xb[j,2]),na.rm=TRUE))]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
202 if (length(corv)>1) {corv=corv[1]}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
203 newval[j]=(vref*xb[j,3]) / corv
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
204 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
205 if((!is.na(newval[j]))&(newval[j]<0)){newval[j]<-0}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
206 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
207 if (detail=="reg") {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
208 liminf=min(xb[,3],na.rm=TRUE);limsup=max(xb[,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
209 firstinj=min(xb[,2],na.rm=TRUE);lastinj=max(xb[,2],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
210 plot(xb[indsp,2],xb[indsp,3],pch=16,main=paste(labion,"batch ",b),ylab="intensity",xlab="injection order",ylim=c(liminf,limsup),xlim=c(firstinj,lastinj))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
211 if(nrow(xb)>(length(indpb)+length(indsp))){points(xb[-c(indpb,indsp),2], xb[-c(indpb,indsp),3],pch=18)}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
212 points(xb[indpb,2], xb[indpb,3],pch=5)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
213 points(reslowess,type="l",col="red")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
214 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
215 ind <- 1
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
216 } else {# if ok_norm != 0 , we perform a correction based on batch pool or sample average
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
217 if((length(which(!is.na(xb[indpb,3])))>0)&(length(which(xb[indpb,3]>0))>0)){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
218 moypool=mean(xb[indpb,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
219 newval = (vref*xb[,3])/moypool
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
220 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
221 moysamp=mean(xb[indsp,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
222 if((!is.na(moysamp))&(moysamp>0)){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
223 cat("Warning: no pool value >0 detected in batch",b,"of ion",labion,": sample mean used as normalisation term.\n")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
224 newval = (vref*xb[,3])/moysamp
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
225 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
226 dev.off()
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
227 stop(paste("\n- - - -\nNo pool nor sample value >0 in batch",b,"of ion",labion,"- correction process aborted.\n- - - -\n"))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
228 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
229 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
230 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
231 newval <- list(norm.ion=newval,norm.ind=ind)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
232 return(newval)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
233 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
234
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
235 normlinear <- function (xb,detail="no",vref=1,b,valneg="none",sm_meta=list(batch="batch", injectionOrder="injectionOrder", sampleType="sampleType",
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
236 sampleTag=list(pool="pool",blank="blank",sample="sample")),min_norm=1){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
237 # Correction function applied to 1 ion in 1 batch.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
238 # Uses a linear regression computed on QC-pools in order to correct samples intensity values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
239 # xb: dataframe with ions in columns and samples in rows; x is a result of concatenation of sample metadata file and ion file
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
240 # detail: level of detail in the outlog file
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
241 # vref: reference value (average of ion)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
242 # b: which batch it is
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
243 # valneg: to determine what to do with generated negative and Inf values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
244 # sm_meta: list of information about sample metadata coding
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
245 # min_norm: minimum value accepted for normalisation term (denominator); should be strictly positive
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
246 indpb = which(xb[[sm_meta$sampleType]] %in% sm_meta$sampleTag$pool)# pools subscripts of current batch
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
247 indsp = which(xb[[sm_meta$sampleType]] %in% sm_meta$sampleTag$sample)# samples of current batch subscripts
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
248 labion=dimnames(xb)[[2]][3]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
249 newval=xb[[3]] # initialisation of corrected values = intial values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
250 ind <- 0 # initialisation of correction indicator
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
251 normTodo=ok_norm(xb[indpb,3],xb[indpb,2], xb[indsp,3],xb[indsp,2],method="linear",normref=vref,valimp=valneg)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
252 if (normTodo==0) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
253 ind <- 1
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
254 reslsfit=lsfit(xb[indpb,2],xb[indpb,3]) # linear regression for QCpools
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
255 reslsfitSample=lsfit(xb[indsp,2],xb[indsp,3]) # linear regression for samples
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
256 ordori=reslsfit$coefficients[1]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
257 pente=reslsfit$coefficients[2]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
258 if (detail=="reg") {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
259 liminf=min(xb[,3],na.rm=TRUE);limsup=max(xb[,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
260 firstinj=min(xb[,2],na.rm=TRUE);lastinj=max(xb[,2],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
261 plot(xb[indsp,2],xb[indsp,3],pch=16,
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
262 main=paste(labion,"batch ",b),ylab="intensity",xlab="injection order",ylim=c(liminf,limsup),xlim=c(firstinj,lastinj))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
263 if(nrow(xb)>(length(indpb)+length(indsp))){points(xb[-c(indpb,indsp),2], xb[-c(indpb,indsp),3],pch=18)}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
264 points(xb[indpb,2], xb[indpb,3],pch=5)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
265 abline(reslsfit)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
266 abline(reslsfitSample,lty=2)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
267 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
268 # correction with rescaling of ion global intensity (vref)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
269 newval = (vref*xb[,3]) / (pente * (xb[,2]) + ordori)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
270 newval[which((pente * (xb[,2]) + ordori)<min_norm)] <- -1 # to handle cases where 0<denominator<1 or negative
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
271 # handling if any negative values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
272 if(length(which((newval==Inf)|(newval<0)))!=0){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
273 toajust <- which((newval==Inf)|(newval<0))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
274 if(valneg=="NA"){ newval[toajust] <- NA
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
275 } else { if(valneg=="0"){ newval[toajust] <- 0
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
276 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
277 mindenom <- (pente * (xb[,2]) + ordori)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
278 mindenom <- min(mindenom[mindenom>=min_norm],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
279 newval[toajust] <- vref * (xb[,3][toajust]) / mindenom
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
280 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
281 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
282 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
283 } else {# if ok_norm != 0 , we perform a correction based on batch pool or sample average
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
284 if((length(which(!is.na(xb[indpb,3])))>0)&(length(which(xb[indpb,3]>0))>0)){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
285 moypool=mean(xb[indpb,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
286 newval = (vref*xb[,3])/moypool
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
287 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
288 moysamp=mean(xb[indsp,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
289 if((!is.na(moysamp))&(moysamp>0)){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
290 cat("Warning: no pool value >0 detected in batch",b,"of ion",labion,": sample mean used as normalisation term.\n")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
291 newval = (vref*xb[,3])/moysamp
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
292 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
293 dev.off()
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
294 stop(paste("\n- - - -\nNo pool nor sample value >0 in batch",b,"of ion",labion,"- correction process aborted.\n- - - -\n"))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
295 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
296 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
297 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
298 newval <- list(norm.ion=newval,norm.ind=ind)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
299 return(newval)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
300 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
301
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
302
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
303 normloess <- function (xb,detail="no",vref=1,b,span=NULL,valneg="none",sm_meta=list(batch="batch", injectionOrder="injectionOrder", sampleType="sampleType",
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
304 sampleTag=list(pool="pool",blank="blank",sample="sample")),min_norm=1){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
305 # Correction function applied to 1 ion in 1 batch.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
306 # Uses a loess regression computed on QC-pools in order to correct samples intensity values.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
307 # xb: dataframe for 1 ion in columns and samples in rows.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
308 # detail: level of detail in the outlog file.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
309 # vref: reference value (average of ion)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
310 # b: batch subscript
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
311 # span: span value for lo(w)ess regression; NULL for default values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
312 # valneg: to determine what to do with generated negative and Inf values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
313 # sm_meta: list of information about sample metadata coding
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
314 # min_norm: minimum value accepted for normalisation term (denominator); should be strictly positive
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
315 indpb = which(xb[[sm_meta$sampleType]] %in% sm_meta$sampleTag$pool) # pools subscripts of current batch
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
316 indsp = which(xb[[sm_meta$sampleType]] %in% sm_meta$sampleTag$sample) # samples of current batch subscripts
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
317 indbt = which(xb[[sm_meta$sampleType]] %in% c(sm_meta$sampleTag$sample,sm_meta$sampleTag$pool))# batch subscripts of samples and QCpools
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
318 labion=dimnames(xb)[[2]][3]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
319 newval=xb[[3]] # initialisation of corrected values = intial values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
320 ind <- 0 # initialisation of correction indicator
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
321 normTodo=ok_norm(xb[indpb,3],xb[indpb,2], xb[indsp,3],xb[indsp,2],method="loess")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
322 if (normTodo==0) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
323 if(length(span)==0){span1<-1}else{span1<-span}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
324 resloess=loess(xb[indpb,3]~xb[indpb,2],span=span1,degree=2,family="gaussian",iterations=4,surface="direct") # loess regression with QCpools
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
325 corv=predict(resloess,newdata=xb[,2])
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
326 if(length(which(corv<min_norm))!=0){ # unconsistant values handling
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
327 toajust <- which(corv<min_norm)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
328 if(valneg=="NA"){ corv[toajust] <- NA
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
329 } else { if(valneg=="0"){ corv[toajust] <- -1
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
330 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
331 mindenom <- min(corv[corv>=min_norm],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
332 corv[toajust] <- mindenom
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
333 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
334 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
335 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
336 newvalps=(vref*xb[indbt,3]) / corv[indbt] # to check if correction generates outlier values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
337 refthresh=max(c(3*(quantile(newvalps,na.rm=TRUE)[4]),1.3*(xb[indbt,3])),na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
338 if(length(which(newvalps>refthresh))>0){ # if outliers
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
339 # in this case no modification of initial value
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
340 newval <- xb[,3]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
341 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
342 newval=(vref*xb[,3]) / corv
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
343 newval[newval<0] <- 0
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
344 ind <- 1 # confirmation of correction
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
345 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
346 if ((detail=="reg")&(ind==1)) { # plot
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
347 liminf=min(xb[,3],na.rm=TRUE);limsup=max(xb[,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
348 firstinj=min(xb[,2],na.rm=TRUE);lastinj=max(xb[,2],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
349 plot(xb[indsp,2],xb[indsp,3],pch=16,main=paste(labion,"batch ",b),ylab="intensity",xlab="injection order",ylim=c(liminf,limsup),xlim=c(firstinj,lastinj))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
350 if(nrow(xb)>(length(indpb)+length(indsp))){points(xb[-c(indpb,indsp),2], xb[-c(indpb,indsp),3],pch=18)}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
351 points(xb[indpb,2], xb[indpb,3],pch=5)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
352 points(cbind(resloess$x,resloess$fitted)[order(resloess$x),],type="l",col="red")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
353 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
354 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
355 if (ind==0) {# if ok_norm != 0 or if correction creates outliers, we perform a correction based on batch pool or sample average
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
356 if((length(which(!is.na(xb[indpb,3])))>0)&(length(which(xb[indpb,3]>0))>0)){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
357 moypool=mean(xb[indpb,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
358 newval = (vref*xb[,3])/moypool
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
359 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
360 moysamp=mean(xb[indsp,3],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
361 if((!is.na(moysamp))&(moysamp>0)){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
362 cat("Warning: no pool value >0 detected in batch",b,"of ion",labion,": sample mean used as normalisation term.\n")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
363 newval = (vref*xb[,3])/moysamp
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
364 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
365 dev.off()
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
366 stop(paste("\n- - - -\nNo pool nor sample value >0 in batch",b,"of ion",labion,"- correction process aborted.\n- - - -\n"))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
367 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
368 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
369 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
370 newval <- list(norm.ion=newval,norm.ind=ind)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
371 return(newval)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
372 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
373
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
374
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
375
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
376 norm_QCpool <- function (x, nbid, outlog, fact, metaion, detail="no", NormMoyPool=FALSE, NormInt=FALSE, method="linear",span="none",valNull="0",
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
377 sm_meta=list(batch="batch", injectionOrder="injectionOrder", sampleType="sampleType",
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
378 sampleTag=list(pool="pool",blank="blank",sample="sample")),min_norm=1) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
379 ### Correction applying linear or lo(w)ess correction function on all ions for every batch of a dataframe.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
380 # x: dataframe with ions in column and samples' metadata
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
381 # nbid: number of sample description columns (id and factors) with at least "batch", "injectionOrder", "sampleType"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
382 # outlog: name of regression plots and PCA pdf file
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
383 # fact: factor to be used as categorical variable for plots
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
384 # metaion: dataframe of ions' metadata
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
385 # detail: level of detail in the outlog file. detail="no" ACP + boxplot of CV before and after correction.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
386 # detail="plot" with plot for all batch before and after correction.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
387 # detail="reg" with added plots with regression lines for all batches.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
388 # NormMoyPool: not used
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
389 # NormInt: not used
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
390 # method: regression method to be used to correct : "linear" or "lowess" or "loess"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
391 # span: span value for lo(w)ess regression; "none" for linear or default values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
392 # valNull: to determine what to do with negatively estimated intensities
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
393 # sm_meta: list of information about sample metadata coding
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
394 # min_norm: minimum value accepted for normalisation term (denominator); should be strictly positive
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
395 indfact=which(dimnames(x)[[2]]==fact)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
396 indtypsamp=which(dimnames(x)[[2]]==sm_meta$sampleType)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
397 indbatch=which(dimnames(x)[[2]]==sm_meta$batch)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
398 indinject=which(dimnames(x)[[2]]==sm_meta$injectionOrder)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
399 lastIon=dim(x)[2]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
400 indpool=which(x[[sm_meta$sampleType]] %in% sm_meta$sampleTag$pool)# QCpools subscripts in all batches
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
401 valref=apply(as.matrix(x[indpool,(nbid+1):(lastIon)]),2,mean,na.rm=TRUE) # reference value for each ion used to still have the same rought size of values
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
402 nbi=lastIon-nbid # number of ions
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
403 nbb=length(levels(x[[sm_meta$batch]])) # Number of batch(es) = number of levels of factor "batch" (can be =1)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
404 Xn=data.frame(x[,c(1:nbid)],matrix(0,nrow=nrow(x),ncol=nbi))# initialisation of the corrected dataframe (=initial dataframe)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
405 dimnames(Xn)=dimnames(x)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
406 cv=data.frame(matrix(NA,nrow=nbi,ncol=2))# initialisation of dataframe containing CV before and after correction
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
407 dimnames(cv)[[2]]=c("avant","apres")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
408 if (detail!="reg" && detail!="plot" && detail!="no") {detail="no"}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
409 pdf(outlog,width=27,height=20)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
410 cat(nbi," ions ",nbb," batch(es) \n")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
411 if (detail=="plot") {if(nbb<6){par(mfrow=c(3,3),ask=F,cex=1.5)}else{par(mfrow=c(4,4),ask=F,cex=1.5)}}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
412 res.ind <- matrix(NA,ncol=nbb,nrow=nbi,dimnames=list(dimnames(x)[[2]][-c(1:nbid)],paste("norm.b",1:nbb,sep="")))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
413 for (p in 1:nbi) {# for each ion
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
414 labion=dimnames(x)[[2]][p+nbid]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
415 pools1=x[indpool,p+nbid]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
416 if(length(which(pools1[!(is.na(pools1))]>0))<2){ # if not enough pools >0 -> no normalisation
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
417 war.note <- paste("Warning: less than 2 pools with values >0 in",labion,"-> no normalisation for this ion.")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
418 cat(war.note,"\n")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
419 Xn[,p+nbid] <- x[,p+nbid]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
420 res.ind[p,] <- rep(0,nbb)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
421 if (detail=="reg" || detail=="plot" ) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
422 par(mfrow=c(2,2),ask=F,cex=1.5)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
423 plot.new()
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
424 legend("center",war.note)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
425 minval=min(x[p+nbid],na.rm=TRUE);maxval=max(x[p+nbid],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
426 plot( x[[sm_meta$injectionOrder]], x[,p+nbid],col=x[[sm_meta$batch]],ylab=labion,ylim=c(minval,maxval),
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
427 main="No correction",xlab="injection order")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
428 points(x[[sm_meta$injectionOrder]][indpool],x[indpool,p+nbid],col="maroon",pch=16,cex=1)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
429 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
430 } else {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
431 if (detail == "reg") {if(nbb<6){par(mfrow=c(3,3),ask=F,cex=1.5)}else{par(mfrow=c(4,4),ask=F,cex=1.5)}}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
432 if (detail == "plot") {par(mfrow=c(2,2),ask=F,cex=1.5)}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
433 cv[p,1]=sd(pools1,na.rm=TRUE)/mean(pools1,na.rm=TRUE)# CV before correction
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
434 for (b in 1:nbb) {# for every batch
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
435 indbt = which(x[[sm_meta$batch]]==(levels(x[[sm_meta$batch]])[b])) # subscripts of all samples
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
436 sub=data.frame(x[(x[[sm_meta$batch]]==levels(x[[sm_meta$batch]])[b]),c(indtypsamp,indinject,p+nbid)])
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
437 if (method=="linear") { res.norm = normlinear(sub,detail,valref[p],b,valNull,sm_meta,min_norm)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
438 } else { if (method=="loess"){ res.norm <- normloess(sub,detail,valref[p],b,span,valNull,sm_meta,min_norm)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
439 } else { if (method=="lowess"){ res.norm <- normlowess(sub,detail,valref[p],b,span,valNull,sm_meta,min_norm)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
440 } else {stop("\n--\nNo valid 'method' argument supplied.\nMust be 'linear','loess' or 'lowess'.\n--\n")}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
441 }}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
442 Xn[indbt,p+nbid] = res.norm[[1]]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
443 res.ind[p,b] <- res.norm[[2]]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
444 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
445 # Post correction CV calculation
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
446 pools2=Xn[indpool,p+nbid]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
447 cv[p,2]=sd(pools2,na.rm=TRUE)/mean(pools2,na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
448 if (detail=="reg" || detail=="plot" ) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
449 # plot before and after correction
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
450 minval=min(cbind(x[p+nbid],Xn[p+nbid]),na.rm=TRUE);maxval=max(cbind(x[p+nbid],Xn[p+nbid]),na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
451 plot( x[[sm_meta$injectionOrder]], x[,p+nbid],col=x[[sm_meta$batch]],ylab=labion,ylim=c(minval,maxval),
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
452 main=paste0("before correction (CV for pools = ",round(cv[p,1],2),")"),xlab="injection order")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
453 points(x[[sm_meta$injectionOrder]][indpool],x[indpool,p+nbid],col="maroon",pch=16,cex=1)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
454 plot(Xn[[sm_meta$injectionOrder]],Xn[,p+nbid],col=x[[sm_meta$batch]],ylab="",ylim=c(minval,maxval),
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
455 main=paste0("after correction (CV for pools = ",round(cv[p,2],2),")"),xlab="injection order")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
456 points(Xn[[sm_meta$injectionOrder]][indpool],Xn[indpool,p+nbid],col="maroon",pch=16,cex=1)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
457 suppressWarnings(plot.design( x[c(indtypsamp,indbatch,indfact,p+nbid)],main="factors effect before correction"))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
458 suppressWarnings(plot.design(Xn[c(indtypsamp,indbatch,indfact,p+nbid)],main="factors effect after correction"))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
459 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
460 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
461 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
462
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
463 if (detail=="reg" || detail=="plot" || detail=="no") {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
464 if (nbi > 3) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
465 # Sum of ions before/after plot
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
466 par(mfrow=c(1,2),ask=F,cex=1.2)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
467 xsum <- rowSums(x[,(nbid+1):lastIon],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
468 Xnsum <- rowSums(Xn[,(nbid+1):lastIon],na.rm=TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
469 plot(x[[sm_meta$injectionOrder]],xsum,col=x[[sm_meta$batch]],ylab="sum of variables' intensities",xlab="injection order",
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
470 ylim=c(min(c(xsum,Xnsum),na.rm=TRUE),max(c(xsum,Xnsum),na.rm=TRUE)),main="Sum of intensities\nBefore correction")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
471 points(x[[sm_meta$injectionOrder]][indpool],xsum[indpool],col="maroon",pch=16,cex=1.2)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
472 plot(x[[sm_meta$injectionOrder]],Xnsum,col=x[[sm_meta$batch]],ylab="sum of variables' intensities",xlab="injection order",
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
473 ylim=c(min(c(xsum,Xnsum),na.rm=TRUE),max(c(xsum,Xnsum),na.rm=TRUE)),main="Sum of intensities\nAfter correction")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
474 points(x[[sm_meta$injectionOrder]][indpool],Xnsum[indpool],col="maroon",pch=16,cex=1.2)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
475 # PCA Plot before/after, normed only and ions plot
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
476 par(mfrow=c(3,4),ask=F,cex=1.2)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
477 acplight(x[,c(indtypsamp,indbatch,indtypsamp,indfact,(nbid+1):lastIon)],"uv",TRUE)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
478 norm.ion <- which(colnames(Xn)%in%(rownames(res.ind)[which(rowSums(res.ind)>=1)]))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
479 acplight(Xn[,c(indtypsamp,indbatch,indtypsamp,indfact,(nbid+1):lastIon)],"uv",TRUE,norm.ion)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
480 if(length(norm.ion)>0){acplight(Xn[,c(indtypsamp,indbatch,indtypsamp,indfact,norm.ion)],"uv",TRUE)}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
481 # Before/after boxplot
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
482 par(mfrow=c(1,2),ask=F,cex=1.2)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
483 cvplot=cv[!is.na(cv[[1]])&!is.na(cv[[2]]),]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
484 if(nrow(cvplot)>0){
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
485 boxplot(cvplot[[1]],ylim=c(min(cvplot),max(cvplot)),main="CV of pools before correction")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
486 boxplot(cvplot[[2]],ylim=c(min(cvplot),max(cvplot)),main="CV of pools after correction")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
487 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
488 dev.off()
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
489 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
490 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
491 if (nbi<=3) {dev.off()}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
492 # transposed matrix is return (format of the initial matrix with ions in rows)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
493 Xr=Xn[,-c(1:nbid)]; dimnames(Xr)[[1]]=Xn[[1]]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
494 Xr=t(Xr) ; Xr <- data.frame(ions=rownames(Xr),Xr)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
495
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
496 res.norm[[1]] <- Xr ; res.norm[[2]] <- data.frame(metaion,res.ind) ; res.norm[[3]] <- x[,c(1:nbid)]
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
497 names(res.norm) <- c("dataMatrix","variableMetadata","sampleMetadata")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
498 return(res.norm)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
499 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
500
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
501
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
502
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
503
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
504
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
505 acplight <- function(ids, scaling="uv", indiv=FALSE,indcol=NULL) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
506 suppressPackageStartupMessages(library(ade4))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
507 suppressPackageStartupMessages(library(pcaMethods))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
508 # Make a PCA and plot scores and loadings.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
509 # First column must contain samples' identifiers.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
510 # Columns 2 to 4 contain factors to colour the plots.
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
511 for (i in 1:3) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
512 idss <- data.frame(ids)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
513 idss[,i+1] <- as.character(idss[,i+1])
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
514 idss[which(is.na(idss[,i+1])),i+1] <- "no_modality"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
515 idss[which(idss[,i+1]=="NA"),i+1] <- "no_modality"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
516 idss[which(idss[,i+1]==""),i+1] <- "no_modality"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
517 classe=as.factor(idss[[i+1]])
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
518 idsample=as.character(idss[[1]])
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
519 colour=1:length(levels(classe))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
520 ions=as.matrix(idss[,5:dim(idss)[2]])
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
521 # Removing ions containing NA (not compatible with standard PCA)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
522 ions=t(na.omit(t(ions)))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
523 if(i==1){if(ncol(ions)!=(ncol(idss)-4)){cat("Note:",(ncol(idss)-4)-ncol(ions),"ions were ignored for PCA display due to NA in intensities.\n")}}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
524 # Scaling choice: "uv","none","pareto"
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
525 object=suppressWarnings(prep(ions, scale=scaling, center=TRUE))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
526 if(i==1){if(length(which(apply(ions,2,var)==0))>0){cat("Warning: there are",length(which(apply(ions,2,var)==0)),"constant ions.\n")}}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
527 # ALGO: nipals,svdImpute, Bayesian, svd, probalistic=F
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
528 result <- pca(object, center=F, method="svd", nPcs=2)
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
529 # ADE4 : to plot samples' ellipsoid for each class
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
530 s.class(result@scores, classe, cpoint = 1,xax=1,yax=2,col=colour,sub=sprintf("Scores - PCs %sx%s",1,2), possub="bottomright")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
531 #s.label(result@loadings,label = ions, cpoint = 0, clabel=0.4, xax=1,yax=2,sub="Loadings",possub="bottomright")
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
532 if(i==1){resulti <- result}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
533 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
534 if(indiv) {
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
535 colour <- rep("darkblue",length(resulti@loadings)) ; if(!is.null(indcol)) {colour[-c(indcol)] <- "red"}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
536 plot(resulti@loadings,col=colour,main="Loadings",xaxt="n",yaxt="n",pch=20,
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
537 xlab=bquote(PC1-R^2==.(resulti@R2[1])),ylab=bquote(PC2 - R^2 == .(resulti@R2[2])))
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
538 abline(h=0,v=0)}
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
539 }
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
540
23314e1192d4 Uploaded
melpetera
parents:
diff changeset
541