%% Explore gene/projection relationship
%% Load data

clearvars
clc
load filt_neurons-bc-somas-proj-axonal-clust_nocompr.mat

%


output_dir=fullfile('figs1','gene_proj');
mkdir(output_dir)

%
%% Look at the average projections of major subclasses


combined_proj=filt_neurons.projmat_norm(filt_neurons.has_proj>0,:);
combined_proj=combined_proj(:,1:12)+combined_proj(:,13:24);
expmat=filt_neurons.expmat(filt_neurons.has_proj>0,:);
subclass=filt_neurons.subclass(filt_neurons.has_proj>0,:);
uniq_subclass=unique(subclass);
uniq_subclass=uniq_subclass(~contains(uniq_subclass,{'Lowqual','NA','gaba'}));



mean_projmat=zeros(numel(uniq_subclass),size(combined_proj,2));
for n=1:numel(uniq_subclass)
    is_subclass=ismember(subclass,uniq_subclass(n));
    mean_projmat(n,:)=mean(combined_proj(is_subclass,:));
end

figure;
imagesc(mean_projmat);
set(gca,'ytick',1:size(mean_projmat,1), ...
    'yticklabel',uniq_subclass, ...
    'xtick',1:size(mean_projmat,2), ...
    'xticklabel',filt_neurons.proj_area(1:12), ...
    'TickLabelInterpreter','none')

%% What about individual neurons?
% Optional task:
% (1) The data look noisy, what are the potential sources of noise? And can you filter them out in analysis?
% (2) It's not obvious which population belongs to which subclass. Remake the plot so that this is more obvious.
% (3) This is a gene-centric approach: we consider transcriptomic clusters as the true grouping, and look at projection differences. Can you group by projections, and look at transcriptomic differences.

[~,idx]=ismember(subclass,uniq_subclass);
[~,I]=sort(idx,'ascend');
figure;
imagesc(combined_proj(I,:)./sum(combined_proj(I,:),2));
set(gca, ...
    'xtick',1:size(mean_projmat,2), ...
    'xticklabel',filt_neurons.proj_area(1:12), ...
    'TickLabelInterpreter','none')




