initial commit; copying from AKG private repo

Clusterless_schematic.m

@@ -0,0 +1,72 @@
+animal='despereaux';
+destdir = '/mnt/stelmo/anna/despereaux/filterframework/';
+d=6;
+e=2;
+t=3; % 2 20 
+marks = loaddatastruct(destdir, animal,'marks',[d e ]);
+pos = loaddatastruct(destdir, animal,'pos',[d e ]);
+
+linposfile = sprintf('%sdecoding/%s_%d_%d_shuffle_0_linearposition_v2.nc', destdir, animal, d, e); 
+linpos = ncread(linposfile,'linpos_flat')+1;
+linvel = ncread(linposfile,'linvel_flat');
+linposts = ncread(linposfile,'time');
+
+immobilefilter = {'ag_get2dstate', '($immobility == 1)','immobility_velocity',4,'immobility_buffer',0};
+immoperiods = kk_evaluatetimefilter(destdir,animal, {immobilefilter}, [d e]);
+immospikes = isExcluded(marks{d}{e}{t}.times, immoperiods{d}{e});    
+% correct for gaps between arms
+nodatarows = histcounts(linpos,[1:max(linpos)+1])==0;
+correction = cumsum(nodatarows);
+linposcorr = double(linpos)'-correction(linpos);
+bounds = .5+find(diff(correction(~nodatarows))>0);
+
+starttime = 3000; %
+endtime = 7500; %
+markinds = marks{d}{e}{t}.times>=starttime & marks{d}{e}{t}.times<=endtime & ~immospikes;
+
+posinds = linposts>=starttime & linposts<=endtime;
+posatmarktimes = linposcorr(lookup(marks{d}{e}{t}.times(markinds),linposts));
+%edges = [0 10 27 44 76 92 108 124 140]; 
+%segatmarktimes = discretize(posatmarktimes,edges);
+%segs = discretize(linpos,edges);
+
+figure; set(gcf,'Position',[86 179 1804 739])
+subplot(1,3,1); hold on
+subsetpos = linposts>=6850 & linposts<=7250;
+scatter(pos{d}{e}.data(subsetpos & posinds & linvel>=4,7),pos{d}{e}.data(subsetpos & posinds & linvel>=4,6),10,linposcorr(subsetpos & posinds & linvel>=4),'filled')
+set(gca,'YDir','reverse'); ylabel('x (cm; rev)'); xlabel('y (cm)');
+title(sprintf('%s d%d e%d t%d',animal,d,e,t)); 
+subplot(1,3,2);  
+scatter(linposts(subsetpos & posinds & linvel>=4)-starttime,linposcorr(subsetpos & posinds & linvel>=4),10,linposcorr(subsetpos & posinds & linvel>=4),'filled') %segs
+ylabel('linearized pos'); xlabel('time (s)'); axis square tight; colorbar
+%scatter(marks{d}{e}{t}.times(markinds),posatmarktimes,10,segatmarktimes)
+
+figure % save separately otherwise too big
+scatter(marks{d}{e}{t}.marks(markinds,1),marks{d}{e}{t}.marks(markinds,3),10,posatmarktimes,'filled'); colorbar %,'MarkerFaceAlpha',.5
+xlabel('ch2 amplitude (uV)'); ylabel('ch3 amplitude (uV)'); ylim([-200 1000]); xlim([-200 800]); 
+
+%% plot spike raster 
+tetinfo = loaddatastruct(destdir, animal,'tetinfo',[d e]);
+tets = evaluatefilter(tetinfo{d}{e},'isequal($area,''ca1'')'); %
+rips = loaddatastruct(destdir, animal,'ca1rippleskons',[d e]);
+immoevents = isExcluded(rips{d}{e}{1}.starttime, immoperiods{d}{e}) & isExcluded(rips{d}{e}{1}.endtime, immoperiods{d}{e});    
+valrips = isExcluded(rips{d}{e}{1}.starttime, immoperiods{d}{e}) & isExcluded(rips{d}{e}{1}.endtime, immoperiods{d}{e}) ;    
+riptsz = [rips{d}{e}{1}.starttime(valrips), rips{d}{e}{1}.endtime(valrips), rips{d}{e}{1}.maxthresh(valrips)];
+
+figure; hold on;
+startend = [5678.892 5682.892];
+ripinds = riptsz(:,1)>startend(1) & riptsz(:,2)<startend(2);
+for t = 1:length(tets)
+    spktimes = marks{d}{e}{tets(t)}.times(marks{d}{e}{tets(t)}.times>=startend(1) & marks{d}{e}{tets(t)}.times<=startend(2));
+    plot([spktimes'; spktimes'],repmat([t-.9;t],1,length(spktimes)),'k')
+end
+patch([riptsz(ripinds,1:2), fliplr(riptsz(ripinds,1:2))],[0 0 t t]','k','FaceAlpha',.1,'EdgeColor','none');
+axis tight; ylabel('tets'); title(sprintf('%s d%d e%d',animal, d, e)); 
+
+
+%% plot velocity
+figure;
+posinds = linposts>=startend(1) & linposts<=startend(2);
+plot(linposts(posinds),linvel(posinds),'k'); 
+patch([riptsz(ripinds,1:2), fliplr(riptsz(ripinds,1:2))],[0 0 75 75]','k','FaceAlpha',.1,'EdgeColor','none');
+axis tight; ylabel('velocity (cm/s)'); title(sprintf('%s d%d e%d',animal, d, e)); ylim([0 75]); set(gca,'ytick',[0 25 50 75])

README.md

@@ -1,2 +1,101 @@
-# Gillespie_Neuron_2021
-analysis code for Gillespie et al, Neuron 2021
+Code for doing analysis and generating figures for Gillespie et al, "Hippocampal replay represents specific past experiences rather than a plan for subsequent choice"
+
+Franklab repos needed:
+
+for data processing: trodes2ff_shared (franklab public) 
+
+for analysis: filterframework_shared (franklab private)
+
+Context:
+Data acquired 2018-2020 using trodes1.6.3 and 30-tet DKR drive version for 4 rats: Jaq, Roquefort, Despereaux, Montague
+Data was extracted using D. Liu's python extractor to run Trodes export functions, then processed and analysed in Matlab 2015b
+
+Decoding results used are "ripdecodesv3.mat" pulled from /filterframework/decoding/animal_day_ep_shuffle_0_posterior_acausalv2_full2state.nc (2-state state space model, decoding all times, 2ms bins, 5cm time bins, uniform and stepwise transmats (equal prob back/forward/stay); decoding done mostly at LLNL)
+
+Figures were saved in as pdf or eps and formatted in Adobe Illustrator
+
+Data in NWB format is available on the DANDI Archive: https://dandiarchive.org/dandiset/000115/draft
+
+saved outputs (f) included:
+
+	- dfs_ripcontent -> ctrl_ripcontent.mat
+	- dfs_ripcontent_ripspeed -> ctrl_ripcontent_ripspeed.mat
+	- dfs_ripcontent_movement -> ctrl_movementquant_full2state_all_withtrialwise.mat
+
+Code to generate figures:
+
+Figure 1: 
+
+    - B,C - plot_behavior_example.m, desp15_2 
+	- D - dfs_ripcontent.m, part 1
+    - E,F - dfs_ripcontent.m, part 2
+
+Figure 2:
+
+	- A,B - dfs_plotripcontent.m, span=full, despereaux session 8
+	- C,D,E - dfs_plotripcontent.m, span=rips, despereaux session 8, tets 25,26,30
+	- F - dfs_ripcontent.m, part 3
+
+Figure 3:
+
+	- A - dfs_ripcontent.m, part 3
+	- B,C - dfs_ripcontent.m, part 4 
+	- D,E - dfs_ripcontent.m, part 5
+
+Figure 4:
+
+	- A,B - dfs_ripcontent.m, part 5
+	- C - dfs_ripcontent.m, part 6
+
+Figure 5:
+
+	- A - dfs_ripcontent.m, part 7
+	- B,C,D - dfs_ripcontent.m, part 6
+	- E - dfs_ripcontent.m, part 8
+	- Correlations (in text, no fig) - dfs_ripcontent.m, part 16
+
+Figure 6: dfs_ripcontent.m, part 9
+
+Figure 7: dfs_ripcontent.m, part 8
+
+Figure 8:
+
+	- A - dfs_ripcontent.m, part 10
+	- B - dfs_ripcontent.m, part 5
+
+Supplementary Fig 1: 
+
+	- A - dfs_ripcontent.m, part 11
+	- C - dfs_ripcontent.m, part 12
+
+Supplementary Fig 2: clusterless_schematic.m
+
+Supplementary Fig 3:
+
+	- A-C - dfs_ripcontent_movement.m, part 2
+	- D-F - dfs_ripcontent_movement.m, part 6
+	- G - dfs_ripcontent.m, part 3
+	- H - dfs_ripcontent.m, part 13
+
+Supplementary Fig 4: dfs_plotripcontent, span=rips, run once per rat (day/ep/rip noted on fig)
+
+Supplementary Fig 5: 
+
+	- first run dfs_ripfeatures.m to calc tetwise features (slow! alternativately, load saved output)
+	- then dfs_ripcontent.m, part 14 to visualize
+
+Supplementary Fig 6: dfs_ripcontent.m, part 5
+
+Supplementary Fig 7: dfs_ripcontent.m, part 6 (commented out)
+
+Supplementary Fig 8: 
+
+	- A-D - dfs_ripcontent_ripspeed.m
+	- E-G - dfs_ripcontent.m, part 15
+
+Control Analyses (not shown in figures):
+
+	- Likelihoods (decode without state space model, 20ms bins): dfs_likcontent.m
+	- Use MUA as event detection instead of SWRs: dfs_muacontent.m
+	- Use 3-state decoder instead of 2-state to exclude stationary events: dfs_ripcontent_3state.m
+	- Use more permissive content threshold: dfs_ripcontent_contentthresh