|
|
楼主 |
发表于 2012-7-9 11:45:15
|
显示全部楼层
haven_t 发表于 2011-8-3 16:49
在Modules 中看到Logicle 和Hyperlog，恕我多口提醒一下，Logicle 转换似乎在美国申请了专利，您可以上美国 ...
haven_t，看到一份有关logicle算法的matlab代码，声明里面说了，虽然stanford大学拥有专利权，但他们不限制非营利或者赢利性应用的使用。
- function [returnHistogramData, returnBins] = logicleHist(FCSdata, varargin)
- %logicleHist Plots data as a histogram on a biexponential X-axis.
- %logicleHist(FCSdata, options)
- % Uses the equation from "A New 'Logicle' Display Method Avoids Deceptive
- % Effects of Logarithmic Scaling for Low Signals and Compensated Data"
- % 4/3/2012: Modified using "Update for the logicle data scale including
- % operationl code implementation."
- % logicleHist(FCSdata) creates a plot using the default scaling parameters
- % [Hist, bins] = logicleHist(FCSdata) creates a plot and returns the data
- %
- % Options:
- % 'W','T','M','A','Span' will adjust the biexponential scaling.
- % ... 'Smooth' true) will use a loess to smooth the data
- % Normalize determines how to scale the y-axis:
- % 1 - length of FCSdata (# of points)
- % 2 - max of histogram
- % 3 - area of histogram
- % 4 - none (don't use this)
- %
- %% Copyright stuff...
- % Copyright (c) 2009, 2011, The Board of Trustees of The Leland Stanford
- % Junior University. All rights reserved. Redistribution and use in source
- % and binary forms, with or without modification, are permitted provided
- % that the following conditions are met:
- %
- % Redistributions of source code must retain the above copyright notice,
- % this list of conditions and the following disclaimer. Redistributions in
- % binary form must reproduce the above copyright notice, this list of
- % conditions and the following disclaimer in the documentation and/or other
- % materials provided with the distribution. Neither the name of the Leland
- % Stanford Junior University nor the names of its contributors may be used
- % to endorse or promote products derived from this software without
- % specific prior written permission.
- % THIS SOFTWARE IS PROVIDED BY THE
- % COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
- % WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
- % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
- % NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY
- % DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- % DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- % OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- % HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
- % STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
- % ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- % POSSIBILITY OF SUCH DAMAGE.
- %
- % The Logicle method is patented under United States Patent 6,954,722.
- % However, Stanford University does not enforce the patent for non-profit
- % academic purposes or for commercial use in the field of flow cytometry.
- %% Parse the input parameters
- p = inputParser; %special element
- %defaults
- DEFAULT_M = 4.5; % number of logarithmic decades
- DEFAULT_T = 2^18; % maximum value (18 bit storage)
- DEFAULT_SPAN = [-200 DEFAULT_T]; %min and max value
- DEFAULT_W = (DEFAULT_M - log10(DEFAULT_T/abs(DEFAULT_SPAN(1))))/2;
- %size of linear range
- DEFAULT_A = 0; %Extra negative range (will start being exponential)
- DEFAULT_NBINS = 250;
- DEFAULT_SMOOTH = true;
- DEFAULT_NORMALIZE = 2;
- DEFAULT_PLOT = true;
- addRequired(p,'FCSdata',@isnumeric);
- addParamValue(p,'W',DEFAULT_W,@isnumeric);
- addParamValue(p,'T',DEFAULT_T,@isnumeric);
- addParamValue(p,'M',DEFAULT_M,@isnumeric);
- addParamValue(p,'A',DEFAULT_A,@isnumeric);
- addParamValue(p,'Span',DEFAULT_SPAN,@isnumeric);
- addParamValue(p,'nBins',DEFAULT_NBINS,@isnumeric);
- addParamValue(p,'Smooth',DEFAULT_SMOOTH,@islogical);
- addParamValue(p,'Normalize',DEFAULT_NORMALIZE,@isscalar);
- addParamValue(p,'Plot',DEFAULT_PLOT,@islogical);
- parse(p,FCSdata, varargin{:});
- W = p.Results.W;
- M = p.Results.M;
- T = p.Results.T;
- A = p.Results.A;
- span = p.Results.Span;
- %note: this is not taking into account A. If you are clever enough to define
- %A, then just define everything.
- if (any(strcmp('Span',p.UsingDefaults))) %span is not specified
- if (any(strcmp('T',p.UsingDefaults))) %T is also undefined
- if (any(strcmp('W',p.UsingDefaults))) %W is undefined
- %nothing is defined, use defaults
- else %span/T is not defined by W is.
- span(1) = T / 10^(W*2+M);
- end
- else %T is defined
- if (any(strcmp('W',p.UsingDefaults))) %W is undefined
- span(2) = T;
- else %T and W are defined, but span isn't
- span = [(-1*T / 10^(M-W*2)), T];
- end
- end
- else %span is defined
- if (any(strcmp('T',p.UsingDefaults))) %T is undefined
- T = span(2);
- end
- if (any(strcmp('W',p.UsingDefaults))) %W is undefined
- W = (M - log10(T/span(1)))/2;
- end
- end
- nBins = p.Results.nBins;
- smoothHistogram = p.Results.Smooth;
- normalize = p.Results.Normalize;
- plotData = p.Results.Plot;
- %Check the inputs
- if T<0
- error('T cannot be negative');
- end
- if (W<0) || (W>M/2)
- error('W must satisfy 0 <= W <= M/2')
- end
- if (-1*W > A) || (A > M-2*W)
- error('A must satisfy -W <= A <= M-2W');
- end
- if (normalize < 1) || (normalize >4)
- error('normalize must be between 1 and 4')
- end
- %% Start the algorithm
- %Convert the parameters to the new biexponential scale
- [a, b, c, d, f, x0, x1, x2] = convertParm(T, M, W, A);
- %Bins are linearly spaced in X domain
- binEdges = linspace(0,1,nBins);
- binCenter = linspace((1/nBins/2),(1-1/nBins/2),(nBins-1));
- %Convert the bins back to the S domain
- %The bins are now biexponentially spaced
- for i = 1:nBins
- sTransform(i) = logicleFun(binEdges(i), a, b, c, d, f, x1, 0);
- end
- %generate the histogram information
- histData = histc(FCSdata,sTransform);
- %the last data point is trash (see histc documentation)
- histData(end) = [];
- %normalize the histogram
- if normalize == 1
- histData = histData./length(FCSdata);
- elseif normalize == 2
- histData = histData./max(histData);
- elseif normalize == 3
- histData = histData./trapz(histData);
- elseif normalize == 4
- %Do nothing
- end
- if smoothHistogram
- %Smooth the histogram
- histData = smooth(binCenter,histData,0.01,'loess');
- end
- if plotData
- %plot the results
- plot(binCenter,histData)
- %Format the figure
- ylim('auto')
- xlabel('Fluorescence (a.u.)');
- ylabel('Count');
- set(gca,'YTick',[])
- yLimits = ylim(gca);
- yLimits(1) = 0;
- ylim(yLimits);
- %% Generate the X tick marks and labels
- options = optimset('Display','off','TolFun',1e-4);
- LINEAR_TICKMARK_SPACING = 50;
- END_LINEAR_SPACING = logicleFun(x0, a, b, c, d, f, x1, 0);
- %linear portion
- sTicks = (floor(span(1)/10)*10):LINEAR_TICKMARK_SPACING: ...
- floor(END_LINEAR_SPACING/10)*10;
- %log portion
- decade = floor(log10(max(sTicks)));
- startLogSpace = 10^decade+ max(sTicks);
- sTicks = [sTicks startLogSpace:10^decade:10^(decade+1)];
- decade = decade+1;
- while max(sTicks)<span(2)
- sTicks = [sTicks 2*10^decade:10^decade:10^(decade+1)];
- decade = decade+1;
- end
- sTicks(sTicks>span(2))=[]; %remove extra sTicks
- for i=1:length(sTicks)
- if (sTicks(i)<=100)
- if(rem(sTicks(i),100)==0)
- xLabel{i}=num2str(sTicks(i));
- else
- xLabel{i}=' ';
- end
- elseif floor(log10(sTicks(i)))==log10(sTicks(i))
- xLabel{i}=num2str(sTicks(i));
- else
- xLabel{i}=' ';
- end
- S=sTicks(i);
- xTickGuess = xtickguess(a,b,c,d,f,x1,S);
- xTicks(i) = fsolve(@(X)logicleFun(X,a,b,c,d,f,x1,S),xTickGuess,options);
- end
- set(gca,'XTick',xTicks);
- set(gca,'XTickLabel',xLabel);
- end
- %% Return the information
- if nargout == 2
- returnHistogramData = histData; returnBins=binCenter;
- end
- end
- function [a, b, c, d, f, x0, x1, x2] = convertParm(T, M, W, A)
- b = (M+A)*log(10);
- w = W / (M+A);
- x2 = A / (M+A);
- x1 = x2 + w;
- x0 = x1 + w;
- d = rtsafe(w,b);
- if (d<0) || (d>b), error('d must satisfy 0 < d < b'); end
- cDivA = exp((b+d)*x0);
- fDivA = -1*( exp(b*x1) - cDivA*exp(-1*d*x1));
- a = T / (exp(b) - cDivA*exp(-1*d) + fDivA);
- c = cDivA * a;
- f = fDivA * a;
- end
- function B = logicleFun(x, a, b, c, d, f, x1, S)
- %logicleFun Computes the logicle transform
- % Uses the equation from "A New 'Logicle' Display Method Avoids Deceptive
- % Effects of Logarithmic Scaling for Low Signals and Compensated Data"
- if nargin == 7 %not necessary to define S
- S = 0;
- elseif nargin == 6 %x1 and S are undefined (lazy)
- S = 0;
- x1 = fsolve(@(x) a*exp(b*x) - c*exp(-1*d*x) + f,0.1);
- end
- if x<x1 %negative regime
- x = x1 + (x1-x);
- B = -1*a*exp(b*x) + c*exp(-1*d*x) - f - S;
- else %positive regime
- B = a*exp(b*x) - c*exp(-1*d*x) + f - S;
- end
- end
- function X = xtickguess(a, b, c, d, f, x1, S)
- %This speads up fsolve
- if S > f %exponential regime
- X = log(S/a) / b;
- else %linear regime
- slope = a*b*exp(b*x1)+c*d*exp(-1*d*x1);
- X = S/slope + x1;
- end
- end
- function d = rtsafe(w,b)
- %Modified from 'Numerical recipes: the art of scientific computing'
- %solves the following equation for d :
- % w = 2ln(d/b) / (b+d)
- %the function
- gFunction = @(d)(w*(b+d) + 2*(log(d)-log(b)));
- derivFunction = @(d)(w+2/d);
- X_ACCURACY = 0.0000001;
- MAX_IT = 100;
- lowerLimit = 0;
- upperLimit = b;
- if ((gFunction(lowerLimit)>0) && (gFunction(upperLimit)>0)) || ...
- ((gFunction(lowerLimit)<0) && (gFunction(upperLimit)<0))
- error('Root must be bracketed');
- end
- if gFunction(lowerLimit)<0
- xLow = lowerLimit;
- xHigh = upperLimit;
- else
- xLow = upperLimit;
- xHigh = lowerLimit;
- end
- root = (lowerLimit + upperLimit)/2;
- dxOld = abs(upperLimit - lowerLimit);
- dx = dxOld;
- g = gFunction(root);
- dg = derivFunction(root);
- for i = 0:MAX_IT
- if (((root-xHigh)*dg-g)*((root-xLow)*dg-g) > 0) || ...
- (abs(2*g) > abs(dxOld*dg))
- %Bisect method
- dxOld = dx;
- dx = (xHigh-xLow)/2;
- root = xLow + dx;
- else
- %Newton method
- dxOld = dx;
- dx = g/dg;
- root = root - dx;
- end
-
- if (abs(dx) < X_ACCURACY)
- d = root;
- return;
- end
-
- g = gFunction(root);
- dg = derivFunction(root);
-
- if (g < 0)
- xLow = root;
- else
- xHigh = root;
- end
- end
- error('Maximum number of interations exceeded in rtsafe');
- end
|
|
楼主: niwanmao
发表于 2012-7-10 21:37:07
