############################
### Schmitt_DLFC_H3K4me3 ###
############################
name.peak <- 'hg19_DLFC_H3K4me3.1e5_peaks.narrowPeak.40Kb'
name.FIRE <- 'Schmitt_FIRE_tissue_DLPFC.txt.gz'  # user can change to file name of FIREs

Bin.Size <- 40000
Bin.Range <- 2e6
Bin.Num <- Bin.Range/Bin.Size * 2 + 1

options(scipen=999)

a <- read.table(name.peak, head=T)
FIRE <- read.table(name.FIRE, header = T)
Binall <- read.csv("Schmitt_primary_cohort_FIREscore.csv")  # A full list of all the bin pairs of all the cell lines and tissues in Schmitt et al. (2016). Users can change the file name

###FIRE data
x <- FIRE[,c(1,2,3)]
### randomly select bins
set.seed(123)
id = sample(1:dim(Binall)[1], dim(x)[1])
y <- Binall[id[order(id)],c(1,2,3)]

x<-unique(x)
y<-unique(y)

dim(a)
dim(x)
dim(y)
colnames(x)<-colnames(a)
colnames(y)<-colnames(a)

a$ID <- paste(a$chr, a$start, sep='_')

outx <- matrix(0, nrow=nrow(x), ncol=Bin.Num)
outy <- matrix(0, nrow=nrow(y), ncol=Bin.Num)

for(id in 1:Bin.Num)
{
  xtmp <- paste0("chr", x[,1], "_", x[,2]+Bin.Size*(id-Bin.Range/Bin.Size-1)) 
  ytmp <- paste0("chr", y[,1], "_", y[,2]+Bin.Size*(id-Bin.Range/Bin.Size-1)) 
  
  outx[,id]<-as.numeric( I( xtmp %in% a$ID ) )
  outy[,id]<-as.numeric( I( ytmp %in% a$ID ) )
  
  if(id%%10==0)
  {
    print( c(id, id/Bin.Num) )
  }
}

out <- matrix(0,nrow=Bin.Num,ncol=3)
out[,1] <- apply(outx, 2, mean)
out[,2] <- apply(outy, 2, mean)
out[,3] <- out[,1]/out[,2]

id2=which(!is.infinite(out[,3]))
# barplot
plot( Bin.Size /1000 *( seq(1,Bin.Num,1)-Bin.Range/Bin.Size-1 ), out[,3], 
      col='red', type='l', ylim=c(1,2), main='', lwd=2,
      xlab='Distance from the FIRE bin (Kb)', ylab='Fold enrichment')


############################
### Schmitt_DLFC_H3K27ac ###
############################
name.peak <- 'hg19_DLFC_H3K27ac.1e5_peaks.narrowPeak.40Kb'
name.FIRE <- 'Schmitt_FIRE_tissue_DLPFC.txt.gz'  # user can change to file name of FIREs

Bin.Size <- 40000
Bin.Range <- 2e6
Bin.Num <- Bin.Range/Bin.Size * 2 + 1

options(scipen=999)

a <- read.table(name.peak, head=T)
FIRE <- read.table(name.FIRE, header = F)
Binall <- read.csv("Schmitt_primary_cohort_FIREscore.csv")  # A full list of all the bin pairs of all the cell lines and tissues in Schmitt et al. (2016). Users can change the file name

###FIRE data
x <- FIRE[,c(1,2,3)]
### randomly select bins
set.seed(123)
id = sample(1:dim(Binall)[1], dim(x)[1])
y <- Binall[id[order(id)],c(1,2,3)]

x<-unique(x)
y<-unique(y)

dim(a)
dim(x)
dim(y)
colnames(x)<-colnames(a)
colnames(y)<-colnames(a)

a$ID <- paste(a$chr, a$start, sep='_')

outx <- matrix(0, nrow=nrow(x), ncol=Bin.Num)
outy <- matrix(0, nrow=nrow(y), ncol=Bin.Num)

for(id in 1:Bin.Num)
{
  xtmp <- paste0("chr", x[,1], "_", x[,2]+Bin.Size*(id-Bin.Range/Bin.Size-1)) 
  ytmp <- paste0("chr", y[,1], "_", y[,2]+Bin.Size*(id-Bin.Range/Bin.Size-1)) 
  
  outx[,id]<-as.numeric( I( xtmp %in% a$ID ) )
  outy[,id]<-as.numeric( I( ytmp %in% a$ID ) )
  
  if(id%%10==0)
  {
    print( c(id, id/Bin.Num) )
  }
}

out <- matrix(0,nrow=Bin.Num,ncol=3)
out[,1] <- apply(outx, 2, mean)
out[,2] <- apply(outy, 2, mean)
out[,3] <- out[,1]/out[,2]

id2=which(!is.infinite(out[,3]))
# barplot
plot( Bin.Size /1000 *( seq(1,Bin.Num,1)-Bin.Range/Bin.Size-1 ), out[,3], 
      col='red', type='l', ylim=c(1,3.5), main='', lwd=2,
      xlab='Distance from the FIRE bin (Kb)', ylab='Fold enrichment')