4DN

library(DT)
library(DOSE)
library(vroom)
library(WGCNA)
library(readr)
library(dplyr)
library(tidyr)
library(qqman)
library(GO.db)
library(caret)
library(fgsea)
library(VplotR)
library(tibble)
library(ChIPQC)
library(GOplot)
library(DESeq2)
library(UpSetR)
library(viridis)
library(ggrepel)
library(forcats)
library(ggplot2)
library(biomaRt)
library(msigdbr)
library(viridis)
library(CEMiTool)
library(pheatmap)
library(tximport)
library(corrplot)
library(GOSemSim)
library(circlize)
library(patchwork)
library(tidyverse)
library(DMRcaller)
library(networkD3)
library(ATACseqQC)
library(hrbrthemes)
library(ChIPseeker)
library(ReactomePA)
library(enrichplot)
library(data.table)
library(multiWGCNA)
library(flashClust)
library(NucleoATACR)
library(rtracklayer)
library(tximportData)
library(directlabels)
library(BiocParallel)
library(org.Mm.eg.db)
library(org.Hs.eg.db)
library(RColorBrewer)
library(ExperimentHub)
library(GenomicRanges)
library(dynamicTreeCut)
library(clusterProfiler)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)

FUNCTION

cluster_p300_viewer = function(num_cluster){

subdata = as.data.frame(UMAP_DIFF_correlation_50_2_1_1[which(UMAP_DIFF_correlation_50_2_1_1$cluster == num_cluster),])

Ets1_Q = c("Q1","Q2","Q3","Q4","Q5","Q6","Q7","Q8","Q9","Q10")

Ets1_Q = as.data.frame(Ets1_Q)
colnames(Ets1_Q) = "Ets1_Q"
Ets1_Q$score = 0

for(i in 1:nrow(Ets1_Q)){
Ets1_Q[i,2] = nrow(subdata[which(subdata$Ets1_Q == Ets1_Q[i,1]),])
    }

subdata = as.data.frame(UMAP_DIFF_correlation_50_2_1_1[which(UMAP_DIFF_correlation_50_2_1_1$cluster == num_cluster),])

p300_Q = c("Q1","Q2","Q3","Q4","Q5","Q6","Q7","Q8","Q9","Q10")

p300_Q = as.data.frame(p300_Q)
colnames(p300_Q) = "p300_Q"
p300_Q$score = 0

for(i in 1:nrow(p300_Q)){
p300_Q[i,2] = nrow(subdata[which(subdata$p300_Q == p300_Q[i,1]),])
    }

ddata= cbind(Ets1_Q,p300_Q)

ddata$label = c("A","B","C","D","E","F","G","H","I","J")

data <- data.frame(
    category=ddata[,5],
    count=ddata[,4])
 
data$fraction <- data$count / sum(data$count)
data$ymax <- cumsum(data$fraction)
data$ymin <- c(0, head(data$ymax, n=-1))
data$labelPosition <- (data$ymax + data$ymin) / 2


data$label <- paste0(data$category, "\n value: ", data$count)
nb.cols <- 19
mycolors <- colorRampPalette(brewer.pal(8, "Pastel2"))(nb.cols)


ggplot(data, aes(ymax=ymax, ymin=ymin, xmax=4, xmin=3, fill=category)) +
geom_rect() +
scale_color_manual(values=c("grey40"))+
scale_fill_manual(values=c("#E0F2F1","#B2DFDB","#80CBC4","#4DB6AC","#26A69A",
                           "#009688","#00897B","#00796B","#00695C","#004D40"))+
coord_polar(theta="y") +
  coord_polar(theta="y") +
  xlim(c(2.5, 4)) +
  theme_void() +
scale_colour_hue(l = 20, c = 100)+
  theme(legend.position = "none")
    }

cluster_Ets1_viewer = function(num_cluster){

subdata = as.data.frame(UMAP_DIFF_correlation_50_2_1_1[which(UMAP_DIFF_correlation_50_2_1_1$cluster == num_cluster),])

Ets1_Q = c("Q1","Q2","Q3","Q4","Q5","Q6","Q7","Q8","Q9","Q10")

Ets1_Q = as.data.frame(Ets1_Q)
colnames(Ets1_Q) = "Ets1_Q"
Ets1_Q$score = 0

for(i in 1:nrow(Ets1_Q)){
Ets1_Q[i,2] = nrow(subdata[which(subdata$Ets1_Q == Ets1_Q[i,1]),])
    }

subdata = as.data.frame(UMAP_DIFF_correlation_50_2_1_1[which(UMAP_DIFF_correlation_50_2_1_1$cluster == num_cluster),])

p300_Q = c("Q1","Q2","Q3","Q4","Q5","Q6","Q7","Q8","Q9","Q10")

p300_Q = as.data.frame(p300_Q)
colnames(p300_Q) = "p300_Q"
p300_Q$score = 0

for(i in 1:nrow(p300_Q)){
p300_Q[i,2] = nrow(subdata[which(subdata$p300_Q == p300_Q[i,1]),])
    }

ddata= cbind(Ets1_Q,p300_Q)

ddata$label = c("A","B","C","D","E","F","G","H","I","J")

options(repr.plot.width = 8, repr.plot.height = 8, repr.plot.res = 1000, repr.plot.pointsize = 10)

data <- data.frame(
    category=ddata[,5],
    count=ddata[,2])
 
data$fraction <- data$count / sum(data$count)
data$ymax <- cumsum(data$fraction)
data$ymin <- c(0, head(data$ymax, n=-1))
data$labelPosition <- (data$ymax + data$ymin) / 2


data$label <- paste0(data$category, "\n value: ", data$count)
nb.cols <- 19
mycolors <- colorRampPalette(brewer.pal(8, "Pastel2"))(nb.cols)


ggplot(data, aes(ymax=ymax, ymin=ymin, xmax=4, xmin=3, fill=category)) +
geom_rect() +
scale_color_manual(values=c("grey40"))+
scale_fill_manual(values=c("#FFF8E1","#FFECB3","#FFE082","#FFD54F","#FFCA28",
                           "#FFC107","#FFB300","#FFA000","#FF8F00","#FF6F00"))+
coord_polar(theta="y") +
  coord_polar(theta="y") +
  xlim(c(2.5, 4)) +
  theme_void() +
scale_colour_hue(l = 20, c = 100)+
  theme(legend.position = "none")
    }

UMAP_viewer = function(inputPATH){

UMAP = read.table(inputPATH,sep = "\t", header=T)

options(repr.plot.width = 8, repr.plot.height = 8, repr.plot.res = 200, repr.plot.pointsize = 30)

color_palette <- c("darkblue", "#B71C1C","#283593", "#3949AB", "#5C6BC0", "#9FA8DA",
                   "#757575", "#FFCDD2", "#E57373", "#E53935")

ggplot(UMAP, aes(X, Y, color = factor(DIFF_Q))) +
  geom_point(size = 1, alpha = 1) +
  scale_color_manual(values = color_palette) +
  theme_classic(base_size = 10) +
  ggtitle(NULL) +
  xlab("UMAP1") +
  ylab("UMAP2") +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5, size = 25, color = "darkblue"),
    axis.text.x = element_text(face = "bold", size = 20),
    axis.text.y = element_text(face = "bold", size = 20),
    axis.title.x = element_text(face = "bold", size = 20, color = "darkblue"),
    axis.title.y = element_text(face = "bold", size = 20, color = "darkblue"),
    legend.position = "none",
    legend.title = element_text(face = "bold", size = 10),
    legend.text = element_text(face = "bold", size = 7)
  )
    }

UMAP_VIZ  = function(X){
    
options(repr.plot.width = 8, repr.plot.height = 8, repr.plot.res = 100, repr.plot.pointsize = 30)
    
    umap_DF = as.data.frame(vroom(X, delim = "\t",show_col_types = FALSE))
    umap_DF = umap_DF[,c(1,2,24628:24637)]

ggplot(umap_DF, aes(x, y)) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q2"), color="#283593", size=1, alpha = 0.5) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q3"), color="#3949AB", size=1, alpha = 0.5) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q4"), color="#5C6BC0", size=1, alpha = 0.5) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q5"), color="#9FA8DA", size=1, alpha = 0.5) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q6"), color="#757575", size=1, alpha = 0.5) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q7"), color="#FFCDD2", size=1, alpha = 0.5) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q8"), color="#E57373", size=1, alpha = 0.5) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q9"), color="#E53935", size=1, alpha = 0.5) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q10"), color="#B71C1C", size=1, alpha = 0.5) +
geom_point(data=subset(umap_DF,DIFF_Q=="Q1"), color="darkblue", size=1, alpha = 0.5) +
theme_classic(base_size = 10) +
# xlim(-10,10)+
# ylim(-10,10)+
ggtitle(NULL)+
xlab("UMAP1") +
ylab("UMAP2") +
theme(plot.title=element_text(face="bold",hjust=0.5,size=25,color = "darkblue"),
      axis.text.x=element_text(face="bold",size=20),
      axis.text.y=element_text(face="bold",size=20), 
      axis.title.x = element_text(face="bold",size = 20,color = "darkblue"),
      axis.title.y = element_text(face="bold",size = 20,color = "darkblue"),
      legend.position="right", 
      legend.title=element_text(face="bold",size=10), 
      legend.text=element_text(face="bold",size=7))
}

p300_to_Ets1_Q  = function(p300_Q){

    options(repr.plot.width = 5, repr.plot.height = 5, repr.plot.res = 200, repr.plot.pointsize = 30)

p300_Q10 = BIN_DATA[which(BIN_DATA$p300_Q == p300_Q),]

A = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q1"),])/nrow(p300_Q10)*100,2)
B = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q2"),])/nrow(p300_Q10)*100,2)
C = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q3"),])/nrow(p300_Q10)*100,2)
D = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q4"),])/nrow(p300_Q10)*100,2)
E = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q5"),])/nrow(p300_Q10)*100,2)
F = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q6"),])/nrow(p300_Q10)*100,2)
G = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q7"),])/nrow(p300_Q10)*100,2)
H = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q8"),])/nrow(p300_Q10)*100,2)
I = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q9"),])/nrow(p300_Q10)*100,2)
J = round(nrow(p300_Q10[which(p300_Q10$Ets1_Q == "Q10"),])/nrow(p300_Q10)*100,2)

data <- data.frame(
    category=c("A","B","C","D","E","F","G","H","I","J"),
    count=c(A,B,C,D,E,F,G,H,I,J)
)
 
data$fraction <- data$count / sum(data$count)
data$ymax <- cumsum(data$fraction)
data$ymin <- c(0, head(data$ymax, n=-1))
data$labelPosition <- (data$ymax + data$ymin) / 2
data$label <- paste0(data$category, "\n value: ", data$count)


ggplot(data, aes(ymax=ymax, ymin=ymin, xmax=4, xmin=3, fill=c("#004D40","#00695C","#00796B","#00897B","#009688",
                                                              "#26A69A","#4DB6AC","#80CBC4","#B2DFDB","#E0F2F1"))) +
geom_rect() +
scale_fill_manual(values=c("#E0F2F1","#B2DFDB","#80CBC4","#4DB6AC","#26A69A",
                           "#009688","#00897B","#00796B","#00695C","#004D40"))+
coord_polar(theta="y") +
xlim(c(2.5, 4)) +
theme_void() +
theme(legend.position = "none")
    }

Ets1_to_p300_Q  = function(Ets1_Q){

    options(repr.plot.width = 5, repr.plot.height = 5, repr.plot.res = 200, repr.plot.pointsize = 30)

Ets1_Q10 = BIN_DATA[which(BIN_DATA$Ets1_Q == Ets1_Q),]

A = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q1"),])/nrow(Ets1_Q10)*100,2)
B = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q2"),])/nrow(Ets1_Q10)*100,2)
C = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q3"),])/nrow(Ets1_Q10)*100,2)
D = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q4"),])/nrow(Ets1_Q10)*100,2)
E = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q5"),])/nrow(Ets1_Q10)*100,2)
F = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q6"),])/nrow(Ets1_Q10)*100,2)
G = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q7"),])/nrow(Ets1_Q10)*100,2)
H = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q8"),])/nrow(Ets1_Q10)*100,2)
I = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q9"),])/nrow(Ets1_Q10)*100,2)
J = round(nrow(Ets1_Q10[which(Ets1_Q10$p300_Q == "Q10"),])/nrow(Ets1_Q10)*100,2)

data <- data.frame(
    category=c("A","B","C","D","E","F","G","H","I","J"),
    count=c(A,B,C,D,E,F,G,H,I,J)
)
 
data$fraction <- data$count / sum(data$count)
data$ymax <- cumsum(data$fraction)
data$ymin <- c(0, head(data$ymax, n=-1))
data$labelPosition <- (data$ymax + data$ymin) / 2
data$label <- paste0(data$category, "\n value: ", data$count)


ggplot(data, aes(ymax=ymax, ymin=ymin, xmax=4, xmin=3, fill=c("#004D40","#00695C","#00796B","#00897B","#009688",
                                                              "#26A69A","#4DB6AC","#80CBC4","#B2DFDB","#E0F2F1"))) +
geom_rect() +
scale_fill_manual(values=c("#FFF8E1","#FFECB3","#FFE082","#FFD54F","#FFCA28",
                           "#FFC107","#FFB300","#FFA000","#FF8F00","#FF6F00"))+
coord_polar(theta="y") +
xlim(c(2.5, 4)) +
theme_void() +
theme(legend.position = "none")
    }

