library(optparse)
library(doParallel)
library(foreach)
library(ggplot2)
library(gridExtra)
library(scales)
library(cowplot)
library(ggpubr)

#parameter list
option_list = list(make_option(c("-g", "--fgdata"), type = "character", default = NULL, help  = "Fusion gene data with FIS scores"),
                   make_option(c("-f","--featurepath"), type = "character", default = NULL, help = "File path to 44 human genomic feautures"),
                   make_option(c("-s","--chrsize"), type = "character", default = NULL, help = "Chromosome size file"),
                   make_option(c("-i","--featureinfo"), type = "character", default = NULL, help = "Feature infomation file"),
                   make_option(c("-o","--output"), type = "character", default = NULL, help = "Output png file"))         

#functions
overlap<-function(fgdata,filename,gene_size){
  options(bedtools.path = "/path/to/bedtools")
  library(bedtoolsr)
  #step1: intersect fgdata with feature
  featurepath<-paste(f2,filename,sep = "")
  outputpath<-paste("../genomecov_result/",filename,sep = "")
  feature<- read.table(featurepath ,header =F, stringsAsFactors = FALSE)
  intersect_result<-bt.intersect(a=fgdata, b=feature)
  
  if (length(intersect_result)==0){
      paste(filename,"has no overlap with fusion gene data",sep = " ")
  }else{
  #step2:map the location into 20kb
       cc<-c(1:nrow(intersect_result))
       for (i in cc){
         if ((intersect_result[i,9]>=0)&(intersect_result[i,10]<=10000)){
           intersect_result[i,12]=intersect_result[i,2]-(intersect_result[i,6]-5000)
           intersect_result[i,13]=intersect_result[i,3]-(intersect_result[i,6]-5000)
         }else if ((intersect_result[i,9]>=10000)&(intersect_result[i,10]<=20000)){
           intersect_result[i,12]=intersect_result[i,2]-(intersect_result[i,8]-15000)
           intersect_result[i,13]=intersect_result[i,3]-(intersect_result[i,8]-15000)
         }
       }
       intersect_result_in20kb<-cbind(intersect_result[,1],intersect_result[,12:13],intersect_result[,3:11],intersect_result[,2:3])
  
       #step3: sort file according to the chr location from low to high
       intersect_result_in20kb_sorted<-bt.sort(i=intersect_result_in20kb)
  
       #step4: genomecov to calculate the overlap depth
       genomecov_result<-bt.genomecov(i=intersect_result_in20kb_sorted, g=gene_size, d = T, split = T)
       write.table(genomecov_result,outputpath,row.names = F,col.names = F,sep = '\t',quote = F)
  }
}

deepthsum<- function(data){
  pos_count<-matrix(0,20000,2)
  pos_count[,1]<-c(1:20000)
  count_pos<-matrix(0,20000,1)
  for (i in 1:20000){
    a<-which(data[,2]==i)
    count_pos[i,1]<-sum(data[a,3]!=0)
    pos_count[i,2]<-sum(data[a,3])/count_pos[i,1]
  }
  pos_count[is.na(pos_count)] <- 0
  pos_count<-as.data.frame(pos_count)
  colnames(pos_count) <- c("x", "y")
  return(pos_count)
}

mytheme1<-theme_bw() +  theme(
  axis.text=element_text(size=6),
  axis.title.x=element_blank(), 
  axis.text.x=element_blank(), 
  axis.ticks.x=element_blank(), 
  axis.title.y=element_blank(), 
  axis.ticks.y=element_blank(),
  axis.text.y = element_text(size = 5),
  plot.margin = unit(c(0, 0.1, 0, 0.5), "cm"),
  panel.background = element_rect(fill = "#c8e5ee",colour = "#c8e5ee"),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank()
)

mytheme2<-theme_bw() +  theme(
  axis.text=element_text(size=6),
  axis.title.x=element_blank(), 
  axis.text.x=element_blank(), 
  axis.ticks.x=element_blank(), 
  axis.title.y=element_blank(), 
  axis.ticks.y=element_blank(),
  axis.text.y = element_text(size = 5),
  plot.margin = unit(c(0, 0.1, 0, 0.5), "cm"),
  panel.background = element_rect(fill = "#c7e1a3",colour = "#c7e1a3"),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank()
)

