#' @description: identify the target genes of TF

library(openxlsx)

library(circlize)

library(dplyr)

library(tibble)

library(data.table)

library(Signac)

library(Seurat)

library(stringr)

library(BSgenome.Hsapiens.UCSC.hg38)

library(chromVARmotifs)

data("human_pwms_v2")

library(ggpubr)

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

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

DefaultAssay(scATAC.data) <- "Peaks"

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

DefaultAssay(scRNA.data) <- "RNA"

celltype <- "Tumor"

cCREs <- readRDS(paste0("6.Co-Accessible/", celltype,"_peak_gene_correlation.rds")) # 4368 cCREs

#gene-cCRE link --- enhancer

cCREs <- cCREs %>% subset(distance_bp/1000>=10 & peak2_type != 'Promoter')

cCREs$target_gene <- cCREs$peak1_nearestGene

cCREs$peak_gene <- paste0(cCREs$Peak2, '_', cCREs$peak1_nearestGene)

cCREs$Peak1 <- gsub("_", "-", cCREs$Peak1)

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

# motif information

motif.info <- data.frame(originName = names(scATAC.data@assays$Peaks@motifs@motif.names), TF = unlist(scATAC.data@assays$Peaks@motifs@motif.names))

rownames(motif.info) <- NULL

#                           originName     TF

#   ENSG00000008196_LINE2_TFAP2B_D_N1 TFAP2B

#      ENSG00000008197_LINE6_TFAP2D_D TFAP2D

#   ENSG00000087510_LINE7_TFAP2C_D_N3 TFAP2C

#     ENSG00000116819_LINE19_TFAP2E_I TFAP2E

#  ENSG00000137203_LINE29_TFAP2A_D_N3 TFAP2A

#  ENSG00000116017_LINE44_ARID3A_D_N1 ARID3A

# subset by celltype

scATAC.sub <- subset(scATAC.data, cellType_low == celltype)

scRNA.sub <- subset(scRNA.data, cellType_low == celltype)

scATAC.sub <- FindTopFeatures(scATAC.sub, min.cutoff=ncol(scATAC.sub)*0.05) # present in what % of cells, 要求95%的细胞都包含peak

scRNA.sub <- FindVariableFeatures(scRNA.sub, nfeatures = 3000)

# use chipseeker annotate each peak in scATAC-seq

peak.info <- readRDS("4.Peak/peak.annotation.simple.ChIPseeker.rds")

# TF target genes

# for each motif, find genes:

use_variable_genes <- F

motif_list <- list()

edge_df.list <- data.frame()

vertex_df.list <- data.frame()

for(i in 1:nrow(motif.info)){

  motif_name <- motif.info[i,2]

  motif_ID <- motif.info[i,1]

  cat("Analysis ", motif_ID, "\n")

  # get list of promoter and enhancer targets of these TFs

  DefaultAssay(scATAC.data) <- "Peaks"

  motif_accessible <- Motifs(scATAC.data)@data[,motif_ID]

  motif_accessible <- names(motif_accessible)[motif_accessible > 0]

  motif_accessible <- motif_accessible[motif_accessible %in% VariableFeatures(scATAC.sub)]

  # promoter locate in 1kb

  motif_accessible_promoters <- motif_accessible[motif_accessible %in% peak.info$peaks[which(peak.info$originType == "Promoter (<=1kb)")]]

  motif_target_genes <- peak.info$SYMBOL[match(motif_accessible_promoters, peak.info$peaks)]

  motif_target_genes <- as.character(motif_target_genes)

  # variable genes only?

  if(use_variable_genes){

    motif_target_genes <- motif_target_genes[motif_target_genes %in% VariableFeatures(scRNA.sub)] %>% as.character %>% unique

  } else{

    motif_target_genes<- motif_target_genes %>% as.character %>% unique

  }

  # enhancer target genes

  motif_accessible_enhancers <- motif_accessible[motif_accessible %in% cCREs$Peak2]

  motif_cCREs <- subset(cCREs, Peak2 %in% motif_accessible_enhancers)

  motif_enhancers_target_genes <- motif_cCREs$target_gene

  # variable genes only?

  if(use_variable_genes){

    motif_enhancers_target_genes <- motif_enhancers_target_genes[motif_enhancers_target_genes %in% VariableFeatures(scRNA.sub)] %>% as.character %>% unique

  }else{

    motif_enhancers_target_genes<- motif_enhancers_target_genes %>% as.character %>% unique

  }

  vertex_df <- data.frame(name = c(motif_name, as.character(unique(c(motif_enhancers_target_genes, motif_target_genes)))))

  # check if there are promoter targets:

  if(length(motif_target_genes) > 0){

    promoter_edge_df <- data.frame(

      from = motif_name,

      to = as.character(motif_target_genes),

      type = 'promoter'

    )

  } else{promoter_edge_df <- data.frame()}

  # check if there are enhancer targets:

  if(length(motif_enhancers_target_genes) > 0){

    enhancer_edge_df <- data.frame(

      from = motif_name,

      to = as.character(motif_enhancers_target_genes),

      type = 'enhancer'

    )

  } else{enhancer_edge_df <- data.frame()}

  #cur_edge_df <- rbind(cur_promoter_edge_df, cur_enhancer_edge_df, cur_repressors_edge_df)

  edge_df <- rbind(promoter_edge_df, enhancer_edge_df)

  edge_df.list <- rbind(edge_df.list, edge_df)

  vertex_df.list <- rbind(vertex_df.list, vertex_df)

}

vertex_df.list <- data.frame(name=na.omit(as.character(unique(vertex_df.list$name))))

vertex_df.list$name <- as.character(vertex_df.list$name)

edge_df.list <- na.omit(edge_df.list)

saveRDS(edge_df.list, file = paste0("7.TF.analysis/targets/", celltype,".All.TF.targets.rds"))