p300_Ets1_toDIFF  = function(Qlevel){

p300_Q10 = BIN_DATA[which(BIN_DATA$p300_Q == Qlevel),]
p300_Ets1_Q10 = p300_Q10[which(p300_Q10$Ets1_Q == Qlevel),]

A = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q1"),])/nrow(p300_Ets1_Q10)*100,2)
B = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q2"),])/nrow(p300_Ets1_Q10)*100,2)
C = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q3"),])/nrow(p300_Ets1_Q10)*100,2)
D = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q4"),])/nrow(p300_Ets1_Q10)*100,2)
E = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q5"),])/nrow(p300_Ets1_Q10)*100,2)
F = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q6"),])/nrow(p300_Ets1_Q10)*100,2)
G = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q7"),])/nrow(p300_Ets1_Q10)*100,2)
H = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q8"),])/nrow(p300_Ets1_Q10)*100,2)
I = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q9"),])/nrow(p300_Ets1_Q10)*100,2)
J = round(nrow(p300_Ets1_Q10[which(p300_Ets1_Q10$DIFF_Q == "Q10"),])/nrow(p300_Ets1_Q10)*100,2)

data <- data.frame(
    category=c("A","B","C","D","E","F","G","H","I","J"),
    count=c(A,B,C,D,E,F,G,H,I,J)
)
 
data$fraction <- data$count / sum(data$count)
data$ymax <- cumsum(data$fraction)
data$ymin <- c(0, head(data$ymax, n=-1))
data$labelPosition <- (data$ymax + data$ymin) / 2
data$label <- paste0(data$category, "\n value: ", data$count)


ggplot(data, aes(ymax=ymax, ymin=ymin, xmax=4, xmin=3, fill=c("#004D40","#00695C","#00796B","#00897B","#009688",
                                                              "#26A69A","#4DB6AC","#80CBC4","#B2DFDB","#E0F2F1"))) +
geom_rect() +
scale_fill_manual(values=c("darkblue","#283593","#3949AB","#5C6BC0","#9FA8DA",      
                           "#E0E0E0","#FFCDD2","#E57373","#E53935","#B71C1C"))+
coord_polar(theta="y") +
xlim(c(2.5, 4)) +
theme_void() +
theme(legend.position = "none")
    }

ATACtoBINvio_0_stat = function(DAY_para,HISTONE_para){
path1 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/AtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path2 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/AtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")
path3 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/BtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path4 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/BtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")

path5 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/partially_AtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path6 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/partially_AtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")
path7 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/partially_BtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path8 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/partially_BtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")

AtoA_DAY0_H3K27ac = read.table(path1,sep = "\t", header=F,fill = TRUE)
AtoB_DAY0_H3K27ac = read.table(path2,sep = "\t", header=F,fill = TRUE)
BtoA_DAY0_H3K27ac = read.table(path3,sep = "\t", header=F,fill = TRUE)
BtoB_DAY0_H3K27ac = read.table(path4,sep = "\t", header=F,fill = TRUE)

partially_AtoA_DAY0_H3K27ac = read.table(path5,sep = "\t", header=F,fill = TRUE)
partially_AtoB_DAY0_H3K27ac = read.table(path6,sep = "\t", header=F,fill = TRUE)
partially_BtoA_DAY0_H3K27ac = read.table(path7,sep = "\t", header=F,fill = TRUE)
partially_BtoB_DAY0_H3K27ac = read.table(path8,sep = "\t", header=F,fill = TRUE)

AtoA_DAY0_H3K27ac_forP = rbind(AtoA_DAY0_H3K27ac,partially_AtoA_DAY0_H3K27ac)
AtoB_DAY0_H3K27ac_forP = rbind(AtoB_DAY0_H3K27ac,partially_AtoB_DAY0_H3K27ac)
BtoA_DAY0_H3K27ac_forP = rbind(BtoA_DAY0_H3K27ac,partially_BtoA_DAY0_H3K27ac)
BtoB_DAY0_H3K27ac_forP = rbind(BtoB_DAY0_H3K27ac,partially_BtoB_DAY0_H3K27ac)

AtoA_DAY0_H3K27ac_forP = AtoA_DAY0_H3K27ac_forP[order(-AtoA_DAY0_H3K27ac_forP[,7]),]
AtoB_DAY0_H3K27ac_forP = AtoB_DAY0_H3K27ac_forP[order(-AtoB_DAY0_H3K27ac_forP[,7]),]
BtoA_DAY0_H3K27ac_forP = BtoA_DAY0_H3K27ac_forP[order(-BtoA_DAY0_H3K27ac_forP[,7]),]
BtoB_DAY0_H3K27ac_forP = BtoB_DAY0_H3K27ac_forP[order(-BtoB_DAY0_H3K27ac_forP[,7]),]

AtoA_DAY0_H3K27ac_fit = AtoA_DAY0_H3K27ac_forP[,5:7]
AtoB_DAY0_H3K27ac_fit = AtoB_DAY0_H3K27ac_forP[,5:7]
BtoA_DAY0_H3K27ac_fit = BtoA_DAY0_H3K27ac_forP[,5:7]
BtoB_DAY0_H3K27ac_fit = BtoB_DAY0_H3K27ac_forP[,5:7]

A_DAY0_H3K27ac_fit = rbind(AtoA_DAY0_H3K27ac_fit, AtoB_DAY0_H3K27ac_fit)
B_DAY0_H3K27ac_fit = rbind(BtoA_DAY0_H3K27ac_fit, BtoB_DAY0_H3K27ac_fit)

A_DAY0_H3K27ac_fit[,1] = "B"
B_DAY0_H3K27ac_fit[,1] = "A"

A_DAY0_H3K27ac_fit = A_DAY0_H3K27ac_fit[order(-A_DAY0_H3K27ac_fit[,3]),]
B_DAY0_H3K27ac_fit = B_DAY0_H3K27ac_fit[order(-B_DAY0_H3K27ac_fit[,3]),]

H3K27ac_plot = rbind(A_DAY0_H3K27ac_fit,B_DAY0_H3K27ac_fit)

colnames(H3K27ac_plot) = c("comp","rawvalue","value")

H3K27ac_plot$z = (H3K27ac_plot$value-mean(H3K27ac_plot$value))/sd(H3K27ac_plot$value)

H3K27ac_plot = H3K27ac_plot[which(H3K27ac_plot$z >= -3),]
H3K27ac_plot = H3K27ac_plot[which(H3K27ac_plot$z <= 3),]

wilcox.test(
    H3K27ac_plot[which(H3K27ac_plot[,1] == "A"),4],
    H3K27ac_plot[which(H3K27ac_plot[,1] == "B"),4], 
    alternative = "two.sided"
)
    }

ATACtoBINvio_4_stat = function(DAY_para,HISTONE_para){
path1 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/AtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path2 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/AtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")
path3 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/BtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path4 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/BtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")

path5 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/partially_AtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path6 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/partially_AtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")
path7 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/partially_BtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path8 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/partially_BtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")

AtoA_DAY0_H3K27ac = read.table(path1,sep = "\t", header=F,fill = TRUE)
AtoB_DAY0_H3K27ac = read.table(path2,sep = "\t", header=F,fill = TRUE)
BtoA_DAY0_H3K27ac = read.table(path3,sep = "\t", header=F,fill = TRUE)
BtoB_DAY0_H3K27ac = read.table(path4,sep = "\t", header=F,fill = TRUE)

partially_AtoA_DAY0_H3K27ac = read.table(path5,sep = "\t", header=F,fill = TRUE)
partially_AtoB_DAY0_H3K27ac = read.table(path6,sep = "\t", header=F,fill = TRUE)
partially_BtoA_DAY0_H3K27ac = read.table(path7,sep = "\t", header=F,fill = TRUE)
partially_BtoB_DAY0_H3K27ac = read.table(path8,sep = "\t", header=F,fill = TRUE)

AtoA_DAY0_H3K27ac_forP = rbind(AtoA_DAY0_H3K27ac,partially_AtoA_DAY0_H3K27ac)
AtoB_DAY0_H3K27ac_forP = rbind(AtoB_DAY0_H3K27ac,partially_AtoB_DAY0_H3K27ac)
BtoA_DAY0_H3K27ac_forP = rbind(BtoA_DAY0_H3K27ac,partially_BtoA_DAY0_H3K27ac)
BtoB_DAY0_H3K27ac_forP = rbind(BtoB_DAY0_H3K27ac,partially_BtoB_DAY0_H3K27ac)

AtoA_DAY0_H3K27ac_forP = AtoA_DAY0_H3K27ac_forP[order(-AtoA_DAY0_H3K27ac_forP[,7]),]
AtoB_DAY0_H3K27ac_forP = AtoB_DAY0_H3K27ac_forP[order(-AtoB_DAY0_H3K27ac_forP[,7]),]
BtoA_DAY0_H3K27ac_forP = BtoA_DAY0_H3K27ac_forP[order(-BtoA_DAY0_H3K27ac_forP[,7]),]
BtoB_DAY0_H3K27ac_forP = BtoB_DAY0_H3K27ac_forP[order(-BtoB_DAY0_H3K27ac_forP[,7]),]

AtoA_DAY0_H3K27ac_fit = AtoA_DAY0_H3K27ac_forP[,5:7]
AtoB_DAY0_H3K27ac_fit = AtoB_DAY0_H3K27ac_forP[,5:7]
BtoA_DAY0_H3K27ac_fit = BtoA_DAY0_H3K27ac_forP[,5:7]
BtoB_DAY0_H3K27ac_fit = BtoB_DAY0_H3K27ac_forP[,5:7]

A_DAY0_H3K27ac_fit = rbind(AtoA_DAY0_H3K27ac_fit, AtoB_DAY0_H3K27ac_fit)
B_DAY0_H3K27ac_fit = rbind(BtoA_DAY0_H3K27ac_fit, BtoB_DAY0_H3K27ac_fit)

A_DAY0_H3K27ac_fit[,1] = "B"
B_DAY0_H3K27ac_fit[,1] = "A"

A_DAY0_H3K27ac_fit = A_DAY0_H3K27ac_fit[order(-A_DAY0_H3K27ac_fit[,3]),]
B_DAY0_H3K27ac_fit = B_DAY0_H3K27ac_fit[order(-B_DAY0_H3K27ac_fit[,3]),]

H3K27ac_plot = rbind(A_DAY0_H3K27ac_fit,B_DAY0_H3K27ac_fit)

colnames(H3K27ac_plot) = c("comp","rawvalue","value")

H3K27ac_plot$z = (H3K27ac_plot$value-mean(H3K27ac_plot$value))/sd(H3K27ac_plot$value)

H3K27ac_plot = H3K27ac_plot[which(H3K27ac_plot$z >= -3),]
H3K27ac_plot = H3K27ac_plot[which(H3K27ac_plot$z <= 3),]

wilcox.test(
    H3K27ac_plot[which(H3K27ac_plot[,1] == "A"),4],
    H3K27ac_plot[which(H3K27ac_plot[,1] == "B"),4], 
    alternative = "two.sided"
)
    }

ATACtoBINvio_4 = function(DAY_para,HISTONE_para){
path1 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/AtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path2 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/AtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")
path3 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/BtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path4 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/BtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")

path5 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/partially_AtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path6 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/partially_AtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")
path7 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/partially_BtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path8 = paste("/data3/psg/4DN_20241106/FiG23/ABcompartment/partially_BtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")

AtoA_DAY0_H3K27ac = read.table(path1,sep = "\t", header=F,fill = TRUE)
AtoB_DAY0_H3K27ac = read.table(path2,sep = "\t", header=F,fill = TRUE)
BtoA_DAY0_H3K27ac = read.table(path3,sep = "\t", header=F,fill = TRUE)
BtoB_DAY0_H3K27ac = read.table(path4,sep = "\t", header=F,fill = TRUE)

partially_AtoA_DAY0_H3K27ac = read.table(path5,sep = "\t", header=F,fill = TRUE)
partially_AtoB_DAY0_H3K27ac = read.table(path6,sep = "\t", header=F,fill = TRUE)
partially_BtoA_DAY0_H3K27ac = read.table(path7,sep = "\t", header=F,fill = TRUE)
partially_BtoB_DAY0_H3K27ac = read.table(path8,sep = "\t", header=F,fill = TRUE)

AtoA_DAY0_H3K27ac_forP = rbind(AtoA_DAY0_H3K27ac,partially_AtoA_DAY0_H3K27ac)
AtoB_DAY0_H3K27ac_forP = rbind(AtoB_DAY0_H3K27ac,partially_AtoB_DAY0_H3K27ac)
BtoA_DAY0_H3K27ac_forP = rbind(BtoA_DAY0_H3K27ac,partially_BtoA_DAY0_H3K27ac)
BtoB_DAY0_H3K27ac_forP = rbind(BtoB_DAY0_H3K27ac,partially_BtoB_DAY0_H3K27ac)

AtoA_DAY0_H3K27ac_forP = AtoA_DAY0_H3K27ac_forP[order(-AtoA_DAY0_H3K27ac_forP[,7]),]
AtoB_DAY0_H3K27ac_forP = AtoB_DAY0_H3K27ac_forP[order(-AtoB_DAY0_H3K27ac_forP[,7]),]
BtoA_DAY0_H3K27ac_forP = BtoA_DAY0_H3K27ac_forP[order(-BtoA_DAY0_H3K27ac_forP[,7]),]
BtoB_DAY0_H3K27ac_forP = BtoB_DAY0_H3K27ac_forP[order(-BtoB_DAY0_H3K27ac_forP[,7]),]

AtoA_DAY0_H3K27ac_fit = AtoA_DAY0_H3K27ac_forP[,5:7]
AtoB_DAY0_H3K27ac_fit = AtoB_DAY0_H3K27ac_forP[,5:7]
BtoA_DAY0_H3K27ac_fit = BtoA_DAY0_H3K27ac_forP[,5:7]
BtoB_DAY0_H3K27ac_fit = BtoB_DAY0_H3K27ac_forP[,5:7]

A_DAY0_H3K27ac_fit = rbind(AtoA_DAY0_H3K27ac_fit, AtoB_DAY0_H3K27ac_fit)
B_DAY0_H3K27ac_fit = rbind(BtoA_DAY0_H3K27ac_fit, BtoB_DAY0_H3K27ac_fit)

A_DAY0_H3K27ac_fit[,1] = "B"
B_DAY0_H3K27ac_fit[,1] = "A"

A_DAY0_H3K27ac_fit = A_DAY0_H3K27ac_fit[order(-A_DAY0_H3K27ac_fit[,3]),]
B_DAY0_H3K27ac_fit = B_DAY0_H3K27ac_fit[order(-B_DAY0_H3K27ac_fit[,3]),]

H3K27ac_plot = rbind(A_DAY0_H3K27ac_fit,B_DAY0_H3K27ac_fit)

colnames(H3K27ac_plot) = c("comp","rawvalue","value")

H3K27ac_plot$z = (H3K27ac_plot$value-mean(H3K27ac_plot$value))/sd(H3K27ac_plot$value)

H3K27ac_plot = H3K27ac_plot[which(H3K27ac_plot$z >= -3),]
H3K27ac_plot = H3K27ac_plot[which(H3K27ac_plot$z <= 3),]

options(repr.plot.width = 3, repr.plot.height = 4, repr.plot.res = 1000, repr.plot.pointsize = 20)

ggplot(H3K27ac_plot, aes(x=factor(comp), y=z, fill=factor(comp)))+
geom_boxplot(fill=c("dodgerblue3","yellow3"),
             size=1,
             color='black',
             width=0.5,
             alpha=0.7)+
geom_violin(adjust=1,
            scale='count',
            fill=c("black"),
            alpha=0.5)+
    guides(fill="none")+
theme_minimal()+
    labs(title=HISTONE_para,
        subtitle = "***\n")+
    xlab("")+
    ylab("enrichment of read count\n(z-score normalization)\n")+
theme(plot.title=element_text(face="bold.italic",hjust=0.5,size=25,color = "black"),
      axis.text.x=element_text(face="bold",size=12,color = "black"),
      axis.text.y=element_text(face="bold",size=10),
      plot.subtitle=element_text(face="bold",hjust=0.5,vjust=-3,size=15,color = "black"),
      axis.title.x = element_text(size=0,color = "black"),
      axis.title.y = element_text(face="bold",size=15,color = "black"))
    }

