%% Explore gene expression
%% Load data
% This data has been matched to projection barcodes, and the cells are clustered to get cell types. These two processes are not covered in this course due to time constraints. For more information, check out https://doi.org/10.1038/s41586-024-07221-6 for clustering and https://doi.org/10.1038/s41593-021-00842-4 for matching projections.
% We normally also register to CCF before doing any spatial analysis. This is not covered in this notebook. Check https://doi.org/10.1038/s41586-024-07221-6 for more information. Currently this is done manually using QuickNII and Visualign (https://quicknii.readthedocs.io/en/latest/Introduction.html), and is time consuming. Automation of this process is underway https://doi.org/10.1038/s41467-024-47883-4 .

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

%
%% Make folders for saving figures
% (Always save figures instead of relying on notebooks, which are not reproducible by themselves)

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

%
%% Plot all cells on umap in transcriptomic space
% Umap is commonly used to "visualize" gene expression space. Not useful for quantitative analysis because it distorts gene expression space in an opaque way, but can be used to get a sense of data initially.
%

% read umap coordinates. These were pre-calculated during clustering so
% that you don't need to recalculate everytime
c=table2cell(readtable('analysis/whole/umap.csv','NumHeaderLines',0));
umap_coor=cell2mat(c);
% find all the cells that passed QC. cells that didn't pass are labeled as
% "filtered" in clusters.
qc=~(contains(filt_neurons.clustid,'filtered'));

is_barcoded = filt_neurons.is_barcoded == 1;

% get all the subclass names
uniq_types=unique(filt_neurons.subclass(qc));
[~,clustidx]=ismember(filt_neurons.subclass,uniq_types);
filt_neurons.clustidx = clustidx;
uniq_clustidx=unique(clustidx);

% get distinguishable colors for plotting
close all;
cmap=distinguishable_colors(numel(uniq_types));
rng(10);
cmap=cmap(randperm(size(cmap,1)),:); % Try changing the rng value to randomize color sequence in different ways to get better visualization
map1=colormap('lines');

uniq_clustidx=unique(clustidx);
uniq_clustidx(uniq_clustidx==0)=[];

% Plot all cells on a umap, colored by subclass
figure('Position',[50 50 800 400]);
subplot(1,2,1)
hold on;
% first plot all the cells, colored by clustidx
scatter(umap_coor(:,1),...
    umap_coor(:,2),...
    3,... % dot size, try playing with this to get better visualization
    cmap(clustidx(qc),:),...
    'filled',...
    'MarkerFaceAlpha',0.1,... % try playing with this to get better visualization
    'MarkerEdgeAlpha',0.1 ... % try playing with this to get better visualization
    );
title('Clusters')
xl=xlim;yl=ylim;pbaspect([range(xl),range(yl),1]);
axis off
subplot(1,2,2)

hold on;
% this part plots "empty" data in the corresponding colors. TThis is a
% shortcut for building the proper legend next
for i=1:numel(uniq_clustidx)
    scatter([],...
        [],...
        3,...
        cmap(i,:),...
        'filled',...
        'MarkerFaceAlpha',1,...
        'MarkerEdgeAlpha',1 ...
        );
end
% add color legends
legend([uniq_types],...
    'Location','eastoutside',... 
    'NumColumns',2,...
    'Interpreter','none')
axis off
% Always save your plots. Notebooks are not reproducible.
exportgraphics(gcf,...
    fullfile(output_dir,['cluster_whole_umap','.tif']),...
    'ContentType','Image',...
    'Resolution', 300 ...
    );

%
%% Plot barcoded cells on the same plot.
% A common concern about Sindbis is toxicity, so we usually start by checking whether barcoded cells still have sufficient gene expression. This should usually be done more quantitatively, but a good starting point is to plot them on the same umap plots to see (1) if they are biased to certain cell types, and (2) if they are low quality (usually this appears as scattered cells outside of the main clusters with low gene counts.

% Where are the barcoded cells?
%% all clusters grey

% get some colors for plotting
cmap1=colormap('lines');
close all
figure('Position',[50 50 800 400]);
subplot(1,2,1)
hold on;
% Plot all the not is_barcoded cells in gray([0.7 0.7 0.7] in rgb).
scatter(umap_coor(~is_barcoded(qc),1),...
    umap_coor(~is_barcoded(qc),2),...
    3,...
    [0.7,0.7,0.7],...
    'filled',...
    'MarkerFaceAlpha',1,...
    'MarkerEdgeAlpha',1 ...
    );
% Add onto barcoded cells in a different color
scatter(umap_coor(is_barcoded(qc),1),...
    umap_coor(is_barcoded(qc),2),...
    3,...
    cmap1(1,:),...
    'filled',...
    'MarkerFaceAlpha',1,...
    'MarkerEdgeAlpha',1 ...
    );
% make x and y the same scale
xl=xlim;yl=ylim;pbaspect([range(xl) range(yl),1]);
axis off
title('BC')
%barcoded positive colored

% make a second plot with no actual plot, but just the legend with the
% proper colors.
subplot(1,2,2)
hold on;
scatter([],...
    [],...
    1,...
    [0.7,0.7,0.7],...
    'filled',...
    'MarkerFaceAlpha',1,...
    'MarkerEdgeAlpha',1 ...
    );
scatter([],...
    [],...
    1,...
    cmap1(1,:),...
    'filled',...
    'MarkerFaceAlpha',1,...
    'MarkerEdgeAlpha',1 ...
    );
legend({'Non-BC','BC'},...
    'Location','eastoutside',...
    'NumColumns',1,...
    'Interpreter','none')
xl=xlim;yl=ylim;pbaspect([range(xl) range(yl) 1])
axis off
    exportgraphics(gcf,fullfile(output_dir,['BC_whole_umap','.tif']),...
        'ContentType','Image',...
        'Resolution', 300 ...
        );

% Based on this, barcoded cells are all scatter across the UMAP, which is consistent with a more or less unbiased labeling by Sindbis.
% This also tells us that there is sufficient transcriptomic information in teh barcoded cells.

    

%