mytheme3<-theme_bw() +  theme(
  axis.text=element_text(size=6),
  axis.title.x=element_blank(), 
  axis.text.x=element_blank(), 
  axis.ticks.x=element_blank(), 
  axis.title.y=element_blank(), 
  axis.ticks.y=element_blank(),
  axis.text.y = element_text(size = 5),
  plot.margin = unit(c(0, 0.1, 0, 0.5), "cm"),
  panel.background = element_rect(fill = "#ffcf9e",colour = "#ffcf9e"),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank()
)

mytheme4<-theme_bw() +  theme(
  axis.text=element_text(size=6),
  axis.title.x=element_blank(), 
  axis.text.x=element_blank(), 
  axis.ticks.x=element_blank(), 
  axis.title.y=element_blank(), 
  axis.ticks.y=element_blank(),
  axis.text.y = element_text(size = 5),
  plot.margin = unit(c(0, 0.1, 0, 0.5), "cm"),
  panel.background = element_rect(fill = "#d4c6e2",colour = "#d4c6e2"),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank()
)

mytheme5<-theme_bw() +  theme(
  axis.text=element_text(size=6),
  axis.title.x=element_blank(), 
  axis.text.x=element_blank(), 
  axis.ticks.x=element_blank(), 
  axis.title.y=element_blank(), 
  axis.ticks.y=element_blank(),
  axis.text.y = element_text(size = 5),
  plot.margin = unit(c(0, 0.1, 0, 0.5), "cm"),
  panel.background = element_rect(fill = "#cfcfcf",colour = "#cfcfcf"),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank()
)

#data processing for visualization
args <- parse_args(OptionParser(option_list=option_list))
f1=args$fgdata
f2=args$featurepath
f3=args$chrsize
dir.create("../genomecov_result")
fisdata <- read.table(f1 ,header =F, stringsAsFactors = FALSE)
fgdata<-as.data.frame(matrix(nrow = 10*nrow(fisdata),ncol = 8))
#select top 1% fusion gene regions
for (i in c(1:nrow(fisdata))){
  data<-as.data.frame(matrix(nrow = 1000,ncol = 9))
  cc<-c(1:1000)
  for (j in cc){
    data[j,1]<-paste(fisdata[i,1],fisdata[i,5],sep = "_")
    data[j,2]<-fisdata[i,2]
    data[j,3]<-fisdata[i,3]
    data[j,4]<-fisdata[i,6]
    data[j,5]<-fisdata[i,7]
    data[j,6]<-20*(j-1)
    data[j,7]<-20*j
    data[j,8]<-paste(fisdata[i,9],fisdata[i,10],sep = "")
    data[j,9]<-fisdata[i,(j+10)]
  }
  data<-data[order(data[,9],decreasing=TRUE)[1:10],]
  fgdata[((10*(i-1)+1):(10*i)),]<-data[,1:8]
}

colnames(fgdata)<-c("fusiongene","Hchr","Hbp","Tchr","Tbp","start","end","sequence")
gene_size<- read.table(f3 ,header =F, stringsAsFactors = FALSE)

#calculate fusion gene's location in whole chromosome
cc<-c(1:nrow(fgdata))
for (i in cc){
  if ((fgdata[i,"start"]>=0)&(fgdata[i,"end"]<=10000)){
    fgdata[i,"chr"]<-fgdata[i,"Hchr"]
    fgdata[i,"newstart"]<-fgdata[i,"Hbp"]-5000+fgdata[i,"start"]
    fgdata[i,"newend"]<-fgdata[i,"Hbp"]-5000+fgdata[i,"end"]
  }else if ((fgdata[i,"start"]>=10000)&(fgdata[i,"end"]<=20000)){
    fgdata[i,"chr"]<-fgdata[i,"Tchr"]
    fgdata[i,"newstart"]<-fgdata[i,"Tbp"]-15000+fgdata[i,"start"]
    fgdata[i,"newend"]<-fgdata[i,"Tbp"]-15000+fgdata[i,"end"]
  }
}
fgdata<-cbind(fgdata[,9:11],fgdata[,1:8])

cl<-makeCluster(10)
registerDoParallel(cl)
file=list.files(f2)
foreach(file=file) %dopar% overlap(fgdata,file,gene_size)
stopCluster(cl)
cat("Data processing done, start visualization")

