```{r}
# Functions
numPts_below_line = function(myVector,slope,x){ # from ROSE
  yPt <- myVector[x]
  b <- yPt-(slope*x)
  xPts <- 1:length(myVector)
  return(sum(myVector<=(xPts*slope+b)))
}

# Libraries
library(ggplot2)
library(ggprism)

# Set working directory
workdir = "./cutll1_superhubs/"

# Load data
setwd(workdir)
peakome_cutll1_ctrl1 = read.csv("peakome_cutll1_ctrl1.csv", stringsAsFactors = F)
```

```{r}
## Identify superhubs

# Prepare data
dat=peakome_cutll1_ctrl1[peakome_cutll1_ctrl1$CUTLL1_ctrl_rep1_HiChIP_H3K27Ac_IA > 0,]
metric="CUTLL1_ctrl_rep1_HiChIP_H3K27Ac_IA"
dat[is.na(dat[,metric]),metric]=0
dat = dat[dat$CUTLL1_ctrl_rep1_ChIP_H3K27Ac > 1,] # H3K27Ac activity > 1 (We want H3K27Ac peaks)
nrow(dat) # 28094
dat = dat[order(dat[,metric],decreasing = T),]
dat$rank=nrow(dat):1

# Compute ROSE-like cutoff
inputVector = dat[,metric]
inputVector = sort(inputVector)
inputVector[inputVector<0]=0 #set those regions with more control than ranking equal to zero
slope = (max(inputVector)-min(inputVector))/length(inputVector) #This is the slope of the line we want to slide. This is the diagonal.
xPt = floor(optimize(numPts_below_line,lower=1,upper=length(inputVector),myVector= inputVector,slope=slope)$minimum) #Find the x-axis point where a line passing through that point has the minimum number of points below it. (ie. tangent)
y_cutoff = inputVector[xPt]
b = y_cutoff-(slope* xPt)
colCut="grey";colCutTan="red"
print(paste0("xPt = ",xPt))
print(paste0("b = ",b))
print(paste0("cutoff = ",y_cutoff))
```

```{r}
# Make Elbow plot
xlim=40000
ylim=500
baseSize=16
force=25
pointSize=0.5
width=5
height=4

plt = ggplot(data = dat,aes(y=dat[,metric], x=rank))+
  geom_point(size=pointSize)+
  ylim(c(0,ylim))+
  xlim(c(0,xlim))+
  xlab("RANK IA (CTRL1)")+
  ylab("IA (CTRL1)")+
  geom_hline(yintercept = y_cutoff,linetype="dotted",color=colCut)+
  geom_abline(slope=slope ,intercept = b,color=colCutTan,size=0.3)+
  scale_color_manual(values=c("black"))+
  theme_prism(base_size = baseSize)+
  theme(legend.position="none")

print(plt)

pdf(paste0("elbow_IA_CUTLL1-ctrl1.pdf"),width=width,height=height)
print(plt)
dev.off()  
```

```{r}
# Save superhubs data
dat_sh = dat[dat[, metric]>=y_cutoff,] # IA > 45 degree cutoff
dat_sh = dat_sh[dat_sh$CUTLL1_ctrl_rep1_HiChIP_H3K27Ac_Hub > 1,] # Hubness > 1 (We want hubs)
nrow(dat_sh) # 2172

write.table(dat_sh[,c(2:4,1,5)],"superhubs_cutll1_ctrl1.bed",row.names = F, col.names=F, sep="\t", quote=F)
```