ATACtoBINvio_0 = function(DAY_para,HISTONE_para){
path1 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/AtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path2 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/AtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")
path3 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/BtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path4 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/BtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")

path5 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/partially_AtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path6 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/partially_AtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")
path7 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/partially_BtoA_",DAY_para,"_",HISTONE_para,".txt",sep="")
path8 = paste("/data3/psg/4DN_20241106/FiG23/",HISTONE_para,"_to_HiC/partially_BtoB_",DAY_para,"_",HISTONE_para,".txt",sep="")

AtoA_DAY0_H3K27ac = read.table(path1,sep = "\t", header=F,fill = TRUE)
AtoB_DAY0_H3K27ac = read.table(path2,sep = "\t", header=F,fill = TRUE)
BtoA_DAY0_H3K27ac = read.table(path3,sep = "\t", header=F,fill = TRUE)
BtoB_DAY0_H3K27ac = read.table(path4,sep = "\t", header=F,fill = TRUE)

partially_AtoA_DAY0_H3K27ac = read.table(path5,sep = "\t", header=F,fill = TRUE)
partially_AtoB_DAY0_H3K27ac = read.table(path6,sep = "\t", header=F,fill = TRUE)
partially_BtoA_DAY0_H3K27ac = read.table(path7,sep = "\t", header=F,fill = TRUE)
partially_BtoB_DAY0_H3K27ac = read.table(path8,sep = "\t", header=F,fill = TRUE)

AtoA_DAY0_H3K27ac_forP = rbind(AtoA_DAY0_H3K27ac,partially_AtoA_DAY0_H3K27ac)
AtoB_DAY0_H3K27ac_forP = rbind(AtoB_DAY0_H3K27ac,partially_AtoB_DAY0_H3K27ac)
BtoA_DAY0_H3K27ac_forP = rbind(BtoA_DAY0_H3K27ac,partially_BtoA_DAY0_H3K27ac)
BtoB_DAY0_H3K27ac_forP = rbind(BtoB_DAY0_H3K27ac,partially_BtoB_DAY0_H3K27ac)

AtoA_DAY0_H3K27ac_forP = AtoA_DAY0_H3K27ac_forP[order(-AtoA_DAY0_H3K27ac_forP[,7]),]
AtoB_DAY0_H3K27ac_forP = AtoB_DAY0_H3K27ac_forP[order(-AtoB_DAY0_H3K27ac_forP[,7]),]
BtoA_DAY0_H3K27ac_forP = BtoA_DAY0_H3K27ac_forP[order(-BtoA_DAY0_H3K27ac_forP[,7]),]
BtoB_DAY0_H3K27ac_forP = BtoB_DAY0_H3K27ac_forP[order(-BtoB_DAY0_H3K27ac_forP[,7]),]

AtoA_DAY0_H3K27ac_fit = AtoA_DAY0_H3K27ac_forP[,5:7]
AtoB_DAY0_H3K27ac_fit = AtoB_DAY0_H3K27ac_forP[,5:7]
BtoA_DAY0_H3K27ac_fit = BtoA_DAY0_H3K27ac_forP[,5:7]
BtoB_DAY0_H3K27ac_fit = BtoB_DAY0_H3K27ac_forP[,5:7]

A_DAY0_H3K27ac_fit = rbind(AtoA_DAY0_H3K27ac_fit, AtoB_DAY0_H3K27ac_fit)
B_DAY0_H3K27ac_fit = rbind(BtoA_DAY0_H3K27ac_fit, BtoB_DAY0_H3K27ac_fit)

A_DAY0_H3K27ac_fit[,1] = "B"
B_DAY0_H3K27ac_fit[,1] = "A"

A_DAY0_H3K27ac_fit = A_DAY0_H3K27ac_fit[order(-A_DAY0_H3K27ac_fit[,3]),]
B_DAY0_H3K27ac_fit = B_DAY0_H3K27ac_fit[order(-B_DAY0_H3K27ac_fit[,3]),]

H3K27ac_plot = rbind(A_DAY0_H3K27ac_fit,B_DAY0_H3K27ac_fit)

colnames(H3K27ac_plot) = c("comp","rawvalue","value")

H3K27ac_plot$z = (H3K27ac_plot$value-mean(H3K27ac_plot$value))/sd(H3K27ac_plot$value)

H3K27ac_plot = H3K27ac_plot[which(H3K27ac_plot$z >= -3),]
H3K27ac_plot = H3K27ac_plot[which(H3K27ac_plot$z <= 3),]

options(repr.plot.width = 3, repr.plot.height = 4, repr.plot.res = 1000, repr.plot.pointsize = 20)

ggplot(H3K27ac_plot, aes(x=factor(comp), y=z, fill=factor(comp)))+
geom_boxplot(fill=c("dodgerblue3","yellow3"),
             size=1,
             color='black',
             width=0.5,
             alpha=0.7)+
geom_violin(adjust=1,
            scale='count',
            fill=c("black"),
            alpha=0.5)+
    guides(fill="none")+
theme_minimal()+
    labs(title=HISTONE_para,
        subtitle = "***\n")+
    xlab("")+
    ylab("enrichment of read count\n(z-score normalization)\n")+
theme(plot.title=element_text(face="bold.italic",hjust=0.5,size=25,color = "black"),
      axis.text.x=element_text(face="bold",size=12,color = "black"),
      axis.text.y=element_text(face="bold",size=10),
      plot.subtitle=element_text(face="bold",hjust=0.5,vjust=-3,size=15,color = "black"),
      axis.title.x = element_text(size=0,color = "black"),
      axis.title.y = element_text(face="bold",size=15,color = "black"))
    }

single_Vplot = function(input){
v <- read_vplot(input)

v$label = c(nrow(v):1)
v = v[order(v$label),-ncol(v)]

data = v

options(repr.plot.width = 2.5, repr.plot.height = 3, repr.plot.res = 1500, repr.plot.pointsize = 10)
conte=data
namen = c(1:nrow(data))
rownames(conte)=as.vector(namen)
PCA=conte
vvv=pheatmap(PCA,
         cluster_rows = F,
         cluster_cols = F,
         show_rownames = F,
         show_colnames = F,
         legend = FALSE,
         border_color = FALSE,
         color = colorRampPalette(c("white","dodgerblue1","dodgerblue2","dodgerblue3","dodgerblue4","darkblue"))(100),
         breaks=seq(from=0,to=2,length.out = 100))
vvv
    }


diff_Vplot = function(DAY0,DAY4){
DAY0 <- read_vplot(DAY0)
DAY4 <- read_vplot(DAY4)

v = DAY4 - DAY0 
v$label = c(nrow(v):1)
v = v[order(v$label),-ncol(v)]

data = v

options(repr.plot.width = 2.5, repr.plot.height = 3, repr.plot.res = 1500, repr.plot.pointsize = 10)
conte=data
namen = c(1:nrow(data))
rownames(conte)=as.vector(namen)
PCA=conte
vvv=pheatmap(PCA,
         cluster_rows = F,
         cluster_cols = F,
         show_rownames = F,
         show_colnames = F,
         legend = FALSE,
         border_color = FALSE,
         color = colorRampPalette(c("darkblue","blue","skyblue","pink","red","darkred"))(100),
         breaks=seq(from=-1,to=1,length.out = 100))
vvv
    }

TOP_BOTTOM_MAT = function(chr_para,ratio,outpath,histone_input){

chr_DF = c("chr1","chr2","chr3","chr4","chr5","chr6","chr7","chr8","chr9","chr10",
           "chr11","chr12","chr13","chr14","chr15","chr16","chr17","chr18","chr19")

len_DF = c(195471971,182113224,160039680,156508116,151834684,149736546,145441459,129401213,124595110,130694993,
           122082543,120129022,120421639,124902244,104043685,98207768,94987271,90702639,61431566)

para_DF = cbind(chr_DF,len_DF)
para_DF = as.data.frame(para_DF)

len = as.numeric(para_DF[which(para_DF[,1] == chr_para),2])

DAY0_WT_H3K27ac_region = read.table(histone_input,sep = "\t", header=F,fill = TRUE)
DAY0_WT_H3K27ac_region = DAY0_WT_H3K27ac_region[,1:3]
colnames(DAY0_WT_H3K27ac_region) = c("chr","starts","ends")

DAY0_WT_H3K27ac_region = DAY0_WT_H3K27ac_region[which(DAY0_WT_H3K27ac_region[,1] == chr_para),]

AtoA_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/AtoA_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
BtoA_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/BtoA_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
AtoB_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/AtoB_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
BtoB_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/BtoB_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)

partially_AtoA_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/partially_AtoA_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
partially_BtoA_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/partially_BtoA_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
partially_AtoB_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/partially_AtoB_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
partially_BtoB_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/partially_BtoB_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)

Acompartment = rbind(AtoA_100kb_bin[,1:3],
                     AtoB_100kb_bin[,1:3],
                     partially_AtoA_100kb_bin[,1:3],
                     partially_AtoB_100kb_bin[,1:3])
Bcompartment = rbind(BtoA_100kb_bin[,1:3],
                     BtoB_100kb_bin[,1:3],
                     partially_BtoA_100kb_bin[,1:3],
                     partially_BtoB_100kb_bin[,1:3])

Acompartment = Acompartment[which(Acompartment[,1] == chr_para),]
Bcompartment = Bcompartment[which(Bcompartment[,1] == chr_para),]

Acompartment[,2] = Acompartment[,2] -1
Bcompartment[,2] = Bcompartment[,2] -1

Acompartment = Acompartment[order(Acompartment[,2]),]
Bcompartment = Bcompartment[order(Bcompartment[,2]),]

Acompartment[,2] = as.integer(Acompartment[,2])
Acompartment[,3] = as.integer(Acompartment[,3])

Bcompartment[,2] = as.integer(Bcompartment[,2])
Bcompartment[,3] = as.integer(Bcompartment[,3])

Acompartment$Asite_label = paste(Acompartment[,1],
                                      Acompartment[,2],
                                      Acompartment[,3])

Acompartment$Bsite_label = paste(Acompartment[,1],
                                      Acompartment[,2],
                                      Acompartment[,3])

Bcompartment$Asite_label = paste(Bcompartment[,1],
                                      Bcompartment[,2],
                                      Bcompartment[,3])

Bcompartment$Bsite_label = paste(Bcompartment[,1],
                                      Bcompartment[,2],
                                      Bcompartment[,3])

start = seq(0,len,by=100000)
start = as.integer(start)
start = as.data.frame(start)
chrname = c(1:nrow(start))

DAY0_WT_corrected_binlist_FULL = cbind(chrname,start)
DAY0_WT_corrected_binlist_FULL[,1] = chr_para
colnames(DAY0_WT_corrected_binlist_FULL) = c("chrA","startA")
DAY0_WT_corrected_binlist_FULL$endA = DAY0_WT_corrected_binlist_FULL[,2]+100000
DAY0_WT_corrected_binlist_FULL[,3] = as.integer(DAY0_WT_corrected_binlist_FULL[,3])

DAY0_WT_corrected_binlist_DB = DAY0_WT_corrected_binlist_FULL
DAY0_WT_corrected_binlist_DB$chrB = 0
DAY0_WT_corrected_binlist_DB$startB = 0
DAY0_WT_corrected_binlist_DB$endB = 0

DAY0_WT_corrected_binlist_DB$chrB = DAY0_WT_corrected_binlist_FULL[1,1]
DAY0_WT_corrected_binlist_DB$startB = DAY0_WT_corrected_binlist_FULL[1,2]
DAY0_WT_corrected_binlist_DB$endB = DAY0_WT_corrected_binlist_FULL[1,3]

for(i in 2:nrow(DAY0_WT_corrected_binlist_FULL)){
    DAY0_WT_corrected_binlist_DB_frag = DAY0_WT_corrected_binlist_FULL
    DAY0_WT_corrected_binlist_DB_frag$chrB = 0
    DAY0_WT_corrected_binlist_DB_frag$startB = 0
    DAY0_WT_corrected_binlist_DB_frag$endB = 0
    
    DAY0_WT_corrected_binlist_DB_frag$chrB = DAY0_WT_corrected_binlist_FULL[i,1]
    DAY0_WT_corrected_binlist_DB_frag$startB = DAY0_WT_corrected_binlist_FULL[i,2]
    DAY0_WT_corrected_binlist_DB_frag$endB = DAY0_WT_corrected_binlist_FULL[i,3]

    DAY0_WT_corrected_binlist_DB = rbind(DAY0_WT_corrected_binlist_DB,DAY0_WT_corrected_binlist_DB_frag)
    }

DAY0_WT_corrected_binAC = DAY0_WT_corrected_binlist_FULL
DAY0_WT_corrected_binAC[,4] = 0

