#' @description peak calling

library(Signac)
library(Seurat)
library(GenomeInfoDb)
library(EnsDb.Hsapiens.v86) #---GRCh38 (hg38)
library(ggplot2)
library(patchwork)
set.seed(101)
library(GenomicRanges)
library(ggpubr)
library(tibble)
library(dplyr)
library(ComplexHeatmap)
library(circlize)
library(openxlsx)

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

scATAC.data <- readRDS("scATAC.data.rds")
Idents(scATAC.data) <- scATAC.data$AnnotatedcellType
DefaultAssay(scATAC.data) <- "ATAC"

####To call peaks on each annotated cell type, we can use the group.by argument
peaks <- CallPeaks(
  object = scATAC.data,
  group.by = "AnnotatedcellType",
  macs2.path = "/home/longzhilin/miniconda3/envs/SingleCell/bin/macs2",
  outdir = "/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210721/scATAC/4.Peak"
)
saveRDS(peaks, "4.Peak/cellType.peak.rds") 

library(future)
plan("multiprocess", workers = 10)
options(future.globals.maxSize = 100000 * 1024^2) 
# remove peaks on nonstandard chromosomes and in genomic blacklist regions
peaks <- keepStandardChromosomes(peaks, pruning.mode = "coarse")
peaks <- subsetByOverlaps(x = peaks, ranges = blacklist_hg38_unified, invert = TRUE)
# quantify counts in each peak
macs2_counts <- FeatureMatrix(
  fragments = Fragments(scATAC.data), # from cellranger fragment result
  features = peaks,
  cells = colnames(scATAC.data)
)
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotation) <- "UCSC"
genome(annotation) <- "hg38"
saveRDS(annotation, file = "4.Peak/annotation.rds")
# create a new assay using the MACS2 peak set and add it to the Seurat object
scATAC.data[["Peaks"]] <- CreateChromatinAssay(
  counts = macs2_counts,
  fragments = Fragments(scATAC.data),
  annotation = annotation,
  genome = "hg38"
)
DefaultAssay(scATAC.data) <- "Peaks"
scATAC.data <- RunTFIDF(scATAC.data)
gene.activities <- GeneActivity(scATAC.data)
# add the gene activity matrix to the Seurat object as a new assay and normalize it
scATAC.data[['Macs2ACTIVITY']] <- CreateAssayObject(counts = gene.activities)
scATAC.data <- NormalizeData(
  object = scATAC.data,
  assay = 'Macs2ACTIVITY',
  normalization.method = 'LogNormalize',
  scale.factor = median(scATAC.data$nCount_Macs2ACTIVITY)
)
saveRDS(scATAC.data, "scATAC.data.pro.rds") #### macs2 calling

pdf("4.Peak/CA9.macs2.pdf")
CoveragePlot(
  object = scATAC.data,
  region = "CA9",
  assay = "Peaks",
  features = "CA9",
  ranges.title = "MACS2",
  expression.assay = "Macs2ACTIVITY",
  annotation = TRUE,
  peaks = F,
  links = F
)
tile_plot <- TilePlot(
  object = scATAC.data,
  region = "CA9"
)
print(tile_plot)
dev.off()

pdf("4.Peak/CA9.pdf")
CoveragePlot(
  object = scATAC.data,
  assay = "ATAC",
  expression.assay = "ACTIVITY",
  region = "CA9",
  features = "CA9",
  annotation = TRUE,
  peaks = TRUE,
  links = TRUE
)
dev.off()

###coverage plot of marker genes
source("/home/longzhilin/Analysis_Code/SingleCell/FindRegion.R")
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
promoter <- getPromoters(TxDb=txdb, upstream=1000, downstream=1000)
plot.cellType <- rev(c("CD4+ T cell", "Treg", "CD8+ T cell", "NK/NKT cell", "B cell", "Macrophage", "Monocyte", "Mast cell", "Endothelium (VCAM1+)", "Endothelium (VCAM1-)", "Mesangial cell", "Tumor"))
scATAC.data$AnnotatedcellType <- factor(scATAC.data$AnnotatedcellType, levels = plot.cellType)
DefaultAssay(scATAC.data) <- "Peaks"
Idents(scATAC.data) <- scATAC.data$AnnotatedcellType

cell.type.markers <- read.table(file = "/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/2.Cluster/AnnotateCellType/cellMarker.txt", header = T, stringsAsFactors = F, sep = "\t")
genes <- cell.type.markers$Gene
genes <- genes[-19]

