%SUMMARY:
%This is an example of how to use the following functions to analyze FACS
%data.  The process is:
% 1. Gate on side-scatter and forward-scatter (FindGate)
% 2. Plot histograms of the timecourse data (FACSTC_PlotHist)
% 3. Find peaks (1 or 2 per timepoint) by fitting a mixture of gaussians 
%   to the logarithm of the intensity data (FACSTimecourse)
% 4. Find Peaks (2 per timepoint) by using a peak-finding function
%   (findpeaks) on the logarithm of the intensity data (
%
%INPUT:
%This is just written as an example script, but can be modified to work for
%another timecourse by changing the variables NameFiles and Time below (see
%details below)
%This script is also not particularly intelligent about how it orders the
%files (so T18 comes after T1), which is why there is some manipulation of
%the P1 and P2 variables below to rearrange everything in the right order.
%The script will either need to be made more intelligent or the user will
%have to take care of those manipulations to modify this for another
%timecourse
%
%OUTPUT:
%The variables for the gating information (FNames and Gates) are saved to
%the variable NameFiles_GateVariables.mat
%
%A series of images are saved to the current directory:
%NameFiles_histograms.fig/.png--histograms of each timepoint all on one plot
%NameFiles_low.fig/.png--the population mean as identified by fitting 1 guassian to
%the logarith of the data
%NameFiles_high.fig/.png--population mean of the induced population as identified by
%fitting 2 gaussians
%NameFiles_peakfind.fi/.png--population peaks as identified using the
%peakfinding algorithm on the logarithm of the intensity data
%
%REQUIRES:
%   FindGate
%   FACSTC_PlotHist
%   FACSTimecourse
%   shadedLogErrorBar
%   FACSTC_PeakFind
%   FitNorms (required by FACSTimecourse
%   
%   fca_readfcs (used by several functions to read in .fcs files.
%   Available from Matlab file exchange: 
%   http://www.mathworks.com/matlabcentral/fileexchange/9608
%
%    Written by Megan McClean, Ph.D.
%               Lewis-Sigler Institute for Integrative Genomics
%               Princeton University
%               120 Carl Icahn
%               Princeton, NJ 08544
%               mmcclean@princeton.edu
%
%    Last revised on February 21, 2012




close all; clear all;

%% Some variables that can be quickly change to modify this for another timecourse
%Name variable so only have to change this to adapt this file to another
%series of data
NameFiles='165'; %Name common to all the files in the timecourse
Time=[0 1  2 3 4 5 6 7]; %Timepoints
Output=1; %Print the output (1) or don't (0) to the screen
column=3; %Need column 4 for the mCherry data, 3 for GFP
Bins=linspace(0,10,100);


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

%% 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 (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');