for(i in 1:nrow(DAY0_WT_corrected_binAC)){
    DAY0_WT_H3K27ac_region_tmp1 = DAY0_WT_H3K27ac_region[which(DAY0_WT_H3K27ac_region[,2] < DAY0_WT_corrected_binAC[i,3]),]
    DAY0_WT_H3K27ac_region_type1 = DAY0_WT_H3K27ac_region_tmp1[which(DAY0_WT_H3K27ac_region_tmp1[,3] > DAY0_WT_corrected_binAC[i,2]),]
    
    overlaped = DAY0_WT_H3K27ac_region_type1
    overlaped = unique(overlaped)

    if(nrow(overlaped)>0){
        DAY0_WT_corrected_binAC[i,4] = nrow(overlaped)
        }
}

colnames(DAY0_WT_corrected_binAC)[4] = "acA"

DAY0_WT_corrected_binAC_AB = DAY0_WT_corrected_binAC
DAY0_WT_corrected_binAC_AB$chrB = 0
DAY0_WT_corrected_binAC_AB$startB = 0
DAY0_WT_corrected_binAC_AB$endB = 0
DAY0_WT_corrected_binAC_AB$acB = 0 

DAY0_WT_corrected_binAC_AB$chrB = DAY0_WT_corrected_binAC[1,1]
DAY0_WT_corrected_binAC_AB$startB = DAY0_WT_corrected_binAC[1,2]
DAY0_WT_corrected_binAC_AB$endB = DAY0_WT_corrected_binAC[1,3]
DAY0_WT_corrected_binAC_AB$acB = DAY0_WT_corrected_binAC[1,4]

for(i in 2:nrow(DAY0_WT_corrected_binAC)){
    DAY0_WT_corrected_binAC_AB_frag = DAY0_WT_corrected_binAC
    DAY0_WT_corrected_binAC_AB_frag$chrB = 0
    DAY0_WT_corrected_binAC_AB_frag$startB = 0
    DAY0_WT_corrected_binAC_AB_frag$endB = 0
    DAY0_WT_corrected_binAC_AB_frag$acB = 0 
    
    DAY0_WT_corrected_binAC_AB_frag$chrB = DAY0_WT_corrected_binAC[i,1]
    DAY0_WT_corrected_binAC_AB_frag$startB = DAY0_WT_corrected_binAC[i,2]
    DAY0_WT_corrected_binAC_AB_frag$endB = DAY0_WT_corrected_binAC[i,3]
    DAY0_WT_corrected_binAC_AB_frag$acB = DAY0_WT_corrected_binAC[i,4]

    DAY0_WT_corrected_binAC_AB = rbind(DAY0_WT_corrected_binAC_AB,DAY0_WT_corrected_binAC_AB_frag)
    }

DAY0_bin_pAC = DAY0_WT_corrected_binAC_AB
DAY0_bin_pAC$pAC = DAY0_bin_pAC[,4]+DAY0_bin_pAC[,8]

DAY0_bin_pAC_topONE_tmp = DAY0_bin_pAC[order(-DAY0_bin_pAC[,9]),]
DAY0_bin_pAC_topONE_tmp[round(nrow(DAY0_bin_pAC)*ratio):nrow(DAY0_bin_pAC_topONE_tmp),9] = 0

DAY0_bin_pAC_topONE_tmp$label = paste(DAY0_bin_pAC_topONE_tmp[,1],
                                            DAY0_bin_pAC_topONE_tmp[,2],
                                            DAY0_bin_pAC_topONE_tmp[,3],
                                            DAY0_bin_pAC_topONE_tmp[,5],
                                            DAY0_bin_pAC_topONE_tmp[,6],
                                            DAY0_bin_pAC_topONE_tmp[,7])

DAY0_bin_pAC_topONE = DAY0_bin_pAC
DAY0_bin_pAC_topONE$label = paste(DAY0_bin_pAC_topONE[,1],
                                        DAY0_bin_pAC_topONE[,2],
                                        DAY0_bin_pAC_topONE[,3],
                                        DAY0_bin_pAC_topONE[,5],
                                        DAY0_bin_pAC_topONE[,6],
                                        DAY0_bin_pAC_topONE[,7])

DAY0_bin_pAC_topONE_tmp = DAY0_bin_pAC_topONE_tmp[,c(10,9)]
DAY0_bin_pAC_topONE = DAY0_bin_pAC_topONE[,-9]
DAY0_bin_pAC_topONE = inner_join(DAY0_bin_pAC_topONE,DAY0_bin_pAC_topONE_tmp,by="label")

DAY0_bin_pAC_bottomONE_tmp = DAY0_bin_pAC[order(-DAY0_bin_pAC[,9]),]
DAY0_bin_pAC_bottomONE_tmp[1:(nrow(DAY0_bin_pAC_bottomONE_tmp)-round(nrow(DAY0_bin_pAC)*ratio)),9] = 0

DAY0_bin_pAC_bottomONE_tmp$label = paste(DAY0_bin_pAC_bottomONE_tmp[,1],
                                            DAY0_bin_pAC_bottomONE_tmp[,2],
                                            DAY0_bin_pAC_bottomONE_tmp[,3],
                                            DAY0_bin_pAC_bottomONE_tmp[,5],
                                            DAY0_bin_pAC_bottomONE_tmp[,6],
                                            DAY0_bin_pAC_bottomONE_tmp[,7])

DAY0_bin_pAC_bottomONE = DAY0_bin_pAC
DAY0_bin_pAC_bottomONE$label = paste(DAY0_bin_pAC_bottomONE[,1],
                                        DAY0_bin_pAC_bottomONE[,2],
                                        DAY0_bin_pAC_bottomONE[,3],
                                        DAY0_bin_pAC_bottomONE[,5],
                                        DAY0_bin_pAC_bottomONE[,6],
                                        DAY0_bin_pAC_bottomONE[,7])

DAY0_bin_pAC_bottomONE_tmp = DAY0_bin_pAC_bottomONE_tmp[,c(10,9)]
DAY0_bin_pAC_bottomONE = DAY0_bin_pAC_bottomONE[,-9]
DAY0_bin_pAC_bottomONE = inner_join(DAY0_bin_pAC_bottomONE,DAY0_bin_pAC_bottomONE_tmp,by="label")
DAY0_bin_pAC_bottomONE[,10] = DAY0_bin_pAC_bottomONE[,10]*(-1)

rowBIN = seq(1:nrow(start))
BINnum = rep(rowBIN, times = nrow(start))
BINnum = as.data.frame(BINnum)

DAY0_bin_pAC_BINnum = cbind(DAY0_bin_pAC,BINnum)

DAY0_bin_pAC_MAT = DAY0_bin_pAC_BINnum[which(DAY0_bin_pAC_BINnum$BINnum == 1),8]
DAY0_bin_pAC_MAT = as.data.frame(DAY0_bin_pAC_MAT)
colnames(DAY0_bin_pAC_MAT) = NULL

for(i in 2:nrow(start)){
    DAY0_bin_pAC_MATfrag = DAY0_bin_pAC_BINnum[which(DAY0_bin_pAC_BINnum$BINnum == i),8]
    DAY0_bin_pAC_MATfrag = as.data.frame(DAY0_bin_pAC_MATfrag)
    colnames(DAY0_bin_pAC_MATfrag) = NULL
    DAY0_bin_pAC_MAT = cbind(DAY0_bin_pAC_MAT,DAY0_bin_pAC_MATfrag)
    }

colnames(DAY0_bin_pAC_MAT) = NULL

DAY0_bin_pAC_MAT = as.matrix(DAY0_bin_pAC_MAT)
pmean <- function(x,y) (x+y)/2
DAY0_bin_pAC_MAT[] <- pmean(DAY0_bin_pAC_MAT, matrix(DAY0_bin_pAC_MAT, nrow(DAY0_bin_pAC_MAT), byrow=TRUE))     

rowBIN = seq(1:nrow(start))
BINnum = rep(rowBIN, times = nrow(start))
BINnum = as.data.frame(BINnum)

DAY0_bin_pAC_topONE_BINnum = cbind(DAY0_bin_pAC_topONE,BINnum)

DAY0_bin_pAC_topONE_MAT = DAY0_bin_pAC_topONE_BINnum[which(DAY0_bin_pAC_topONE_BINnum$BINnum == 1),10]
DAY0_bin_pAC_topONE_MAT = as.data.frame(DAY0_bin_pAC_topONE_MAT)
colnames(DAY0_bin_pAC_topONE_MAT) = NULL

for(i in 2:nrow(start)){
    DAY0_bin_pAC_topONE_MATfrag = DAY0_bin_pAC_topONE_BINnum[which(DAY0_bin_pAC_topONE_BINnum$BINnum == i),10]
    DAY0_bin_pAC_topONE_MATfrag = as.data.frame(DAY0_bin_pAC_topONE_MATfrag)
    colnames(DAY0_bin_pAC_topONE_MATfrag) = NULL
    DAY0_bin_pAC_topONE_MAT = cbind(DAY0_bin_pAC_topONE_MAT,DAY0_bin_pAC_topONE_MATfrag)
    }

colnames(DAY0_bin_pAC_topONE_MAT) = NULL

DAY0_bin_pAC_topONE_MAT = as.matrix(DAY0_bin_pAC_topONE_MAT)
pmean <- function(x,y) (x+y)/2
DAY0_bin_pAC_topONE_MAT[] <- pmean(DAY0_bin_pAC_topONE_MAT, matrix(DAY0_bin_pAC_topONE_MAT, nrow(DAY0_bin_pAC_topONE_MAT), byrow=TRUE))

rowBIN = seq(1:nrow(start))
BINnum = rep(rowBIN, times = nrow(start))
BINnum = as.data.frame(BINnum)

DAY0_bin_pAC_bottomONE_BINnum = cbind(DAY0_bin_pAC_bottomONE,BINnum)

DAY0_bin_pAC_bottomONE_MAT = DAY0_bin_pAC_bottomONE_BINnum[which(DAY0_bin_pAC_bottomONE_BINnum$BINnum == 1),10]
DAY0_bin_pAC_bottomONE_MAT = as.data.frame(DAY0_bin_pAC_bottomONE_MAT)
colnames(DAY0_bin_pAC_bottomONE_MAT) = NULL

for(i in 2:nrow(start)){
    DAY0_bin_pAC_bottomONE_MATfrag = DAY0_bin_pAC_bottomONE_BINnum[which(DAY0_bin_pAC_bottomONE_BINnum$BINnum == i),10]
    DAY0_bin_pAC_bottomONE_MATfrag = as.data.frame(DAY0_bin_pAC_bottomONE_MATfrag)
    colnames(DAY0_bin_pAC_bottomONE_MATfrag) = NULL
    DAY0_bin_pAC_bottomONE_MAT = cbind(DAY0_bin_pAC_bottomONE_MAT,DAY0_bin_pAC_bottomONE_MATfrag)
    }

colnames(DAY0_bin_pAC_bottomONE_MAT) = NULL

DAY0_bin_pAC_bottomONE_MAT = as.matrix(DAY0_bin_pAC_bottomONE_MAT)
pmean <- function(x,y) (x+y)/2
DAY0_bin_pAC_bottomONE_MAT[] <- pmean(DAY0_bin_pAC_bottomONE_MAT, matrix(DAY0_bin_pAC_bottomONE_MAT, nrow(DAY0_bin_pAC_bottomONE_MAT), byrow=TRUE)) 


DAY0_bin_pAC_TEN_MAT = DAY0_bin_pAC_topONE_MAT+DAY0_bin_pAC_bottomONE_MAT

DAY0_bin_pAC_bottomTEN = DAY0_bin_pAC_bottomONE_BINnum[which(DAY0_bin_pAC_bottomONE_BINnum$pAC < 0),c(1,2,3,5,6,7)]

DAY0_bin_pAC_bottomTEN$Asite_label = paste(DAY0_bin_pAC_bottomTEN[,1],
                                                DAY0_bin_pAC_bottomTEN[,2],
                                                DAY0_bin_pAC_bottomTEN[,3])

DAY0_bin_pAC_bottomTEN$Bsite_label = paste(DAY0_bin_pAC_bottomTEN[,4],
                                                DAY0_bin_pAC_bottomTEN[,5],
                                                DAY0_bin_pAC_bottomTEN[,6])

DAY0_bin_pAC_bottomTEN_Asite = DAY0_bin_pAC_bottomTEN

DAY0_bin_pAC_bottomTEN_Asite$Asite = "NONE"

for(i in 1:nrow(DAY0_bin_pAC_bottomTEN_Asite)){
    Acompartment_corr = Acompartment[which(Acompartment[,4] == DAY0_bin_pAC_bottomTEN_Asite[i,7]),]
    Bcompartment_corr = Bcompartment[which(Bcompartment[,4] == DAY0_bin_pAC_bottomTEN_Asite[i,7]),]

    if(nrow(Acompartment_corr)>0){
        DAY0_bin_pAC_bottomTEN_Asite[i,9] = "A"
        }

    if(nrow(Bcompartment_corr)>0){
        DAY0_bin_pAC_bottomTEN_Asite[i,9] = "B"
        }
}

DAY0_bin_pAC_bottomTEN_FULL = DAY0_bin_pAC_bottomTEN_Asite

DAY0_bin_pAC_bottomTEN_FULL$Bsite = "NONE"

for(i in 1:nrow(DAY0_bin_pAC_bottomTEN_FULL)){
    Acompartment_corr = Acompartment[which(Acompartment[,4] == DAY0_bin_pAC_bottomTEN_FULL[i,8]),]
    Bcompartment_corr = Bcompartment[which(Bcompartment[,4] == DAY0_bin_pAC_bottomTEN_FULL[i,8]),]

    if(nrow(Acompartment_corr)>0){
        DAY0_bin_pAC_bottomTEN_FULL[i,10] = "A"
        }

    if(nrow(Bcompartment_corr)>0){
        DAY0_bin_pAC_bottomTEN_FULL[i,10] = "B"
        }
}

DAY0_bin_pAC_topTEN = DAY0_bin_pAC_topONE_BINnum[which(DAY0_bin_pAC_topONE_BINnum$pAC > 0),c(1,2,3,5,6,7)]

DAY0_bin_pAC_topTEN$Asite_label = paste(DAY0_bin_pAC_topTEN[,1],
                                                DAY0_bin_pAC_topTEN[,2],
                                                DAY0_bin_pAC_topTEN[,3])

DAY0_bin_pAC_topTEN$Bsite_label = paste(DAY0_bin_pAC_topTEN[,4],
                                                DAY0_bin_pAC_topTEN[,5],
                                                DAY0_bin_pAC_topTEN[,6])

