%quick QC check for new batches of data
clearvars
clc

F = load('filt_neurons.mat');
filt_neurons = F.filt_neurons;

output_dir=fullfile('figs','QC');
mkdir(output_dir)

%% 
figure;imagesc(log(1+filt_neurons.expmat*10))

sgtitle('All genes')
exportgraphics(gcf,fullfile(output_dir,['Gene_expression','.pdf']),'ContentType','Image','Resolution',300);
%close all;
%%
%cells

genes_per_cell = sum(filt_neurons.expmat>0,2);
avg_genes_per_cell = mean(genes_per_cell)

counts_per_cell = sum(filt_neurons.expmat,2);
avg_counts_per_cell = mean(counts_per_cell)

figure;
hold on;
plot(genes_per_cell,counts_per_cell,'k.')
line([avg_genes_per_cell avg_genes_per_cell],[0 400],'linewidth',2)
line([0 100],[avg_counts_per_cell avg_counts_per_cell],'linewidth',2)
%prettyFig();
ylabel('counts per cell')
xlabel('genes per cell')

all_gpc =  nonzeros(genes_per_cell);
all_cpc = nonzeros(counts_per_cell);
qc = find(all_gpc>=5 & all_cpc >= 20);

qc_gpc = all_gpc(qc);
median_qc_genes_per_cell = median(qc_gpc)
qc_cpc = all_cpc(qc);
median_qc_counts_per_cell = median(qc_cpc)

fraction_cells_passed_qc = length(qc_cpc)/length(filt_neurons.slice)

sgtitle('QC cells')
exportgraphics(gcf,fullfile(output_dir,['QC_genespercell','.pdf']),'ContentType','Image','Resolution',300);