Delete 'MatlabScripts/k_Fig6ThreeRewards.m'

MatlabScripts/k_Fig6ThreeRewards.m

@@ -1,339 +0,0 @@
-%plotting 3rewards data
-
-clear all;
-load ('R_3R.mat');
-load ('RAW.mat');
-
-%get parameters
-BinBase=R_3R.Param.BinBase;
-BinDura=R_3R.Param.BinDura;
-bins=R_3R.Param.bins;
-binint=R_3R.Param.binint;
-binstart=R_3R.Param.binstart;
-PStatBins=0.01; %using more stringent cutoff to reduce pre-delivery noise
-
-%which bins bound the area examined for reward-selectivity? (in seconds)
-Time1=0.4; %seconds
-Time2=3; %seconds
-Bin1=(((Time1-BinDura(2)/2)-binstart)/binint); %convert to bin name
-Bin2=(((Time2-BinDura(2)/2)-binstart)/binint); %convert to bin name
-
-%sorting bin -- which bin the neurons' activity is sorted on for heatmap(in seconds)
-SortBinTime=1.1; %seconds
-SortBin=round((((SortBinTime-BinDura(2)/2)-binstart)/binint)); %convert to bin name
-
-sucrose=[1  0.6  0.1];
-maltodextrin=[.9  0.3  .9];
-water=[0.00  0.75  0.75];
-total=[0.3  0.1  0.8];
-exc=[0 113/255 188/255];
-inh=[240/255 0 50/255];
-
-%% conduct lick analysis
-global Dura Tm BSIZE Tbin
-path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\3RLick_paper.xls';
-
-%Main settings
-BSIZE=0.01; %Do not change
-Dura=[-22 20]; Tm=Dura(1):BSIZE:Dura(2);
-Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
-MinNumTrials=5;
-Lick=[];Lick.Ninfo={};LL=0;Nlick=0;
-
-%Smoothing
-Smoothing=1; %0 for raw and 1 for smoothing
-SmoothTYPE='lowess'; %can change this between lowess and rlowess (more robust, ignores outliers more)
-SmoothSPAN=50; %percentage of total data points
-if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
-
-% List of events to analyze and analysis windows EXTRACTED from excel file
-[~,Erefnames]=xlsread(path,'Windows','a3:a8'); % cell that contains the event names
-
-%Finds the total number of sessions
-for i=1:length(RAW)
-    if strcmp('TH',RAW(i).Type(1:2))
-        Nlick=Nlick+1;
-        Lick.Linfo(i,1)=RAW(i).Ninfo(1,1);
-    end
-end
-Lick.Erefnames= Erefnames;
-
-%preallocating the result matrix
-for k=1:length(Erefnames)
-    Lick.Ev(k).PSTHraw(1:Nlick,1:length(Tm))=NaN(Nlick,length(Tm));
-    Lick.Ev(k).BW(1:Nlick,1)=NaN;
-    Lick.Ev(k).NumberTrials(1:Nlick,1)=NaN;
-end
-
-for i=1:length(RAW) %loops through sessions
-    if strcmp('TH',RAW(i).Type(1:2))
-        LL=LL+1; %lick session counter
-        for k=1:length(Erefnames) %loops thorough the events
-            EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
-            LickInd=strcmp('Licks',RAW(i).Einfo(:,2)); %find the event id number from RAW
-            if sum(EvInd)==0
-                fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
-            end
-
-            Lick.Ev(k).NumberTrials(LL,1)=length(RAW(i).Erast{EvInd});
-            if  ~isempty(EvInd) && Lick.Ev(k).NumberTrials(LL,1)>MinNumTrials %avoid analyzing sessions where that do not have enough trials
-                [PSR1,N1]=MakePSR04(RAW(i).Erast(LickInd),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
-                if ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
-                    %forcing 100ms bin size to keep it consistent across
-                    %sessions (no reason is should be different for licks)
-                    [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Erast{LickInd},1),1),{2,0,0.1});%these values force bin size to be 100ms
-                    PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
-
-                    %------------- Fills the R.Ev(k) fields --------------
-                    Lick.Ev(k).BW(LL,1)=BW1;
-                    Lick.Ev(k).PSTHraw(LL,1:length(Tm))=PTH1;
-                end
-            end
-        end
-    end
-end
-
-Xaxis=[-2 12];
-Ishow=find(Tm>=Xaxis(1) & Tm<=Xaxis(2));
-time1=Tm(Ishow);
-c=[-1000 7500];ClimE=sign(c).*abs(c).^(1/4);%ColorMapExc
-colormap(jet);
-sucrose=[1  0.6  0.1];
-maltodextrin=[.9  0.3  .9];
-water=[0.00  0.75  0.75];
-total=[0.3  0.1  0.8];
-exc=[0 113/255 188/255];
-inh=[240/255 0 50/255];
-
-Ev1=strcmp('RD1', Lick.Erefnames);
-Ev2=strcmp('RD2', Lick.Erefnames);
-Ev3=strcmp('RD3', Lick.Erefnames);
-
-psth1=nanmean(Lick.Ev(Ev1).PSTHraw(:,Ishow),1);
-sem1=nanste(Lick.Ev(Ev1).PSTHraw(:,Ishow),1); %calculate standard error of the mean
-up1=psth1+sem1;
-down1=psth1-sem1;
-
-psthE=nanmean(Lick.Ev(Ev2).PSTHraw(:,Ishow),1);
-semE=nanste(Lick.Ev(Ev2).PSTHraw(:,Ishow),1); %calculate standard error of the mean
-upE=psthE+semE;
-downE=psthE-semE;
-
-psth3=nanmean(Lick.Ev(Ev3).PSTHraw(:,Ishow),1);
-sem3=nanste(Lick.Ev(Ev3).PSTHraw(:,Ishow),1); %calculate standard error of the mean
-up3=psth3+sem3;
-down3=psth3-sem3;
-
-%plotting
-subplot(2,3,1);
-hold on;
-plot(time1,psth1,'Color',sucrose,'linewidth',1);
-plot(time1,psthE,'Color',maltodextrin,'linewidth',1);
-plot(time1,psth3,'Color',water,'linewidth',1);
-patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
-patch([time1,time1(end:-1:1)],[upE,downE(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
-patch([time1,time1(end:-1:1)],[up3,down3(end:-1:1)],water,'EdgeColor','none');alpha(0.5);
-title('Mean lick rate');
-%plot([-2 10],[0 0],':','color','k');
-plot([0 0],[-2 8],':','color','k','linewidth',0.75);
-axis([-2 12 0 7.5]);
-xlabel('Seconds from reward delivery');
-ylabel('Licks/s');
-legend('Sucrose trials','Maltodextrin trials','Water trials');
-
-
-%% bins analysis
-
-%first and last bin aren't included because you can't compare to the previous/subsequent bin
-%this axis plots the bins on their centers
-xaxis=linspace(binstart+binint+BinDura(2)/2,binstart+(bins-2)*binint+BinDura(2)/2,bins-2);
-
-for i=2:(bins-1) %no including first or last bin because can't compare to previous/subsequent bin
-    
-    %finds out whether firing is stronger (high excitation or lower inhibition) for 1 or 2
-    for k=1:length(R_3R.Ninfo) %runs through neurons
-        if R_3R.BinStatRwrd{i,1}(k).IntSig < PStatBins %if neuron is significant for this bin
-            if R_3R.BinStatRwrd{i-1,1}(k).IntSig < PStatBins || R_3R.BinStatRwrd{i+1,1}(k).IntSig < PStatBins %either previous or subsequent bin must be significant
-                if R_3R.BinStatRwrd{i,1}(k).R1Mean >= R_3R.BinStatRwrd{i,1}(k).R2Mean && R_3R.BinStatRwrd{i,1}(k).R1Mean >= R_3R.BinStatRwrd{i,1}(k).R3Mean %if firing is greater on sucrose trials than maltodextrin and water. if there's a tie, make it sucrose (this is highly unlikely)
-                    R_3R.BinRewPref{i,1}(k,1)=1; %sucrose preferring
-                elseif R_3R.BinStatRwrd{i,1}(k).R2Mean > R_3R.BinStatRwrd{i,1}(k).R1Mean && R_3R.BinStatRwrd{i,1}(k).R2Mean >= R_3R.BinStatRwrd{i,1}(k).R3Mean %if mal is greater than suc and water
-                    R_3R.BinRewPref{i,1}(k,1)=2; %maltodextrin preferring
-                else
-                    R_3R.BinRewPref{i,1}(k,1)=3; %water preferring
-                end
-            else
-                R_3R.BinRewPref{i,1}(k,1)=0; %if not significant in 2 consecutive bins
-            end
-        else
-            R_3R.BinRewPref{i,1}(k,1)=0; %if no sig reward modulation
-        end
-        
-    end
-    %find how many NAc neurons have significant reward modulation in each bin
-    NN1perBin(i,1)=sum(R_3R.BinRewPref{i,1}(:,1)==1); %sucrose pref
-    NN2perBin(i,1)=sum(R_3R.BinRewPref{i,1}(:,1)==2); %malto pref
-    NN3perBin(i,1)=sum(R_3R.BinRewPref{i,1}(:,1)==3); %malto pref
-    NNperBin(i,1)=sum(R_3R.BinRewPref{i,1}(:,1)>0); %any
-    %normalize to number of neurons in population
-    NN1norm=NN1perBin./length(R_3R.Ninfo);
-    NN2norm=NN2perBin./length(R_3R.Ninfo);
-    NN3norm=NN3perBin./length(R_3R.Ninfo);
-    NNnorm=NNperBin./length(R_3R.Ninfo);
-    
-end
-
-
-%plotting number of significantly modulated neurons across time
-%NAc
-subplot(2,3,2);
-hold on;
-plot(xaxis,NNnorm(2:bins-1),'Color', total,'linewidth',1.5);
-plot(xaxis,NN1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
-plot(xaxis,NN2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
-plot(xaxis,NN3norm(2:bins-1),'Color',water,'linewidth',1.5);
-plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
-plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
-axis([xaxis(1) xaxis(end) 0 0.85]);
-legend('Total','Suc > mal & wat','Mal > suc & wat','Water > suc & mal','Location','northeast')
-ylabel('Fraction of population');
-xlabel('Seconds from RD');
-title('Reward-selective neurons');
-
-
-%% plotting sucrose-selective neurons
-
-%color map
-[magma,inferno,plasma,viridis]=colormaps;
-colormap(plasma);
-c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
-
-%events we're looking at
-RD1=strcmp('RD1', R_3R.Erefnames);
-RD2=strcmp('RD2', R_3R.Erefnames);
-RD3=strcmp('RD3', R_3R.Erefnames);
-
-%setting up parameters
-Xaxis=[-2 5];
-inttime=find(R_3R.Param.Tm>=Xaxis(1) & R_3R.Param.Tm<=Xaxis(2));
-paramtime=R_3R.Param.Tm(inttime);
-
-%find all neurons with greater firing for sucrose
-for i = 1:(Bin2-Bin1+1)
-    %the added +1 is necessary because bin 20 is the 21st entry in the matrix
-    Pref1(:,i)=R_3R.BinRewPref{Bin1+i}==1; %get neurons that have greater firing for sucrose in any of the bins bounded above
-    Resp11(:,i)=Pref1(:,i)&cat(1,R_3R.BinStatRwrd{Bin1+i,1}.SucRespDir)==1; %get neurons with excitation to sucrose
-    Resp12(:,i)=Pref1(:,i)&cat(1,R_3R.BinStatRwrd{Bin1+i,1}.MalRespDir)==1;%get neurons with inhibition to maltodextrin
-    Resp13(:,i)=Pref1(:,i)&cat(1,R_3R.BinStatRwrd{Bin1+i,1}.WatRespDir)==-1;%get neurons with inhibition to maltodextrin
-end
-
-Sel=sum(Pref1,2)>0;  %all neurons selective in any bin
-Sel1=sum(Resp11,2)>0; %all selective neurons sucrose excited in any bin
-Sel3=sum(Resp12,2)>0; %all selective neurons malto inhibited in any bin
-Sel6=sum(Resp13,2)>0; %all selective neurons malto inhibited in any bin
-
-subplot(2,4,5); %heatmap of suc preferring neurons' response to sucrose
-Neurons=R_3R.Ev(RD1).PSTHz(Sel,inttime); %get the firing rates of neurons of interest
-SucResp=cat(1,R_3R.BinStatRwrd{SortBin+1,1}.R1Mean); %sucrose responses
-TMP=SucResp(Sel); %find the magnitude of the excitations for water for this binTMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
-[~,SORTimg]=sort(TMP);
-Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-imagesc(paramtime,[1,sum(Sel,1)],Neurons,ClimE);
-title(['Sucrose trials']);
-xlabel('Seconds from RD');
-ylabel('Individual neurons');
-hold on;
-plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
-
-
-subplot(2,4,6); %heatmap of suc preferring neurons' response to maltodextrin
-Neurons=R_3R.Ev(RD2).PSTHz(Sel,inttime); %get the firing rates of neurons of interest
-Neurons=Neurons(SORTimg,:); %sort the neurons same as before
-imagesc(paramtime,[1,sum(Sel,1)],Neurons,ClimE);
-title(['Maltodextrin trials']);
-xlabel('Seconds from RD');
-hold on;
-plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
-
-subplot(2,4,7); %heatmap of suc preferring neurons' response to water
-Neurons=R_3R.Ev(RD3).PSTHz(Sel,inttime); %get the firing rates of neurons of interest
-Neurons=Neurons(SORTimg,:); %sort the neurons same as before
-imagesc(paramtime,[1,sum(Sel,1)],Neurons,ClimE);
-title(['Water trials']);
-xlabel('Seconds from RD');
-hold on;
-plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
-
-%plot suc preferring neurons to suc
-psthE=nanmean(R_3R.Ev(RD1).PSTHz(Sel,inttime),1); 
-semE=nanste(R_3R.Ev(RD1).PSTHz(Sel,inttime),1); %calculate standard error of the mean
-upE=psthE+semE;
-downE=psthE-semE;
-
-%plot suc preferring neurons to malt
-psth2=nanmean(R_3R.Ev(RD2).PSTHz(Sel,inttime),1); 
-sem2=nanste(R_3R.Ev(RD2).PSTHz(Sel,inttime),1); %calculate standard error of the mean
-up2=psth2+sem2;
-down2=psth2-sem2;
-
-%plot suc preferring neurons to water
-psth3=nanmean(R_3R.Ev(RD3).PSTHz(Sel,inttime),1); 
-sem3=nanste(R_3R.Ev(RD3).PSTHz(Sel,inttime),1); %calculate standard error of the mean
-up3=psth3+sem3;
-down3=psth3-sem3;
-
-%plotting
-subplot(2,3,3);
-hold on;
-plot(paramtime,psthE,'Color',sucrose,'linewidth',1);
-plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
-plot(paramtime,psth3,'Color',water,'linewidth',1);
-patch([paramtime,paramtime(end:-1:1)],[upE,downE(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
-patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
-patch([paramtime,paramtime(end:-1:1)],[up3,down3(end:-1:1)],water,'EdgeColor','none');alpha(0.5);
-legend('Suc trials','Mal trials','Wat trials');
-plot([-2 5],[0 0],':','color','k','linewidth',0.75);
-plot([0 0],[-2 8],':','color','k','linewidth',0.75);
-axis([-2 5 -2 4.7]);
-ylabel('Mean firing (z-score)');
-title(['Suc>mal&wat (n=' num2str(sum(Sel)) ' of ' num2str(length(Sel)) ')'])
-xlabel('Seconds from RD');
-
-G = sum(Sel1&Sel3&Sel6);
-F = sum(Sel3&Sel6);
-E = sum(Sel1&Sel6);
-D = sum(Sel1&Sel3);
-A = sum(Sel1);
-B = sum(Sel3);
-C = sum(Sel6);
-
-%vertical venn diagram
-x = [0 0 1 1];
-y1 = [C-E C-E+A C-E+A C-E];
-y2 = [C-F C-F+B C-F+B C-F];
-y3 = [0 C C 0];
-
-subplot(2,35,64);
-hold on;
-s = patch(x,y1,sucrose);
-m = patch(x,y2,maltodextrin);
-w = patch(x,y3,water);
-alpha(s,0.7);
-alpha(w,0.5);
-alpha(m,0.5);
-set(gca,'xtick',[]);
-ylabel('Distribution of neurons');
-axis([0 1 0 C-E+A]);
-
-
-%% stats on reward firing averaged together
-%for simplicity, just looking at average activity in bin centered at 1sec
-%because I already have that data collected
-
-SucResp=cat(1,R_3R.BinStatRwrd{SortBin+1,1}.R1Mean); %suc responses
-MalResp=cat(1,R_3R.BinStatRwrd{SortBin+1,1}.R2Mean); %mal responses
-WatResp=cat(1,R_3R.BinStatRwrd{SortBin+1,1}.R3Mean); %wat responses
-
-[~,R_3R.RewRespStat{1,1},R_3R.RewRespStat{1,2}]=anovan(cat(1,SucResp(Sel),MalResp(Sel),WatResp(Sel)),cat(1,zeros(sum(Sel),1),ones(sum(Sel),1),2*ones(sum(Sel),1)),'varnames',{'Reward'},'Display','off');
-R_3R.RewRespStat{1,3}=multcompare(R_3R.RewRespStat{1,2},'display','off');
-
-save('R_3R.mat','R_3R');