DAY0_bin_pAC_topTEN_Asite = DAY0_bin_pAC_topTEN

DAY0_bin_pAC_topTEN_Asite$Asite = "NONE"

for(i in 1:nrow(DAY0_bin_pAC_topTEN_Asite)){
    Acompartment_corr = Acompartment[which(Acompartment[,4] == DAY0_bin_pAC_topTEN_Asite[i,7]),]
    Bcompartment_corr = Bcompartment[which(Bcompartment[,4] == DAY0_bin_pAC_topTEN_Asite[i,7]),]

    if(nrow(Acompartment_corr)>0){
        DAY0_bin_pAC_topTEN_Asite[i,9] = "A"
        }

    if(nrow(Bcompartment_corr)>0){
        DAY0_bin_pAC_topTEN_Asite[i,9] = "B"
        }
}

DAY0_bin_pAC_topTEN_FULL = DAY0_bin_pAC_topTEN_Asite

DAY0_bin_pAC_topTEN_FULL$Bsite = "NONE"

for(i in 1:nrow(DAY0_bin_pAC_topTEN_FULL)){
    Acompartment_corr = Acompartment[which(Acompartment[,4] == DAY0_bin_pAC_topTEN_FULL[i,8]),]
    Bcompartment_corr = Bcompartment[which(Bcompartment[,4] == DAY0_bin_pAC_topTEN_FULL[i,8]),]

    if(nrow(Acompartment_corr)>0){
        DAY0_bin_pAC_topTEN_FULL[i,10] = "A"
        }

    if(nrow(Bcompartment_corr)>0){
        DAY0_bin_pAC_topTEN_FULL[i,10] = "B"
        }
}
    
write.table(DAY0_bin_pAC_TEN_MAT, file = outpath ,col.names=T, row.names=F, quote=F, sep="\t")
    }

GENE_PLOT = function(x){

DESeq2_normalized = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG1/3_DEGs/count/DESeq2_normalized.txt",sep = "\t", header=T)
DESeq2_normalized = DESeq2_normalized[,1:9]

DEG_FULL_P = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG1/3_DEGs/DEG_result/condition_DAY4_vs_DAY0_0.05_1.5/DEG_FULL.txt",sep = "\t", header=T,fill = TRUE)


geneDATA = DESeq2_normalized[which(DESeq2_normalized[,1] == x),]

DAY0_mean = (geneDATA[,2] + geneDATA[,3] + geneDATA[,4] + geneDATA[,5])/4
DAY0_sd = sd(c(geneDATA[,2],geneDATA[,3],geneDATA[,4], geneDATA[,5]))
DAY0_se = DAY0_sd/(4^(1/2))

DAY4_mean = (geneDATA[,6] + geneDATA[,7] + geneDATA[,8] + geneDATA[,9])/4
DAY4_sd = sd(c(geneDATA[,6],geneDATA[,7],geneDATA[,8], geneDATA[,9]))
DAY4_se = DAY4_sd/(4^(1/2))

DAY0_DATA = "DAY0"
DAY0_DATA = as.data.frame(DAY0_DATA)
DAY0_DATA$DAY = "DAY0"
DAY0_DATA$mean = DAY0_mean
DAY0_DATA$sd = DAY0_sd
DAY0_DATA$se = DAY0_se
colnames(DAY0_DATA) = c("Condition","DAY","mean","sd","se")  

DAY4_DATA = "DAY4"
DAY4_DATA = as.data.frame(DAY4_DATA)
DAY4_DATA$DAY = "DAY4"
DAY4_DATA$mean = DAY4_mean
DAY4_DATA$sd = DAY4_sd
DAY4_DATA$se = DAY4_se
colnames(DAY4_DATA) = c("Condition","DAY","mean","sd","se")      

RNAplotDATA = rbind(DAY0_DATA,
                    DAY4_DATA)

options(repr.plot.width = 2, repr.plot.height = 3, repr.plot.res = 300, repr.plot.pointsize = 40)

p<-ggplot(RNAplotDATA,aes(x=DAY,y=mean,fill=Condition))
p2=p+geom_bar(stat="identity",position="dodge",width=0.6)+
scale_fill_manual(values = c("gray25","gray25"))+
theme_classic(base_size = 10) +
ggtitle(label = geneDATA[,1],
        subtitle = paste("p = ",round(DEG_FULL_P[which(DEG_FULL_P[,1] == x),6],4),
                         sep="")) +
ylab("Normalized gene expression\n") +
#scale_x_discrete(limits=c("control", "exosome", "exosome\nRunx3")) + 
theme(plot.title=element_text(face="bold.italic",hjust=0.5,size=20,color = "black"),
      plot.subtitle=element_text(face="bold.italic",hjust=0.5,size=10,color = "black"),
      axis.text.x=element_text(face="bold",size=13,color = "black"),
      axis.text.y=element_text(face="bold",size=12,color = "black"),
      axis.title.x = element_text(face="bold",size = 0,color = "darkblue"),
      axis.title.y = element_text(face="bold",size = 12,color = "black"),
      legend.title=element_text(face="bold",size=10), 
      legend.text=element_text(size=8),
      legend.position = "none") +
geom_errorbar(aes(ymin=mean-se,ymax=mean+se),
              position=position_dodge(0.9),width=0.2)
p2
    }

ValcanoMAKER = function(DEG_FULL_path,xrange,yrange,FC,adjP,labeltext){

    n_xrange = xrange*(-1)
    adjP = as.numeric(adjP)
    
    DEG_FULL = read.table(DEG_FULL_path,sep = "\t", header=T,fill = TRUE)
    DEG_FULL$label = 0
    
    DEG_Fcut = DEG_FULL[which(DEG_FULL[,7] <= adjP),]
    
    DEG_UP = DEG_Fcut[which(DEG_Fcut[,3] >= log2(FC)),]
    DEG_DOWN = DEG_Fcut[which(DEG_Fcut[,3] <= -log2(FC)),]

    for(i in 1:nrow(DEG_UP)){
        DEG_FULL[which(DEG_FULL[,1] == DEG_UP[i,1]),8] = 2
        }

    for(i in 1:nrow(DEG_DOWN)){
        DEG_FULL[which(DEG_FULL[,1] == DEG_DOWN[i,1]),8] = 1
        }
    
    options(repr.plot.width = 5, repr.plot.height = 5, repr.plot.res = 300, repr.plot.pointsize = 10)
    
    ggplot(DEG_FULL,aes(x=log2FoldChange, y= -log2(padj), group=label)) +
    geom_point(data=subset(DEG_FULL,label==2),color="red", size=1) +
    geom_point(data=subset(DEG_FULL,label==1),color="blue", size=1) +
    geom_point(data=subset(DEG_FULL,label==0),color="grey25", size=1) +
    xlab("Log2FoldChange") + 
    ylab("P.adjust\n(Log2 Transformation)\n")+
    ggtitle(label = labeltext,
            subtitle = paste("FC :",FC,"/","P.adjust :",adjP,"\n\n")) +
    xlim(n_xrange,xrange)+
    theme_classic(base_size = 10) +
    scale_y_continuous(expand = c(0, 0),limit = c(0,yrange) ) +
    geom_vline(xintercept=c(-log2(FC),log2(FC)), color = "black",size = 1)+
    geom_hline(yintercept=-log2(adjP), color = "black",size = 1)+
    theme(plot.title=element_text(face="bold",hjust=0.5,size=25,color = "black"),
          plot.subtitle=element_text(face="bold",hjust=0.5,size=15,color = "black"),
          axis.text.x=element_text(face="bold",size=12),
          axis.text.y=element_text(face="bold",size=12),
          axis.title.x = element_text(face="bold",size = 15),
          axis.title.y = element_text(face="bold",size = 15),
          legend.title=element_text(face="bold",size=15), 
          legend.text=element_text(face="bold",size=15))
    }

mm10_RE = function(cate){
path = paste("/data3/psg/NGS_2025/4DN/HiC_data/FiG1/DBs/mm10_RE_",cate,"_tmp.txt",sep="")
mm10_RE_Histone_tmp = read.table(path, sep = "\t", header=F,fill = TRUE)
mm10_RE_Histone_tmp$ID = paste(mm10_RE_Histone_tmp[,1],mm10_RE_Histone_tmp[,2],mm10_RE_Histone_tmp[,3],mm10_RE_Histone_tmp[,4])

mm10_RE_Histone = mm10_RE_Histone_tmp[,c(1:4,15)]
mm10_RE_Histone = unique(mm10_RE_Histone)

mm10_RE_Histone$histone = 0
colnames(mm10_RE_Histone)[6] = cate

for(i in 1:nrow(mm10_RE_Histone)){mm10_RE_Histone[i,6] = nrow(mm10_RE_Histone_tmp[which(mm10_RE_Histone_tmp$ID == mm10_RE_Histone[i,5]),])}

path = paste("/data3/psg/NGS_2025/4DN/HiC_data/FiG1/DBs/mm10_RE_",cate,".txt",sep="")
    
write.table(mm10_RE_Histone, file = path,col.names=T, row.names=F, quote=F, sep="\t")
    }

Histone_filter_frag = function(dir,min,max){

main_path = paste("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/",dir,"/",sep="")

DAY0_WT_H3K27ac = read.table( paste(main_path,"DAY0_WT_H3K27ac.bedpe",sep=""),sep = "\t", header=F,fill = TRUE)
DAY0_WT_H3K27ac$para = DAY0_WT_H3K27ac[,6] - DAY0_WT_H3K27ac[,2]
DAY0_WT_H3K27ac = DAY0_WT_H3K27ac[which(DAY0_WT_H3K27ac[,11] >= min),]
DAY0_WT_H3K27ac = DAY0_WT_H3K27ac[which(DAY0_WT_H3K27ac[,11] <= max),]
DAY0_WT_H3K27ac = DAY0_WT_H3K27ac[,c(1,2,6,7)]
DAY0_WT_H3K27ac = DAY0_WT_H3K27ac[grepl("chr",DAY0_WT_H3K27ac[,1]), ]

write.table(DAY0_WT_H3K27ac,
            file = paste(main_path,"DAY0_WT_H3K27ac_150_fit.bedpe",sep=""),
            col.names=F, row.names=F, quote=F, sep="\t")

DAY0_WT_H3K4me1 = read.table( paste(main_path,"DAY0_WT_H3K4me1.bedpe",sep=""),sep = "\t", header=F,fill = TRUE)
DAY0_WT_H3K4me1$para = DAY0_WT_H3K4me1[,6] - DAY0_WT_H3K4me1[,2]
DAY0_WT_H3K4me1 = DAY0_WT_H3K4me1[which(DAY0_WT_H3K4me1[,11] >= min),]
DAY0_WT_H3K4me1 = DAY0_WT_H3K4me1[which(DAY0_WT_H3K4me1[,11] <= max),]
DAY0_WT_H3K4me1 = DAY0_WT_H3K4me1[,c(1,2,6,7)]
DAY0_WT_H3K4me1 = DAY0_WT_H3K4me1[grepl("chr",DAY0_WT_H3K4me1[,1]), ]

write.table(DAY0_WT_H3K4me1,
            file = paste(main_path,"DAY0_WT_H3K4me1_150_fit.bedpe",sep=""),
            col.names=F, row.names=F, quote=F, sep="\t")


DAY0_WT_H3K4me3 = read.table( paste(main_path,"DAY0_WT_H3K4me3.bedpe",sep=""),sep = "\t", header=F,fill = TRUE)
DAY0_WT_H3K4me3$para = DAY0_WT_H3K4me3[,6] - DAY0_WT_H3K4me3[,2]
DAY0_WT_H3K4me3 = DAY0_WT_H3K4me3[which(DAY0_WT_H3K4me3[,11] >= min),]
DAY0_WT_H3K4me3 = DAY0_WT_H3K4me3[which(DAY0_WT_H3K4me3[,11] <= max),]
DAY0_WT_H3K4me3 = DAY0_WT_H3K4me3[,c(1,2,6,7)]
DAY0_WT_H3K4me3 = DAY0_WT_H3K4me3[grepl("chr",DAY0_WT_H3K4me3[,1]), ]

write.table(DAY0_WT_H3K4me3,
            file = paste(main_path,"DAY0_WT_H3K4me3_150_fit.bedpe",sep=""),
            col.names=F, row.names=F, quote=F, sep="\t")


DAY4_WT_H3K27ac = read.table( paste(main_path,"DAY4_WT_H3K27ac.bedpe",sep=""),sep = "\t", header=F,fill = TRUE)
DAY4_WT_H3K27ac$para = DAY4_WT_H3K27ac[,6] - DAY4_WT_H3K27ac[,2]
DAY4_WT_H3K27ac = DAY4_WT_H3K27ac[which(DAY4_WT_H3K27ac[,11] >= min),]
DAY4_WT_H3K27ac = DAY4_WT_H3K27ac[which(DAY4_WT_H3K27ac[,11] <= max),]
DAY4_WT_H3K27ac = DAY4_WT_H3K27ac[,c(1,2,6,7)]
DAY4_WT_H3K27ac = DAY4_WT_H3K27ac[grepl("chr",DAY4_WT_H3K27ac[,1]), ]

write.table(DAY4_WT_H3K27ac,
            file = paste(main_path,"DAY4_WT_H3K27ac_150_fit.bedpe",sep=""),
            col.names=F, row.names=F, quote=F, sep="\t")

DAY4_WT_H3K4me1 = read.table( paste(main_path,"DAY4_WT_H3K4me1.bedpe",sep=""),sep = "\t", header=F,fill = TRUE)
DAY4_WT_H3K4me1$para = DAY4_WT_H3K4me1[,6] - DAY4_WT_H3K4me1[,2]
DAY4_WT_H3K4me1 = DAY4_WT_H3K4me1[which(DAY4_WT_H3K4me1[,11] >= min),]
DAY4_WT_H3K4me1 = DAY4_WT_H3K4me1[which(DAY4_WT_H3K4me1[,11] <= max),]
DAY4_WT_H3K4me1 = DAY4_WT_H3K4me1[,c(1,2,6,7)]
DAY4_WT_H3K4me1 = DAY4_WT_H3K4me1[grepl("chr",DAY4_WT_H3K4me1[,1]), ]

write.table(DAY4_WT_H3K4me1,
            file = paste(main_path,"DAY4_WT_H3K4me1_150_fit.bedpe",sep=""),
            col.names=F, row.names=F, quote=F, sep="\t")


DAY4_WT_H3K4me3 = read.table( paste(main_path,"DAY4_WT_H3K4me3.bedpe",sep=""),sep = "\t", header=F,fill = TRUE)
DAY4_WT_H3K4me3$para = DAY4_WT_H3K4me3[,6] - DAY4_WT_H3K4me3[,2]
DAY4_WT_H3K4me3 = DAY4_WT_H3K4me3[which(DAY4_WT_H3K4me3[,11] >= min),]
DAY4_WT_H3K4me3 = DAY4_WT_H3K4me3[which(DAY4_WT_H3K4me3[,11] <= max),]
DAY4_WT_H3K4me3 = DAY4_WT_H3K4me3[,c(1,2,6,7)]
DAY4_WT_H3K4me3 = DAY4_WT_H3K4me3[grepl("chr",DAY4_WT_H3K4me3[,1]), ]

write.table(DAY4_WT_H3K4me3,
            file = paste(main_path,"DAY4_WT_H3K4me3_150_fit.bedpe",sep=""),
            col.names=F, row.names=F, quote=F, sep="\t")
    }

