#' @description: this script will annotate cicero ccan with the ChIPSeeker database to determine what genomic

# regions are linked by the connections and create circos plots

library(ChIPseeker)

library(TxDb.Hsapiens.UCSC.hg38.knownGene)

library(openxlsx)

library(circlize)

library(dplyr)

library(tibble)

library(data.table)

library(Signac)

library(Seurat)

library(stringr)

library(RColorBrewer)

library(cicero)

source(file = "/home/longzhilin/Analysis_Code/Plot_colorPaletters.R")

source(file = "/home/longzhilin/Analysis_Code/SingleCell/user.CoveragePlot.R")

setwd("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scATAC")

# read in cell-type-specific CCAN and filter by coaccess threshold > 0.2

dar_files <- list.files("6.Co-Accessible/ccans", recursive = TRUE, pattern = "ciceroConns.*.csv$", full.names = TRUE)

idents <- gsub(".*ciceroConns.", "", dar_files)

idents <- gsub(".csv", "", idents)

list.dar <- lapply(dar_files, function(file_path) {

  fread(file_path) %>%

    dplyr::filter(coaccess > 0.2) %>%

    dplyr::select(Peak1, Peak2)

})

names(list.dar) <- idents

# count cis-regulatory number of coaccess > 0.2

peak.num <- sapply(list.dar, function(x){return(nrow(x))})

names(peak.num) <- idents

# convert the DAR to GRanges objects to annotate

list.dar.peak1.gr <- lapply(seq(list.dar), function(x) {

  df <- list.dar[[x]]

  gr <- StringToGRanges(df$Peak1, sep = c("_","_"))

  return(gr)

})

names(list.dar.peak1.gr) <- idents

list.dar.peak2.gr <- lapply(seq(list.dar), function(x) {

  df <- list.dar[[x]]

  gr <- StringToGRanges(df$Peak2, sep = c("_","_"))

  return(gr)

})

names(list.dar.peak2.gr) <- idents

# annotate the list of GRanges DAR for each cell type with the peak location

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

list.peak1.Anno <- lapply(list.dar.peak1.gr, annotatePeak, TxDb = txdb,

                        tssRegion = c(-3000, 3000), annoDb="org.Hs.eg.db", verbose = FALSE)

saveRDS(list.peak1.Anno, file = "6.Co-Accessible/ccans/list.peak1.Anno.rds")                        

list.peak1.loc <- lapply(seq(list.peak1.Anno), function(x) {

  loc <- list.peak1.Anno[[x]]@anno$annotation 

  loc <- str_split(loc, pattern = " ", simplify = TRUE)[,1]

  loc <- str_replace(loc, pattern="3'", replacement="3p_UTR")

  loc <- str_replace(loc, pattern="5'", replacement="5p_UTR")

  loc <- str_replace(loc, pattern="Distal", replacement="Intergenic")

  return(loc)

})

list.peak2.Anno <- lapply(list.dar.peak2.gr, annotatePeak, TxDb = txdb,

                          tssRegion = c(-3000, 3000), annoDb="org.Hs.eg.db", verbose = FALSE)

saveRDS(list.peak2.Anno, file = "6.Co-Accessible/ccans/list.peak2.Anno.rds")

list.peak2.loc <- lapply(seq(list.peak2.Anno), function(x) {

  loc <- list.peak2.Anno[[x]]@anno$annotation 

  loc <- str_split(loc, pattern = " ", simplify = TRUE)[,1]

  loc <- str_replace(loc, pattern="3'", replacement="3p_UTR")

  loc <- str_replace(loc, pattern="5'", replacement="5p_UTR")

  loc <- str_replace(loc, pattern="Distal", replacement="Intergenic")

  return(loc)

})

# create df with peak1 and peak2 locations

list.peaks.loc.df <- lapply(seq(idents), function(x) {

  df <- cbind(list.peak1.loc[[x]], list.peak2.loc[[x]]) %>% as.data.frame()

  colnames(df) <- c("Peak1","Peak2")

  counts <- dplyr::count(df, Peak1, Peak2) %>% as.data.frame()

})

names(list.peaks.loc.df) <- idents

# generate circos plots of predicted chromatin-chromatin interactions and save to plots directory

Plot_Cicero_Anno <- function(ident, list.peaks.loc.df) {

  clusterID <- ident

  toplot <- as.data.frame(list.peaks.loc.df[[ident]])

  colnames(toplot) <- c("Peak1","Peak2","n")

  print(clusterID)

  # convert to an adjacency list with a value indicating the number of connections between

  # each of the unique genomic location pairs

  unique_combos <- !duplicated(t(apply(toplot, 1, sort)))

  toplot <- toplot[unique_combos, ]

  toplot <- toplot[order(toplot$n), ]

  # toplot$n <- log2(toplot$n)

  # grid.col = c('3p_UTR'="red", '5p_UTR'="black", Distal_Intergenic="blue",

  # Downstream="grey", Exon="purple", Intron="orange", Promoter="green")

  grid.col = brewer.pal(7, "Paired")

  circos.clear()

  pdf(paste0("6.Co-Accessible/ccans/plots/ccans.annotation.",clusterID,".pdf"), width=10, height=5)

  par(mar=c(0.5,0.5,0.5,0.5))

  circos.par(gap.after = 5)

  chordDiagram(toplot,

             grid.col = grid.col,

             annotationTrack = "grid",

             preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(toplot))))))

  # we go back to the first track and customize sector labels

  circos.track(track.index = 1, panel.fun = function(x, y) {

    circos.text(CELL_META$xcenter,

                CELL_META$ylim[1],

                CELL_META$sector.index,

                facing = "clockwise",

                niceFacing = TRUE,

                adj = c(0, 0.5))

   }, bg.border = NA) # here set bg.border to NA is important

  dev.off()

  }

