Add configuration options for customizing processing / reporting:

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,35 +299,40 @@
             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
     if ~isempty(force_peak), first_peak = force_peak-peakOffset; end
     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

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;
 

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,13 +22,20 @@
 GMMG = class(GMMG,'GMMGating',Filter());
 
 if nargin==0, return; end;
+if nargin<2, prefilters = {}; end;
 
 % Model is a gmdistribution
 
 file = ensureDataFile(file);
 
 [~, 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;
 

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

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';

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);

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