ccf_phase_plot2_V3.R

ccf_phase_plot2<-function(file,gal4,namecolor){
  #rm(list=ls())
  #rm(res)
  source('~/Dropbox/R/by_date2_V2.R')
  if(gal4=="tim"){cell=c(1:5)}
  if(gal4=="pdf"){cell=c(2,3)}
  if(gal4=="pdfrB"){cell=c(1,6,3)}
  if(gal4=="pdf80"){cell=c(4,5,1)}
  if(gal4=="pdfrB_pdf80"){cell=c(1,6)}
  if(gal4=="929"){cell=c(1,2)}
  if(gal4=="645"){cell=c(4,5,1,3)}
  if(gal4=="681"){cell=c(4,5,1,2)}
  
  if(gal4=="TH"){cell=c(1,2,4,5)}
  #if(gal4=="TH"){cell=c(1,2,3,4,5)}
  if(gal4=="THC"){cell=c(1,2)}
  if(gal4=="THD"){cell=c(4,5)}
  #if(gal4=="TH-PPM1"){cell=c(4)}
  if(is.numeric(gal4)){
    cell=seq(1,gal4)
    #cell=gal4
  }
  
  by_date2(file,gal4)
  
  
  
  setwd(sprintf("~/Dropbox/lab/data sheet/%s/date",file))
  name=list.files()
  cells=length(name)
  
  load(sprintf("~/Dropbox/lab/data sheet/%s/%s_celltype.RData",file,file))
  lgd=c("LNd","lLNv","sLNv","DN3","DN1","fifth_sLNv")
  a=list()
  ####################
  par(mai=c(0.1,0.1,0.1,0.1))
  par(cex=2)
  ####################
  
  for (k in cell){
    list=get(lgd[k])
    if(is.list(list)){
      ff=matrix(,144,length(list))
      for (i in 1:length(list)){
        ff[,i]=list[[i]]
      }
    }
    else {ff=list}  
    
    #cells=dim(ff)[2]
    mean=colMeans(ff,na.rm=TRUE)
    ma=matrix(,cells,cells)
    lag=ma
    for (i in 1:cells){
      for (j in 1:cells){
        ac=ccf(ff[,i],ff[,j],lag.max=144,na.action=na.pass,plot=FALSE)
        ma[i,j]=max(ac[[1]])
        if(!is.na(ac[[1]])){
          lag[i,j]=(ac[[4]][which.max(ac[[1]])])*10/60
        }
        
      }
    }
    for(i in 1:cells){
      if(!is.na(ff[,i])){
        lag[,i]=lag[,i]+which.max(ff[,i])*(24/144) 
      }
    }
    
    source('~/Dropbox/R/phaseplot.R')
    DA=0
    if(grepl("TH",gal4)){
      DA=1 
    }
    if(!missing(namecolor)){
      if (exists('res')) {phaseplot2(lag,k,res,DA=DA,namecolor=namecolor)}
      else{res=phaseplot2(lag,k,DA=DA,namecolor=namecolor)}
      }
    else{  
      if (exists('res')) {  
        phaseplot2(lag,k,res,DA=DA)
      }
      else {res=phaseplot2(lag,k,DA=DA)} }
    
    a[length(a)+1]=list(phaseplot3(lag))

  }
  

  #if(grepl("TH",gal4)){ccf_phase_plot2_eb(file,5,res)}
  
  
  
  #title(main="Calcium rhythms",cex=1)
  #x=circular(seq(mean(a[[1]],na.rm=T), mean(a[[3]],na.rm=T),-0.01),units="hours",template="clock24")
  #lines(x,rep(-0.5,length(x)),plot.info=res,lwd=2)
  
  
  ##############summary#########
  name=meanphase=varphase=n=rayleigh.test.p=matrix(,6,1)
  sum=data.frame(name,meanphase,varphase,n,rayleigh.test.p)
  psi=matrix(,6,cells)
  for(i in 1:length(a)){
    k=cell[i]
    sum$name[k]=lgd[k]
    sum$meanphase[k]=signif(mean(a[[i]],na.rm=T)[[1]],2)
    sum$varphase[k]=signif(var(a[[i]],na.rm=T),2)
    sum$n[k]=length(a[[i]])-sum(is.na(a[[i]]))
    sum$rayleigh.test.p[k]=signif(rayleigh.test(a[[i]], mu = NULL)$p.value,2) 
    
    #psi[k,]=a[[i]]-a[[3]]
  }
  write.csv(sum,file=sprintf("~/Dropbox/lab/data sheet/%s/summary.csv",file))
  
  rownames(psi)=lgd
  
  for (i in 1:6){
    for (j in 1:cells){
      if(!is.na(psi[i,j])){
        if(psi[i,j]<0){psi[i,j]=24+psi[i,j]}   
      }    
    }
  }
  write.csv(psi,file=sprintf("~/Dropbox/lab/data sheet/%s/psi.csv",file))
  
  #print(mean(a[[1]],na.rm=T)-mean(a[[3]],na.rm=T))
  # E=a[[1]]-mean(a[[3]],na.rm=T)
  # M=-a[[3]]+mean(a[[1]],na.rm=T)
  #ME=a[[1]]-a[[3]]
  #write.csv(list(M,E,ME),file=sprintf("~/Dropbox/lab/data sheet/%s/ME.csv",file))
  
  #return(sum$meanphase)
  return(a)
}