ABcomp_to_Histone = function(ChIP,ABcomp){

DAY0_path = paste("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/",ChIP,"/DAY0_",ABcomp,"_",sep="")
DAY0_partially_path = paste("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/",ChIP,"/DAY0_partially_",ABcomp,"_",sep="")
DAY4_path = paste("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/",ChIP,"/DAY4_",ABcomp,"_",sep="")
DAY4_partially_path = paste("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/",ChIP,"/DAY4_partially_",ABcomp,"_",sep="")
AB_path = paste("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/AB_status_202505/",ABcomp,"_bulk_sorted.txt",sep="")

DAY0_ABcomp_H3K27ac_tmp = read.table( paste(DAY0_path,"H3K27ac.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
DAY0_partially_ABcomp_H3K27ac_tmp = read.table( paste(DAY0_partially_path,"H3K27ac.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
DAY4_ABcomp_H3K27ac_tmp = read.table( paste(DAY4_path,"H3K27ac.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
DAY4_partially_ABcomp_H3K27ac_tmp = read.table( paste(DAY4_partially_path,"H3K27ac.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
AB_comp = read.table(AB_path,sep = "\t", header=F,fill = TRUE)

DAY0_ABcomp_H3K27ac_tmp = rbind(DAY0_ABcomp_H3K27ac_tmp,DAY0_partially_ABcomp_H3K27ac_tmp)
DAY0_ABcomp_H3K27ac_tmp$ID = paste(DAY0_ABcomp_H3K27ac_tmp[,1],DAY0_ABcomp_H3K27ac_tmp[,2],DAY0_ABcomp_H3K27ac_tmp[,3])
DAY4_ABcomp_H3K27ac_tmp = rbind(DAY4_ABcomp_H3K27ac_tmp,DAY4_partially_ABcomp_H3K27ac_tmp)
DAY4_ABcomp_H3K27ac_tmp$ID = paste(DAY4_ABcomp_H3K27ac_tmp[,1],DAY4_ABcomp_H3K27ac_tmp[,2],DAY4_ABcomp_H3K27ac_tmp[,3])
AB_comp_H3K27ac = AB_comp
AB_comp_H3K27ac$ID = paste(AB_comp_H3K27ac[,1],AB_comp_H3K27ac[,2],AB_comp_H3K27ac[,3])

AB_comp_H3K27ac$DAY0_H3K27ac = 0
AB_comp_H3K27ac$DAY4_H3K27ac = 0
for(i in 1:nrow(AB_comp_H3K27ac)){AB_comp_H3K27ac[i,7] = nrow(DAY0_ABcomp_H3K27ac_tmp[which(DAY0_ABcomp_H3K27ac_tmp$ID == AB_comp_H3K27ac[i,6]),])}
for(i in 1:nrow(AB_comp_H3K27ac)){AB_comp_H3K27ac[i,8] = nrow(DAY4_ABcomp_H3K27ac_tmp[which(DAY4_ABcomp_H3K27ac_tmp$ID == AB_comp_H3K27ac[i,6]),])}

AB_comp_H3K27ac$len = (AB_comp_H3K27ac[,3]-AB_comp_H3K27ac[,2])/100000
AB_comp_H3K27ac$DAY0_H3K27ac_ratio = AB_comp_H3K27ac[,7]/AB_comp_H3K27ac[,9]
AB_comp_H3K27ac$DAY4_H3K27ac_ratio = AB_comp_H3K27ac[,8]/AB_comp_H3K27ac[,9]

AB_comp_H3K27ac = AB_comp_H3K27ac[,c(1,2,3,4,5,6,10,11)]

DAY0_ABcomp_H3K4me1_tmp = read.table( paste(DAY0_path,"H3K4me1.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
DAY0_partially_ABcomp_H3K4me1_tmp = read.table( paste(DAY0_partially_path,"H3K4me1.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
DAY4_ABcomp_H3K4me1_tmp = read.table( paste(DAY4_path,"H3K4me1.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
DAY4_partially_ABcomp_H3K4me1_tmp = read.table( paste(DAY4_partially_path,"H3K4me1.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
AB_comp = read.table(AB_path,sep = "\t", header=F,fill = TRUE)

DAY0_ABcomp_H3K4me1_tmp = rbind(DAY0_ABcomp_H3K4me1_tmp,DAY0_partially_ABcomp_H3K4me1_tmp)
DAY0_ABcomp_H3K4me1_tmp$ID = paste(DAY0_ABcomp_H3K4me1_tmp[,1],DAY0_ABcomp_H3K4me1_tmp[,2],DAY0_ABcomp_H3K4me1_tmp[,3])
DAY4_ABcomp_H3K4me1_tmp = rbind(DAY4_ABcomp_H3K4me1_tmp,DAY4_partially_ABcomp_H3K4me1_tmp)
DAY4_ABcomp_H3K4me1_tmp$ID = paste(DAY4_ABcomp_H3K4me1_tmp[,1],DAY4_ABcomp_H3K4me1_tmp[,2],DAY4_ABcomp_H3K4me1_tmp[,3])
AB_comp_H3K4me1 = AB_comp
AB_comp_H3K4me1$ID = paste(AB_comp_H3K4me1[,1],AB_comp_H3K4me1[,2],AB_comp_H3K4me1[,3])

AB_comp_H3K4me1$DAY0_H3K4me1 = 0
AB_comp_H3K4me1$DAY4_H3K4me1 = 0
for(i in 1:nrow(AB_comp_H3K4me1)){AB_comp_H3K4me1[i,7] = nrow(DAY0_ABcomp_H3K4me1_tmp[which(DAY0_ABcomp_H3K4me1_tmp$ID == AB_comp_H3K4me1[i,6]),])}
for(i in 1:nrow(AB_comp_H3K4me1)){AB_comp_H3K4me1[i,8] = nrow(DAY4_ABcomp_H3K4me1_tmp[which(DAY4_ABcomp_H3K4me1_tmp$ID == AB_comp_H3K4me1[i,6]),])}

AB_comp_H3K4me1$len = (AB_comp_H3K4me1[,3]-AB_comp_H3K4me1[,2])/100000
AB_comp_H3K4me1$DAY0_H3K4me1_ratio = AB_comp_H3K4me1[,7]/AB_comp_H3K4me1[,9]
AB_comp_H3K4me1$DAY4_H3K4me1_ratio = AB_comp_H3K4me1[,8]/AB_comp_H3K4me1[,9]

AB_comp_H3K4me1 = AB_comp_H3K4me1[,c(1,2,3,4,5,6,10,11)]

DAY0_ABcomp_H3K4me3_tmp = read.table( paste(DAY0_path,"H3K4me3.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
DAY0_partially_ABcomp_H3K4me3_tmp = read.table( paste(DAY0_partially_path,"H3K4me3.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
DAY4_ABcomp_H3K4me3_tmp = read.table( paste(DAY4_path,"H3K4me3.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
DAY4_partially_ABcomp_H3K4me3_tmp = read.table( paste(DAY4_partially_path,"H3K4me3.txt",sep="") ,sep = "\t", header=F,fill = TRUE)
AB_comp = read.table(AB_path,sep = "\t", header=F,fill = TRUE)

DAY0_ABcomp_H3K4me3_tmp = rbind(DAY0_ABcomp_H3K4me3_tmp,DAY0_partially_ABcomp_H3K4me3_tmp)
DAY0_ABcomp_H3K4me3_tmp$ID = paste(DAY0_ABcomp_H3K4me3_tmp[,1],DAY0_ABcomp_H3K4me3_tmp[,2],DAY0_ABcomp_H3K4me3_tmp[,3])
DAY4_ABcomp_H3K4me3_tmp = rbind(DAY4_ABcomp_H3K4me3_tmp,DAY4_partially_ABcomp_H3K4me3_tmp)
DAY4_ABcomp_H3K4me3_tmp$ID = paste(DAY4_ABcomp_H3K4me3_tmp[,1],DAY4_ABcomp_H3K4me3_tmp[,2],DAY4_ABcomp_H3K4me3_tmp[,3])
AB_comp_H3K4me3 = AB_comp
AB_comp_H3K4me3$ID = paste(AB_comp_H3K4me3[,1],AB_comp_H3K4me3[,2],AB_comp_H3K4me3[,3])

AB_comp_H3K4me3$DAY0_H3K4me3 = 0
AB_comp_H3K4me3$DAY4_H3K4me3 = 0
for(i in 1:nrow(AB_comp_H3K4me3)){AB_comp_H3K4me3[i,7] = nrow(DAY0_ABcomp_H3K4me3_tmp[which(DAY0_ABcomp_H3K4me3_tmp$ID == AB_comp_H3K4me3[i,6]),])}
for(i in 1:nrow(AB_comp_H3K4me3)){AB_comp_H3K4me3[i,8] = nrow(DAY4_ABcomp_H3K4me3_tmp[which(DAY4_ABcomp_H3K4me3_tmp$ID == AB_comp_H3K4me3[i,6]),])}

AB_comp_H3K4me3$len = (AB_comp_H3K4me3[,3]-AB_comp_H3K4me3[,2])/100000
AB_comp_H3K4me3$DAY0_H3K4me3_ratio = AB_comp_H3K4me3[,7]/AB_comp_H3K4me3[,9]
AB_comp_H3K4me3$DAY4_H3K4me3_ratio = AB_comp_H3K4me3[,8]/AB_comp_H3K4me3[,9]

AB_comp_H3K4me3 = AB_comp_H3K4me3[,c(1,2,3,4,5,6,10,11)]

ABcomp_HistoneLevel = inner_join(AB_comp_H3K27ac,AB_comp_H3K4me1[,6:8],by = "ID")
ABcomp_HistoneLevel = inner_join(ABcomp_HistoneLevel,AB_comp_H3K4me3[,6:8],by = "ID")

out_path = paste("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/",ChIP,"/",ABcomp,"_",sep="")
write.table(ABcomp_HistoneLevel,file = paste(out_path,"HistoneLevel.txt",sep=""),col.names=F, row.names=F, quote=F, sep="\t")
    }

his_MAP = function(chrom,len_chrom,histone_input,output_path){

chr_para = chrom
len = len_chrom

DAY0_WT_H3K27ac_region = read.table(histone_input,sep = "\t", header=F,fill = TRUE)
DAY0_WT_H3K27ac_region = DAY0_WT_H3K27ac_region[,1:3]
colnames(DAY0_WT_H3K27ac_region) = c("chr","starts","ends")

DAY0_WT_H3K27ac_region = DAY0_WT_H3K27ac_region[which(DAY0_WT_H3K27ac_region[,1] == chr_para),]

AtoA_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/AtoA_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
BtoA_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/BtoA_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
AtoB_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/AtoB_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
BtoB_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/BtoB_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)

partially_AtoA_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/partially_AtoA_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
partially_BtoA_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/partially_BtoA_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
partially_AtoB_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/partially_AtoB_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)
partially_BtoB_100kb_bin = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG23/ABcompartment/100kb_bin/partially_BtoB_100kb_bin.txt",sep = "\t", header=F,fill = TRUE)

Acompartment = rbind(AtoA_100kb_bin[,1:3],
                     AtoB_100kb_bin[,1:3],
                     partially_AtoA_100kb_bin[,1:3],
                     partially_AtoB_100kb_bin[,1:3])
Bcompartment = rbind(BtoA_100kb_bin[,1:3],
                     BtoB_100kb_bin[,1:3],
                     partially_BtoA_100kb_bin[,1:3],
                     partially_BtoB_100kb_bin[,1:3])

Acompartment = Acompartment[which(Acompartment[,1] == chr_para),]
Bcompartment = Bcompartment[which(Bcompartment[,1] == chr_para),]

Acompartment[,2] = Acompartment[,2] -1
Bcompartment[,2] = Bcompartment[,2] -1

Acompartment = Acompartment[order(Acompartment[,2]),]
Bcompartment = Bcompartment[order(Bcompartment[,2]),]

Acompartment[,2] = as.integer(Acompartment[,2])
Acompartment[,3] = as.integer(Acompartment[,3])

Bcompartment[,2] = as.integer(Bcompartment[,2])
Bcompartment[,3] = as.integer(Bcompartment[,3])

Acompartment$Asite_label = paste(Acompartment[,1],
                                      Acompartment[,2],
                                      Acompartment[,3])

Acompartment$Bsite_label = paste(Acompartment[,1],
                                      Acompartment[,2],
                                      Acompartment[,3])

Bcompartment$Asite_label = paste(Bcompartment[,1],
                                      Bcompartment[,2],
                                      Bcompartment[,3])

Bcompartment$Bsite_label = paste(Bcompartment[,1],
                                      Bcompartment[,2],
                                      Bcompartment[,3])

start = seq(0,len,by=100000)
start = as.integer(start)
start = as.data.frame(start)
chrname = c(1:nrow(start))

DAY0_WT_corrected_binlist_FULL = cbind(chrname,start)
DAY0_WT_corrected_binlist_FULL[,1] = chr_para
colnames(DAY0_WT_corrected_binlist_FULL) = c("chrA","startA")
DAY0_WT_corrected_binlist_FULL$endA = DAY0_WT_corrected_binlist_FULL[,2]+100000
DAY0_WT_corrected_binlist_FULL[,3] = as.integer(DAY0_WT_corrected_binlist_FULL[,3])

DAY0_WT_corrected_binlist_DB = DAY0_WT_corrected_binlist_FULL
DAY0_WT_corrected_binlist_DB$chrB = 0
DAY0_WT_corrected_binlist_DB$startB = 0
DAY0_WT_corrected_binlist_DB$endB = 0

DAY0_WT_corrected_binlist_DB$chrB = DAY0_WT_corrected_binlist_FULL[1,1]
DAY0_WT_corrected_binlist_DB$startB = DAY0_WT_corrected_binlist_FULL[1,2]
DAY0_WT_corrected_binlist_DB$endB = DAY0_WT_corrected_binlist_FULL[1,3]

for(i in 2:nrow(DAY0_WT_corrected_binlist_FULL)){
    DAY0_WT_corrected_binlist_DB_frag = DAY0_WT_corrected_binlist_FULL
    DAY0_WT_corrected_binlist_DB_frag$chrB = 0
    DAY0_WT_corrected_binlist_DB_frag$startB = 0
    DAY0_WT_corrected_binlist_DB_frag$endB = 0
    
    DAY0_WT_corrected_binlist_DB_frag$chrB = DAY0_WT_corrected_binlist_FULL[i,1]
    DAY0_WT_corrected_binlist_DB_frag$startB = DAY0_WT_corrected_binlist_FULL[i,2]
    DAY0_WT_corrected_binlist_DB_frag$endB = DAY0_WT_corrected_binlist_FULL[i,3]

    DAY0_WT_corrected_binlist_DB = rbind(DAY0_WT_corrected_binlist_DB,DAY0_WT_corrected_binlist_DB_frag)
    }

DAY0_WT_corrected_binAC = DAY0_WT_corrected_binlist_FULL
DAY0_WT_corrected_binAC[,4] = 0

for(i in 1:nrow(DAY0_WT_corrected_binAC)){
    DAY0_WT_H3K27ac_region_tmp1 = DAY0_WT_H3K27ac_region[which(DAY0_WT_H3K27ac_region[,2] < DAY0_WT_corrected_binAC[i,3]),]
    DAY0_WT_H3K27ac_region_type1 = DAY0_WT_H3K27ac_region_tmp1[which(DAY0_WT_H3K27ac_region_tmp1[,3] > DAY0_WT_corrected_binAC[i,2]),]
    
    overlaped = DAY0_WT_H3K27ac_region_type1
    overlaped = unique(overlaped)

    if(nrow(overlaped)>0){
        DAY0_WT_corrected_binAC[i,4] = nrow(overlaped)
        }
}

colnames(DAY0_WT_corrected_binAC)[4] = "acA"

DAY0_WT_corrected_binAC_AB = DAY0_WT_corrected_binAC
DAY0_WT_corrected_binAC_AB$chrB = 0
DAY0_WT_corrected_binAC_AB$startB = 0
DAY0_WT_corrected_binAC_AB$endB = 0
DAY0_WT_corrected_binAC_AB$acB = 0 

DAY0_WT_corrected_binAC_AB$chrB = DAY0_WT_corrected_binAC[1,1]
DAY0_WT_corrected_binAC_AB$startB = DAY0_WT_corrected_binAC[1,2]
DAY0_WT_corrected_binAC_AB$endB = DAY0_WT_corrected_binAC[1,3]
DAY0_WT_corrected_binAC_AB$acB = DAY0_WT_corrected_binAC[1,4]

for(i in 2:nrow(DAY0_WT_corrected_binAC)){
    DAY0_WT_corrected_binAC_AB_frag = DAY0_WT_corrected_binAC
    DAY0_WT_corrected_binAC_AB_frag$chrB = 0
    DAY0_WT_corrected_binAC_AB_frag$startB = 0
    DAY0_WT_corrected_binAC_AB_frag$endB = 0
    DAY0_WT_corrected_binAC_AB_frag$acB = 0 
    
    DAY0_WT_corrected_binAC_AB_frag$chrB = DAY0_WT_corrected_binAC[i,1]
    DAY0_WT_corrected_binAC_AB_frag$startB = DAY0_WT_corrected_binAC[i,2]
    DAY0_WT_corrected_binAC_AB_frag$endB = DAY0_WT_corrected_binAC[i,3]
    DAY0_WT_corrected_binAC_AB_frag$acB = DAY0_WT_corrected_binAC[i,4]

    DAY0_WT_corrected_binAC_AB = rbind(DAY0_WT_corrected_binAC_AB,DAY0_WT_corrected_binAC_AB_frag)
    }

DAY0_bin_pAC = DAY0_WT_corrected_binAC_AB
DAY0_bin_pAC$pAC = DAY0_bin_pAC[,4]+DAY0_bin_pAC[,8]

DAY0_bin_pAC_topONE_tmp = DAY0_bin_pAC[order(-DAY0_bin_pAC[,9]),]
DAY0_bin_pAC_topONE_tmp[round(nrow(DAY0_bin_pAC)*0.1):nrow(DAY0_bin_pAC_topONE_tmp),9] = 0

DAY0_bin_pAC_topONE_tmp$label = paste(DAY0_bin_pAC_topONE_tmp[,1],
                                            DAY0_bin_pAC_topONE_tmp[,2],
                                            DAY0_bin_pAC_topONE_tmp[,3],
                                            DAY0_bin_pAC_topONE_tmp[,5],
                                            DAY0_bin_pAC_topONE_tmp[,6],
                                            DAY0_bin_pAC_topONE_tmp[,7])

DAY0_bin_pAC_topONE = DAY0_bin_pAC
DAY0_bin_pAC_topONE$label = paste(DAY0_bin_pAC_topONE[,1],
                                        DAY0_bin_pAC_topONE[,2],
                                        DAY0_bin_pAC_topONE[,3],
                                        DAY0_bin_pAC_topONE[,5],
                                        DAY0_bin_pAC_topONE[,6],
                                        DAY0_bin_pAC_topONE[,7])

DAY0_bin_pAC_topONE_tmp = DAY0_bin_pAC_topONE_tmp[,c(10,9)]
DAY0_bin_pAC_topONE = DAY0_bin_pAC_topONE[,-9]
DAY0_bin_pAC_topONE = inner_join(DAY0_bin_pAC_topONE,DAY0_bin_pAC_topONE_tmp,by="label")

DAY0_bin_pAC_bottomONE_tmp = DAY0_bin_pAC[order(-DAY0_bin_pAC[,9]),]
DAY0_bin_pAC_bottomONE_tmp[1:(nrow(DAY0_bin_pAC_bottomONE_tmp)-round(nrow(DAY0_bin_pAC)*0.1)),9] = 0

DAY0_bin_pAC_bottomONE_tmp$label = paste(DAY0_bin_pAC_bottomONE_tmp[,1],
                                            DAY0_bin_pAC_bottomONE_tmp[,2],
                                            DAY0_bin_pAC_bottomONE_tmp[,3],
                                            DAY0_bin_pAC_bottomONE_tmp[,5],
                                            DAY0_bin_pAC_bottomONE_tmp[,6],
                                            DAY0_bin_pAC_bottomONE_tmp[,7])

DAY0_bin_pAC_bottomONE = DAY0_bin_pAC
DAY0_bin_pAC_bottomONE$label = paste(DAY0_bin_pAC_bottomONE[,1],
                                        DAY0_bin_pAC_bottomONE[,2],
                                        DAY0_bin_pAC_bottomONE[,3],
                                        DAY0_bin_pAC_bottomONE[,5],
                                        DAY0_bin_pAC_bottomONE[,6],
                                        DAY0_bin_pAC_bottomONE[,7])

DAY0_bin_pAC_bottomONE_tmp = DAY0_bin_pAC_bottomONE_tmp[,c(10,9)]
DAY0_bin_pAC_bottomONE = DAY0_bin_pAC_bottomONE[,-9]
DAY0_bin_pAC_bottomONE = inner_join(DAY0_bin_pAC_bottomONE,DAY0_bin_pAC_bottomONE_tmp,by="label")
DAY0_bin_pAC_bottomONE[,10] = DAY0_bin_pAC_bottomONE[,10]*(-1)

rowBIN = seq(1:nrow(start))
BINnum = rep(rowBIN, times = nrow(start))
BINnum = as.data.frame(BINnum)

DAY0_bin_pAC_BINnum = cbind(DAY0_bin_pAC,BINnum)

DAY0_bin_pAC_MAT = DAY0_bin_pAC_BINnum[which(DAY0_bin_pAC_BINnum$BINnum == 1),8]
DAY0_bin_pAC_MAT = as.data.frame(DAY0_bin_pAC_MAT)
colnames(DAY0_bin_pAC_MAT) = NULL

for(i in 2:nrow(start)){
    DAY0_bin_pAC_MATfrag = DAY0_bin_pAC_BINnum[which(DAY0_bin_pAC_BINnum$BINnum == i),8]
    DAY0_bin_pAC_MATfrag = as.data.frame(DAY0_bin_pAC_MATfrag)
    colnames(DAY0_bin_pAC_MATfrag) = NULL
    DAY0_bin_pAC_MAT = cbind(DAY0_bin_pAC_MAT,DAY0_bin_pAC_MATfrag)
    }

colnames(DAY0_bin_pAC_MAT) = NULL

DAY0_bin_pAC_MAT = as.matrix(DAY0_bin_pAC_MAT)
pmean <- function(x,y) (x+y)/2
DAY0_bin_pAC_MAT[] <- pmean(DAY0_bin_pAC_MAT, matrix(DAY0_bin_pAC_MAT, nrow(DAY0_bin_pAC_MAT), byrow=TRUE))     

    DAY0_bin_pAC_MAT = as.data.frame(DAY0_bin_pAC_MAT)

    write.table(DAY0_bin_pAC_MAT,
            file = output_path, 
            col.names=F, row.names=F, quote=F,sep="\t")
    }

GENE_PLOT = function(x){

DESeq2_normalized = read.table("/data3/psg/4DN_20241106/FiG1/3_DEGs/count/DESeq2_normalized.txt",sep = "\t", header=T)
DESeq2_normalized = DESeq2_normalized[,1:9]

DEG_FULL_P = read.table("/data3/psg/4DN_20241106/FiG1/3_DEGs/DEG_result/condition_DAY4_vs_DAY0_0.05_1.5/DEG_FULL.txt",sep = "\t", header=T,fill = TRUE)

geneDATA = DESeq2_normalized[which(DESeq2_normalized[,1] == x),]

DAY0_mean = (geneDATA[,2] + geneDATA[,3] + geneDATA[,4] + geneDATA[,5])/4
DAY0_sd = sd(c(geneDATA[,2],geneDATA[,3],geneDATA[,4], geneDATA[,5]))
DAY0_se = DAY0_sd/(4^(1/2))

DAY4_mean = (geneDATA[,6] + geneDATA[,7] + geneDATA[,8] + geneDATA[,9])/4
DAY4_sd = sd(c(geneDATA[,6],geneDATA[,7],geneDATA[,8], geneDATA[,9]))
DAY4_se = DAY4_sd/(4^(1/2))

DAY0_DATA = "DAY0"
DAY0_DATA = as.data.frame(DAY0_DATA)
DAY0_DATA$DAY = "DAY0"
DAY0_DATA$mean = DAY0_mean
DAY0_DATA$sd = DAY0_sd
DAY0_DATA$se = DAY0_se
colnames(DAY0_DATA) = c("Condition","DAY","mean","sd","se")  

DAY4_DATA = "DAY4"
DAY4_DATA = as.data.frame(DAY4_DATA)
DAY4_DATA$DAY = "DAY4"
DAY4_DATA$mean = DAY4_mean
DAY4_DATA$sd = DAY4_sd
DAY4_DATA$se = DAY4_se
colnames(DAY4_DATA) = c("Condition","DAY","mean","sd","se")      

RNAplotDATA = rbind(DAY0_DATA,DAY4_DATA)

geneDATA = DEG_FULL_P[which(DEG_FULL_P[,1] == x),]

    if(geneDATA[1,6] <= 0.05){
        if(geneDATA[1,6] > 0){
            anno = "*"}
        }
    if(geneDATA[1,6] <= 0.001){
        if(geneDATA[1,6] > 0.05){
            anno = "**"}
        }
    if(geneDATA[1,6] <= 0.0001){
            anno = "***"}

options(repr.plot.width = 2, repr.plot.height = 3, repr.plot.res = 1000, repr.plot.pointsize = 40)

p<-ggplot(RNAplotDATA,aes(x=DAY,y=mean,fill=Condition))
p2=p+geom_bar(stat="identity",position="dodge",width=0.6)+
scale_fill_manual(values = c("grey30","grey30"))+
theme_classic(base_size = 10) +
ggtitle(label = geneDATA[,1],
        subtitle = anno) +
ylab("Normalized gene expression\n") +
#scale_x_discrete(limits=c("control", "exosome", "exosome\nRunx3")) + 
theme(plot.title=element_text(face="bold.italic",hjust=0.5,size=20,color = "black"),
      plot.subtitle=element_text(face="bold.italic",hjust=0.5,size=10,color = "black"),
      axis.text.x=element_text(face="bold",size=13,color = "black"),
      axis.text.y=element_text(face="bold",size=12,color = "black"),
      axis.title.x = element_text(face="bold",size = 0,color = "darkblue"),
      axis.title.y = element_text(face="bold",size = 12,color = "black"),
      legend.title=element_text(face="bold",size=10), 
      legend.text=element_text(size=8),
      legend.position = "none") +
geom_errorbar(aes(ymin=mean-se,ymax=mean+se),
              position=position_dodge(0.9),width=0.2)
p2
    }

GENE_PLOT_AtoB = function(x){

DESeq2_normalized = read.table("/data3/psg/4DN_20241106/FiG1/3_DEGs/count/DESeq2_normalized.txt",sep = "\t", header=T)
DESeq2_normalized = DESeq2_normalized[,1:9]

DEG_FULL_P = read.table("/data3/psg/4DN_20241106/FiG1/3_DEGs/DEG_result/condition_DAY4_vs_DAY0_0.05_1.5/DEG_FULL.txt",sep = "\t", header=T,fill = TRUE)

geneDATA = DESeq2_normalized[which(DESeq2_normalized[,1] == x),]

DAY0_mean = (geneDATA[,2] + geneDATA[,3] + geneDATA[,4] + geneDATA[,5])/4
DAY0_sd = sd(c(geneDATA[,2],geneDATA[,3],geneDATA[,4], geneDATA[,5]))
DAY0_se = DAY0_sd/(4^(1/2))

DAY4_mean = (geneDATA[,6] + geneDATA[,7] + geneDATA[,8] + geneDATA[,9])/4
DAY4_sd = sd(c(geneDATA[,6],geneDATA[,7],geneDATA[,8], geneDATA[,9]))
DAY4_se = DAY4_sd/(4^(1/2))

DAY0_DATA = "DAY0"
DAY0_DATA = as.data.frame(DAY0_DATA)
DAY0_DATA$DAY = "DAY0"
DAY0_DATA$mean = DAY0_mean
DAY0_DATA$sd = DAY0_sd
DAY0_DATA$se = DAY0_se
colnames(DAY0_DATA) = c("Condition","DAY","mean","sd","se")  

DAY4_DATA = "DAY4"
DAY4_DATA = as.data.frame(DAY4_DATA)
DAY4_DATA$DAY = "DAY4"
DAY4_DATA$mean = DAY4_mean
DAY4_DATA$sd = DAY4_sd
DAY4_DATA$se = DAY4_se
colnames(DAY4_DATA) = c("Condition","DAY","mean","sd","se")      

RNAplotDATA = rbind(DAY0_DATA,DAY4_DATA)

geneDATA = DEG_FULL_P[which(DEG_FULL_P[,1] == x),]

    if(geneDATA[1,6] <= 0.05){
        if(geneDATA[1,6] > 0){
            anno = "*"}
        }
    if(geneDATA[1,6] <= 0.001){
        if(geneDATA[1,6] > 0.05){
            anno = "**"}
        }
    if(geneDATA[1,6] <= 0.0001){
            anno = "***"}

options(repr.plot.width = 1.5, repr.plot.height = 3, repr.plot.res = 1000, repr.plot.pointsize = 40)

p<-ggplot(RNAplotDATA,aes(x=DAY,y=mean,fill=Condition))
p2=p+geom_bar(stat="identity",position="dodge",width=0.6)+
scale_fill_manual(values = c("skyblue3","skyblue3"))+
theme_classic(base_size = 10) +
ggtitle(label = geneDATA[,1],
        subtitle = anno) +
ylab("Normalized gene expression\n") +
#scale_x_discrete(limits=c("control", "exosome", "exosome\nRunx3")) + 
theme(plot.title=element_text(face="bold.italic",hjust=0.5,size=20,color = "black"),
      plot.subtitle=element_text(face="bold.italic",hjust=0.5,size=10,color = "black"),
      axis.text.x=element_text(face="bold",size=13,color = "black"),
      axis.text.y=element_text(face="bold",size=12,color = "black"),
      axis.title.x = element_text(face="bold",size = 0,color = "darkblue"),
      axis.title.y = element_text(face="bold",size = 12,color = "black"),
      legend.title=element_text(face="bold",size=10), 
      legend.text=element_text(size=8),
      legend.position = "none") +
geom_errorbar(aes(ymin=mean-se,ymax=mean+se),
              position=position_dodge(0.9),width=0.2)
p2
    }

GENE_PLOT_BtoA = function(x){

DESeq2_normalized = read.table("/data3/psg/4DN_20241106/FiG1/3_DEGs/count/DESeq2_normalized.txt",sep = "\t", header=T)
DESeq2_normalized = DESeq2_normalized[,1:9]

DEG_FULL_P = read.table("/data3/psg/4DN_20241106/FiG1/3_DEGs/DEG_result/condition_DAY4_vs_DAY0_0.05_1.5/DEG_FULL.txt",sep = "\t", header=T,fill = TRUE)

geneDATA = DESeq2_normalized[which(DESeq2_normalized[,1] == x),]

DAY0_mean = (geneDATA[,2] + geneDATA[,3] + geneDATA[,4] + geneDATA[,5])/4
DAY0_sd = sd(c(geneDATA[,2],geneDATA[,3],geneDATA[,4], geneDATA[,5]))
DAY0_se = DAY0_sd/(4^(1/2))

DAY4_mean = (geneDATA[,6] + geneDATA[,7] + geneDATA[,8] + geneDATA[,9])/4
DAY4_sd = sd(c(geneDATA[,6],geneDATA[,7],geneDATA[,8], geneDATA[,9]))
DAY4_se = DAY4_sd/(4^(1/2))

DAY0_DATA = "DAY0"
DAY0_DATA = as.data.frame(DAY0_DATA)
DAY0_DATA$DAY = "DAY0"
DAY0_DATA$mean = DAY0_mean
DAY0_DATA$sd = DAY0_sd
DAY0_DATA$se = DAY0_se
colnames(DAY0_DATA) = c("Condition","DAY","mean","sd","se")  

DAY4_DATA = "DAY4"
DAY4_DATA = as.data.frame(DAY4_DATA)
DAY4_DATA$DAY = "DAY4"
DAY4_DATA$mean = DAY4_mean
DAY4_DATA$sd = DAY4_sd
DAY4_DATA$se = DAY4_se
colnames(DAY4_DATA) = c("Condition","DAY","mean","sd","se")      

RNAplotDATA = rbind(DAY0_DATA,DAY4_DATA)

geneDATA = DEG_FULL_P[which(DEG_FULL_P[,1] == x),]

    if(geneDATA[1,6] <= 0.05){
        if(geneDATA[1,6] > 0){
            anno = "*"}
        }
    if(geneDATA[1,6] <= 0.001){
        if(geneDATA[1,6] > 0.05){
            anno = "**"}
        }
    if(geneDATA[1,6] <= 0.0001){
            anno = "***"}

options(repr.plot.width = 2.5, repr.plot.height = 2.5, repr.plot.res = 1000, repr.plot.pointsize = 40)

p<-ggplot(RNAplotDATA,aes(x=DAY,y=mean,fill=Condition))
p2=p+geom_bar(stat="identity",position="dodge",width=0.6)+
scale_fill_manual(values = c("gold3","gold3"))+
theme_classic(base_size = 10) +
ggtitle(label = geneDATA[,1],
        subtitle = anno) +
ylab("Normalized gene expression\n") +
#scale_x_discrete(limits=c("control", "exosome", "exosome\nRunx3")) + 
theme(plot.title=element_text(face="bold.italic",hjust=0.5,size=20,color = "black"),
      plot.subtitle=element_text(face="bold.italic",hjust=0.5,size=10,color = "black"),
      axis.text.x=element_text(face="bold",size=13,color = "black"),
      axis.text.y=element_text(face="bold",size=12,color = "black"),
      axis.title.x = element_text(face="bold",size = 0,color = "darkblue"),
      axis.title.y = element_text(face="bold",size = 12,color = "black"),
      legend.title=element_text(face="bold",size=10), 
      legend.text=element_text(size=8),
      legend.position = "none") +
geom_errorbar(aes(ymin=mean-se,ymax=mean+se),
              position=position_dodge(0.9),width=0.2)
p2
    }

GENE_PLOT_AtoB = function(x){

DESeq2_normalized = read.table("/data3/psg/4DN_20241106/FiG1/3_DEGs/count/DESeq2_normalized.txt",sep = "\t", header=T)
DESeq2_normalized = DESeq2_normalized[,1:9]

DEG_FULL_P = read.table("/data3/psg/4DN_20241106/FiG1/3_DEGs/DEG_result/condition_DAY4_vs_DAY0_0.05_1.5/DEG_FULL.txt",sep = "\t", header=T,fill = TRUE)

geneDATA = DESeq2_normalized[which(DESeq2_normalized[,1] == x),]

DAY0_mean = (geneDATA[,2] + geneDATA[,3] + geneDATA[,4] + geneDATA[,5])/4
DAY0_sd = sd(c(geneDATA[,2],geneDATA[,3],geneDATA[,4], geneDATA[,5]))
DAY0_se = DAY0_sd/(4^(1/2))

DAY4_mean = (geneDATA[,6] + geneDATA[,7] + geneDATA[,8] + geneDATA[,9])/4
DAY4_sd = sd(c(geneDATA[,6],geneDATA[,7],geneDATA[,8], geneDATA[,9]))
DAY4_se = DAY4_sd/(4^(1/2))

DAY0_DATA = "DAY0"
DAY0_DATA = as.data.frame(DAY0_DATA)
DAY0_DATA$DAY = "DAY0"
DAY0_DATA$mean = DAY0_mean
DAY0_DATA$sd = DAY0_sd
DAY0_DATA$se = DAY0_se
colnames(DAY0_DATA) = c("Condition","DAY","mean","sd","se")  

DAY4_DATA = "DAY4"
DAY4_DATA = as.data.frame(DAY4_DATA)
DAY4_DATA$DAY = "DAY4"
DAY4_DATA$mean = DAY4_mean
DAY4_DATA$sd = DAY4_sd
DAY4_DATA$se = DAY4_se
colnames(DAY4_DATA) = c("Condition","DAY","mean","sd","se")      

RNAplotDATA = rbind(DAY0_DATA,DAY4_DATA)

geneDATA = DEG_FULL_P[which(DEG_FULL_P[,1] == x),]

    if(geneDATA[1,6] <= 0.05){
        if(geneDATA[1,6] > 0){
            anno = "*"}
        }
    if(geneDATA[1,6] <= 0.001){
        if(geneDATA[1,6] > 0.05){
            anno = "**"}
        }
    if(geneDATA[1,6] <= 0.0001){
            anno = "***"}

options(repr.plot.width = 2.5, repr.plot.height = 2.5, repr.plot.res = 1000, repr.plot.pointsize = 40)

p<-ggplot(RNAplotDATA,aes(x=DAY,y=mean,fill=Condition))
p2=p+geom_bar(stat="identity",position="dodge",width=0.6)+
scale_fill_manual(values = c("skyblue3","skyblue3"))+
theme_classic(base_size = 10) +
ggtitle(label = geneDATA[,1],
        subtitle = anno) +
ylab("Normalized gene expression\n") +
#scale_x_discrete(limits=c("control", "exosome", "exosome\nRunx3")) + 
theme(plot.title=element_text(face="bold.italic",hjust=0.5,size=20,color = "black"),
      plot.subtitle=element_text(face="bold.italic",hjust=0.5,size=10,color = "black"),
      axis.text.x=element_text(face="bold",size=13,color = "black"),
      axis.text.y=element_text(face="bold",size=12,color = "black"),
      axis.title.x = element_text(face="bold",size = 0,color = "darkblue"),
      axis.title.y = element_text(face="bold",size = 12,color = "black"),
      legend.title=element_text(face="bold",size=10), 
      legend.text=element_text(size=8),
      legend.position = "none") +
geom_errorbar(aes(ymin=mean-se,ymax=mean+se),
              position=position_dodge(0.9),width=0.2)
p2
    }

ORA_finder = function(inputPATH,outputPATH){
    
DF = read.table(inputPATH,sep = "\t", header=T,fill = TRUE)

ego <- enrichGO(gene          = DF[,1],
                OrgDb         = org.Mm.eg.db,
                keyType       = 'SYMBOL',
                ont           = "BP")
egoplotdata1 <- mutate(ego, richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))
egoplotdata2 <- mutate(egoplotdata1, BgRatio_num = as.numeric(sub("/\\d+", "", BgRatio))/1)
egoplotdata <- mutate(egoplotdata2, GeneRatio_num = as.numeric(sub("/\\d+", "", GeneRatio))/1)
egoplotdata = as.data.frame(egoplotdata)

write.table(egoplotdata,
            file = outputPATH,
            col.names=T, row.names=F, quote=F,sep="\t")
    }

Fig1

B

RAWDATA = read.table("/data3/psg/NGS_2025/4DN/HiC_data/FiG1/3_DEGs/count/DESeq2_normalized.txt",sep = "\t", header=T,fill = TRUE)

condition = c("DAY0","DAY0","DAY0","DAY0",
              "DAY4","DAY4","DAY4","DAY4",
              "DAY7","DAY7","DAY7","DAY7",
              "DAY10","DAY10","DAY10","DAY10")

# label = c("DAY0","","","",
#           "DAY4","","","",
#           "","","DAY7","",
#           "","","","DAY10")

label = c("","","","",
          "","","","",
          "","","","",
          "","","","")

#################################################################################################################
#################################################################################################################
#################################################################################################################

    
predata = t(RAWDATA)
predata=as.data.frame(predata)
predata_contents = predata
predata_contents = predata_contents
predata_contents = predata_contents[-1,]
predata_colname = as.vector(predata[1,])
colnames(predata_contents)=predata_colname
predata=predata_contents
data=as.data.frame(predata)

for(i in 1:ncol(data)){
    data[,i]=as.numeric(data[,i])
}

datapre = data
dataset=datapre[ , which(apply(datapre, 2, var) != 0)]
pca <- prcomp(dataset, center = TRUE,scale. = TRUE)
pca_perc=round(100*pca$sdev^2/sum(pca$sdev^2),1)


df_pca_data=data.frame(PC1 = pca$x[,1], 
                       PC2 = pca$x[,2], 
                       PC3 = pca$x[,3], 
                       sample = rownames(dataset), 
                       condition=condition)

options(repr.plot.width = 5, repr.plot.height = 5, repr.plot.res = 1000, repr.plot.pointsize = 20)

ggplot(df_pca_data,aes(PC1,PC2,color = condition))+
geom_point(size=5)+
labs(x=paste0("PC1 (",pca_perc[1],")"), 
     y=paste0("PC2 (",pca_perc[2],")"))+
ylim(-130,130)+
xlim(-130,130)+
theme_classic(base_size = 10) +
scale_color_manual(values = c("pink3","purple","skyblue3","tan3")) +
theme(plot.title=element_text(face="bold",hjust=0.5,size=25,color = "darkblue"),
      axis.text.x=element_text(face="bold",size=13),
      axis.text.y=element_text(face="bold",size=13),
      axis.title.x = element_text(face="bold",size = 15),
      axis.title.y = element_text(face="bold",size = 15),
      legend.title=element_text(face="bold",size=0,color = "darkblue"),
      legend.text=element_text(face="bold",size=8))+
guides(color="none")