new_LD187_preprocessing.m

% ATK 200115 new_LD187_preprocessing for use with suite2p
% Needs constrained_foopsi https://github.com/epnev/constrained-foopsi
% Needs HarveyLab helper functions (abfload) https://github.com/HarveyLab/helperFunctions
% Needs matlab signal processing toolbox

%% make sessionList_all variable for each mouse and save 
masterPath = '/n/groups/htem/temcagt/datasets/ppc/2P_data/';
mouse = 'LD187';

sessionList = dir('/n/groups/htem/temcagt/datasets/ppc/2P_data/scanimage/LD187/LD187*');
sessionList_all = cell(size(sessionList,1),1);
for s = 1:size(sessionList,1)
    sessionList_all{s,1} = sessionList(s,1).name;
end

%% Select Session and lineup

%%%%%%%%%
%session = sessionList_all{36};
session = 'LD187_141212';
disp(session)
%%%%%%%%%

%% line up 
lineUpSession(session) 
% Exports VirminCombined to ppc/2P_data/code_workspace/LD187/virmen/FILE.mat

%% Load suite2p Data
disp(['Loading suite2p data for ' session]);

% Load suite2p output
suite2pPath = fullfile(masterPath,'scanimage',mouse,session,'suite2p');
suite2pOutput = fullfile(suite2pPath,'combined/Fall.mat');
s2p = load(suite2pOutput);
numCells = size(s2p.F,1);
neucoeff = 0.7;
s2p.Fsub = s2p.F - neucoeff * s2p.Fneu;
disp(['Suite2p video data has ' num2str(size(s2p.F,2)) ' frames']);

%% Deconvolution
disp(session)
tic
[c,b,c1,g,sn,sp] = run_constrained_foopsi(double(s2p.Fsub));
toc

output_dir_deconv = fullfile(masterPath, 'code_workspace',mouse,'deconvData');   
if ~exist(output_dir_deconv,'dir')
    mkdir(output_dir_deconv)
end
save(fullfile(output_dir_deconv,session),'c','b','c1','g','sn','sp');

%% ATK 201030 for loading from ext HD
%%%%%%%%%%%%%%%%%%%%%%%%%
% From this point on you don't need access to /n/groups/htem
masterPath = "/Volumes/Aaron's PPC/ppc/2P_data/";
%masterPath = "/media/atk13/Aaron's 5TB HD/ppc/2P_data/";
mouse = 'LD187';
session = 'LD187_141206';
disp(session)

% Load aligned virmen data
linedUpPath = fullfile(masterPath,'code_workspace',mouse,'virmen',[session '.mat']);
vData = load(linedUpPath);
disp(['Virmen Data has ' num2str(size(vData.VirmenCombined,2)) ' frames']);

% Load deconv data
disp(['Loading deconv data for ' session]);
%disp(session)
output_dir_deconv = fullfile(masterPath, 'code_workspace',mouse,'deconvData');   
load(fullfile(output_dir_deconv,session),'c','b','c1','g','sn','sp');

% Calculate SNR metrics
disp(['Calculating SNR for ' session]);
snr = nan([],length(sn));
for i = 1:length(sn)
    snr(i) = std(c(i,:))^2/sn(i)^2;
end

% Parse virmen trials
disp(['Saving trial aligned data for ' session]);
% Hack for 141210 %% trialAlignedData = parseVirmenTrials(vData.VirmenCombined(:,6001:16000), sp);
trialAlignedData = parseVirmenTrials(vData.VirmenCombined, sp);
output_dir = fullfile(masterPath, 'code_workspace',mouse,'syncedData');   
if ~exist(output_dir,'dir')
    mkdir(output_dir)
end
trialAlignedData.SNR_raw = snr;
save(fullfile(output_dir,session),'trialAlignedData','-v7.3');

