Code covered by the BSD License  

Highlights from
Flow Cytometry Data Reader and Visualization

image thumbnail

Flow Cytometry Data Reader and Visualization

by

Robert Henson (view profile)

 

06 Sep 2005 (Updated )

Flow Cytometry Data Reader and Visualization Example

fcsread(filename)
function [data,paramVals,textHeader] = fcsread(filename)
% FCSREAD read a FACS FCS format file.
%
%   DATA = FCSREAD(FILENAME) reads an FCS format file FILENAME, returning
%   the data in the file as an array with the data for each parameter
%   stored in the columns of the array. 
%
%   [DATA, PARAMS] = FCSREAD(FILENAME) returns information about the
%   parameters, such as the name, range and amplification factors in a
%   structure.
%
%   [DATA, PARAMS, HEADER] = FCSREAD(FILENAME) returns the TEXT field
%   header information in an Nx2 cell array with the field names in the
%   first column and the values in the second column.
%
%   Example:
%
%       % Read a sample FCS file
%       [data,params] = fcsread('SampleFACS.fcs');
%
%       % Display a histogram of the first parameter values
%       hist(data(:,1),max(data(:,1))-min(data(:,1)));
%       title(params(1).LongName);
%
%       % Create a scatter plot the data using a log scale
%       numDecades = params(2).Amplification(1);
%       theRange = params(2).Range;
%       semilogy(data(:,1),10.^(data(:,2)/(theRange/numDecades)),'+');
%       xlabel(params(1).Name); ylabel(params(2).Name);
%

%   Copyright 2005 The MathWorks, Inc.
%   $Revision:  $  $Date: 2005/3/24 20:41:57 $

% FCS format specified here:
% http://www.isac-net.org/resources/fcs3.htm

% open the file...

correctEndian = false;
endianCount = 1;
endianOptions = {'l','b'};
while correctEndian == false
    try
        fid = fopen(filename,'rb',endianOptions{endianCount});
        if fid == -1
            if ~exist(filename)
                error('fcsread:CannotFindFile','%s does not appear to exist.',filename)
            end
        end
    catch
        if ~ischar(filename)
            error('fcsread:InvalidInput','Input must be a character array')
        end

        if endianCount == 3
            error('fcsread:CannotDetermineEndian','%s does not appear to be a valid FCS file.',filename)
        end
        rethrow(lasterr);
    end


    % read the header
    fileHeader = readHeader(fid,filename);

    % Check endian-ness
    correctEndian = checkByteOrder(fileHeader.textData,endianCount);
    if ~correctEndian
        fclose(fid);
        endianCount = endianCount+1;
    end
end

% get the number of parameters

textData = strrep(fileHeader.textData,char(0),'');
textCell = strread(textData,'%s','delimiter','\\');
if isempty(textCell{1})
    textCell(1) = [];
end
textHeader = reshape(textCell,2,numel(textCell)/2)';
numValues = lookupNumericData(textCell,'$TOT');

% get details for each parameter
[numParams, paramVals] = getParamVals(textCell);

% read the data
if any(diff([paramVals.Bits]))
    error('fcsread:MixedIntegersNotSupported',...
        'Data in %s is stored in an unsupported format.',filename)
end
checkInts = strcmpi('I',lookupTextData(textCell,'$DATATYPE'));
if ~checkInts
    error('fcsread:NonIntegersNotSupported',...
        'Data in %s is stored in an unsupported format.',filename)
end
checkMode = strcmpi('L',lookupTextData(textCell,'$MODE'));
if ~checkMode
    error('fcsread:NonListNotSupported',...
        'Data in %s is stored in an unsupported format.',filename)
end

intSize = paramVals(1).Bits;
dataFormat = sprintf('uint%d',intSize);

% for now we only support one size of uints in a binary list.
fseek(fid,fileHeader.dataStart,-1);
data = fread(fid,numValues*numParams,dataFormat);
fclose(fid);
% Check we read the right amount of data
if numel(data) ~= numValues*numParams
    error('fcsread:DataSizeMismatch',...
        'Error reading the correct number of data values.')
end

% reshape the data into columns
data = reshape(data,numParams,numValues)';

% check ranges
maxParams = max(data);
paramRanges = [paramVals.Range];
if any(maxParams > paramRanges)
    error('fcsread:DataRangeError',...
        'Error reading the correct number of data values.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [numParams, paramVals] = getParamVals(textCell);
numParams = lookupNumericData(textCell,'$PAR');


for count = 1:numParams
    paramVals(count).Name = lookupTextData(textCell,sprintf('$P%dN',count));
    paramVals(count).Range = lookupNumericData(textCell,sprintf('$P%dR',count));
    paramVals(count).Bits= lookupNumericData(textCell,sprintf('$P%dB',count));
    paramVals(count).Amplification = lookupNumericData(textCell,sprintf('$P%dE',count));
    paramVals(count).LongName= lookupTextData(textCell,sprintf('$P%dS',count));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function header = readHeader(fid,filename)
% Extract the header information

% we could check here for
header.ver = fread(fid,6,'*char')';
if isempty(strfind(header.ver,'FCS'))
    error('fcsread:NoFCSTag','%s does not appear to be a valid FCS file.',filename)
end
header.skip = fread(fid,4,'*char');
header.textStart = str2num(fread(fid,8,'*char')');
header.textEnd = str2num(fread(fid,8,'*char')');
header.dataStart = str2num(fread(fid,8,'*char')');
header.dataEnd = str2num(fread(fid,8,'*char')');
header.analysisStart = str2num(fread(fid,8,'*char')');
header.analysisEnd = str2num(fread(fid,8,'*char')');
fseek(fid,header.textStart,-1);
header.textData = fread(fid,header.textEnd -header.textStart+1,'*char')';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function checkEndian = checkByteOrder(textData,endianCount);
% check the endianness of the file.
if endianCount == 1
    checkEndian = any(strfind(textData,'1,2,3,4'));
else
    checkEndian = any(strfind(textData,'4,3,2,1'));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val = lookupNumericData(textCell,fieldName)

val = find(strcmpi(textCell,fieldName));
if ~isempty(val)
    val = str2num(textCell{val+1});
else
    val = NaN;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val = lookupTextData(textCell,fieldName)

val = find(strcmpi(textCell,fieldName));
if ~isempty(val)
    val = textCell{val+1};
else
    val = '';
end

Contact us