123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223 |
- %decodes trial identity from spiking across bins for each individual neuron
- clear all
- load('RAWpsvpALL.mat')
- load('RAWvpppsALL.mat')
- load('RAWvppaALL.mat')
- load('RAWvpppaALL.mat')
- if ~exist('RPPSall'),load('RPPSall.mat'); end
- if ~exist('RPSall'),load('RPSall.mat'); end
- if ~exist('RvppaALL'),load('RvppaALL.mat'); end
- if ~exist('RvpppaALL'),load('RvpppaALL.mat'); end
- for PavInst=1; %DS is 1, Pav is 2
- if PavInst==1
- RAW=RAWvpppsALL;
- selection=RPPSall.Structure==10 & RPPSall.CSPlusRatio>=.5 & (RPPSall.CSPlusRatio./(RPPSall.CSPlusRatio+RPPSall.CSMinusRatio))>=.7;
- else if PavInst==2
- RAW=RAWpsvpALL;
- selection=RPSall.Structure==10 & RPSall.CSPlusRatio>=.5 & (RPSall.CSPlusRatio./(RPSall.CSPlusRatio+RPSall.CSMinusRatio))>=.7;
- else if PavInst==3
- RAW=RAWvpppaALL;
- selection=RvpppaALL.Structure==10 & RvpppaALL.CSPlusRatio>=.5 & (RvpppaALL.CSPlusRatio./(RvpppaALL.CSPlusRatio+RvpppaALL.CSMinusRatio))>=.7;
- else if PavInst==4
- RAW=RAWvppaALL;
- selection=RvppaALL.Structure==10 & RvppaALL.CSPlusRatio>=.5 & (RvppaALL.CSPlusRatio./(RvppaALL.CSPlusRatio+RvppaALL.CSMinusRatio))>=.7;
- end
- end
- end
- end
-
-
-
-
-
- TotalNeurons = 0;
-
-
- NumNeurons=[1 5 10 25 50 100]; %matrix of how many neurons to use on each iteration
- repetitions=20; %how many times to run the analysis
- folds = 5; %number of times cross-validated:
- Dura=[0 0.3]; %size of bin used for decoding:
- bins = 1; %number of bins
- binint = 0.3; %spacing of bins
- binstart = 0; %start time of first bin relative to event
- shuffs = 1; %number of shuffled models created
-
-
-
- for i=1:length(NumNeurons)
- PoolMdlAcc{i,1}=NaN(folds,bins);
- PoolMdlAccShuff{i,1}=NaN(folds,bins);
- end
-
- %total number of neurons in all sessions and events per session
- %for i=1:length(RAW)
- TotalNeurons=sum(selection);
- % Ev1perSess(i)=length(RAW(i).Erast{Ev1,1});
- % Ev2perSess(i)=length(RAW(i).Erast{Ev2,1});
- %end
-
- %figure out how many trials of each event can be used
- % Ev1s=min(Ev1perSess);
- % Ev2s=min(Ev2perSess);
-
- %Need to keep number of trials used in VP and NAc constant, so have to pick
- %minimum number across all sessions in both regions. Right now with
- %existing sessions it's 20 Ev1s and 20 Ev2s
- Ev1s=27;
- Ev2s=27;
-
- for v = 1:length(NumNeurons)
- for u=1:repetitions
- %pick which neurons to use
- SetupSel=cat(1,ones(NumNeurons(v),1),zeros(TotalNeurons-NumNeurons(v),1));
- Sel=(SetupSel(randperm(length(SetupSel)))==1);
-
- %setup the event identity matrix for decoding
- DecodeRs(1:Ev1s,1)=1;
- DecodeRs((Ev1s+1):(Ev1s+Ev2s),1)=2;
-
- for l=1:bins
- DecodeSpikes=NaN(1,1);
- AllSpikes=NaN(1,1);
- LowHz=zeros(1,NumNeurons(v));
- NN = 1;
- AllNN=1;
- for i=1:length(RAW) %loops through sessions
- for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
- if selection(AllNN)==1
- Ev1Spikes=NaN(1,1);
- Ev2Spikes=NaN(1,1);
-
- %events being compared
- if (PavInst==1 | PavInst==2)
- Ev1=strcmp('CSPlus1', RAW(i).Einfo(:,2));
- Ev2=strcmp('CSMinus1', RAW(i).Einfo(:,2));
- else
- Ev1=strcmp('CSPlus', RAW(i).Einfo(:,2));
- Ev2=strcmp('CSMinus', RAW(i).Einfo(:,2));
- end
-
-
- %get number of spikes for this neuron in this bin for all
- %Ev1 trials
- [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev1},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
- for m=1:length(PSR1)
- Ev1Spikes(m,1)=sum(PSR1{1,m}>(binstart));
- end
-
- %pick which trials get used for decoding
- SetupTrials=cat(1,ones(Ev1s,1),zeros(length(PSR1)-Ev1s,1));
- Trials=(SetupTrials(randperm(length(SetupTrials)))==1);
-
- %put spikes from those trials in a matrix
- AllSpikes(Trials,NN)=Ev1Spikes(Trials,1);
-
- %get all the spikes from reward 1p2 trials
- [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev2},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
- for n=1:length(PSR2)
- Ev2Spikes(n,1)=sum(PSR2{1,n}>(binstart));
- end
-
- %pick which trials get used for decoding
- SetupTrials2=cat(1,ones(Ev2s,1),zeros(length(PSR2)-Ev2s,1));
- Trials2=(SetupTrials2(randperm(length(SetupTrials2)))==1);
-
- %put spikes from those trials in a matrix
- AllSpikes((Ev1s+1):(Ev1s+Ev2s),NN)=Ev2Spikes(Trials2,1);
- NN=NN+1; %neuron counter
- AllNN=AllNN+1;
- else AllNN=AllNN+1;
- end
- end
- end
-
-
-
- %Setup spikes matrix for decoding
- DecodeSpikes=AllSpikes(:,Sel);
-
- %find neurons with too few spikes
-
- for t=1:NumNeurons(v)
- if sum(DecodeSpikes(1:Ev1s,t))<3 || sum(DecodeSpikes((1+Ev1s):(Ev1s+Ev2s),t))<3
- LowHz(1,t)=1;
- end
- end
-
- %remove neurons with too few spikes
- % if sum(LowHz)>0
- % Sel2=LowHz<1;
- % DecodeSpikes=DecodeSpikes(:,Sel2);
- % end
-
- %set up validation result matrices
- CVacc = NaN(folds,1);
- CVaccSh = NaN(folds,1);
-
- %normal model
- Partitions = cvpartition(DecodeRs,'KFold',folds);
- for r = 1:folds
-
- SVMModel = fitcdiscr(DecodeSpikes(Partitions.training(r),:),DecodeRs(Partitions.training(r)),'discrimType','diaglinear');
- [prediction,score] = predict(SVMModel,DecodeSpikes(Partitions.test(r),:));
- actual = DecodeRs(Partitions.test(r));
- correct = prediction - actual;
- CVacc(r) = sum(correct==0) / length(correct);
-
- end
-
- %shuffled model
- for q=1:shuffs
- DecodeRsSh=DecodeRs(randperm(length(DecodeRs)));
- PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
- for s = 1:folds
- SVMModelSh = fitcdiscr(DecodeSpikes(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)),'discrimType','diaglinear');
- [predictionSh, score] = predict(SVMModelSh,DecodeSpikes(PartitionsSh.test(s),:));
- actualSh = DecodeRsSh(PartitionsSh.test(s));
- correctSh = predictionSh - actualSh;
- CVaccSh(s) = sum(correctSh==0) / length(correctSh);
- end
- AccShuff(q,1) = nanmean(CVaccSh);
- end
-
- PoolMdlAcc{v,1}(u,l) = nanmean(CVacc);
- PoolMdlAccShuff{v,1}(u,l) = nanmean(AccShuff);
-
- fprintf(['Condition #' num2str(v) ', Rep #' num2str(u) ', Bin #' num2str(l) '\n']);
- end %bins
- end
- end
-
-
-
- if PavInst==1
- PoolMdlAccDSSucrose=PoolMdlAcc;
- PoolMdlAccDSSucroseShuff=PoolMdlAccShuff;
- save('PoolMdlAccDSSucrose.mat','PoolMdlAccDSSucrose')
- save('PoolMdlAccDSSucroseShuff.mat','PoolMdlAccDSSucroseShuff')
- else if PavInst==2
- PoolMdlAccPavSucrose=PoolMdlAcc;
- PoolMdlAccPavSucroseShuff=PoolMdlAccShuff;
- save('PoolMdlAccPavSucrose.mat','PoolMdlAccPavSucrose')
- save('PoolMdlAccPavSucroseShuff.mat','PoolMdlAccPavSucroseShuff')
- else if PavInst==3
- PoolMdlAccDSAlcohol=PoolMdlAcc;
- PoolMdlAccDSAlcoholShuff=PoolMdlAccShuff;
- save('PoolMdlAccDSAlcohol.mat','PoolMdlAccDSAlcohol')
- save('PoolMdlAccDSAlcoholShuff.mat','PoolMdlAccDSAlcoholShuff')
- else if PavInst==4
- PoolMdlAccPavAlcohol=PoolMdlAcc;
- PoolMdlAccPavAlcoholShuff=PoolMdlAccShuff;
- save('PoolMdlAccPavAlcohol.mat','PoolMdlAccPavAlcohol')
- save('PoolMdlAccPavAlcoholShuff.mat','PoolMdlAccPavAlcoholShuff')
-
- end
- end
- end
- end
-
- end
|