genes <- c("CD8A", "CD4", "GNLY", "MS4A1", "CD163", "S100A12", "TPSAB1", "PECAM1", "PDGFRB", "CA9")
pdf("2.Cluster/AnnotateCellType/cellType.coverage.plot.origin2.pdf", height = unit(3, "inches"))
res <- sapply(genes, function(x){
  cat(x, "...\n")
  regions <- FindRegion(object = scATAC.data, region = x, assay = "Peaks", extend.upstream = 1000, extend.downstream = 1000)
  idx <- data.frame(findOverlaps(regions, promoter))
  if(nrow(idx)>0){
    p <- CoveragePlot(
      object = scATAC.data,
      region = x,
      ranges.title = "MACS2",
      links = F,
      peaks = T,
      extend.upstream = 1000,
      extend.downstream = 1000,
      region.highlight = promoter[idx[,2],])
  }else{
    p <- CoveragePlot(
      object = scATAC.data,
      region = x,
      ranges.title = "MACS2",
      links = F,
      extend.upstream = 1000,
      extend.downstream = 1000,
      peaks = T)
  }
  print(p)
  return(regions)
})
dev.off()

pdf("2.Cluster/AnnotateCellType/cellType.coverage.plot2.pdf", width = unit(3, "inches"), height = unit(3, "inches"))
res <- sapply(genes, function(x){
  cat(x, "...\n")
  regions <- FindRegion(object = scATAC.data, region = x, assay = "Peaks", extend.upstream = 1000, extend.downstream = 1000)
  idx <- data.frame(findOverlaps(regions, promoter))
  if(nrow(idx)>0){
    p <- CoveragePlot(
      object = scATAC.data,
      region = x,
      ranges.title = "MACS2",
      links = F,
      peaks = F,
      extend.upstream = 1000,
      extend.downstream = 1000,
      region.highlight = promoter[idx[,2],])
  }else{
    p <- CoveragePlot(
      object = scATAC.data,
      region = x,
      ranges.title = "MACS2",
      links = F,
      extend.upstream = 1000,
      extend.downstream = 1000,
      peaks = F)
  }
  print(p)
  return(regions)
})
dev.off()

############################# identify differentially accessible chromatin regions between celltypes
DefaultAssay(scATAC.data) <- "Peaks"
Idents(scATAC.data) <- scATAC.data$AnnotatedcellType
idents <- as.character(levels(scATAC.data))
cellType.DARs <- FindAllMarkers(scATAC.data, 
                                test.use = 'LR',
                                logfc.threshold=0, 
                                min.pct = 0.05, # often necessary to lower the min.pct threshold
                                latent.vars = "peak_region_fragments")
cf <- ClosestFeature(scATAC.data, regions = rownames(cellType.DARs)) # Find the closest feature to a given set of genomic regions
cellType.DARs <- cbind(cellType.DARs, gene=cf$gene_name, gene_biotype = cf$gene_biotype, type = cf$type, distance=cf$distance)
colnames(cellType.DARs)[6:7] <- c("cellType", "genomicRegion")
saveFormat <- lapply(idents, function(x){
  index <- which(cellType.DARs$cellType == x)
  DARs <- cellType.DARs[index,]
  DARs.up <- DARs %>% filter(avg_log2FC>0) %>% arrange(desc(avg_log2FC))
  DARs.down <- DARs %>% filter(avg_log2FC<0) %>% arrange(avg_log2FC)
  DARs <- rbind(DARs.up, DARs.down)
  return(DARs)
})
write.xlsx(saveFormat, file = "4.Peak/celltype.all.DARs.xlsx", sheetName = idents, rowNames = F)
saveRDS(cellType.DARs, file = "4.Peak/cellType.all.DARs.rds")

#require logfc.threshold >= 0.25 & p_val_adj < 0.05
cellType.sig.pos.DARs <- cellType.DARs %>% filter(avg_log2FC >=0.25 & p_val_adj < 0.05) %>% arrange(desc(avg_log2FC)) # 31925 peaks
saveFormat <- lapply(idents, function(x){
  index <- which(cellType.sig.pos.DARs$cellType == x)
  DARs <- cellType.sig.pos.DARs[index,]
  DARs <- DARs %>% arrange(desc(avg_log2FC))
  return(DARs)
})
names(saveFormat) <- idents
write.xlsx(saveFormat, file = "4.Peak/celltype.sig.pos.DARs.xlsx", sheetName = idents, rowNames = F)
saveRDS(cellType.sig.pos.DARs, file = "4.Peak/cellType.sig.pos.DARs.rds")

