diff --git a/code/@ColorModel/beads_to_ERF_model.m b/code/@ColorModel/beads_to_ERF_model.m index a42d4e0..acbdecc 100644 --- a/code/@ColorModel/beads_to_ERF_model.m +++ b/code/@ColorModel/beads_to_ERF_model.m @@ -19,7 +19,7 @@ % exception, as described in the file LICENSE in the TASBE analytics % package distribution's top directory. -function [UT CM] = beads_to_ERF_model(CM, beadfile) +function [UT, CM] = beads_to_ERF_model(CM, beadfile) ERF_channel = CM.ERF_channel; if(isUnprocessed(ERF_channel)) TASBESession.error('TASBE:Beads','ERFUnprocessed','ERF channel %s cannot be unprocessed',getName(ERF_Channel)); @@ -29,9 +29,12 @@ visiblePlots = TASBEConfig.get('beads.visiblePlots'); plotPath = TASBEConfig.get('beads.plotPath'); plotSize = TASBEConfig.get('beads.plotSize'); +beadERFs = TASBEConfig.getexact('beads.beadERFs',[]); beadModel = TASBEConfig.get('beads.beadModel'); beadChannel = TASBEConfig.get('beads.beadChannel'); beadBatch = TASBEConfig.getexact('beads.beadBatch',[]); +fitThreshold = TASBEConfig.get('beads.fitThreshold'); +widthThreshold = TASBEConfig.get('beads.widthThreshold'); force_peak = TASBEConfig.getexact('beads.forceFirstPeak',[]); if ~isempty(force_peak) @@ -48,10 +51,20 @@ erfChannelName=getName(ERF_channel); i_ERF = find(CM,ERF_channel); -[PeakERFs,units,actualBatch] = get_bead_peaks(beadModel,beadChannel,beadBatch); -CM.standardUnits = units; -if ~strcmp(units,'MEFL') - TASBESession.warn('TASBE:Beads','NonMEFL','MEFL units are recommended, rather than %s',units); +if numel(beadERFs)>0 + units = beadERFs{1}; + PeakERFs = beadERFs{2}; + if numel(beadERFs)>2 + actualBatch = beadERFs{3}; + else + actualBatch = 'Unspecified Lot'; + end +else + [PeakERFs,units,actualBatch] = get_bead_peaks(beadModel,beadChannel,beadBatch); + CM.standardUnits = units; + if ~strcmp(units,'MEFL') + TASBESession.warn('TASBE:Beads','NonMEFL','MEFL units are recommended, rather than %s',units); + end end ERF_channel_idx = indexof(CM.Channels,CM.ERF_channel); CM.Channels{ERF_channel_idx} = setUnits(CM.Channels{ERF_channel_idx},units); @@ -147,6 +160,11 @@ segment_peak_means(n_peaks) = mean(segment_data(which)); % arithmetic q. beads are primarily measurement noise peak_means(n_peaks) = mean(bead_data(which)); peak_counts(n_peaks) = sum(which); + % check for problematic width + peak_width = max(bead_data(which))/min(bead_data(which)); + if peak_width>widthThreshold + TASBESession.warn('TASBE:Beads','WidePeak','Unexpectedly wide peak: %.2d-fold range',peak_width); + end end end end @@ -215,7 +233,7 @@ if ~issorted(alt_peak_maximas) TASBESession.warn('TASBE:Beads','PeaksNotAscending','Peaks are not in ascending order. May need to adjust rangeMin or rangeMax for %s.',clean_for_latex(getPrintName(CM.Channels{i}))); if alt_peak_maximas(alt_n_peaks) < max(alt_peak_maximas) %if the last peak is lower than the maximum peak - TASBESession.warn('TASBE:Beads','PotentialBeadClump','Last peak may consist of beads stuck together. May need to adjust rangeMax or peakThresholdfor %s.',clean_for_latex(getPrintName(CM.Channels{i}))); + TASBESession.warn('TASBE:Beads','PotentialBeadClump','Last peak may consist of beads stuck together. May need to adjust rangeMax or peakThreshold for %s.',clean_for_latex(getPrintName(CM.Channels{i}))); end end % Check to make sure that the true peak was properly identified @@ -233,7 +251,7 @@ if makePlots graph_max = max(alt_range_bin_counts); h = figure('PaperPosition',[1 1 plotSize]); - if(~visiblePlots), set(h,'visible','off'); end; + if(~visiblePlots), set(h,'visible','off'); end semilogx(range_bin_centers,alt_range_bin_counts,'b-'); hold on; for j=1:alt_n_peaks semilogx([alt_peak_means(j) alt_peak_means(j)],[0 graph_max],'r-'); @@ -281,7 +299,9 @@ TASBESession.succeed('TASBE:Beads','PeakIdentification','Matched multiple peaks in reasonable range'); end else % 2 peaks - TASBESession.warn('TASBE:Beads','PeakIdentification','Only two bead peaks found, assuming brightest two'); + if numQuantifiedPeaks>2 % Only warn that we're using two if we were expecting more than two + TASBESession.warn('TASBE:Beads','PeakIdentification','Only 2 bead peaks found out of %i, assuming brightest two',numQuantifiedPeaks); + end [poly,S] = polyfit(log10(peak_means),log10(quantifiedPeakERFs(end-1:end)),1); fit_error = S.normr; model = poly; first_peak = numQuantifiedPeaks-1; end @@ -289,27 +309,30 @@ constrained_fit = mean(log10(quantifiedPeakERFs((1:n_peaks)+first_peak-1)) - log10(peak_means)); cf_error = mean(10.^abs(log10((quantifiedPeakERFs((1:n_peaks)+first_peak-1)./peak_means) / 10.^constrained_fit))); % Final fit_error should be close to zero / 1-fold - if(cf_error>1.05), - TASBESession.warn('TASBE:Beads','PeakFitQuality','Bead calibration may be incorrect: fit more than 5 percent off: error = %.2d',cf_error); + if(cf_error>(1+fitThreshold)) + TASBESession.warn('TASBE:Beads','PeakFitQuality','Bead calibration may be incorrect: fit more than %i percent off: error = %.2d',fitThreshold*100,cf_error); else TASBESession.succeed('TASBE:Beads','PeakFitQuality','Bead fit quality acceptable: error = %.2d',cf_error); - end; + end %if(abs(model(1)-1)>0.05), TASBESession.warn('TASBE:Beads','Bead calibration probably incorrect: fit more than 5 percent off: slope = %.2d',model(1)); end; + fit_error = cf_error; % report on the constrained fit k_ERF = 10^constrained_fit; elseif(n_peaks==1) % 1 peak - TASBESession.warn('TASBE:Beads','PeakIdentification','Only one bead peak found, assuming brightest'); - TASBESession.skip('TASBE:Beads','PeakFitQuality','Fit quality irrelevant for single peak'); + if numQuantifiedPeaks>1 % Only warn that we're using one if we were expecting more than one + TASBESession.warn('TASBE:Beads','PeakIdentification','Only 1 bead peak found out of %i, assuming brightest',numQuantifiedPeaks); + TASBESession.skip('TASBE:Beads','PeakFitQuality','Fit quality irrelevant for single peak'); + end fit_error = 0; first_peak = numQuantifiedPeaks; if ~isempty(force_peak), first_peak = force_peak-peakOffset; end k_ERF = quantifiedPeakERFs(first_peak)/peak_means; else % n_peaks = 0 TASBESession.warn('TASBE:Beads','PeakIdentification','Bead calibration failed: found no bead peaks; using single dummy peak'); - TASBESession.skip('TASBE:Beads','PeakFitQuality','Fit quality irrelevant for single peak'); + TASBESession.skip('TASBE:Beads','PeakFitQuality','Fit quality irrelevant for dummy peak'); k_ERF = 1; fit_error = Inf; first_peak = NaN; CM.standardUnits = 'arbitrary units'; -end; +end % Plot fitted channel if makePlots diff --git a/code/@ColorModel/find.m b/code/@ColorModel/find.m index 07c8cbd..52a4d63 100644 --- a/code/@ColorModel/find.m +++ b/code/@ColorModel/find.m @@ -11,25 +11,38 @@ function index = find(CM, channel) +if ischar(channel) % if we were given a string, just match on name + channel_name = channel; +else % otherwise, assume it's a channel object + channel_name = getName(channel); +end + foundset=[]; nameeqset=[]; for i=1:numel(CM.Channels) - if(CM.Channels{i} == channel) - foundset(end+1) = i; - if(strcmp(getName(CM.Channels{i}),getName(channel))) + if ischar(channel) % if we were given a string, just match on name (either given or print) + if(strcmp(getName(CM.Channels{i}),channel) || strcmp(getPrintName(CM.Channels{i}),channel)) + foundset(end+1) = i; nameeqset(end+1) = i; end - end; + else % otherwise, assume it's a channel object + if(CM.Channels{i} == channel) + foundset(end+1) = i; + if(strcmp(getName(CM.Channels{i}),channel_name)) + nameeqset(end+1) = i; + end + end; + end end; if(isempty(foundset)), - TASBESession.error('TASBE:ColorModel','MissingChannel','Unable to find channel %s',getName(channel)); + TASBESession.error('TASBE:ColorModel','MissingChannel','Unable to find channel %s',channel_name); elseif numel(foundset)==1 index = foundset(1); elseif numel(nameeqset)==1 % more than one match, but only one matches name precisely % Turning this warning off, since it gets called too much and isn't very useful - %TASBESession.warn('TASBE:ColorModel','DisambiguateChannel','Multiple channels match %s, discriminating by name',getName(channel)); + %TASBESession.warn('TASBE:ColorModel','DisambiguateChannel','Multiple channels match %s, discriminating by name',channel_name); index = nameeqset(1); else - TASBESession.error('TASBE:ColorModel','MultipleChannels','Multiple channels match %s, and cannot discriminate by name',getName(channel)); + TASBESession.error('TASBE:ColorModel','MultipleChannels','Multiple channels match %s, and cannot discriminate by name',channel_name); end; \ No newline at end of file diff --git a/code/@GMMGating/GMMGating.m b/code/@GMMGating/GMMGating.m index 6cd08e5..c995efa 100644 --- a/code/@GMMGating/GMMGating.m +++ b/code/@GMMGating/GMMGating.m @@ -10,7 +10,7 @@ % exception, as described in the file LICENSE in the TASBE analytics % package distribution's top directory. -function varargout = GMMGating(file) +function varargout = GMMGating(file, prefilters) GMMG.selected_components = []; GMMG.channel_names = {}; @@ -22,6 +22,7 @@ GMMG = class(GMMG,'GMMGating',Filter()); if nargin==0, return; end; +if nargin<2, prefilters = {}; end; % Model is a gmdistribution @@ -29,6 +30,12 @@ [~, fcshdr, rawfcs] = fca_read(file); +% apply supplied prefilters +for i=1:numel(prefilters) + rawfcs = applyFilter(prefilters{i},fcshdr,rawfcs); +end + + channel_names = TASBEConfig.get('gating.channelNames'); n_channels = numel(channel_names); @@ -179,8 +186,8 @@ if(~visiblePlots), set(h,'visible','off'); end; %smoothhist2D([channel_data{1} channel_data{2}],10,[200, 200],[],type,range,largeOutliers); smoothhist2D([channel_data(:,i) channel_data(:,i+1)],5,[500, 500],[],type,range,largeOutliers); - xlabel([clean_for_latex(channel_names{i}) ' a.u.']); - ylabel([clean_for_latex(channel_names{i+1}) ' a.u.']); + xlabel(['log_{10} ' clean_for_latex(channel_names{i}) ' a.u.']); + ylabel(['log_{10} ' clean_for_latex(channel_names{i+1}) ' a.u.']); title('2D Gaussian Gate Fit'); hold on; diff --git a/code/@RangeFilter/applyFilter.m b/code/@RangeFilter/applyFilter.m index a3a370b..9067611 100644 --- a/code/@RangeFilter/applyFilter.m +++ b/code/@RangeFilter/applyFilter.m @@ -35,7 +35,7 @@ end else % assume it's a cell array of channel objects: for j=1:numel(fcshdr) - if(strcmp(RF.channels{i},getName(fcshdr{j}))) + if(strcmp(RF.channels{i},getName(fcshdr{j})) || strcmp(RF.channels{i},getPrintName(fcshdr{j}))) channeldata=rawfcs(:,j); found=true; continue; end diff --git a/code/TASBEConfig.m b/code/TASBEConfig.m index 93c8af2..8312e5b 100644 --- a/code/TASBEConfig.m +++ b/code/TASBEConfig.m @@ -25,6 +25,8 @@ % Generic flow data analysis s.flow = struct(); doc.flow = struct(); doc.flow.about = 'General settings for flow cytometry data analysis'; + doc.flow.ignoreChannelParameters = 'If true, ignore "$PnS" parameter entries in FCS headers'; + s.flow.ignoreChannelParameters = 0; doc.flow.maxEvents = 'Drop events above this count to avoid memory issues'; s.flow.maxEvents = 1e6; doc.flow.rangeMin = 'bin minimum (log10 scale)'; @@ -249,6 +251,8 @@ s.beads.rangeMax = 7; doc.beads.binIncrement = 'Resolution of histogram bins used for finding bead peaks (log10 scale)'; s.beads.binIncrement = 0.02; + doc.beads.beadERFs= 'When set, override model information and use supplied {unit, ERF} values instead'; + s.beads.beadERFs = []; doc.beads.beadModel = 'Model of beads that are being used. Should match an option in BeadCatalog.xlsx'; s.beads.beadModel = 'SpheroTech RCP-30-5A'; doc.beads.beadChannel = 'Laser/filter channel being used, defaults to FITC (MEFL). Should match an option in BeadCatalog.xlsx'; @@ -257,6 +261,10 @@ s.beads.beadBatch = []; doc.beads.forceFirstPeak = 'If set to N, lowest observed bead peak is forced to be interpreted as Nth peak'; s.beads.forceFirstPeak = []; + doc.beads.fitThreshold = 'If fit error is at or above this threshold, report a warning'; + s.beads.fitThreshold = 0.05; + doc.beads.widthThreshold = 'If at least one peak width is at or above this threshold (fold), report a warning'; + s.beads.widthThreshold = 1.5; doc.beads.plot = 'When true, make diagnostic plots while computing bead unit calibration'; s.beads.plot = []; defaults('beads.plot') = 'calibration.plot'; diff --git a/code/fca_readfcs.m b/code/fca_readfcs.m index 7cd7d39..57fabfd 100644 --- a/code/fca_readfcs.m +++ b/code/fca_readfcs.m @@ -238,7 +238,11 @@ % $PnS Name used for parameter n. replacename = get_mnemonic_value(['$P',num2str(i),'S'],fcsheader_main, mnemonic_separator); % Replace name if nickname is present and not just whitespace - if(numel(replacename)>0 && ~isempty(find(~isspace(replacename), 1))), fcshdr.par(i).name = replacename; end; + if(numel(replacename)>0 && ~isempty(find(~isspace(replacename), 1))), + if ~TASBEConfig.get('flow.ignoreChannelParameters'), + fcshdr.par(i).name = replacename; + end + end; % FCS 3.1-only parameters: if(strcmp(fcsheader_type,'FCS3.1')), fcshdr.par(i).calibration = get_mnemonic_value(['$P',num2str(i),'CALIBRATION'],fcsheader_main, mnemonic_separator); diff --git a/code/plot_plusminus_comparison.m b/code/plot_plusminus_comparison.m index b598d55..25aedc0 100644 --- a/code/plot_plusminus_comparison.m +++ b/code/plot_plusminus_comparison.m @@ -134,7 +134,7 @@ function plot_plusminus_comparison(all_pm_results, batch_names) legend(lines, legendentries,'Location','Best'); if(TASBEConfig.isSet('OutputSettings.FixedInputAxis')), xlim(TASBEConfig.get('OutputSettings.FixedInputAxis')); end if(TASBEConfig.isSet('OutputSettings.FixedOutputAxis')), ylim(TASBEConfig.get('OutputSettings.FixedOutputAxis')); end -title(['Raw All ',clean_for_latex(stemName),' Transfer Curves']); +title(['All ',clean_for_latex(stemName),' Transfer Curves']); if(strcmp(clean_for_latex(stemName), ' ') || strcmp(clean_for_latex(stemName), '')) outputfig(h,'all-mean',directory); else