-
Notifications
You must be signed in to change notification settings - Fork 0
/
dscatter.m
205 lines (190 loc) · 6.13 KB
/
dscatter.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
function hAxes = dscatter(X,Y, varargin)
% DSCATTER creates a scatter plot coloured by density.
%
% DSCATTER(X,Y) creates a scatterplot of X and Y at the locations
% specified by the vectors X and Y (which must be the same size), colored
% by the density of the points.
%
% DSCATTER(...,'MARKER',M) allows you to set the marker for the
% scatter plot. Default is 's', square.
%
% DSCATTER(...,'MSIZE',MS) allows you to set the marker size for the
% scatter plot. Default is 10.
%
% DSCATTER(...,'FILLED',false) sets the markers in the scatter plot to be
% outline. The default is to use filled markers.
%
% DSCATTER(...,'PLOTTYPE',TYPE) allows you to create other ways of
% plotting the scatter data. Options are "surf','mesh' and 'contour'.
% These create surf, mesh and contour plots colored by density of the
% scatter data.
%
% DSCATTER(...,'BINS',[NX,NY]) allows you to set the number of bins used
% for the 2D histogram used to estimate the density. The default is to
% use the number of unique values in X and Y up to a maximum of 200.
%
% DSCATTER(...,'SMOOTHING',LAMBDA) allows you to set the smoothing factor
% used by the density estimator. The default value is 20 which roughly
% means that the smoothing is over 20 bins around a given point.
%
% DSCATTER(...,'LOGY',true) uses a log scale for the yaxis.
%
% Examples:
%
% [data, params] = fcsread('SampleFACS');
% dscatter(data(:,1),10.^(data(:,2)/256),'log',1)
% % Add contours
% hold on
% dscatter(data(:,1),10.^(data(:,2)/256),'log',1,'plottype','contour')
% hold off
% xlabel(params(1).LongName); ylabel(params(2).LongName);
%
% See also FCSREAD, SCATTER.
% Copyright 2003-2004 The MathWorks, Inc.
% $Revision: $ $Date: $
% Reference:
% Paul H. C. Eilers and Jelle J. Goeman
% Enhancing scatterplots with smoothed densities
% Bioinformatics, Mar 2004; 20: 623 - 628.
lambda = [];
nbins = [];
plottype = 'scatter';
contourFlag = false;
msize = 10;
marker = 's';
logy = false;
filled = true;
if nargin > 2
if rem(nargin,2) == 1
error('Bioinfo:IncorrectNumberOfArguments',...
'Incorrect number of arguments to %s.',mfilename);
end
okargs = {'smoothing','bins','plottype','logy','marker','msize','filled'};
for j=1:2:nargin-2
pname = varargin{j};
pval = varargin{j+1};
k = strmatch(lower(pname), okargs); %#ok
if isempty(k)
error('Bioinfo:UnknownParameterName',...
'Unknown parameter name: %s.',pname);
elseif length(k)>1
error('Bioinfo:AmbiguousParameterName',...
'Ambiguous parameter name: %s.',pname);
else
switch(k)
case 1 % smoothing factor
if isnumeric(pval)
lambda = pval;
else
error('Bioinfo:InvalidScoringMatrix','Invalid smoothing parameter.');
end
case 2
if isscalar(pval)
nbins = [ pval pval];
else
nbins = pval;
end
case 3
plottype = pval;
case 4
logy = pval;
Y = log10(Y);
%logx
X=log10(X);
case 5
contourFlag = pval;
case 6
marker = pval;
case 7
msize = pval;
case 8
filled = pval;
end
end
end
end
minx = min(X,[],1);
maxx = max(X,[],1);
miny = min(Y,[],1);
maxy = max(Y,[],1);
if isempty(nbins)
nbins = [min(numel(unique(X)),200) ,min(numel(unique(Y)),200) ];
end
if isempty(lambda)
lambda = 20;
end
edges1 = linspace(minx, maxx, nbins(1)+1);
ctrs1 = edges1(1:end-1) + .5*diff(edges1);
edges1 = [-Inf edges1(2:end-1) Inf];
edges2 = linspace(miny, maxy, nbins(2)+1);
ctrs2 = edges2(1:end-1) + .5*diff(edges2);
edges2 = [-Inf edges2(2:end-1) Inf];
[n,p] = size(X);
bin = zeros(n,2);
% Reverse the columns to put the first column of X along the horizontal
% axis, the second along the vertical.
[dum,bin(:,2)] = histc(X,edges1);
[dum,bin(:,1)] = histc(Y,edges2);
H = accumarray(bin,1,nbins([2 1])) ./ n;
G = smooth1D(H,nbins(2)/lambda);
F = smooth1D(G',nbins(1)/lambda)';
% = filter2D(H,lambda);
if logy
ctrs2 = 10.^ctrs2;
Y = 10.^Y;
%logX
ctrs1=10.^ctrs1;
X=10.^X;
end
okTypes = {'surf','mesh','contour','image','scatter'};
k = strmatch(lower(plottype), okTypes); %#ok
if isempty(k)
error('dscatter:UnknownPlotType',...
'Unknown plot type: %s.',plottype);
elseif length(k)>1
error('dscatter:AmbiguousPlotType',...
'Ambiguous plot type: %s.',plottype);
else
switch(k)
case 1 %'surf'
h = surf(ctrs1,ctrs2,F,'edgealpha',0);
case 2 % 'mesh'
h = mesh(ctrs1,ctrs2,F);
case 3 %'contour'
[dummy, h] =contour(ctrs1,ctrs2,F);
case 4 %'image'
nc = 256;
F = F./max(F(:));
colormap(repmat(linspace(1,0,nc)',1,3));
h =image(ctrs1,ctrs2,floor(nc.*F) + 1);
case 5 %'scatter'
F = F./max(F(:));
ind = sub2ind(size(F),bin(:,1),bin(:,2));
col = F(ind);
if filled
h = scatter(X,Y,msize,col,marker,'filled');
else
h = scatter(X,Y,msize,col,marker);
end
end
end
if logy
set(gca,'yscale','log');
set(gca,'xscale','log');
end
if nargout > 0
hAxes = get(h,'parent');
end
%%%% This method is quicker for symmetric data.
% function Z = filter2D(Y,bw)
% z = -1:(1/bw):1;
% k = .75 * (1 - z.^2);
% k = k ./ sum(k);
% Z = filter2(k'*k,Y);
function Z = smooth1D(Y,lambda)
[m,n] = size(Y);
E = eye(m);
D1 = diff(E,1);
D2 = diff(D1,1);
P = lambda.^2 .* D2'*D2 + 2.*lambda .* D1'*D1;
Z = (E + P) \ Y;