%% Explore projection data
%% Load data

clearvars
clc
load filt_neurons-bc-somas-proj-axonal-clust_nocompr.mat
output_dir=fullfile('figs1','projmat');
mkdir(output_dir);

%
%% Visualize projection matrix
% When doing quantitative analysis, we usually use the spike-in normalized projection matrix (projmat_norm)
% First plot all neurons with projections

%%
close all;
figure('Position',[50 50 400 250]);
imagesc(filt_neurons.projmat_norm(filt_neurons.has_proj>0,:)./ ...
    sum(filt_neurons.projmat_norm(filt_neurons.has_proj>0,:),2));
set(gca, ...
    'xtick',1:numel(filt_neurons.proj_area),...
    'xticklabel',filt_neurons.proj_area);
colorbar;
ylabel('Neurons')
xlabel('Areas');
title('Spike-in corrected barcode count, normalized by row')
exportgraphics(gcf,fullfile(output_dir,'rownorm_mat.pdf'))

% You can also visualize raw barcode counts, which is more useful when binarizing data (does a cell project to an area or not?)

%%
figure('Position',[50 50 400 250]);
imagesc(log10(1+filt_neurons.projmat(filt_neurons.has_proj>0,:)));colorbar;
set(gca, ...
    'xtick',1:numel(filt_neurons.proj_area),...
    'xticklabel',filt_neurons.proj_area);
colorbar;
ylabel('Neurons')
xlabel('Areas');
title('Log10(1+Raw barcode count)')
exportgraphics(gcf,fullfile(output_dir,'raw_mat.pdf'))

% Alternatively, we can also work with log scaled spike-in corrected projection data, without row normalization. This is useful sometimes for picking up pouplations with weak projections, because the log transform semi-binarizes the data.


figure('Position',[50 50 400 250]);
imagesc(log10(1+1000*filt_neurons.projmat_norm(filt_neurons.has_proj>0,:)));colorbar;
set(gca, ...
    'xtick',1:numel(filt_neurons.proj_area),...
    'xticklabel',filt_neurons.proj_area);
colorbar;
ylabel('Neurons')
xlabel('Areas');
title('Log10(1+1000*corrected barcode count)')
exportgraphics(gcf,fullfile(output_dir,'log_mat.pdf'))

%
%% What's the distribution of barcode count per area?
% Plot histograms for each area separately.

proj_thresh=1; % try changing this.
for n=1:numel(filt_neurons.proj_area)
    figure('Position',[50 50 300 250]);
    has_proj=filt_neurons.projmat(:,n)>=proj_thresh;
    histogram(filt_neurons.projmat(has_proj,n),0:5:100); %try changing the bins
    ylabel('# of neurons')
    xlabel('max bc count per neuron')
    title(sprintf('%u_%s',n,filt_neurons.proj_area{n}), ...
        'Interpreter','none');
    exportgraphics(gcf,fullfile(output_dir,sprintf('%u_%s.pdf',n,filt_neurons.proj_area{n})));
end

%
%% Clustering projections
% In many cases, it's convenient to explore the data by clustering neurons into types. This is frequently done in transcriptomics, and similarly it can be done for projections. Many method exists for this. We'll try a simple k-means clustering on the row-normalized data, but other methods usually work better.
% Optional task: Try playing with the number of clusters (k), other clustering methods (hierarchical, louvaine, etc.), or clustering after other processing (e.g. PCA on projmat_norm);

% how many clusters do we want?

k=10;

rng(10);% for reproducibility
norm_proj=filt_neurons.projmat(filt_neurons.has_proj>0,:);
norm_proj=norm_proj./sum(norm_proj,2);

% clustering
idx=kmeans(norm_proj,k);

% sort cells by cluster identity, then plot them
[~,I]=sort(idx,'ascend');
figure('Position',[50 50 400 250]);
imagesc(norm_proj(I,:));colorbar;
set(gca, ...
    'xtick',1:numel(filt_neurons.proj_area),...
    'xticklabel',filt_neurons.proj_area);
colorbar;
ylabel('Neurons')
xlabel('Areas');
title('rownorm counts, sorted by cluster')
exportgraphics(gcf,fullfile(output_dir,'clutsered.pdf'))


%
% Alternatively, think of neurons as sum of parts that are composed of sets of projections. We can find these "proto-projections" by nmf. This is similar to looking for co-expression in transcriptomics.
% Optional task: try other dimensionality reduction methods, PCA, sparse-nmf, etc. Or try changing the number of components in nmf.

rng(100)
k=6;

% because the projection matrix contains data from two brains with similar
% areas, first combine them
combined_proj=filt_neurons.projmat_norm(filt_neurons.has_proj>0,:);
combined_proj=combined_proj(:,1:12)+combined_proj(:,13:24);




rng(10);% for reproducibility

[W,H] = nnmf(combined_proj,k);
figure;
imagesc(H);
set(gca, ...
    'xtick',1:numel(filt_neurons.proj_area),...
    'xticklabel',filt_neurons.proj_area);
xlabel('Areas');
ylabel('Components');

% Can also do this on the log scaled data. Sometimes, it is easier to pick up sets of projections because it is more binarized. In this dataset, doesn't make a big difference.

k=6
combined_proj=log10(1+1000*filt_neurons.projmat_norm(filt_neurons.has_proj>0,:));
combined_proj=combined_proj(:,1:12)+combined_proj(:,13:24);


rng(10);% for reproducibility

[W,H] = nnmf(combined_proj,k);
figure;
imagesc(H);
set(gca, ...
    'xtick',1:numel(filt_neurons.proj_area),...
    'xticklabel',filt_neurons.proj_area);
xlabel('Areas');
ylabel('Components');