library(data.table)
library(stats)

samples=c("CKp25-NueN+veGH2AX-ve","CKp25-allGH2AX+ve",
          "CKp25-Con","CKp25-Step1","CKp25-Step2",
          "DMSO_Replicate1","DMSO_Replicate2","ETP5uM2h-Recovery36h_Replicate1","ETP5uM2h-Recovery36h_Replicate2","ETP5uM12h-Recovery36h_Replicate1","ETP5uM12h-Recovery36h_Replicate2",
          "X-rayControl_Replicate1","X-rayControl_Replicate2","X-rayRecovery24h_Replicate1","X-rayRecovery24h_Replicate2")
           

h=data.frame(dist=seq(0, 200e6-10000, 10000))
totalreads=NULL

###read mate-pair files and create histogram the paired read distances###
for(sample in samples){
  cat("Sample: ", sample,"\n")
  dat=fread(paste0("./allMatepairs/",sample,"_matepair.txt"), header=T) #read files
  totalreads=c(totalreads,nrow(dat))
  datb=subset(dat, CHR1==CHR2)
  datb$dist=abs(datb$POS1-datb$POS2)                                    # Absolute distances between paired reads
  dath=hist(datb$dist, breaks=seq(0,200e6,10000))

h=cbind(h,dath$counts)

cat("Total reads: ", sum(dath$counts),"\n\n\n")
}

names(h)=c("dist",samples)
totalreads=data.frame(samples,totalreads)

##convert to percentages###
for (i in 2:ncol(h)) h[,i] = h[,i]/totalreads$totalreads[i-1] * 100

write.table(h,"Mate-pair percentage reads vs distance.txt", row.names=F, sep="\t", quote=F)


###############Mate pair distance plots#############

h=read.delim("Mate-pair percentage reads vs distance.txt", header=T)
h=h[,c(1:8,11:16,19:20)]

##subset distances##
hp=subset(h,h$dist>=200000 & h$dist<1.5e6)

##Loess smooth##
L=NULL
for(i in 2:ncol(hp)){
  l=loess(hp[,i]~hp$dist,span=0.3)
  L=cbind(L,l$fitted)
}

h.lp=cbind(hp$dist,L)
h.lp=as.data.frame(h.lp)

names(h.lp)=names(hp)


##################plots#####################

########Figure 1
myxlim=c(200000,1.5e6); myylim=c(0,0.0045)
plot(hp$'CKp25-Con'~hp$dist,ylim=myylim,xlim=myxlim,type="p",pch=19,cex=0.5)
lines(h.lp$'CKp25-Con'~h.lp$dist,lwd=4)

myred <- rgb(1, 0, 0, alpha = 0.7)
points(hp$'CKp25-Step1'~hp$dist,pch=19,cex=0.5,col=myred)
lines(h.lp$'CKp25-Step1'~h.lp$dist,col=myred,lwd=4)

points(hp$'CKp25-Step2'~hp$dist,pch=19,cex=0.5,col="red4")
lines(h.lp$'CKp25-Step2'~h.lp$dist,col="red4",lwd=4)


##Etoposide##Figure 1
myxlim=c(200000,1.5e6); myylim=c(0,0.0045)
plot((hp$'DMSO_Replicate1'+ hp$'DMSO_Replicate2')/2 ~hp$dist,ylim=myylim,xlim=myxlim,type="p",pch=19,cex=0.5)
lines((h.lp$'DMSO_Replicate1'+ h.lp$'DMSO_Replicate2')/2 ~hp$dist,lwd=4)

lightblue <- rgb(173/255, 216/255, 230/255, alpha = 0.7)
points((hp$'ETP5uM2h-Recovery36h_Replicate1'+hp$'ETP5uM2h-Recovery36h_Replicate2')/2 ~hp$dist,pch=19,cex=0.5, col=lightblue)
lines((h.lp$'ETP5uM2h-Recovery36h_Replicate1'+h.lp$'ETP5uM2h-Recovery36h_Replicate2')/2 ~h.lp$dist,col=lightblue,lwd=4)

points((hp$'ETP5uM12h-Recovery36h_Replicate1'+hp$'ETP5uM12h-Recovery36h_Replicate2')/2 ~hp$dist,pch=19,cex=0.5, col="blue")
lines((h.lp$'ETP5uM12h-Recovery36h_Replicate1'+h.lp$'ETP5uM12h-Recovery36h_Replicate2')/2 ~h.lp$dist,col="blue",lwd=4)




##X-ray and all yH2AX Supplementary###

##CKp25###Supplementary Figure 2##
myxlim=c(200000,1.5e6); myylim=c(0,0.0015)

plot(hp$'CKp25-NueN+veGH2AX-ve'~hp$dist,ylim=myylim,xlim=myxlim,type="p",pch=19,cex=0.5)
lines(h.lp$'CKp25-NueN+veGH2AX-ve'~h.lp$dist,lwd=4)

points(hp$'CKp25-allGH2AX+ve'~hp$dist,pch=19,cex=0.5,col="red")
lines(h.lp$'CKp25-allGH2AX+ve'~h.lp$dist,col="red",lwd=4)


###X-ray####

plot( (hp$'X-rayControl_Replicate1'+hp$'X-rayControl_Replicate2')/2 ~ hp$dist,ylim=myylim,xlim=myxlim,type="p",pch=19,cex=0.5)
lines((h.lp$'X-rayControl_Replicate1' + h.lp$'X-rayControl_Replicate2')/2 ~h.lp$dist,lwd=4)

points((hp$'X-rayRecovery24h_Replicate1'+hp$'X-rayRecovery24h_Replicate2')/2 ~hp$dist,pch=19,cex=0.5, col="blue")
lines((h.lp$'X-rayRecovery24h_Replicate1'+h.lp$'X-rayRecovery24h_Replicate2')/2 ~h.lp$dist,col="blue",lwd=4)