# save circlize plots to file

lapply(idents, function(ident) {

  Plot_Cicero_Anno(ident, list.peaks.loc.df=list.peaks.loc.df)

  })

##############################identify specific cis-regulatory elements

# focus on the tumor or Macrophage cell 

cell.type <- "Tumor"

#Computing average expression & accessibility

scRNA <- readRDS("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/data.merge.pro.rds")

DefaultAssay(scRNA) <- "RNA"

Idents(scRNA) <- scRNA$cellType_low

scATAC <- readRDS("scATAC.data.pro.rds")

DefaultAssay(scATAC) <- "Peaks"

Idents(scATAC) <- scATAC$AnnotatedcellType

average_exp_all <- AverageExpression(scRNA, assays = "RNA")

average_acc_all <- AverageExpression(scATAC, assays = "Peaks")

# extract the same cell type between scRNA and scATAC

average_exp_all$RNA <- average_exp_all$RNA[, colnames(average_acc_all$Peaks)]

####require one of the peaks overlap a gene's promoter (1kb)

peak1.anno <- list.peak1.Anno[[cell.type]]@anno

peak2.anno <- list.peak2.Anno[[cell.type]]@anno

conns <- fread(paste0("6.Co-Accessible/ccans/ciceroConns.", cell.type,".csv")) # 250302

conns <- conns %>% dplyr::filter(coaccess > 0.2)

conns <- conns %>% mutate(peak1_type = peak1.anno$annotation, peak1_distanceToTSS = peak1.anno$distanceToTSS, peak1_nearestGene = peak1.anno$SYMBOL,

                                      peak2_type = peak2.anno$annotation, peak2_distanceToTSS = peak2.anno$distanceToTSS, peak2_nearestGene = peak2.anno$SYMBOL)

bg <- na.omit(unique(c(conns$peak1_nearestGene, conns$peak2_nearestGene))) # 作为background gene：17567

saveRDS(bg, file = paste0("6.Co-Accessible/bg.", cell.type,".rds")) # coaccess > 0.2

#select the peaks overlap the gene's promoter in conns

cCREs <- conns %>% dplyr::filter(str_detect(peak1_type, pattern = "Promoter \\(<=1kb\\)"))

cCREs <- cCREs[!is.na(cCREs$peak1_nearestGene),]

conns_list <- group_split(cCREs, peak1_nearestGene) #基于peak1_nearestGene分组，通常按字母顺序排序

names(conns_list) <- sort(unique(cCREs$peak1_nearestGene)) 

cTargetGenes <- names(conns_list)[names(conns_list) %in% rownames(average_exp_all$RNA)] # 11453 genes

# loop over all genes to compute correlation between accessibility & expression

df <- data.frame()

corr_list <- lapply(cTargetGenes, function(x){

  cat("calculated gene:", x, "\n")

  # subset by cur_gene

  cur_conns <- conns_list[[x]]

  cur_conns <- cur_conns[!(cur_conns$Peak2 %in% cur_conns$Peak1),]

  cur_conns <- subset(cur_conns, coaccess >= 0.2)

  # skip this gene if there are no co-accessible connections

  if(nrow(cur_conns) == 0){return(data.frame())}

  # get average exp and acc for this gene and peaks that are co-accessible

  cur_conns$Peak2 <- gsub("_", "-", cur_conns$Peak2)

  average_acc <- average_acc_all$Peaks[as.character(cur_conns$Peak2),, drop = F]

  average_exp <- average_exp_all$RNA[x,]

  # correlation between expression and accessibility:

  cor_mat <- apply(average_acc, 1, function(x){

    correlation <- cor.test(as.numeric(average_exp), as.numeric(x), method='pearson')

    data.frame("pval"=correlation$p.value, "pcc"=correlation$estimate)

  })

  # collapse individual correlation dfs, add info

  cor_df <- Reduce(rbind, cor_mat)

  cur_conns$pcc <- cor_df$pcc

  cur_conns$pval <- cor_df$pval

  return(cur_conns)

})

# combine into one df and remove incomplete entries:

df <- Reduce(rbind, corr_list)

df <- df[!is.na(df$pcc),]

# compute FDR:

df$FDR <- p.adjust(df$pval, method='fdr')

df <- df %>% dplyr::filter(FDR < 0.05) %>% arrange(desc(pcc))

df$peak2_type <- gsub(" .*", "", df$peak2_type)

df$peak2_type <- gsub("'", "_UTR", df$peak2_type)

# distance between peak and target gene

peak1_ranges <- StringToGRanges(gsub("_", "-", df$Peak1))

peak2_ranges <- StringToGRanges(gsub("_", "-", df$Peak2))

df$distance_bp <- abs(start(peak1_ranges) - start(peak2_ranges))

# write df to file:

write.table(df, file = paste0("6.Co-Accessible/", cell.type,"_peak_gene_correlation.csv"), sep = ",", quote=FALSE, row.names=F)

saveRDS(df, file = paste0("6.Co-Accessible/", cell.type,"_peak_gene_correlation.rds"))