%% From here down is obselete %%%%%%%%%%%%%%%%%
%% Calc Mutual Information with trial type (not used but dev'd here)
% requires https://github.com/nmtimme/Neuroscience-Information-Theory-Toolbox
% ATK 210616
syncedDataDir = fullfile(masterPath, 'code_workspace',mouse,'syncedData');   
syncedData = open(fullfile(syncedDataDir,[session '.mat']));
DataRaster = cell(1);
% Format trialAlignedData for info theory toolbox
trialTypes = syncedData.trialAlignedData.trialType;
trialTypes(1,1) = nan; %exclude first trial (has nans in ITI) 

% Only include trial types 2 and 3 (R and L)
% 1 - Ca Data
DataRaster{1} = permute(syncedData.trialAlignedData.CaData(:,ismember(trialTypes,[2,3]),:),[1,3,2]);
% 2 - trial types
DataRaster{2}(1,:,:) = repmat(trialTypes(:,ismember(trialTypes, [2,3]))', [1,76])';

% Binary binning (active or not)
numBins=2;
numROIs = length(DataRaster{1});
StatesRaster{1} = DataRaster{1}>0;
StatesRaster{2} = DataRaster{2};
% Uniform Binning
%{
numBins = 4; % bin into 20 uniform width bins
numROIs = length(DataRaster{1});
MethodAssign = cell(numROIs, 4);
for i = 1:numROIs
    MethodAssign(i,:) = {1, i,'UniWB',{numBins}}; 
end
StatesRaster = data2states(DataRaster, MethodAssign);
%}
tic 
% Calculate MI with trial type
disp('Calculating MI with trial type')
choiceMI = nan(numROIs, 76);
for id = 1:numROIs
    for t = 1:76
        VariableIDs = {1,id,t;2,1,1};
        choiceMI(id, t) = instinfo(StatesRaster, 'PairMI', VariableIDs);
    end
end
[choiceMI_max, choiceMI_max_idx] = max(choiceMI, [], 2);

% Calc significance of maxMI by shuffling trial types
disp('Calculating shuffle distributions')
numShuf = 500;
numTrials = size(StatesRaster{2},3);
choiceMI_thresh = nan(numROIs,1);
for id = 1:numROIs
    disp(['Calc shuf for ROI ' num2str(id)])
    choiceMI_shuf = nan(numShuf,1);
    for shuf_idx = 1:numShuf
        % shuffle trial types vector
        StatesRasterShuf = StatesRaster;
        StatesRasterShuf{2} = StatesRasterShuf{2}(:,:,randperm(numTrials));
        t = choiceMI_max_idx(id);
        choiceMI_shuf(shuf_idx) = instinfo(StatesRasterShuf,'PairMI', {1,id,t;2,1,1});            
    end
    choiceMI_thresh(id) = prctile(choiceMI_shuf, 99.94);
end
maxMI_sig = choiceMI_max > choiceMI_thresh;
toc
% Plot significant peaks as function of timepoint
frac_sig = nan(76,1);
for i = 1:76
    frac_sig(i) = sum(maxMI_sig(find(choiceMI_max_idx==i)))/sum(choiceMI_max_idx==i);
    num_sig(i) = sum(maxMI_sig(find(choiceMI_max_idx==i)));
end
disp([num2str(sum(maxMI_sig>0)/numROIs) ' of ROIs are sig'])
figure; plot(num_sig); xlabel('trial timepoint'); ylabel('# ROIs with sig MI with trial type');title(session);
xline(13.5,'k');xline(26.5,'k');xline(38.5,'k--');xline(51.5,'k');xline(64.5,'k--');
figure; plot(frac_sig); xlabel('trial timepoint'); ylabel('frac of ROIs with sig MI with trial type');title(session);
xline(13.5,'k');xline(26.5,'k');xline(38.5,'k--');xline(51.5,'k');xline(64.5,'k--');

%% Example MI with trial type
Method = 'PairMI';
MI = zeros(76,1);
id = 6; %101; %44; %2

for t = 1:76
    VariableIDs = {1,id,t;2,1,1};
    MI(t) = instinfo(StatesRaster, Method, VariableIDs);
end
figure; hold on;set(gca,'fontsize',16)
%plot(1:76,nanmean(DataRaster{1}(id,1:76,DataRaster{2}(1,50,:)==2),3),'LineWidth',2);
%plot(1:76,nanmean(DataRaster{1}(id,1:76,DataRaster{2}(1,50,:)==3),3),'LineWidth',2);
plot(1:76,nanmean(StatesRaster{1}(id,1:76,StatesRaster{2}(1,50,:)==2),3),'LineWidth',2);
plot(1:76,nanmean(StatesRaster{1}(id,1:76,StatesRaster{2}(1,50,:)==3),3),'LineWidth',2);
ylabel('Spike Rate (arb)')
yyaxis right
ylabel('Mutual Information (bits)')
plot(1:76,MI(1:76),'LineWidth',2); 
legend({'Left','Right','MI'})
title(['Cell ID: ' num2str(id) '; numBins: ' num2str(numBins)])

%% Example pairwise MI
Method = 'PairMI';
MI = zeros(t,1);
id1 = 88; %44; %2
id2 = 3;%15 ; %12

for t = 1:76
    VariableIDs = {id1,t;id2,t};
    MI(t) = instinfo(StatesRaster, Method, VariableIDs);
end
figure; hold on;set(gca,'fontsize',16)
plot(1:76,nanmean(DataRaster{1}(id1,1:76,:),3),'LineWidth',2);
plot(1:76,nanmean(DataRaster{1}(id2,1:76,:),3),'LineWidth',2);
ylabel('Spike Rate (arb)')
yyaxis right
ylabel('Mutual Information (bits)')
plot(1:76,MI(1:76),'LineWidth',2); 
legend({'id1','id2','MI'})

%% Trial average MI
Method = 'PairMI';
n_rois = 100;
tic
MI_avg = nan(n_rois,n_rois);
MI_max = nan(n_rois,n_rois);
for i_src = 1:n_rois
    for j_tgt = 1:n_rois
        for t = 1:76
            MI = zeros(t,1);
            VariableIDs = {i_src,t;j_tgt,t};
            MI(t) = instinfo(StatesRaster, Method, VariableIDs);
            MI_avg(i_src,j_tgt) = mean(MI);
            MI_max(i_src,j_tgt) = max(MI);
        end
    end
    MI_avg(i_src,i_src) = nan;
    MI_max(i_src,i_src) = nan;
end
toc

figure; hold on; set(gca,'fontsize',16); 
histogram(MI_max(:),50)
xlabel('MI trial average (bits)');
ylabel('# neuron pairs');
%% Compare with linear corr

figure; hold on; set(gca, 'fontsize',16);
corr_trial_avg = trialAlignedData.corr_Ca_concat(1:n_rois,1:n_rois);
plot(corr_trial_avg(:), MI_max(:), '.');
xlabel('linear correlation');
ylabel('trial avg mutual information (bits)');

%% Calc Mutual Information - single trial
% ATK 210616
output_dir = fullfile(masterPath, 'code_workspace',mouse,'syncedData');   
open(fullfile(output_dir,[session '.mat']));
DataRaster = cell(1);
DataRaster{1}(:,:,1) = trialAlignedData.Ca_concat;
%exclude first trial that has NaNs
MethodAssign = cell(size(DataRaster{1},1), 3);
for i = 1:size(DataRaster{1},1)
    % bin into 20 uniform width bins
    MethodAssign(i,:) = {i,'UniWB',{20}}; 
end
[StatesRaster,MethodResults] = data2states(DataRaster, MethodAssign);

%% Example pairwise MI
Method = 'PairMI';
MI = zeros(t,1);
id1 = 76;%88; %44; %2
id2 = 1;%3;%15 ; %12

VariableIDs = {id1,1;id2,1}; % only relative timing important here
MI = instinfo(StatesRaster, Method, VariableIDs);
figure; hold on;set(gca,'fontsize',16)
plot(DataRaster{1}(id1,:),'LineWidth',2);
plot(DataRaster{1}(id2,:),'LineWidth',2);
%% MI array
Method = 'PairMI';
n_rois = 100;
tic
MI_st = nan(n_rois,n_rois);
for i_src = 1:n_rois
    for j_tgt = 1:n_rois
        VariableIDs = {i_src,1;j_tgt,1};
        MI_st(i_src, j_tgt) = instinfo(StatesRaster, Method, VariableIDs);
    end
    MI_st(i_src,i_src) = nan;
end
   
toc

figure; hold on; set(gca,'fontsize',16); 
histogram(MI_st(:),50)
xlabel('MI trial average (bits)');
ylabel('# neuron pairs');

%% Compare with linear corr

figure; hold on; set(gca, 'fontsize',16);
corr_trial_avg = trialAlignedData.corr_Ca_concat(1:n_rois,1:n_rois);
plot(corr_trial_avg(:), MI_st(:), '.');

xlabel('linear correlation');
ylabel('single trial mutual information (bits)');

%% Plot example traces
figure; hold on;
frames = 1:length(sp);
dt = 1/5.3; % 5.3 Hz sampling rate
time = frames * dt;
for i = 1:30%size(mySignal,1)
    plot(time,sp(i,:)+100*i);
    
end
xlabel('time (sec)');
ylabel('spike rate');

%% Plot example activity as test

s2p_cid = 8;
matlab_cid = s2p_cid+1;

% Define colors
cueEarlyColor = [0 161 75]/255; 
cueLateColor = [255 222 23]/255;
delayEarlyColor = [237 28 36]/255;
delayTurnColor = [127 63 152]/255;
turnITIcolor = [33 64 154]/255;
% Define blocks
cueBlockEarly = 14:26; % 14 is running onset + 12 frames after
cueBlockLate = 27:38; % 12 frames before cue offset (frame 39)
delayBlockEarly = 39:51; % 39 is cue offset + 12 frames after
delayTurnBlock = 52:64; % 12 frames before end of trial (turn a certain amt)
turnITIblock = 65:76; % Trial end (reward given) an 12 frames after (dark, ITI)

figure;
tiledlayout(1,2)

ax1 = nexttile;
%subplot(1,2,1); 
hold on;
plot(cueBlockEarly,trialAlignedData.bR_trials.Ca_trialMean(matlab_cid,cueBlockEarly),'color',cueEarlyColor);
plot(cueBlockLate,trialAlignedData.bR_trials.Ca_trialMean(matlab_cid,cueBlockLate),'color',cueLateColor);
plot(delayBlockEarly,trialAlignedData.bR_trials.Ca_trialMean(matlab_cid,delayBlockEarly),'color',delayEarlyColor);
plot(delayTurnBlock,trialAlignedData.bR_trials.Ca_trialMean(matlab_cid,delayTurnBlock),'color',delayTurnColor);
plot(turnITIblock,trialAlignedData.bR_trials.Ca_trialMean(matlab_cid,turnITIblock),'color',turnITIcolor);
title('right');

ax2 = nexttile;
%subplot(1,2,2); 
hold on;
plot(cueBlockLate,trialAlignedData.wL_trials.Ca_trialMean(matlab_cid,cueBlockLate),'color',cueLateColor);
plot(delayBlockEarly,trialAlignedData.wL_trials.Ca_trialMean(matlab_cid,delayBlockEarly),'color',delayEarlyColor);
plot(delayTurnBlock,trialAlignedData.wL_trials.Ca_trialMean(matlab_cid,delayTurnBlock),'color',delayTurnColor);
plot(turnITIblock,trialAlignedData.wL_trials.Ca_trialMean(matlab_cid,turnITIblock),'color',turnITIcolor);
title('left');
linkaxes([ax1 ax2],'xy')

%% compare suite2p and constrained_foopsi deconv
cell_idx = find(s2p.iscell(:,1)==1);

figure;
num_ex = 3;
start = 100;
num_pts = 400;
tsec = (1:num_pts)/5.6;
for i = start:start+num_ex-1
    disp(cell_idx(i))
    subplot(num_ex,3,3*(i-start)+1);
    plot(tsec,s2p.Fsub(cell_idx(i),1:num_pts));
    title('raw');
    subplot(num_ex,3,3*(i-start) + 2);
    plot(tsec,s2p.spks(cell_idx(i),1:num_pts));
    title('suite2p');
    subplot(num_ex,3,3*(i-start)+3);    
    plot(tsec,sp(cell_idx(i),1:num_pts));
    title('const foopsi');
end


%% Plot SNR 
figure; plot(snr); xlabel('mask ID'); ylabel('SNR');
figure; histogram(snr, 100); xlabel('SNR'); ylabel('# of masks');  xlim([0,150]);
%% Plot example data
cid =674;
dF = s2p.Fsub;
sorted_dF = sort(dF(cid,:));
baseline_adj = mode(round(100*sorted_dF(.05*size(sorted_dF,2):...
        .95*size(sorted_dF,2)))/100);
dF_zeroed = dF(cid,:) - baseline_adj;

figure; hold on; 
plot(dF_zeroed);
%plot(c(cid,:));
title(['SNR: '  num2str(snr(cid))])