#visualization
pos_count<-matrix(0,20000,2)
pos_count[,1]<-c(1:20000)
data<-pos_count
data<-as.data.frame(data)
names(data) <- c("x", "y")
p_empty_blue<-ggplot(data, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill='#c8e5ee',color='#c8e5ee')  +mytheme1 +ylim(0,500)
p_empty_green<-ggplot(data, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill='#c7e1a3',color='#c7e1a3')  +mytheme2 +ylim(0,500)
p_empty_orange<-ggplot(data, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill='#ffcf9e',color='#ffcf9e')  +mytheme3 +ylim(0,500)
p_empty_orange2<-ggplot(data, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill='#ffcf9e',color='#ffcf9e')  +scale_x_continuous(breaks=c(0, 5000, 10000,15000,20000), labels=c("0", "5kb", "10kb","15kb","20kb"))+theme_bw()+ylim(0,500)+theme(axis.text=element_text(size=6),axis.title.x=element_blank(), axis.ticks.x=element_blank(), axis.title.y=element_blank(), axis.ticks.y=element_blank(),axis.text.y = element_text(size = 5),plot.margin = unit(c(0, 0.1, 0, 0.5), "cm"),panel.background = element_rect(fill = "#ffcf9e",colour = "#ffcf9e"),panel.grid.major = element_blank(),panel.grid.minor = element_blank())
p_empty_purple<-ggplot(data, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill='#d4c6e2',color='#d4c6e2')  +mytheme4 +ylim(0,500)
p_empty_grey<-ggplot(data, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill='#cfcfcf',color='#cfcfcf')  +mytheme5 +ylim(0,500)
p_empty_grey2<-ggplot(data, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill='#cfcfcf',color='#cfcfcf')  +scale_x_continuous(breaks=c(0, 5000, 10000,15000,20000), labels=c("0", "5kb", "10kb","15kb","20kb"))+theme_bw()+ylim(0,500)+theme(axis.text=element_text(size=6),axis.title.x=element_blank(), axis.ticks.x=element_blank(), axis.title.y=element_blank(), axis.ticks.y=element_blank(),axis.text.y = element_text(size = 5),plot.margin = unit(c(0, 0.1, 0, 0.5), "cm"),panel.background = element_rect(fill = "#cfcfcf",colour = "#cfcfcf"),panel.grid.major =element_blank(),panel.grid.minor = element_blank())

dat<-read.delim(args$featureinfo,stringsAsFactors=F)
file=list.files(f2)
fd<-as.data.frame(file)
fd$mark<-0

y_col<-c()
y_col[1:5]<-"black"     #gene_expression 5
y_col[6:20]<-"purple"     #chromatain_status 15 
y_col[21:25]<-"orange"      #SV types 5
y_col[26:38]<-"darkgreen"     #repeats 13 
y_col[39:44]<-"blue"     #viruses 6

dat$Feature = factor(dat$Feature, levels=c("TAD_boundary","ReplicationTiming","Promoters","Methylation","CPGisland","15_Quies","14_ReprPCWk","13_ReprPC","12_EnhBiv","11_BivFlnk","10_TssBiv","9_Het","8_ZNF_Rpts","7_Enh","6_EnhG","5_TxWk","4_Tx","3_TxFlnk","2_TssAFlnk","1_TssA","INV","INS","DUP","DEL","CNV","Z-DNA motifs","Mirror repeats","MIR repeats","Microsatellites","Low_complexity A/T rich regions","L2 repeats","L1 repeats","Inverted repeats","G-Quadruplex forming repeats","DNA transposons","Directed repeats","A-Phased repeats","Alu repeats","HTLV-1","HPV","HIV","HBV","EBV","AAV2"))

p1<-ggplot(dat, aes(x = Feature,y = Percentage))+
  geom_bar(stat ="identity",width = 0.6,position ="stack",fill='white',color='white')+ coord_flip()+
  theme(legend.position = "none",axis.text.x = element_blank(), axis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.y = element_text(size=11,colour = y_col), axis.title.x = element_blank(), axis.title.y = element_blank())+
  theme(panel.background = element_rect(fill = "white",colour = "white"),panel.grid.major = element_line(colour = "white",size = 0.03,linetype = "solid"),panel.grid.minor = element_line(colour = "white",size = 0.01,linetype = "dashed"))