#plot--- differentially accessible chromatin regions heatmap
sig.region <- cellType.sig.pos.DARs %>% select(genomicRegion) %>% distinct() 
#average fragment of each peak in each cell type
sig.region.mean <- AverageExpression(scATAC.data, features = sig.region$genomicRegion, assays = "Peaks")
sig.region.mean.scale <- scale(t(sig.region.mean$Peaks))
pdf("4.Peak/cellType.sig.pos.DAR.pdf")
Heatmap(sig.region.mean.scale, name = "z-score", show_column_dend = F, show_row_dend = F, show_column_names = F, row_names_gp = gpar(fontsize = 10), width = unit(10, "cm"), height = unit(8, "cm"))
dev.off()

##plot---DAR distribution
cellType.sig.pos.DARs.ratio <- as.data.frame(table(cellType.sig.pos.DARs$cellType))
cellType.sig.pos.DARs.ratio$Type <- rep("Lymphoid", nrow(cellType.sig.pos.DARs.ratio))
cellType.sig.pos.DARs.ratio$Type[c(4, 6, 12)] <- "Myeloid"
cellType.sig.pos.DARs.ratio$Type[c(3, 7, 8)] <- "Other"
cellType.sig.pos.DARs.ratio$Type[c(9)] <- "Tumor"
cellType.sig.pos.DARs.ratio$Type <- factor(cellType.sig.pos.DARs.ratio$Type, levels = c("Lymphoid", "Myeloid", "Tumor", "Other"))
pdf("4.Peak/cellType.sig.pos.DAR.ratio.pdf")
ggbarplot(cellType.sig.pos.DARs.ratio, x="Var1", y="Freq", fill = "Type", color = "Type",
          sort.by.groups=FALSE, sort.val = "desc", palette = colors <- c("#00A087", "#4DBBD5", "#E64B35", "#3C5488"),#不按组排序
          label = T, xlab = "", ylab = "Number of DAR") + rotate_x_text(60)
dev.off()

############################# cell tpye differentially accessible chromatin and genes
## load DEGs
scRNA.DEGs <- readRDS("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/2.Cluster/AnnotateCellType/cellType.sig.pos.DEGs.rds")
scRNA.DEGs <- scRNA.DEGs %>% filter(avg_log2FC >= 0.25 & p_val_adj < 0.05)
compared.idents <- as.character(levels(scATAC.data))

# calculate the intersection gene between scRNA-seq and scATAC-seq in same cell type
overlap.gene.list <- sapply(compared.idents, function(x){
  idx <- which(scRNA.DEGs$cluster == x)
  DEGs <- scRNA.DEGs$gene[idx]

  idx <- which(cellType.sig.pos.DARs$cellType == x)
  DARs <- unique(cellType.sig.pos.DARs$gene[idx])

  overlap <- intersect(DEGs, DARs)
  return(overlap)
})
saveRDS(overlap.gene.list, file = "4.Peak/DEG.DAR.overlap.genes.rds")

# calculate the differential ration
overlap.ratio <- sapply(names(overlap.gene.list), function(x){
  genes <- overlap.gene.list[[x]]

  #scRNA-seq
  idx <- which(scRNA.DEGs$cluster == x)
  DEGs <- length(scRNA.DEGs$gene[idx])
  DEGs.with.DARs <- length(genes)
  Prop.DEGs.with.DARs <- DEGs.with.DARs/DEGs

  #scATAC-seq
  idx <- which(cellType.sig.pos.DARs$cellType == x)
  DARs <- length(cellType.sig.pos.DARs$genomicRegion[idx])
  DARs.near.DEGs <- length(which(cellType.sig.pos.DARs$gene[idx] %in% genes))
  Prop.DARs.near.DEGs <- DARs.near.DEGs/DARs
  return(c(DEGs, DEGs.with.DARs, Prop.DEGs.with.DARs, DARs, DARs.near.DEGs, Prop.DARs.near.DEGs))
})
overlap.ratio <- t(overlap.ratio)
colnames(overlap.ratio) <- c("DEGs", "DEGs with DARs", "Prop DEGs with DARs", "DARs", "DARs near DEGs", "Prop DARs near DEGs")

#calculate the min max mean sd
res <- apply(overlap.ratio, 2, function(x){
  return(c(min(x), max(x), mean(x), sd(x)))
})
rownames(res) <- c("min", "max", "mean", "sd")
write.xlsx(list(overlap.ratio, res), file = "4.Peak/overlap.ratio.xlsx", sheetName = c("overlap of DEGs and DARs", "Statistics"), rowNames = T)