123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527 |
- %% generate figure 6 and figure 6 - supplement figure 1
- % fourrier analysis - hierarchical clustering - Random forest approach on
- % early training dataset
- % statistics for the early training dataset
- clc;clear;close all;
- load ('Rearlytraining_light.mat');
- load ('Celltype_earlyTraining.mat');
- load('TRN_DT5_earlytrain.mat')
- Xevent=[-0.25 0.25];
- Yaxis=[-5 10];Xaxis=[1 357];
- Ishow=find(Tm>=Xevent(1) & Tm<=Xevent(2));
- time=[1:1:357];
- BrainRegion=[10 20];
- % ----------- COLORS ---------
- myblue=[0 113/255 188/255];
- mypurple=[200/255 0 255/255];
- mydarkblue=[0 0/255 255/255];
- TickSize=[0.015 0.02];
- Clim=[-5 5];
- Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
- DSconcat_ACQ_allsession=[];Coord_allsession=[];Ninfo_allsession=[];DSmeanz_ACQ_allsession=[];
- TRN_allsession=[];EXCINH_allsession=[];sessionID_allsession=[];
- for i=1:length(Ses)
-
- %normalization based on max
- for k=1:length(Erefnames)
- Ev(k).PSTH_norm1=normalize(Ses(i).Ev(k).PSTHraw,Ses(i).Ev(k).PSTHrawBL,1);
- end
-
- DSconcat_ACQ=[];DSmeanz_ACQ=[];
- %concat with normalized data
- DSconcat_ACQ = Ev(7).PSTH_norm1(:,Ishow);
- DSmeanz_ACQ = Ses(i).Ev(7).Meanz(:,1);
- for k=1:5
- DSconcat_ACQ = cat(2,DSconcat_ACQ,Ev(k).PSTH_norm1(:,Ishow));
- DSmeanz_ACQ = cat(2,DSmeanz_ACQ,Ses(i).Ev(k).MeanzPRE(:,1));
- DSmeanz_ACQ = cat(2,DSmeanz_ACQ,Ses(i).Ev(k).Meanz(:,1));
- end
- DSconcat_ACQ = cat(2,DSconcat_ACQ,Ev(6).PSTH_norm1(:,Ishow));
- DSmeanz_ACQ = cat(2,DSmeanz_ACQ,Ses(i).Ev(6).MeanzPRE(:,1));
-
- DSconcat_ACQ_allsession=cat(1,DSconcat_ACQ_allsession,DSconcat_ACQ);
- Coord_allsession=cat(1,Coord_allsession,Ses(i).Coord(:,4));
- Ninfo_allsession=cat(1,Ninfo_allsession,Ses(i).Ninfo);
- DSmeanz_ACQ_allsession=cat(1,DSmeanz_ACQ_allsession,DSmeanz_ACQ);
- TRN_allsession=cat(1,TRN_allsession,sesTRN(i).TRN(:,1));
- sessionID_allsession=cat(1,sessionID_allsession,repmat(i,length(Ses(i).TRN(:,1)),1));
- if i>3
- EXCINH_allsession=cat(1,EXCINH_allsession,sesTRN(i).TRN(:,2));
- end
- end
- s2={'YVG04','YVG05','YVG08','YVJ11','YVJ15'};
- for j=1:5
- for i=1:length(Ninfo_allsession)
- RatID(i,j)=strncmp(Ninfo_allsession{i,1},s2(j),5);
- end
- end
- RatID_ACQ=sum(RatID,2);
- % concat with Zscore
- DSconcatZ_ACQ_allsession=[];DSevent_ACQ_allsession=[];
- for i=1:length(Ses)
- DSconcat_ACQ=[];DSevent_ACQ=[];
-
- DSconcat_ACQ=Ses(i).Ev(7).PSTHz(:,Ishow);
- DSevent_ACQ=Ses(i).Ev(7).Meanz(:,1);
- for k=1:6 %loop through event
- DSconcat_ACQ=cat(2,DSconcat_ACQ,Ses(i).Ev(k).PSTHz(:,Ishow));
- DSevent_ACQ = cat(2,DSevent_ACQ,Ses(i).Ev(k).Meanz(:,1));
- end
- DSconcatZ_ACQ_allsession=cat(1,DSconcatZ_ACQ_allsession,DSconcat_ACQ);
- DSevent_ACQ_allsession=cat(1,DSevent_ACQ_allsession,DSevent_ACQ);
- end
- selection=Coord_allsession~=0 & Celltype(:,2)==1 & ~isnan(mean(DSconcat_ACQ_allsession,2));
- DSconcat_ACQshort=DSconcat_ACQ_allsession(selection,:);
- DSevent_ACQshort=DSevent_ACQ_allsession(selection,:);
- Coord_short=Coord_allsession(selection,:);
- Ninfo_short=Ninfo_allsession(selection,:);
- %% Clustering on phasic and non phasic neurons with Fourier analysis
- close all;
- f_low=1;f_high=4;
- [X_all,f_all,P_all,Nt_all]=myfft(DSconcat_ACQshort',0.01,1);
- total_power=sum(P_all);
- P_all_norm = bsxfun(@rdivide,P_all, total_power);
- %figure;plot(f_all,P_all_norm);set(gca,'ylim',[0 0.6]);xlabel('Hz');ylabel('Normalized power');
- %title('Normalized FFT of all neurons');
- a=find(f_all>f_low);b=find(f_all<f_high);
- interval_to_sum=[a(1) b(end)];
- lower_window_sum_all=sum(P_all_norm(1:a(1),:));
- upper_window_sum_all=sum(P_all_norm(a(1):b(end),:));
- ratio_all=lower_window_sum_all./upper_window_sum_all;
- close all
- X_ratio=[[lower_window_sum_all]' [upper_window_sum_all]']; %#ok<NBRAK>
- eva_Ch_ratio = evalclusters(X_ratio,'linkage','CalinskiHarabasz','KList',[1:10]);
- phasicID=eva_Ch_ratio.OptimalY;
- Z = linkage(X_ratio,'ward'); % X is the data set, neuron in rows, features in columns
- c = cluster(Z,'Maxclust',3); % determine the number cluster here)
- D = pdist(X_ratio);
- leafOrder = optimalleaforder(Z,D);
- cutoff = median([Z(end-2,3) Z(end-1,3)]);
- %
- % figure
- % subplot(9,4,[29 33])
- % dendrogram(Z,'ColorThreshold',cutoff, 'Reorder',leafOrder);
- %% Figure Fourier transform
- c_ratio=eva_Ch_ratio.OptimalY;
- for i=1:length(c_ratio)
- if c_ratio(i,1)==1
- color(i,:)=[255/255 0/255 0/255]; %phasic red
- elseif c_ratio(i,1)==2
- color(i,:)=[0 255/255 0/255]; %sustained green
- end
- end
- figure
- subplot(5,3,1)
- scatter(lower_window_sum_all,upper_window_sum_all,[],color)
- hold on
- title(strcat('hierarchical clustering'))
- xlabel('power in low freq') % x-axis label
- ylabel('power in int freq')% y-axis label
- legend('phasic','non-phasic')
- hold off
- %% visualization classes
- for i=1:2
- ind=find(phasicID==i);
- t=num2str(i);
- Clim=[-5 5];
- DSconcat_ACQ_class=[];SORTimg=[];MaxValDLS=[];MeanValDLS=[]; MaxTimeDLS=[];TMP=[];NNid=[];
- DSconcat_ACQ_class=DSconcat_ACQshort(ind,:);
-
- for NN=1:size(DSconcat_ACQ_class,1)
- MeanValDLS(NN,1)=mean(DSconcat_ACQ_class(NN,52:306),2);
- [MaxValDLS(NN,1),MaxInd]=max(DSconcat_ACQ_class(NN,:));
- MaxTimeDLS(NN,1)=time(time(1)+MaxInd-1);
- if i==2
- MaxTimePhasic=MaxTimeDLS;
- end
- end
- if i==1
- TMP=MeanValDLS;
- subplot(5,3,[4 7 10 13])
- elseif i==2
- TMP=MeanValDLS;
- subplot(5,3,[5 8 11 14])
- end
- TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
- [~,SORTimg]=sort(TMP);
- %NNid=cell2mat(Ninfo(phasic_neuron_selection,3));NNid2=NNid(SORTimg);
- DSconcat_ACQ_class=DSconcat_ACQ_class(SORTimg,:);
-
-
- imagesc(time,[1 size(DSconcat_ACQ_class,1)],DSconcat_ACQ_class);%colorbar;axis tight;
- colormap(jet);
- hold on
- plot([51 51],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
- plot([102 102],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
- plot([153 153],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
- plot([204 204],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
- plot([255 255],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
- plot([306 306],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
- set(gca,'XTick',[0:25.5:357]);
- set(gca,'xticklabel',Eventnames)
- hold off
- if i==1
- ylabel('Non-Phasic neurons')
- else
- ylabel('Phasic neurons')
- end
- end
- ClassACQ(1:length(Celltype),1)=0;
- ClassACQ(selection,1)=phasicID;
- %% Run random forest for phasic neurons
- load('ClassHC.mat' );
- XTrainingDATA=ClassPHAS(:,2:4);
- YTrainingDATA=ClassPHAS(:,1);
- LIwin=[26:51];%post lever insertion window
- LPwin=[52:280];%pre-1stLP to pre-lastLP window
- PEwin=[307:332];% pre-PE window
- binLI=mean(DSconcat_ACQ_allsession(:,LIwin),2);
- binLP=mean(DSconcat_ACQ_allsession(:,LPwin),2);
- binPE=mean(DSconcat_ACQ_allsession(:,PEwin),2);
- binnedDS=cat(2,binLI,binLP,binPE);
- XTestingDATA=binnedDS(ClassACQ(:,1)==2,:);
- model = classRF_train(XTrainingDATA, YTrainingDATA); % train Random Forest model
- % prepare the data for FOR loop
- X = {XTrainingDATA, XTestingDATA};
- Y = {YTrainingDATA};
- classifier = {'Raw Data','Random Forest'};
- mc = ['r','g']; msym = ['o','o']; msiz = [3,15];
- for i = 1 : 2 % iterate through training and test data
- x = X{i}; y = Y;
- for j = 1 : 2 % iterate through different classifiers
-
- switch(j) % y_hat = running models ? model predictions : y
- case 1 % raw data
- y_hat = y;
- case 2 % Random Forest
- y_hat = classRF_predict(x ,model); % Random Forest predictions
- if i==2
- Class_Acq(:,1)=y_hat;
- end
- end
- end
- end
- ClassACQ(1:length(Celltype),2)=0;
- ClassACQ(ClassACQ(:,1)==2,2)=Class_Acq;
- %% Random forest for non phasic neurons
- XTrainingDATA_NP=ClassNONPHASIC(:,2:4);
- YTrainingDATA_NP=ClassNONPHASIC(:,1);
- XTestingDATA_NP=binnedDS(ClassACQ(:,1)==1,:);
- model = classRF_train(XTrainingDATA_NP, YTrainingDATA_NP); % train Random Forest model
- % prepare the data for FOR loop
- X = {XTrainingDATA_NP, XTestingDATA_NP};
- Y = {YTrainingDATA_NP};
- classifier = {'Raw Data','Random Forest'};
- mc = ['r','g']; msym = ['o','o']; msiz = [3,15];
- for i = 1 : 2 % iterate through training and test data
- x = X{i}; y = Y;
- for j = 1 : 2 % iterate through different classifiers
-
- switch(j) % y_hat = running models ? model predictions : y
- case 1 % raw data
- y_hat = y;
- case 2 % Random Forest
- y_hat = classRF_predict(x ,model); % Random Forest predictions
- if i==2
- Class_Acq_NP(:,1)=y_hat;
- end
- end
- end
- end
- ClassACQ(ClassACQ(:,1)==1,2)=Class_Acq_NP+3;
- %% selection of neurons per session per class
- % 1-3 is for phasic and 4-5 for non phasic
- for i=1:max(Celltype(:,1))
- InDSsel=Ses(i).Coord(:,4)~=0;
- Ses(i).Class(InDSsel,1)=ClassACQ(Celltype(:,1)==i,2);
- Ses(i).Celltype(InDSsel,1)=Celltype(Celltype(:,1)==i,2);
- end
- for k=1:max(ClassACQ(:,2))
- DLSconcat_ACQ_3ds(k).DLS=[];
- DMSconcat_ACQ_3ds(k).DMS=[];
- for i=4:13
- DLSselection=Ses(i).Coord(:,4)==10 & Ses(i).Class(:,1)==k & Ses(i).Celltype(:,1)==1;
- DMSselection=Ses(i).Coord(:,4)==20 & Ses(i).Class(:,1)==k & Ses(i).Celltype(:,1)==1;
-
- DLSconcat_ACQ=[];DMSconcat_ACQ=[];
- DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(7).PSTHz(DLSselection,Ishow));
- DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(7).PSTHz(DMSselection,Ishow));
- for j=1:6
- DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(j).PSTHz(DLSselection,Ishow));
- DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(j).PSTHz(DMSselection,Ishow));
- end
- MeanDLS(k).img(i-3,:)=nanmean(DLSconcat_ACQ,1);
- MeanDMS(k).img(i-3,:)=nanmean(DMSconcat_ACQ,1);
- if i>10
- DLSconcat_ACQ_3ds(k).DLS=cat(1,DLSconcat_ACQ_3ds(k).DLS,DLSconcat_ACQ);
- DMSconcat_ACQ_3ds(k).DMS=cat(1,DMSconcat_ACQ_3ds(k).DMS,DMSconcat_ACQ);
-
- DLSsem(k,:)=nanste(DLSconcat_ACQ(:,time),1);
- DMSsem(k,:)=nanste(DMSconcat_ACQ(:,time),1);
- end
- end
- end
- %% Proportion of neurons in phasic vs non pHasic classes
- for j=1:2
- for i=4:13
- DLSselectionP=Ses(i).Coord(:,4)==BrainRegion(1) & (Ses(i).Class(:,1)<4) & Ses(i).Celltype(:,1)==1 ;
- DMSselectionP=Ses(i).Coord(:,4)==BrainRegion(2) & (Ses(i).Class(:,1)<4) & Ses(i).Celltype(:,1)==1 ;
- if j==1;DSselection=logical(DLSselectionP); else DSselection=logical(DMSselectionP);end
-
- DS_Nb_tot_phasic(i-3,j)=sum(DSselection);
- PercentPhasic(i-3,1+(j-1)*2)=100*sum(DSselection)/sum(Ses(i).Coord(:,4)==BrainRegion(j) & Ses(i).Celltype(:,1)==1);
- PercentPhasic(i-3,2+(j-1)*2)=100-PercentPhasic(i-3,1+(j-1)*2);
- end
- end
- subplot(5,3,2) % plot DLS phasic vs non phasic
- bp1=barh(PercentPhasic(:,1:2),1,'stacked');
- axis([0 100 0.5 10.5])
- xlabel('% DLS neurons')
- ylabel('sessions')
- subplot(5,3,3) % plot DMS phasic vs non phasic
- bp2=barh(PercentPhasic(:,3:4),1,'stacked');
- axis([0 100 0.5 10.5])
- xlabel('% DMS neurons')
- ylabel('sessions')
- %% heat maps classes of neurons
- Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
- titlegraphDLS={'DLS - Start','DLS - Stop', 'DLS - middle','DLS - INH','DLS - EXC'};
- titlegraphDMS={'DMS - Start','DMS - Stop', 'DMS - middle','DMS - INH','DMS - EXC'};
- %titlePSTH={'Inhibited neurons', 'Excited neurons'};
- BrainRegion=[10 20];
- figure
- nplot=0;
- for k=1:max(ClassACQ(:,2)) % loop thru non phasic classes
- if k<4 || k==5
- nplot=nplot+1;
- elseif k==4
- nplot=nplot+10;
- end
- subplot(6,4,2+(nplot-1))
- imagesc(time,[1 size(MeanDLS(k).img,1)],MeanDLS(k).img,Clim); %colorbar;axis tight;
- colormap(jet);
- %colormap(plasma);
- set(gca,'XTick',[0:25.5:357]);
- set(gca,'xticklabel',Eventnames)
- title(titlegraphDLS(k));
- ylabel('sessions');
- hold on
- plot([51 51],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
- plot([102 102],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
- plot([153 153],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
- plot([204 204],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
- plot([255 255],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
- plot([306 306],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
- hold off
-
-
- subplot(6,4,6+(nplot-1))
- imagesc(time,[1 size(MeanDMS(k).img,1)],MeanDMS(k).img,Clim); %colorbar;axis tight;
- colormap(jet);
- set(gca,'XTick',[0:25.5:357]);
- set(gca,'xticklabel',Eventnames)
- title(titlegraphDMS(k));
- ylabel('sessions');
- hold on
- plot([51 51],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
- plot([102 102],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
- plot([153 153],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
- plot([204 204],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
- plot([255 255],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
- plot([306 306],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
- hold off
-
- subplot(6,4,10+(nplot-1))
-
- DLSpsth=nanmean(DLSconcat_ACQ_3ds(k).DLS,1);
- DLS_sem=nanste(DLSconcat_ACQ_3ds(k).DLS,1);
- DLSup=DLSpsth+DLS_sem;
- DLSdown=DLSpsth-DLS_sem;
-
- DMSpsth=nanmean(DMSconcat_ACQ_3ds(k).DMS,1);
- DMS_sem=nanste(DMSconcat_ACQ_3ds(k).DMS,1);
- DMSup=DMSpsth+DMS_sem;
- DMSdown=DMSpsth-DMS_sem;
-
- hold on;
- patch([time,time(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
- plot(time,DLSpsth,'Color',myblue,'linewidth',1.5);
- patch([time,time(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
- plot(time,DMSpsth,'r','linewidth',1.5);
-
- yaxis=[round(min(DMSdown)-0.5,0) round(max(DMSup)+0.5,0)];
- plot([0 357], [0 0],'LineStyle',':','Color','k');
- plot([51 51],yaxis,'LineStyle',':','Color','k')
- plot([102 102],yaxis,'LineStyle',':','Color','k')
- plot([153 153],yaxis,'LineStyle',':','Color','k')
- plot([204 204],yaxis,'LineStyle',':','Color','k')
- plot([255 255],yaxis,'LineStyle',':','Color','k')
- plot([306 306],yaxis,'LineStyle',':','Color','k')
- axis([Xaxis,yaxis]);
- set(gca,'XTick',0:25.5:357);
- set(gca,'xticklabel',Eventnames)
- ylabel('z-score');
- hold off
- end
- %% Proportion of neurons
- for k=1:3
- for i=4:13
- DLSselection=Ses(i).Coord(:,4)==BrainRegion(1) & (Ses(i).Class(:,1)==k) & Ses(i).Celltype(:,1)==1;
- DMSselection=Ses(i).Coord(:,4)==BrainRegion(2) & (Ses(i).Class(:,1)==k) & Ses(i).Celltype(:,1)==1;
-
- DLS_prop_Phasic(i-3,k)=100*sum(DLSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(1) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)<4 & Ses(i).Class(:,1)~=0));
- DMS_prop_Phasic(i-3,k)=100*sum(DMSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(2) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)<4 & Ses(i).Class(:,1)~=0));
-
- DLSselectionNumber(i-3,k)=sum(DLSselection);
- DMSselectionNumber(i-3,k)=sum(DMSselection);
- end
- end
- subplot(6,4,1)
- b1=barh(DLS_prop_Phasic,1,'stacked');
- axis([0 100 0.5 10.5])
- xlabel('% DLS neurons')
- ylabel('sessions')
- subplot(6,4,5)
- b2=barh(DMS_prop_Phasic,1,'stacked');
- axis([0 100 0.5 10.5])
- xlabel('% DMS neurons')
- ylabel('sessions')
- for k=1:2
- for i=4:13
- DLSselection=Ses(i).Coord(:,4)==BrainRegion(1) & (Ses(i).Class(:,1)==k+3) & Ses(i).Celltype(:,1)==1;
- DMSselection=Ses(i).Coord(:,4)==BrainRegion(2) & (Ses(i).Class(:,1)==k+3) & Ses(i).Celltype(:,1)==1;
-
- DLS_prop_nonPhasic(i-3,k)=100*sum(DLSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(1) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)>3 & Ses(i).Class(:,1)~=0));
- DMS_prop_nonPhasic(i-3,k)=100*sum(DMSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(2) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)>3 & Ses(i).Class(:,1)~=0));
-
- DLSselectionNumber(i-3,k)=sum(DLSselection);
- DMSselectionNumber(i-3,k)=sum(DMSselection);
- end
- end
- subplot(6,4,13)
- b1=barh(DLS_prop_nonPhasic,1,'stacked');
- axis([0 100 0.5 10.5])
- xlabel('% DLS neurons')
- ylabel('sessions')
- subplot(6,4,17)
- b2=barh(DMS_prop_nonPhasic,1,'stacked');
- axis([0 100 0.5 10.5])
- xlabel('% DMS neurons')
- ylabel('sessions')
- %% %% make table for statistics
- selectDT5Session=sessionID_allsession(:,1)>3;
- table_ACQ=cat(2,sessionID_allsession(selectDT5Session,1),Coord_allsession(selectDT5Session,1),Celltype(selectDT5Session,2),ClassACQ(selectDT5Session,1),ClassACQ(selectDT5Session,2),DSmeanz_ACQ_allsession(selectDT5Session,:),TRN_allsession(selectDT5Session,1),EXCINH_allsession(:,1), RatID_ACQ(selectDT5Session,1));
- save('table_ACQ.mat','table_ACQ')
- nbevent=[1:12];
- selectionSTAT=table_ACQ(:,2)~=0 & table_ACQ(:,3)==1 & table_ACQ(:,18)~=0;
- sessionCat=categorical(table_ACQ(selectionSTAT,1));
- tableSTAT=array2table(table_ACQ(selectionSTAT,1:20),'VariableNames',{'session','region','celltype','phasicID','class','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12','TRN','EXCINH','habit'});
- tableSTAT.session=categorical(tableSTAT.session);
- tableSTAT.region=categorical(tableSTAT.region);
- tableSTAT.habit=categorical(tableSTAT.habit);
-
- rm = fitrm(tableSTAT,'win1-win12~region*session','WithinDesign',nbevent');
- ranovatbl_event = ranova(rm);
- Btw_ranovatbl_event = anova(rm);
- SphericityAssumption=mauchly(rm);
- posthoc = multcompare(rm,'session', 'by', 'region');
-
-
- rmh = fitrm(tableSTAT,'win1-win12~region*session*habit','WithinDesign',nbevent');
- ranovatbl_event_habit = ranova(rmh);
- Btw_ranovatbl_event_habit = anova(rmh);
- SphericityAssumption_habit=mauchly(rmh);
- posthoc_habit = multcompare(rmh,'region', 'by', 'habit');
-
- % for i=1:max(table_ACQ(:,5))
- for i=1:3 %loop thru phasic neurons, problem, not enough neurons in non phasic classes
- selectionSTAT1=table_ACQ(:,2)~=0 & table_ACQ(:,3)==1 & table_ACQ(:,5)==i;
-
- tableSTAT1=array2table(table_ACQ(selectionSTAT1,1:20),'VariableNames',{'session','region','celltype','phasicID','class','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12','TRN','EXCINH','habit'});
- tableSTAT1.session=categorical(tableSTAT1.session);
- tableSTAT1.region=categorical(tableSTAT1.region);
- tableSTAT1.habit=categorical(tableSTAT1.habit);
-
- rm = fitrm(tableSTAT1,'win1-win12~region*session','WithinDesign',nbevent');
- stat(i).ranovatbl_event = ranova(rm);
- stat(i).Btw_ranovatbl_event = anova(rm);
-
- [~,stat(i).LIevent,~]=anovan(table_ACQ(selectionSTAT1,6),{table_ACQ(selectionSTAT1,2)},'varnames',{'region'}, 'display','off');
- [~,stat(i).LRevent,~]=anovan(table_ACQ(selectionSTAT1,16),{table_ACQ(selectionSTAT1,2)},'varnames',{'region'}, 'display','off');
- [~,stat(i).PEevent,~]=anovan(table_ACQ(selectionSTAT1,17),{table_ACQ(selectionSTAT1,2)},'varnames',{'region'}, 'display','off');
- rmh = fitrm(tableSTAT1,'win1-win12~region*session*habit','WithinDesign',nbevent');
- stat(i).ranovatbl_event_habit = ranova(rmh);
- stat(i).Btw_ranovatbl_habit = anova(rmh);
- stat(i).posthoc_habit = multcompare(rmh,'region', 'by', 'habit');
- end
- for i=1:10
- selectionSTAT2=table_ACQ(:,2)~=0 & table_ACQ(:,3)==1 & table_ACQ(:,1)==i+3 & table_ACQ(:,4)==1;
- [chi(i).tbl1,chi(i).chi2,chi(i).p1,chi(i).labels] = crosstab(table_ACQ(selectionSTAT2,2),table_ACQ(selectionSTAT2,5));
- end
-
- save('Rearlytraining_light.mat','Erefnames','Ses','Tm')
|