McClean: Matlab Code for Analyzing FACS Data

Summary

This is a collection of Matlab functions for plotting flow cytometry data.

All of the code requires the function fca_readfcs written by Laszlo Balkay available from the Mathworks file exchange at: FCS data reader

Information about flow cytometry at Princeton can be found at Flow Cytometry Resource Facility

Example

The following example code is for visualization of a time series of GFP fluorescence (expressed from pGAL1-GFP in yeast cells) collected using flow cytometry. The full m-file along with the individual functions can be found here:

Example: Visualizing a time-series of GFP fluorescence measured with flow cytometry

Variables to pass to the individual functions:

%% Some variables that can be quickly changed to modify this for another timecourse
NameFiles='165'; %Name common to all the files in the timecourse (For this timecourse the files are named '_165T0.fcs','_165T1.fcs',....'_165T7.fcs')
Time=[0 1  2 3 4 5 6 7]; %Timepoints
Output=0; %Print the output (1) or don't (0) to the screen from the various functions
column=3; %Need column 4 for the mCherry data, 3 for GFP, this may be different for your FACS experiment
Bins=linspace(0,10,100); %Bins to use for plotting histograms

User gates on the side- and forward-scatter

%% 1. User visually defines the gates based on side-scatter and forward-scatter using the mouse to map out a region
[Gates,FNames]=FindGate(strcat('*',NameFiles,'*','fcs'), Output);
save(strcat(NameFiles,'_GateVariables'),'Gates','FNames'); %Save gate info so analysis can be reproduced later

Plot histograms of the

%% 2. Histograms of the timecourses are plotted (FACSSTC also saves a plot of the histograms)
 [LOGM]=FACSTC_PlotHist(FNames, Gates, column, Bins,strcat('h',NameFiles,'_histograms'));
 pause;
 close all;

%% 3. Find peaks by fitting logarithm of intensities to normal distributions (this assumes that the intensity data is originally lognormal distributed) (FACSTimecourse)


%%Fit to normals or two normals (the logarithm of the data) and then plot
%%the centers of the distributions as a function of time
[P1 P2]=FACSTimecourse(FNames, column,  Gates,Output); %P1 contains the 1 gaussian fit (means and standard deviation) P2 contains the 2 guassian fit (mean in col 1, std deviation in col2)

%%Adjust the first returned by FACSTimecourse so that they are in the
%%correct order.  For this particular dataset, because of the ways the
%%files were names, _168 T18 comes directly after _168T1 not at the end
%%where it should.
%  P1=[P1; P1(3,:)];
%  P2=[P2; P2(3,:)];
% 
%  P1(3,:)=[];
%  P2(3,:)=[];

 figure; hold;
 H=shadedLogErrorBar(Time, exp(P1(:,1)),exp(P1(:,2)),{'bo-','LineWidth',2});
 title('One Gaussian Fit','FontSize',14);
 xlabel('Time (hours)','FontSize',14); ylabel('Intensity [a.u.]','FontSize',14);
 H=gcf;
 saveas(H,strcat(NameFiles,'_low'),'fig')
 saveas(H,strcat(NameFiles,'_low'),'png')

figure; hold;
 H=shadedLogErrorBar(Time, max(exp(P2(:,2)),exp(P2(:,3))),max(exp(P2(:,4)),exp(P2(:,5))),{'ro-','LineWidth',2});
title('Two Gaussian Fit','FontSize',14);
xlabel('Time (hours)','FontSize',14); ylabel('Intensity [a.u.]','FontSize',14);
H=gcf;
 saveas(H,strcat(NameFiles,'_high'),'fig')
 saveas(H,strcat(NameFiles,'_high'),'png')
pause; close all;

 %% 4. Find peaks using a different function (FACSTC_PeakFind) based on Matlab's findpeaks
%Find the midpoints of populations by using a different function based on
%Matlab's findpeaks function:
[PeakLocations]=FACSTC_PeakFind(FNames, column,linspace(0,10,100),Gates,Output);

P3=PeakLocations;
close all;

figure; hold;
plot(Time,exp(P3(:,1)),'ko-','LineWidth',2);
plot(Time,exp(P3(:,2)),'bo-','LineWidth',2);
title('Induced and Uninduced Population Intensity Modes','FontSize',14);
xlabel('Time (hours)','FontSize',14); ylabel('Intensity [a.u.]','FontSize',14);
h=gcf;
saveas(h,strcat(NameFiles,'_peakfind'),'fig')
saveas(h,strcat(NameFiles,'_peakfind'),'png')
pause; close all;

%%
%%Save the fitted peak centers to be used later for plotting
P1_low=P1;
P2_high=P2;
P3_peaks=P3;
save(strcat(NameFiles,'_Peaks'),'P1_low','P2_high','P3_peaks');

Code

function h=EBplotyy(x,y)
%INPUT:
%x-independent variable (often Time)
%y-dependent variable
%OUTPUT:
%h-handle of the errorbar graphics object
%
%Save this code in an m-file named EBplotyy.m

s=nanstd(y);
h=errorbar(x,nanmedian(y),nanstd(y));  %NOTE: You might want to change nanmedian to nanmean, nanstd to standard error, etc depending on what you want to plot

Notes

Please feel free to post comments, questions, or improvements to this protocol. Happy to have your input!

• Megan N McClean 17:27, 30 January 2012 (EDT): Obviously this is nothing fancy, but the code has come in handy for me a number of times, so I thought I would stick it on the wiki in case it is useful to anyone else in the lab.

References

Mathworks Online Help: plotyy

Matlab Newsreader: Plotyy with errorbar

Contact

• Megan N McClean 14:01, 30 January 2012 (EDT)

or instead, discuss this protocol.