outpath<-args$output
path<-"../genomecov_result/"
setwd(path)
  for (j in dir(path)){
    deepthdata<-read.table(j,header=F,sep="\t")
    result<-deepthsum(deepthdata)
    if (j=='AAV2.txt'){
      p_AAV2<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme1
    }else if(fd[which(fd$file=='AAV2.txt'),2]==1){
      p_AAV2<-p_AAV2
    }else{
      p_AAV2<-p_empty_blue
    }
    if (j=='EBV.txt'){
      p_EBV<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme1
    }else if(fd[which(fd$file=='EBV.txt'),2]==1){
      p_EBV<-p_EBV
    }else{
      p_EBV<-p_empty_blue
    }
    if (j=='HBV.txt'){
      p_HBV<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme1
    }else if(fd[which(fd$file=='HBV.txt'),2]==1){
      p_HBV<-p_HBV
    }else{
      p_HBV<-p_empty_blue
    }
    if (j=='HIV.txt'){
      p_HIV<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme1
    }else if(fd[which(fd$file=='HIV.txt'),2]==1){
      p_HIV<-p_HIV
    }else{
      p_HIV<-p_empty_blue
    }
    if (j=='HPV.txt'){
      p_HPV<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme1
    }else if(fd[which(fd$file=='HPV.txt'),2]==1){
      p_HPV<-p_HPV
    }else{
      p_HPV<-p_empty_blue
    }
    if (j=='HTLV_1.txt'){
      p_HTLV_1<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme1
    }else if(fd[which(fd$file=='HTLV_1.txt'),2]==1){
      p_HTLV_1<-p_HTLV_1
    }else{
      p_HTLV_1<-p_empty_blue
    }
    if (j=='Alu.txt'){
      p_Alu<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='Alu.txt'),2]==1){
      p_Alu<-p_Alu
    }else{
      p_Alu<-p_empty_green
    }
    if (j=='A-Phased_repeats.txt'){
      p_A_Phased<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='A-Phased_repeats.txt'),2]==1){
      p_A_Phased<-p_A_Phased
    }else{
      p_A_Phased<-p_empty_green
    }
    if (j=='Direct_repeats.txt'){
      p_Direct<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='Direct_repeats.txt'),2]==1){
      p_Direct<-p_Direct
    }else{
      p_Direct<-p_empty_green
    }
    if (j=='DNA.txt'){
      p_DNA<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='DNA.txt'),2]==1){
      p_DNA<-p_DNA
    }else{
      p_DNA<-p_empty_green
    }
    if (j=='G-Quadruplex.txt'){
      p_G_Quadruplex<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='G-Quadruplex.txt'),2]==1){
      p_G_Quadruplex<-p_G_Quadruplex
    }else{
      p_G_Quadruplex<-p_empty_green
    }
    if (j=='Inverted_repeats.txt'){
      p_Inverted<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='Inverted_repeats.txt'),2]==1){
      p_Inverted<-p_Inverted
    }else{
      p_Inverted<-p_empty_green
    }
    if (j=='L1.txt'){
      p_L1<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='L1.txt'),2]==1){
      p_L1<-p_L1
    }else{
      p_L1<-p_empty_green
    }
    if (j=='L2.txt'){
      p_L2<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='L2.txt'),2]==1){
      p_L2<-p_L2
    }else{
      p_L2<-p_empty_green
    }
    if (j=='Low_complexity.txt'){
      p_Low_complexity<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='Low_complexity.txt'),2]==1){
      p_Low_complexity<-p_Low_complexity
    }else{
      p_Low_complexity<-p_empty_green
    }
    if (j=='Microsatellites.txt'){
      p_Microsatellite<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='Microsatellites.txt'),2]==1){
      p_Microsatellite<-p_Microsatellite
    }else{
      p_Microsatellite<-p_empty_green
    }
    if (j=='MIR.txt'){
      p_MIR<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='MIR.txt'),2]==1){
      p_MIR<-p_MIR
    }else{
      p_MIR<-p_empty_green
    }
    if (j=='Mirror_repeats.txt'){
      p_Mirror<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='Mirror_repeats.txt'),2]==1){
      p_Mirror<-p_Mirror
    }else{
      p_Mirror<-p_empty_green
    }
    if (j=='Z_DNA.txt'){
      p_Z_DNA<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme2
    }else if(fd[which(fd$file=='Z_DNA.txt'),2]==1){
      p_Z_DNA<-p_Z_DNA
    }else{
      p_Z_DNA<-p_empty_green
    }
    if (j=='CNV.txt'){
      p_CNV<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme3
    }else if(fd[which(fd$file=='CNV.txt'),2]==1){
      p_CNV<-p_CNV
    }else{
      p_CNV<-p_empty_orange
    }
    if (j=='DEL.txt'){
      p_DEL<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme3
    }else if(fd[which(fd$file=='DEL.txt'),2]==1){
      p_DEL<-p_DEL
    }else{
      p_DEL<-p_empty_orange
    }
    if (j=='DUP.txt'){
      p_DUP<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red') +scale_x_continuous(breaks=c(0, 5000, 10000,15000,20000), labels=c("0", "5kb", "10kb","15kb","20kb")) + theme_bw() + theme(axis.text=element_text(size=6),axis.title.x=element_blank(), axis.ticks.x=element_blank(), axis.title.y=element_blank(), axis.ticks.y=element_blank(),axis.text.y = element_text(size = 5)) + theme(plot.margin = unit(c(0, 0.1, 0, 0.5), "cm"))+ theme(panel.background = element_rect(fill = "#ffcf9e",colour = "#ffcf9e"),panel.grid.major = element_blank(),panel.grid.minor = element_blank())
    }else if(fd[which(fd$file=='DUP.txt'),2]==1){
      p_DUP<-p_DUP
    }else{
      p_DUP<-p_empty_orange2
    }
    if (j=='INS.txt'){
      p_INS<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme3
    }else if(fd[which(fd$file=='INS.txt'),2]==1){
      p_INS<-p_INS
    }else{
      p_INS<-p_empty_orange
    }
    if (j=='INV.txt'){
      p_INV<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme3
    }else if(fd[which(fd$file=='INV.txt'),2]==1){
      p_INV<-p_INV
    }else{
      p_INV<-p_empty_orange
    }
    if (j=='1_TssA.txt'){
      p_1_TssA<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='1_TssA.txt'),2]==1){
      p_1_TssA<-p_1_TssA
    }else{
      p_1_TssA<-p_empty_purple
    }
    if (j=='2_TssAFlnk.txt'){
      p_2_TssAFlnk<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='2_TssAFlnk.txt'),2]==1){
      p_2_TssAFlnk<-p_2_TssAFlnk
    }else{
      p_2_TssAFlnk<-p_empty_purple
    }
    if (j=='3_TxFlnk.txt'){
      p_3_TxFlnk<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='3_TxFlnk.txt'),2]==1){
      p_3_TxFlnk<-p_3_TxFlnk
    }else{
      p_3_TxFlnk<-p_empty_purple
    }
    if (j=='4_Tx.txt'){
      p_4_Tx<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='4_Tx.txt'),2]==1){
      p_4_Tx<-p_4_Tx
    }else{
      p_4_Tx<-p_empty_purple
    }
    if (j=='5_TxWk.txt'){
      p_5_TxWk<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='5_TxWk.txt'),2]==1){
      p_5_TxWk<-p_5_TxWk
    }else{
      p_5_TxWk<-p_empty_purple
    }
    if (j=='6_EnhG.txt'){
      p_6_EnhG<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='6_EnhG.txt'),2]==1){
      p_6_EnhG<-p_6_EnhG
    }else{
      p_6_EnhG<-p_empty_purple
    }
    if (j=='7_Enh.txt'){
      p_7_Enh<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='7_Enh.txt'),2]==1){
      p_7_Enh<-p_7_Enh
    }else{
      p_7_Enh<-p_empty_purple
    }
    if (j=='8_ZNF_Rpts.txt'){
      p_8_ZNFRpts<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='8_ZNF_Rpts.txt'),2]==1){
      p_8_ZNFRpts<-p_8_ZNFRpts
    }else{
      p_8_ZNFRpts<-p_empty_purple
    }
    if (j=='9_Het.txt'){
      p_9_Het<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='9_Het.txt'),2]==1){
      p_9_Het<-p_9_Het
    }else{
      p_9_Het<-p_empty_purple
    }
    if (j=='10_TssBiv.txt'){
      p_10_TssBiv<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='10_TssBiv.txt'),2]==1){
      p_10_TssBiv<-p_10_TssBiv
    }else{
      p_10_TssBiv<-p_empty_purple
    }
    if (j=='11_BivFlnk.txt'){
      p_11_BivFlnk<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='11_BivFlnk.txt'),2]==1){
      p_11_BivFlnk<-p_11_BivFlnk
    }else{
      p_11_BivFlnk<-p_empty_purple
    }
    if (j=='12_EnhBiv.txt'){
      p_12_EnhBiv<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='12_EnhBiv.txt'),2]==1){
      p_12_EnhBiv<-p_12_EnhBiv
    }else{
      p_12_EnhBiv<-p_empty_purple
    }
    if (j=='13_ReprPC.txt'){
      p_13_ReprPC<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='13_ReprPC.txt'),2]==1){
      p_13_ReprPC<-p_13_ReprPC
    }else{
      p_13_ReprPC<-p_empty_purple
    }
    if (j=='14_ReprPCWk.txt'){
      p_14_ReprPCWk<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='14_ReprPCWk.txt'),2]==1){
      p_14_ReprPCWk<-p_14_ReprPCWk
    }else{
      p_14_ReprPCWk<-p_empty_purple
    }
    if (j=='15_Quies.txt'){
      p_15_Quies<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme4
    }else if(fd[which(fd$file=='15_Quies.txt'),2]==1){
      p_15_Quies<-p_15_Quies
    }else{
      p_15_Quies<-p_empty_purple
    }
    if (j=='CPGisland.txt'){
      p_CPGisland<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme5
    }else if(fd[which(fd$file=='CPGisland.txt'),2]==1){
      p_CPGisland<-p_CPGisland
    }else{
      p_CPGisland<-p_empty_grey
    }
    if (j=='Methylation.txt'){
      p_Methylation<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme5
    }else if(fd[which(fd$file=='Methylation.txt'),2]==1){
      p_Methylation<-p_Methylation
    }else{
      p_Methylation<-p_empty_grey
    }
    if (j=='promoters.txt'){
      p_promoters<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme5
    }else if(fd[which(fd$file=='promoters.txt'),2]==1){
      p_promoters<-p_promoters
    }else{
      p_promoters<-p_empty_grey
    }
    if (j=='ReplicationTiming_Avg.txt'){
      p_ReplicationTiming<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red')  +mytheme5
    }else if(fd[which(fd$file=='ReplicationTiming_Avg.txt'),2]==1){
      p_ReplicationTiming<-p_ReplicationTiming
    }else{
      p_ReplicationTiming<-p_empty_grey
    }
    if (j=='TAD_boundary.txt'){
      p_TAD_boundary<-ggplot(result, mapping = aes(x = x, y = y)) +geom_bar(stat = 'identity', position = 'identity',fill = 'red', colour = 'red') +scale_x_continuous(breaks=c(0, 5000, 10000,15000,20000), labels=c("0", "5kb", "10kb","15kb","20kb")) + theme_bw() + theme(axis.text=element_text(size=6),axis.title.x=element_blank(), axis.ticks.x=element_blank(), axis.title.y=element_blank(), axis.ticks.y=element_blank(),axis.text.y = element_text(size = 5)) + theme(plot.margin = unit(c(0, 0.1, 0, 0.5), "cm"))+ theme(panel.background = element_rect(fill = "#cfcfcf",colour = "#cfcfcf"),panel.grid.major = element_blank(),panel.grid.minor = element_blank())
    }else if(fd[which(fd$file=='TAD_boundary.txt'),2]==1){
      p_TAD_boundary<-p_TAD_boundary
    }else{
      p_TAD_boundary<-p_empty_grey2
    }
    fd[which(fd$file==j),2]=1
  }
  p2 <- ggarrange(p_AAV2,p_EBV,p_HBV,p_HIV,p_HPV,p_HTLV_1,p_Alu,p_A_Phased,p_Direct,p_DNA,p_G_Quadruplex,p_Inverted,p_L1,p_L2,p_Low_complexity,p_Microsatellite,p_MIR,p_Mirror,p_Z_DNA,p_CNV,p_DEL,p_DUP,nrow = 22, labels = c(" ","","","","","","","","","","","","","","","","","","","","","",""),align = "v")
  p3 <- ggarrange(p_INS,p_INV,p_1_TssA,p_2_TssAFlnk,p_3_TxFlnk,p_4_Tx,p_5_TxWk,p_6_EnhG,p_7_Enh,p_8_ZNFRpts,p_9_Het,p_10_TssBiv,p_11_BivFlnk,p_12_EnhBiv,p_13_ReprPC,p_14_ReprPCWk,p_15_Quies,p_CPGisland,p_Methylation,p_promoters,p_ReplicationTiming,p_TAD_boundary, nrow = 22, labels = c(" ","","","","","","","","","","","","","","","","","","","","","",""),align = "v")
  p4<-ggarrange(p1,p2,p3,ncol=3)
  p<-annotate_figure(p4, top = text_grob("Overlap between the top 1% FIS regions and 44 human \n genomic features across 20Kb fusion DNA sequence",hjust = 0.18, size = 12))
  ggsave(p,file=outpath,dpi = 90)