J_generateFigure10_s3_shuffle_RRFR.m 8.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194
  1. %% generate Figure 10 figure supplement 3
  2. clear all;clc;
  3. tic
  4. global Dura Baseline Tm Tbase BSIZE Tbin
  5. SAVE_FLAG=1;
  6. BSIZE=0.05;
  7. Dura=[-25 60]; Tm=-1:BSIZE:60;
  8. Baseline=[-2 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %used to calculate Z-scores
  9. Tbin=-1:0.005:1; %window used to determine the optimal binsize
  10. PStat=0.05;
  11. prewin=[Dura(1) -20];
  12. postwin=[0 .5];
  13. EventWin=[-0.5 0.5];
  14. postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
  15. Slopebounds=[-1.5:0.25:1.5];
  16. R2Bounds=[0:0.05:0.5];
  17. PvalBounds=[0:0.01:0.5];
  18. Struct=[10 20];
  19. XI=1;
  20. RunRegress=0;
  21. SAVEFIG=1;
  22. %R=[];R.Ninfo={};
  23. NN=0;
  24. pre=[]; post=[];
  25. if ~exist('RAW'), load ('RAWextendedtraining.mat'); RAW=RAW; end
  26. if ~exist('Ev'), load('Rextendedtraining_light.mat'); end
  27. load('C.mat');
  28. load('Celltype_extendedTraining.mat')
  29. % 1 for inclusion of omission in analysis, max time=60s / 0 for exclusion omission
  30. %Smoothing
  31. Smoothing=1; %0 for raw and 1 for smoothing
  32. SmoothTYPE='lowess';
  33. SmoothSPAN=5;
  34. if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
  35. %%
  36. % PREALLOCATING
  37. for i=1:length(RAW)
  38. EvInd=strmatch('First_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
  39. Rind=strmatch('ReliableLP_DT',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
  40. if sum(EvInd)~=0 && sum(Rind)~=0
  41. DS(i)=length(RAW(i).Erast{EvInd});
  42. Neur(i)=length(RAW(i).Nrast);
  43. end
  44. end
  45. MaxTrials=max(DS); MaxNeur=sum(Neur); MaxWin=length(postwin);
  46. Cshuffle.preRRRR=NaN(MaxNeur,MaxTrials); % (neurons, trials)
  47. Cshuffle.SEQ=NaN(MaxNeur,MaxTrials); % (neurons, trials)
  48. Cshuffle.postwinRR=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
  49. for i=1:length(RAW) %loops through sessions
  50. EvInd=strmatch('First_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
  51. Rind=strmatch('ReliableLP_DT',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
  52. if sum(EvInd)~=0 && sum(Rind)~=0
  53. Dcell2=MakePSR05(RAW(i).Erast(Rind),RAW(i).Erast{EvInd},[-1 60],{2}); % makes trial by trail
  54. for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
  55. NN=NN+1;
  56. [PSRNeur,Nneur]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},[-1 60],{1});% collapsed raster
  57. [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
  58. [PTHneur,BW1,~]=MakePTH07(PSRNeur,repmat(N2, size(RAW(i).Nrast{j},1),1),{2,0,BSIZE});
  59. PTHneur=smooth(PTHneur,SmoothSPAN,SmoothTYPE)';
  60. Cshuffle.PSTHraw(NN,1:length(Tm))=PTHneur;
  61. Search=find(Tm>=EventWin(1) & Tm<=EventWin(2));
  62. [Cshuffle.MaxVal(NN,1),MaxInd]=max(Cshuffle.PSTHraw(NN,Search));
  63. Cshuffle.MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
  64. for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
  65. SequenceDur=RAW(i).Eint{1,6}(m)-RAW(i).Eint{1,5}(m);
  66. Cshuffle.RR(NN,m)=length(Dcell2{1,m})/SequenceDur;
  67. Cshuffle.preRR(NN,m)=sum(PSR2{1,m}<prewin(2) & PSR2{1,m}>prewin(1))/(abs(prewin(2)-prewin(1))); %neurons, trials
  68. Cshuffle.SEQ(NN,m)=sum(PSR2{1,m}<SequenceDur & PSR2{1,m}>0)/SequenceDur;
  69. end %m loop
  70. fprintf('Neuron ID # %d\n',NN);
  71. end %j loop
  72. end %if
  73. end %i loop
  74. %% Runs the regression analysis on XX number of iterations of shuffled data
  75. XX=1000;
  76. for n=1:XX
  77. for NN=1:MaxNeur
  78. Cshuffle.SEQShuffled(NN,:)=shuffle(Cshuffle.SEQ(NN,:));
  79. RegX=[Cshuffle.SEQShuffled(NN,:)', ones(size(Cshuffle.SEQShuffled,2),1)];
  80. [Cshuffle.SLOPEpostshuffled(NN,n),Cshuffle.STATSpostshuffled(NN,n)]= corr(Cshuffle.RR(NN,:)',Cshuffle.SEQShuffled(NN,:)','rows','pairwise','type','Spearman');
  81. end
  82. fprintf('Shuffled Iteration # %d\n',n);
  83. save('Cshuffle.mat', 'Cshuffle');
  84. end
  85. %% Calculates number of significant correlations and the mean correlation coefficient for each
  86. selection1=Coord(:,4)~=0;
  87. for q=1:length(Cshuffle.SLOPEpostshuffled(1,:))
  88. Cshuffle.SLOPEmeanshuffled(q)=nanmean(Cshuffle.SLOPEpostshuffled(selection1,q));
  89. Cshuffle.SIGshuffled(q)=sum(Cshuffle.STATSpostshuffled(selection1,q)<0.05);
  90. Cshuffle.NEGSIGshuffled(q)=sum(Cshuffle.STATSpostshuffled(selection1,q)<0.05 & Cshuffle.SLOPEpostshuffled(selection1,q)<0);
  91. end
  92. %% Calculates correlation coefficients and significant values for real data
  93. for NN=1:MaxNeur
  94. RegX=[Cshuffle.SEQ(NN,:)', ones(size(Cshuffle.SEQ,2),1)];
  95. [Cshuffle.SLOPEpre,Cshuffle.STATSpre]=corr(Cshuffle.preRR(NN,:)',Cshuffle.SEQ(NN,:)','rows','pairwise','type','Spearman');
  96. [Cshuffle.SLOPEpost(NN,:),Cshuffle.STATSpost(NN,:)]= corr(Cshuffle.RR(NN,:)',Cshuffle.SEQ(NN,:)','rows','pairwise','type','Spearman');
  97. [Cshuffle.SLOPEpostpre(NN,:),Cshuffle.STATSpostpre(NN,:)]= corr((Cshuffle.RR(NN,:)-Cshuffle.preRR(NN,:))',Cshuffle.SEQ(NN,:)','rows','pairwise','type','Spearman');
  98. fprintf('CORREL Neuron ID # %d\n',NN);
  99. end
  100. save('Cshuffle.mat', 'Cshuffle');
  101. %% Figure of distributions of means correlation coefficients and number of significant neurons
  102. selection1=Coord(:,4)==10 & Celltype(:,1)==1;
  103. selection2=Coord(:,4)==20 & Celltype(:,1)==1;
  104. figure;
  105. % Plot of distribution of mean correlation coefficients for all of the
  106. % shuffled data
  107. subplot(2,2,3);
  108. histogram(Cshuffle.SLOPEmeanshuffled,100);
  109. hold on;
  110. MEANsel1=nanmean(Cshuffle.SLOPEpost(selection1,1));
  111. MEANsel2=nanmean(Cshuffle.SLOPEpost(selection2,1));
  112. plot([MEANsel1 MEANsel1], [0 50],'Color','b');
  113. plot([MEANsel2 MEANsel2], [0 50],'Color','r');
  114. xlabel('Mean correlation coefficient') % x-axis label
  115. ylabel('# of iterations')% y-axis label
  116. %Plot of distribution of number of significantly correlated neurons across
  117. %all the shuffled data sets
  118. subplot(2,2,4);
  119. histogram((Cshuffle.SIGshuffled*100/849),40);
  120. xlabel('% of significantly correlated unit') % x-axis label
  121. ylabel('# of iterations')% y-axis label
  122. hold on;
  123. SigNUM1=100*sum(Cshuffle.STATSpost(selection1,1)<0.05)/(sum(selection1));
  124. SigNUM2=100*sum(Cshuffle.STATSpost(selection2,1)<0.05)/(sum(selection2));
  125. plot([SigNUM1 SigNUM1], [0 100],'Color','b');
  126. plot([SigNUM2 SigNUM2], [0 100],'Color','r');
  127. % Plot of correlation coefficients from example shuffled data set
  128. subplot(2,2,1);
  129. sel1=Coord(:,4)~=0;
  130. sel2=sel1 & Cshuffle.STATSpostshuffled(:,1)<PStat;
  131. xmin=floor(min(Cshuffle.SLOPEpostshuffled(sel1,1)));
  132. xmax=ceil(max(Cshuffle.SLOPEpostshuffled(sel1,1)));
  133. step=0.05;
  134. xcenters=xmin:step:xmax;
  135. Nsel1=hist(Cshuffle.SLOPEpostshuffled(sel1,1),xcenters);
  136. Nsel2=hist(Cshuffle.SLOPEpostshuffled(sel2,1),xcenters);
  137. hold on;
  138. bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
  139. bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
  140. plot([0 0], [0 30],'LineStyle',':','Color','k');
  141. MEANsel1=mean(Cshuffle.SLOPEpostshuffled(sel1,1));
  142. plot([MEANsel1 MEANsel1], [0 10],'Color','r');
  143. xlabel('% correlation coefficient shuffle data') % x-axis label
  144. ylabel('# of neurons')% y-axis label
  145. hold on;
  146. % %Plot of correlation coefficients with significance marked based on
  147. % %bootstrapping
  148. subplot(2,2,2);
  149. sel1=Coord(:,4)~=0;
  150. sel2= sel1 & Cshuffle.STATSpost(:,1)<PStat;
  151. xmin=floor(min(Cshuffle.SLOPEpost(sel1,1)));
  152. xmax=ceil(max(Cshuffle.SLOPEpost(sel1,1)));
  153. step=0.05;
  154. xcenters=xmin:step:xmax;
  155. Nsel1=hist(Cshuffle.SLOPEpost(sel1,1),xcenters);
  156. Nsel2=hist(Cshuffle.SLOPEpost(sel2,1),xcenters);
  157. hold on;
  158. bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
  159. bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
  160. plot([0 0], [0 30],'LineStyle',':','Color','k');
  161. MEANsel1=mean(Cshuffle.SLOPEpost(sel2,1));
  162. plot([MEANsel1 MEANsel1], [0 10],'Color','r');
  163. xlabel('% correlation coefficient real data') % x-axis label
  164. ylabel('# of neurons')% y-axis label