-
Notifications
You must be signed in to change notification settings - Fork 1
/
M2S_calculateResiduals.m
354 lines (303 loc) · 19.5 KB
/
M2S_calculateResiduals.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
% [Residuals_X,Residuals_trendline] = M2S_calculateResiduals(refSet,targetSet,Xr_connIdx,Xt_connIdx, nrNeighbors, neighMethod, plotOrNot)
% This function is part of M2S toolbox to match metabolomics features across untargeted datasets.
%
% It calculates the difference (residuals) between each match and the
% inter-dataset shift of two datasets, in each dimension.
%
% INPUT:
% refSet,targetSet: nx3 matrices of MZ,RT,FI that have been matched
% (non-unique entries, containing clusters of multiple matches).
% Xr_connIdx,Xt_connIdx: single-column indices of features in the initial
% datasets that have been matched to each other.
% nrNeighbors: a percentage of the total of features in the refSet (e.g. 0.05) or a specific number (e.g. 25).
% neighMethod:method to calculate inter-dataset shifts can be "cross" (uses RT, MZ, FI)
% or (default) "circle" (uses only RT, MZ). "Circle" does not subtract the shift from FI.
% It has been updated 12Oct2023 to use a robust loess to calculate the
% inter-dataset shift trend in both methods (though without FI in "circle")
% Update 14Nov2023: if nrNeighbours=1 it uses the actual variable, rather
% than its closest neighbour.
% plotOrNot: 0 - no plots; 1 - default plots; 2 - extra plots
% NOTE: % The following call makes automatic choice of neighbour number:
% [Residuals_X,Residuals_trendline] = M2S_findNeighTrendsResiduals(ref_MatchSet, target_MatchSet, refSet, targetSet, Xr_connIdx_inMatchSet, Xr_connIdx)
%
% OUTPUT:
% Residuals_X: the distances between each match and the inter-dataset shift
% Residuals_trendline: the inter-dataset shift at each reference
%
% Rui Climaco Pinto, 2021
% Imperial College London
function [Residuals_X,Residuals_trendline] = M2S_calculateResiduals_v2 ...
(refSet,targetSet,Xr_connIdx,Xt_connIdx, nrNeighbors, neighMethod, pctPointsLoess, plotOrNot)
% Run the function to define the ref and target sets from which to choose
% neighbours. These only contain single matches, not multiple ones (in clusters).
if plotOrNot == 2
[ref_MatchSet,target_MatchSet,Xr_connIdx_inMatchSet,Xt_connIdx_inMatchSet] = ...
M2S_createMatchSets(refSet,targetSet,Xr_connIdx,Xt_connIdx,1);
else
[ref_MatchSet,target_MatchSet,Xr_connIdx_inMatchSet,Xt_connIdx_inMatchSet] = ...
M2S_createMatchSets(refSet,targetSet,Xr_connIdx,Xt_connIdx,0);
end
if nargin == 4
nrNeighbors = 10 + ceil(0.01*size(ref_MatchSet,1));% at least 10 neighbours plus 1% of the size of the reference set
neighMethod = 'circle';
pctPointsLoess = 0;
plotOrNot=1;
elseif nargin == 5
neighMethod = 'circle';
pctPointsLoess = 0;
plotOrNot=1;
elseif nargin == 6
plotOrNot=1;
pctPointsLoess = 0;
elseif nargin == 7
plotOrNot=1;
end
if nrNeighbors<1 % it is a fraction of the size of ref_MatchSet
nrNeighbors = round(nrNeighbors*size(ref_MatchSet,1));
end
% Transform FImed TO LOG10
ref_MatchSet(:,3) = log10(ref_MatchSet(:,3));
target_MatchSet(:,3) = log10(target_MatchSet(:,3));
refSet(:,3) = log10(refSet(:,3));
targetSet(:,3) = log10(targetSet(:,3));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. FIND THE NEIGHBOURS OF EACH FEATURE AND DISTANCES TARGET to REF (TR)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CROSS-TYPE NEIGHB0URS (uses RT, MZ and FI)
if strcmp(neighMethod,'cross')
% This is a function to find neighbours distances.
% It calculates RTdist_neighbours and MZdist_neighbours with neighbours in
% the RT and in the MZ domain separately, resulting in a "cross" shape.
wb1 = waitbar(0,'Calculating distances across sets for the neighbours of each feature','Name','Info');
idx_neighbors_Ref_RT=M2S_findNeighbours(ref_MatchSet,Xr_connIdx_inMatchSet,refSet,Xr_connIdx,nrNeighbors,1);
waitbar(0.33,wb1,'Calculating distances across sets for the neighbours of each feature','Name','Info');
idx_neighbors_Ref_MZ=M2S_findNeighbours(ref_MatchSet,Xr_connIdx_inMatchSet,refSet,Xr_connIdx,nrNeighbors,2);
waitbar(0.66,wb1,'Calculating distances across sets for the neighbours of each feature','Name','Info');
idx_neighbors_Ref_FI=M2S_findNeighbours(ref_MatchSet,Xr_connIdx_inMatchSet,refSet,Xr_connIdx,nrNeighbors,3);
waitbar(0.99,wb1,'Calculating distances across sets for the neighbours of each feature','Name','Info');
% calculate median RTdist of neighborsTR (Target-Reference)
RTdist_neighborsTR=[];
MZdist_neighborsTR=[];
FIdist_neighborsTR=[];
for columnNr = 1:nrNeighbors
RTdist_neighborsTR(:,columnNr) = target_MatchSet(idx_neighbors_Ref_RT(:,columnNr),1) - ref_MatchSet(idx_neighbors_Ref_RT(:,columnNr),1);
MZdist_neighborsTR(:,columnNr) = target_MatchSet(idx_neighbors_Ref_MZ(:,columnNr),2) - ref_MatchSet(idx_neighbors_Ref_MZ(:,columnNr),2);
FIdist_neighborsTR(:,columnNr) = target_MatchSet(idx_neighbors_Ref_FI(:,columnNr),3) - ref_MatchSet(idx_neighbors_Ref_FI(:,columnNr),3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% THIS PART OVERRIDES THE PREVIOUS CALCULATION
%% If nrNeighbors == 1, we use the actual feature, no neighbours
if nrNeighbors == 1
RTdist_neighborsTR = targetSet(:,1) - refSet(:,1);
MZdist_neighborsTR = targetSet(:,2) - refSet(:,2);
FIdist_neighborsTR = targetSet(:,3) - refSet(:,3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
waitbar(1,wb1,'Done!','Name','Info');
pause(1)
close(wb1);
%% CIRCLE-TYPE NEIGHB0URS (uses only RT and MZ)
elseif strcmp(neighMethod,'circle')
% This is a function to find neighbours distances.
% It calculates RTdis_neighbours and MZdist_neighbours with neighbours in
% the RT and in the MZ domain together, resulting in a "circle" shape.
% FI does not enter the calculation.
wb1 = waitbar(0,'Calculating distances across sets for the neighbours of each feature','Name','Info');
idx_neighbors_Ref=M2S_findNeighbours(ref_MatchSet,Xr_connIdx_inMatchSet,refSet,Xr_connIdx,nrNeighbors,[1,2]);% only neighbours in RT + MZ
% calculate median RTdist of neighborsTR (Target-Reference)
RTdist_neighborsTR = [];
MZdist_neighborsTR = [];
FIdist_neighborsTR = [];
for columnNr = 1:nrNeighbors
waitbar(columnNr/nrNeighbors,wb1,'Calculating distances across sets for the neighbours of each feature','Name','Info');
RTdist_neighborsTR(:,columnNr) = target_MatchSet(idx_neighbors_Ref(:,columnNr),1) - ref_MatchSet(idx_neighbors_Ref(:,columnNr),1);
MZdist_neighborsTR(:,columnNr) = target_MatchSet(idx_neighbors_Ref(:,columnNr),2) - ref_MatchSet(idx_neighbors_Ref(:,columnNr),2);
%FIdist_neighborsTR(:,columnNr) = target_MatchSet(idx_neighbors_Ref(:,columnNr),3) - ref_MatchSet(idx_neighbors_Ref(:,columnNr),3);
end
FIdist_neighborsTR = zeros(size(MZdist_neighborsTR));% defined as zero
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% THIS PART OVERRIDES THE PREVIOUS CALCULATION
%% If nrNeighbors == 1, we use the actual feature, no neighbours
if nrNeighbors == 1
RTdist_neighborsTR = targetSet(:,1) - refSet(:,1);
MZdist_neighborsTR = targetSet(:,2) - refSet(:,2);
FIdist_neighborsTR = zeros(size(MZdist_neighborsTR));% defined as zero
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
waitbar(1,wb1,'Done!','Name','Info');
pause(1)
close(wb1);
end
%% 2. CALCULATE THE NEIGHBOUR DISTANCES IN EACH DIMENSION
% For each feature find the median difference of distances of its
% neighbours from reference to target (SHIFTS TRENDLINE). Normalize.
% NOTICE that in "circle" method the median_FIdist_neighborsTR is zero (does not correct FI_dist_TR)
median_RTdist_neighborsTR = nanmedian(RTdist_neighborsTR,2);% REPRESENTS TRENDLINE OF RT(target-ref)
median_MZdist_neighborsTR = nanmedian(MZdist_neighborsTR,2);% REPRESENTS TRENDLINE OF MZ(target-ref)
median_FIdist_neighborsTR = nanmedian(FIdist_neighborsTR,2);% REPRESENTS TRENDLINE OF FI(target-ref)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% THIS PART OVERRIDES THE PREVIOUS CALCULATION
%% If nrNeighbors == 1, we use the actual feature, no neighbours
if nrNeighbors == 1
median_RTdist_neighborsTR = RTdist_neighborsTR;
median_MZdist_neighborsTR = MZdist_neighborsTR;
median_FIdist_neighborsTR = FIdist_neighborsTR;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if strcmp(neighMethod,'cross')
% pctPointsLoess is a percentage ]0-0.75[ of the number of points in the
% reference matched dataset (refSet), to define span of loess.
% e.g. use pctPointsLoess = 0.3
if pctPointsLoess > 0 && pctPointsLoess <0.999
%% Use robust loess to calculate trend (only for "cross" method)
% nrPointsLoess = 0.1 * length(unique(refSet(:,1)));
nrPointsLoess = round(pctPointsLoess * length(unique(refSet(:,1))));
median_RTdist_neighborsTR_v2 = smooth(refSet(:,1),refSet(:,1)+median_RTdist_neighborsTR,nrPointsLoess,'rloess');% to do loess
median_MZdist_neighborsTR_v2 = smooth(refSet(:,2),refSet(:,2)+median_MZdist_neighborsTR,nrPointsLoess,'rloess');% to do loess
if strcmp(neighMethod,'cross')% method 'circle' does not calculate FI
median_FIdist_neighborsTR_v2 = smooth(refSet(:,3),refSet(:,3)+median_FIdist_neighborsTR,nrPointsLoess,'rloess');% to do loess
else
median_FIdist_neighborsTR_v2 = refSet(:,3);
end
median_RTdist_neighborsTR = median_RTdist_neighborsTR_v2 - refSet(:,1);
median_MZdist_neighborsTR = median_MZdist_neighborsTR_v2 - refSet(:,2);
median_FIdist_neighborsTR = median_FIdist_neighborsTR_v2 - refSet(:,3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get RTdist etc of actual feature pairs
RT_dist_TR = targetSet(:,1) - refSet(:,1);% FEATURES TARGET - REF
MZ_dist_TR = targetSet(:,2) - refSet(:,2);% FEATURES TARGET - REF
FI_dist_TR = targetSet(:,3) - refSet(:,3);% FEATURES TARGET - REF
% Calculate deltaRTdist (RESIDUALS) between each median_RTdist_neighborsTR and RTdist (the deltaRT to neighbors)
RTdist_medNeigh_to_TRpair = RT_dist_TR - median_RTdist_neighborsTR; % IMPORTANT ONE
MZdist_medNeigh_to_TRpair = MZ_dist_TR - median_MZdist_neighborsTR; % IMPORTANT ONE
FIdist_medNeigh_to_TRpair = FI_dist_TR - median_FIdist_neighborsTR; % IMPORTANT ONE
Residuals_trendline = [median_RTdist_neighborsTR,median_MZdist_neighborsTR,median_FIdist_neighborsTR];
Residuals_X = [RTdist_medNeigh_to_TRpair,MZdist_medNeigh_to_TRpair,FIdist_medNeigh_to_TRpair];
% THESE ARE THE RESIDUALS USED TO CALCULATE THE PENALISATION SCORES
% NOTE: Residuals_X are the "RTMZFIdist_medNeigh_to_TRpair"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PLOTS ******************************************************
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
greyColor = [0.5 0.5 0.5];
if plotOrNot >= 1
featureNr = round(size(refSet,1)/2); % only used for the plot of neighbours of a feature
if strcmp(neighMethod,'cross')
if plotOrNot==2
%% Figure to check the neighbours of one of the features
temp_idx_neighbours_Ref_RT = (idx_neighbors_Ref_RT(featureNr,:))';
temp_idx_neighbours_Ref_MZ = (idx_neighbors_Ref_MZ(featureNr,:))';
temp_idx_neighbours_Ref_FI = (idx_neighbors_Ref_FI(featureNr,:))';
M2S_figureH(0.65,0.35);
set(gcf,'Name',['Example of neighbours for feature number ',num2str(featureNr)])
subplot(1,4,1)
plot(ref_MatchSet(:,1),ref_MatchSet(:,2),'.k'), hold on
plot(ref_MatchSet(temp_idx_neighbours_Ref_RT,1),ref_MatchSet(temp_idx_neighbours_Ref_RT,2),'.','Color','r','MarkerSize',14)
plot(ref_MatchSet(temp_idx_neighbours_Ref_MZ,1),ref_MatchSet(temp_idx_neighbours_Ref_MZ,2),'.','Color','b','MarkerSize',14)
plot(refSet(featureNr,1),refSet(featureNr,2),'x','Color',greyColor,'MarkerSize',14,'LineWidth',3)
plot(nanmedian(ref_MatchSet(temp_idx_neighbours_Ref_RT,1)),nanmedian(ref_MatchSet(temp_idx_neighbours_Ref_MZ,2)),'o','Color',greyColor,'MarkerSize',12,'LineWidth',3)
xlabel('RT reference (Minutes)'), ylabel('MZ reference (m/z units)'), grid on, axis tight
subplot(1,4,2)
plot(ref_MatchSet(:,1),target_MatchSet(:,1)-ref_MatchSet(:,1),'.k'), hold on
plot(ref_MatchSet(temp_idx_neighbours_Ref_RT,1),target_MatchSet(temp_idx_neighbours_Ref_RT,1)-ref_MatchSet(temp_idx_neighbours_Ref_RT,1),'.','Color','r','MarkerSize',14)
plot(refSet(featureNr,1),targetSet(featureNr,1)-refSet(featureNr,1),'x','Color',greyColor,'MarkerSize',14,'LineWidth',3)
plot(refSet(featureNr,1),nanmedian(target_MatchSet(temp_idx_neighbours_Ref_RT,1)-ref_MatchSet(temp_idx_neighbours_Ref_RT,1)),'o','Color',greyColor,'MarkerSize',12,'LineWidth',3)
xlabel('RT reference (Minutes)'), ylabel('RTdist (Minutes)'), grid on, axis tight
subplot(1,4,3)
plot(ref_MatchSet(:,2),target_MatchSet(:,2)-ref_MatchSet(:,2),'.k'), hold on
plot(ref_MatchSet(temp_idx_neighbours_Ref_MZ,2),target_MatchSet(temp_idx_neighbours_Ref_MZ,2)-ref_MatchSet(temp_idx_neighbours_Ref_MZ,2),'.','Color','b','MarkerSize',14)
plot(refSet(featureNr,2),targetSet(featureNr,2)-refSet(featureNr,2),'x','Color',greyColor,'MarkerSize',14,'LineWidth',3)
plot(refSet(featureNr,2),nanmedian(target_MatchSet(temp_idx_neighbours_Ref_MZ,2)-ref_MatchSet(temp_idx_neighbours_Ref_MZ,2)),'o','Color',greyColor,'MarkerSize',12,'LineWidth',3)
xlabel('MZ reference (m/z units)'), ylabel('MZdist (m/z units)'), grid on, axis tight
subplot(1,4,4)
plot(ref_MatchSet(:,3),target_MatchSet(:,3)-ref_MatchSet(:,3),'.k'), hold on
plot(ref_MatchSet(temp_idx_neighbours_Ref_FI,3),target_MatchSet(temp_idx_neighbours_Ref_FI,3)-ref_MatchSet(temp_idx_neighbours_Ref_FI,3),'.','Color','b','MarkerSize',14)
plot(refSet(featureNr,3),targetSet(featureNr,3)-refSet(featureNr,3),'x','Color',greyColor,'MarkerSize',14,'LineWidth',3)
plot(refSet(featureNr,3),nanmedian(target_MatchSet(temp_idx_neighbours_Ref_FI,3)-ref_MatchSet(temp_idx_neighbours_Ref_FI,3)),'o','Color',greyColor,'MarkerSize',12,'LineWidth',3)
xlabel('log10FI reference (log10FI units)'), ylabel('log10FIdist (log10FI units)'), grid on, axis tight
drawnow
end
%% Figure with the trends for RT, MZ, FI
M2S_figureH(0.65,0.35); set(gcf,'Name','Trends for RT, MZ and log10FI')
subplot(1,3,1),
plot(refSet(:,1),targetSet(:,1) - refSet(:,1),'.k'), axis tight, grid on
y_lim=ylim; x_lim = xlim; hold on; plot(refSet(:,1),median_RTdist_neighborsTR,'or','MarkerSize',4), grid on
ylim(y_lim); xlim(x_lim); xlabel('RT reference (Minutes)'), ylabel('RTdist (Minutes)')
subplot(1,3,2)
plot(refSet(:,2),targetSet(:,2) - refSet(:,2),'.k'), axis tight, grid on
y_lim=ylim; x_lim = xlim; hold on; plot(refSet(:,2),median_MZdist_neighborsTR,'or','MarkerSize',4), grid on
ylim(y_lim); xlim(x_lim); xlabel('MZ reference (m/z units)'), ylabel('MZdist (m/z units)')
subplot(1,3,3)
plot(refSet(:,3),targetSet(:,3) - refSet(:,3),'.k'), axis tight, grid on
y_lim=ylim; x_lim = xlim; hold on; plot(refSet(:,3),median_FIdist_neighborsTR,'or','MarkerSize',4), grid on
ylim(y_lim); xlim(x_lim); xlabel('log10FI reference'), ylabel('log10FIdist')
drawnow
%% Figure with the Residuals_X
M2S_figureH(0.65,0.35); set(gcf,'Name','Residuals for RT, MZ and log10FI')
subplot(1,3,1),
plot(refSet(:,1),Residuals_X(:,1),'.k'), axis tight, hold on, xlim1 = xlim; grid on
%plot(xlim',[0;0],'-k')
xlabel('RT of reference feature'), ylabel('RT residuals')
subplot(1,3,2),plot(refSet(:,2),Residuals_X(:,2),'.k'), axis tight, hold on, xlim1 = xlim; grid on
%plot(xlim',[0;0],'-k')
xlabel('MZ of reference feature'), ylabel('MZ residuals')
subplot(1,3,3),plot(refSet(:,3),Residuals_X(:,3),'.k'), axis tight, hold on, xlim1 = xlim; grid on
%plot(xlim',[0;0],'-k')
xlabel('Log10FI of reference feature'), ylabel('Log10FI residuals')
elseif strcmp(neighMethod,'circle')
if plotOrNot==2
% Figure to check the neighbours of one of the features
temp_idx_neighbours_Ref = (idx_neighbors_Ref(featureNr,:))';
M2S_figureH(0.65,0.35); set(gcf,'Name',['Example of neighbours for feature number ',num2str(featureNr)])
subplot(1,3,1)
plot(ref_MatchSet(:,1),ref_MatchSet(:,2),'.k'), hold on
plot(refSet(featureNr,1),refSet(featureNr,2),'ok')
plot(ref_MatchSet(temp_idx_neighbours_Ref,1),ref_MatchSet(temp_idx_neighbours_Ref,2),'.','Color','r','MarkerSize',14)
plot(refSet(featureNr,1),refSet(featureNr,2),'x','Color',greyColor,'MarkerSize',14,'LineWidth',3)
plot(nanmedian(ref_MatchSet(temp_idx_neighbours_Ref,1)),nanmedian(ref_MatchSet(temp_idx_neighbours_Ref,2)),'o','Color',greyColor,'MarkerSize',12,'LineWidth',3)
xlabel('RT reference (Minutes)'), ylabel('MZ reference (m/z units)'), grid on, axis tight
subplot(1,3,2)
plot(ref_MatchSet(:,1),target_MatchSet(:,1)-ref_MatchSet(:,1),'.k'), hold on
plot(refSet(featureNr,1),targetSet(featureNr,1)-refSet(featureNr,1),'ok')
plot(ref_MatchSet(temp_idx_neighbours_Ref,1),target_MatchSet(temp_idx_neighbours_Ref,1)-ref_MatchSet(temp_idx_neighbours_Ref,1),'.','Color','r','MarkerSize',14)
plot(refSet(featureNr,1),targetSet(featureNr,1)-refSet(featureNr,1),'x','Color',greyColor,'MarkerSize',14,'LineWidth',3)
plot(refSet(featureNr,1),nanmedian(target_MatchSet(temp_idx_neighbours_Ref,1)-ref_MatchSet(temp_idx_neighbours_Ref,1)),'o','Color',greyColor,'MarkerSize',12,'LineWidth',3)
xlabel('RT reference (Minutes)'), ylabel('RTdist (Minutes)'), grid on, axis tight
subplot(1,3,3)
plot(ref_MatchSet(:,2),target_MatchSet(:,2)-ref_MatchSet(:,2),'.k'), hold on
plot(refSet(featureNr,2),targetSet(featureNr,2)-refSet(featureNr,2),'ok')
plot(ref_MatchSet(temp_idx_neighbours_Ref,2),target_MatchSet(temp_idx_neighbours_Ref,2)-ref_MatchSet(temp_idx_neighbours_Ref,2),'.','Color','r','MarkerSize',14)
plot(refSet(featureNr,2),targetSet(featureNr,2)-refSet(featureNr,2),'x','Color',greyColor,'MarkerSize',14,'LineWidth',3)
plot(refSet(featureNr,2),nanmedian(target_MatchSet(temp_idx_neighbours_Ref,2)-ref_MatchSet(temp_idx_neighbours_Ref,2)),'o','Color',greyColor,'MarkerSize',12,'LineWidth',3)
xlabel('MZ reference (m/z units)'), ylabel('MZdist (m/z units)'), grid on, axis tight
drawnow
end
% Figure with the trends for RT and MZ
M2S_figureH(0.65,0.35); set(gcf,'Name','Trends for RT and MZ')
subplot(1,2,1),plot(refSet(:,1),targetSet(:,1) - refSet(:,1),'.k'), axis tight, grid on
y_lim=ylim; x_lim = xlim;
%subplot(2,1,2),
hold on
plot(refSet(:,1),median_RTdist_neighborsTR,'or','MarkerSize',4), grid on
ylim(y_lim); xlim(x_lim);
xlabel('RT reference (Minutes)'), ylabel('RTdist (Minutes)')
subplot(1,2,2),plot(refSet(:,2),targetSet(:,2) - refSet(:,2),'.k'), axis tight, grid on
y_lim=ylim; x_lim = xlim;
%subplot(2,1,2),
hold on
plot(refSet(:,2),median_MZdist_neighborsTR,'or','MarkerSize',4), grid on
ylim(y_lim); xlim(x_lim);
xlabel('MZ reference (m/z units)'), ylabel('MZdist (m/z units)')
% Figure with the Residuals_X
M2S_figureH(0.65,0.35); set(gcf,'Name','Residuals for RT and MZ')
subplot(1,3,1),plot(refSet(:,1),Residuals_X(:,1),'.k'), axis tight, hold on, xlim1 = xlim; grid on
%plot(xlim',[0;0],'-k')
xlabel('RT of reference feature'), ylabel('RT residuals (minutes)')
subplot(1,3,2),plot(refSet(:,2),Residuals_X(:,2),'.k'), axis tight, hold on, xlim1 = xlim; grid on
%plot(xlim',[0;0],'-k')
xlabel('MZ of reference feature'), ylabel('MZ residuals (m/z units)')
subplot(1,3,3),plot(refSet(:,3),Residuals_X(:,3),'.k'), axis tight, hold on, xlim1 = xlim; grid on
%plot(xlim',[0;0],'-k')
xlabel('Log10FI of reference feature'), ylabel('Log10FI residuals')
end
end