% get the names of all fcs files in the current folder
all_filenames = getfilenames;
all_filenames = all_filenames(isInListEnd(all_filenames,'.fcs'));

% give each file a short name (the following two lines needs to be customed for your own set of FCS files)
file_annot = cell(length(all_filenames));
file_annot = [{'file1'}];  
% if you have multiple FCS files in this folder, this looks like:
%       file_annot = [{'file1'},{'file2'},{'file3'}] 
% NOTE: length of "file_annot" variable needs to be consistent with the length of "all_filenames", 


% % define markers used for both downsampling and tree construction
% % if the fcs files are generated by overlapping panel, must use their names (the names have to be exact).
% % if all fcs files share the same panel, can use channel index
used_markers = [{'marker1'}, ...
                {'marker2'}]';


% define data transformation and downsampling parameters
arcsinh_cofactor=5;                 % use arcsinh to transform the data. If you don't want to perform this transformation, set the value to be zero
is_normalize=0;                     % whether or not normalize each channel to zero-mean unit-variance. 0 means no normalization, 1 means normalize after arcsinh transformation. 
downsampling_scaling_factor=5;      % this parameter is used for estimating the local density of each cell. SPADE is not sensitive to this parameter. You can set it to be anywhere between 5 and 50
exclude_prctile = 1;                % exclude the lower 1 prctile sparse cells, please refer to the method section of the manuscript
target_prctile = 20000;             % if this variable is <100, it is the target perctile (refer to method section in paper). Sometimes, it is desirable to define the desired number of cells that survive downsampling. If this variable is > 100, this is the number of remaining cells after downsampling each file. 
num_cells = 50000;                  % max number of cells that feed into the clustering step. This is the max number of cells that the SPADE algorithm is comfortable with, please leave it as is for now. 
num_target_clusters = 100;          % number of desired cell clusters. If the value of this parameter is x, the number of bubbles in the final tree will be between x-x/2 and x+x/2, that is 100~300.
FileInds_used_in_SPADE = [1];       % Define which files are used to construct tree. If you set it to be empty or [1 2 ... number_of_files], all fcs files will be used to construct the tree
% (we used this parameter setting for processing the 30 fcs files from the human bone marrow data discussed in the paper)

% you can ignore this part 
PooledDownsampledDataFilename='PooledData.mat';
FlowSPADE_input_filename = 'DataForSPADE.mat';
FlowSPADE_output_filename = getfilenames; 
ind = setdiff(isInList(FlowSPADE_output_filename,[FlowSPADE_input_filename(1:end-4),'_']),isInList(FlowSPADE_output_filename,['_node_median'])); 
if ~isempty(ind)
    if length(ind)==1
        FlowSPADE_output_filename = FlowSPADE_output_filename{ind}; 
    else
        FlowSPADE_output_filename = FlowSPADE_output_filename{ind(2)}; 
    end
else 
    FlowSPADE_output_filename=[]; 
end
FlowSPADE_output_upsample_filename = getfilenames; 
ind = intersect(isInList(FlowSPADE_output_upsample_filename,[FlowSPADE_input_filename(1:end-4),'_']),isInList(FlowSPADE_output_upsample_filename,['_node_median'])); 
if ~isempty(ind), FlowSPADE_output_upsample_filename = FlowSPADE_output_upsample_filename{ind}; else, FlowSPADE_output_upsample_filename=[]; end


% % for each fcs file, compute the local density of each cell
% % save the data and the local density in an intermediate .mat file
% % one intermediate .mat file is created for each fcs file
% % (20 min per file if use parallel computing toolbox; if not, ~60 mins)
for i=1:length(all_filenames)
    filename = all_filenames{i};
    if exist([all_filenames{i}(1:end-4),'.mat'])==2
        continue;
    end
    tic 
    % % if you have parallel computing toolbox installed, use line 60, which is faster; otherwise, use line 61
    % parfor_1_read_individual_fcs_get_downsampling_info(filename,arcsinh_cofactor,downsampling_scaling_factor,used_markers,is_normalize);
    for_1_read_individual_fcs_get_downsampling_info(filename,arcsinh_cofactor,downsampling_scaling_factor,used_markers,is_normalize);
    toc
end


% % This script actually does the downsampling.
% % each fcs file is downsampled separately, 
% % all the downsampled cells are pooled into a single file
% % this step generates "PooledData.mat"
% % although named as parfor_..., this function does not need parallel computing toolbox
parfor_2_prepare_PooledDownsampledData(all_filenames, target_prctile, exclude_prctile,PooledDownsampledDataFilename)


% % 
% % prepare input file for the clustering step
% % if the total number of pooled cells is too large, we uniformly downsample the data further, because the current clustering algorithm cannot handle too many cells without memory problem or running forever. We limit the total number to be 50000 
% % although named as parfor_..., this function does not need parallel computing toolbox
parfor_3_Prepare_DataForSPADE(num_cells,PooledDownsampledDataFilename,FlowSPADE_input_filename,FileInds_used_in_SPADE)


% % clustering and MST construction
% % this step generate a .mat file that has a very long 
% % although named as parfor_..., this function does not need parallel computing toolbox
[FlowSPADE_output_filename] = parfor_4_Process_DataForSPADE_replace_GUI(FlowSPADE_input_filename,num_target_clusters);


% % this function draws pictures and save them into a subfolder
load(FlowSPADE_input_filename,'marker_names','used_markers')
draw_marker_ind = 1:length(marker_names);
DrawTrees_SPADEoutputFile(FlowSPADE_output_filename,draw_marker_ind, 1, 0, 'JET');
% % the above line generates jpg files
% % other functions for drawing results are as follows:
% %     DrawTrees_SPADEoutputFile_NoColorEPS    % generate eps, no colors
% %     DrawTrees_SPADEoutputFileEPS            % generate eps
% %     DrawTrees_SPADEoutputFilePDF            % generate pdf


% % upsampling (<20 min per file)
% % if you have parallel computing toolbox installed, use the following commented line. Otherwise, use the line started with for_...
% % This file generates the .mat file with the longest name. 
% [FlowSPADE_output_upsample_filename] = parfor_5_upsampling_v2(FlowSPADE_output_filename,all_filenames);
[FlowSPADE_output_upsample_filename] = for_5_upsampling_v2(FlowSPADE_output_filename,all_filenames);


% % After upsampling is performed, additional functions for drawing trees can be used:
% %     DrawTrees_SPADEoutputUpsample           % generate jpg, draw one tree for one marker and one fcs file. Color of a node is defined by the average intensity of the selected marker across the cells in the selected fcs file which also belongs to this node.
% %     DrawTrees_TwoSample_CellFreqDiff        % generate jpg, used when comparing two fcs files that corresponds to two tissues. Draw one tree for one marker and two files. Color of a node is defined by the difference between percentages of cells in this node for the two files.
% %     DrawTrees_UnstimStimRatio               % generate jpg, used when comparing two fcs files, stimulated vs unstimulated, as shown in our paper. Draw one tree for one marker and two files. Color of a node is defined by the ratio between the average intensities of the selected marker in the two selected files. 
% %     DrawTrees_UnstimStimRatioEPS            % generate eps, function is same as the above line.
% %     DrawTrees_UnstimStimRatioPDF            % generate pdf, function is same as the above line.
