#' @description: identify the tumor-specific TFs

library(Signac)

library(Seurat)

library(BSgenome.Hsapiens.UCSC.hg38)

library(patchwork)

set.seed(101)

library(ggpubr)

library(openxlsx)

library(ComplexHeatmap)

library(circlize)

library(tidyverse)

library(matrixStats)

library(future)

plan("multiprocess", workers = 10)

options(future.globals.maxSize = 50000 * 1024^2) #50G

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

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

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

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

motif.info$originName <- gsub("_", "-", motif.info$originName)

#####基于chromVAR数据分析

cellType.motifs.chromVAR <- readRDS("5.Motif/cellType.motifs.chromVAR.human_pwms_v2.rds")

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

###################1.plot motifs deviation score heatmap and number ratioplot

top5 <- sapply(cellType.motifs.chromVAR, function(x){

    x <- arrange(x, desc(avg_log2FC))

    return(x$gene[1:3])

})

top5 <- unique(as.character(top5))

top5 <- c(top5, "EOMES", "TBX10")

####Plot --- cellType specific motif heatmap

sig.motifs <- cellType.motifs.chromVAR %>% 

                                bind_rows %>%

                                filter(avg_log2FC > 1 & p_val_adj < 0.05) %>%

                                select(gene) %>%

                                distinct()

idx <- match(sig.motifs$gene, motif.info$TF)

sig.motifNames <- motif.info$originName[idx]

motifs.avgExp <- AverageExpression(scATAC.data, features = sig.motifNames, assays = "chromvar")

motifs.avgExp <- motifs.avgExp$chromvar

rownames(motifs.avgExp) <- sig.motifs$gene

zScore <- function(x){(x - mean(x)) /sd(x)}

motifs.avgExp.scale <- apply(motifs.avgExp, 1, zScore) %>% t() # row: TF; column: cell type

mark.idx <- match(top5, rownames(motifs.avgExp.scale))

pdf("5.Motif/Analysis/chromVAR.cellType.sig.motifs.heatmap.pdf")

ha <- rowAnnotation(link = anno_mark(at = mark.idx, labels = top5, link_width = unit(2, "mm"), labels_gp = gpar(fontsize = 5), padding = unit(1, "mm")))

Heatmap(motifs.avgExp.scale, name = "Deviation score", 

        width = unit(6, "cm"), height = unit(8, "cm"), right_annotation = ha,

        row_names_gp = gpar(fontsize = 8), column_names_gp = gpar(fontsize = 8),

        show_row_dend = F, show_column_dend = F, show_column_names = T, show_row_names = F, 

        heatmap_legend_param = list(labels_gp = gpar(fontsize = 8), by_row = T))

Heatmap(motifs.avgExp.scale, name = "Deviation score", 

        width = unit(6, "cm"), height = unit(8, "cm"), left_annotation = ha, 

        row_names_gp = gpar(fontsize = 8), column_names_gp = gpar(fontsize = 8),

        show_row_dend = F, show_column_dend = F, show_column_names = T, show_row_names = F, 

        heatmap_legend_param = list(labels_gp = gpar(fontsize = 8), by_row = T))      

dev.off()

####Plot --- top5 cell-type specific TFs

cellType.sig.motifs <- lapply(names(cellType.motifs.chromVAR), function(x){

    sig <- cellType.motifs.chromVAR[[x]] %>% filter(avg_log2FC > 1 & p_val_adj < 0.05) %>% mutate(celltype = x)

    return(sig)

})

names(cellType.sig.motifs) <- names(cellType.motifs.chromVAR)

sig.motifs.num <- cellType.sig.motifs %>% bind_rows()

sig.motifs.num <- as.data.frame(table(sig.motifs.num$celltype))

pdf("5.Motif/Analysis/chromVAR.cellType.sig.motifs.barplot.pdf")

ggbarplot(sig.motifs.num, x="Var1", y="Freq", fill = "Var1", color = "Var1",

          sort.by.groups=FALSE, sort.val = "desc", #不按组排序

          label = T, xlab = "", ylab = "Cell Number") + theme(legend.position="none") + rotate_x_text(30)

dev.off()

###################2.Screen cell-specific enriched motifs

source(file = "/home/longzhilin/Analysis_Code/code/analysis.diff.survival.TCGA.R")

DESeq2.normalized_counts <- readRDS("/data/active_data/lzl/RenalTumor-20200713/Data/TCGA/KIRC/Result/DESeq2.normalized_counts.rds")

DESeq2.normalized_counts <- log2(DESeq2.normalized_counts+1)

DESeq2.result <- readRDS("/data/active_data/lzl/RenalTumor-20200713/Data/TCGA/KIRC/Result/DESeq2.result.rds")

clin.data <- readRDS("/data/active_data/lzl/RenalTumor-20200713/Data/TCGA/KIRC/Result/clin.data.rds")

# extreact the avg_log2FC and FDR value

order.TFs <- cellType.motifs.chromVAR[[1]]$gene

motifs.FC <- sapply(cellType.motifs.chromVAR, function(x){

  index <- match(order.TFs, x$gene)

  return(round(x$avg_log2FC[index], 2))

})

rownames(motifs.FC) <- order.TFs

motifs.fdr <- sapply(cellType.motifs.chromVAR, function(x){

  index <- match(order.TFs, x$gene)

  return(x$p_val_adj[index])

})

rownames(motifs.fdr) <- order.TFs

a <- motifs.fdr

a[which(a==0)] <- 2

motifs.fdr[which(motifs.fdr==0)] <- min(a)*0.001 ###4.940656e-324, Multiply the minimum value by 0.01, instead of 0

motifs.fdr <- round(-log10(motifs.fdr), 2)

#3.deviation score

all.motifs.avgExp <- AverageExpression(scATAC.data, features = motif.info$originName, assays = "chromvar")

all.motifs.avgExp <- all.motifs.avgExp$chromvar

rownames(all.motifs.avgExp) <- motif.info$TF

################################## screening Strategy

#1.deviation score: sd > median(sd) + 4*mad

#2.Tumor cells: fdr> 0.0001 & avg_log2FC>=4

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

avg.sd <- rowSds(all.motifs.avgExp)

names(avg.sd) <- rownames(all.motifs.avgExp)

sd.mad <- DoubleMAD(avg.sd)

mad.threshold <- median(avg.sd) + 4*sd.mad[2] #Take the right

#Only in cell type: FDR < 0.0001 & log2FC>=4; log2FC<1 in other cell types

sig.label <- sapply(order.TFs, function(x){

     if(avg.sd[x] > mad.threshold){

        fdr <- motifs.fdr[x,]

        FC <- motifs.FC[x,]

        sig.fdr <- which(fdr > -log10(0.0001))

        sig.fc <- which(FC >= 4)  

        others <- which(FC >= 1)

        sig.idx <- intersect(sig.fdr, sig.fc)

        common.len <- length(sig.idx)

        if(common.len==1 & length(others)==1){

        #if(common.len==1){

            return("specific")

        }else{

            return("common")

        }

    }else{

        return("low variation")

    }

})

cellType.high.specific.TF <- names(which(sig.label=="specific"))

heatmapEM.fdr <- motifs.fdr[cellType.high.specific.TF,]

heatmapEM.FC <- motifs.FC[cellType.high.specific.TF,]

cellType.sd <- avg.sd[cellType.high.specific.TF]

write.xlsx(list(FDR = heatmapEM.fdr, FC = heatmapEM.FC), file = "5.Motif/Analysis/cellType.specific.TFs.xlsx", sheetName = c("FDR", "Log2FC"), rowNames = T)

#### Tumor specific TFs

idx <- which(colnames(heatmapEM.FC) == "Tumor")

index <- which(heatmapEM.FC[, idx] >= 4)

tumor.specific.TFs <- data.frame(Name = rownames(heatmapEM.fdr)[index], FDR = heatmapEM.fdr[index, idx], avg_log2FC = heatmapEM.FC[index, idx], sd = cellType.sd[index]) # 49

heatmapEM.fdr.tumor <- heatmapEM.fdr[index,]

heatmapEM.FC.tumor <- heatmapEM.FC[index,]

saveRDS(tumor.specific.TFs, file = "5.Motif/Analysis/tumor.specific.TFs.rds")

## Filter the TF result

source(file = "/home/longzhilin/Analysis_Code/Visualization/Plot.EnhancedVolcano.R")

##plot --- heatmap

gene.labels <- c("HNF1A", "HNF1B", "HNF4G", "HNF4A", 

                 "HOXC5", "OTP", "ISL1", "VENTX", 

                 "HOXB5", "HOXB4", "HOXA4", "HOXB7", "HOXB3", "HOXB2",

                 "BARX2", "POU6F2", "JUNB", "JUN", "JUND", "BATF")

mark.idx <- match(gene.labels, rownames(heatmapEM.FC.tumor))

pdf("5.Motif/Analysis/tumor.specific.TFs.heatmap.pdf")

#Plot.EnhancedVolcano(TF.DESeq2, x = "log2FoldChange", y = "padj", select.num = nrow(TF.DESeq2), drawConnectors = F)

ha = rowAnnotation(FDR = anno_barplot(heatmapEM.fdr.tumor[,idx], border = F, gp = gpar(fill = "#e5e4e2"), direction = "reverse"))

col_fun <- colorRamp2(c(min(heatmapEM.FC.tumor), 0, max(heatmapEM.FC.tumor)), c("blue", "white", "red"))

ht2 <- Heatmap(heatmapEM.FC.tumor, left_annotation = ha, col = col_fun, 

               show_row_dend = F, show_column_dend = F, 

               column_names_gp = gpar(fontsize = 8), row_names_gp = gpar(fontsize = 8), 

               name = "avg_log2FC", width = unit(7.5, "cm"), height = unit(10, "cm"), 

               show_heatmap_legend = T, 

               heatmap_legend_param = list(col_fun = col_fun, title = "avg_log2FC", direction = "horizontal", title_position = "topcenter")) +

rowAnnotation(link = anno_mark(at = mark.idx, labels = gene.labels, link_width = unit(3, "mm"), labels_gp = gpar(fontsize = 8), padding = unit(1, "mm")))

draw(ht2, heatmap_legend_side = "top")

dev.off()

# survival in TCGA-KIRC data

write.xlsx(tumor.specific.TFs, file = "5.Motif/Analysis/tumor.specific.TFs.xlsx", sheetName = c("Tumor specific TFs"), rowNames = T)

pdf("5.Motif/Analysis/tumor.specific.TFs.survival.pdf")

tumor.TFs.signature.res <- analysis.diff.survival.TCGA(interest.gene = tumor.specific.TFs$Name, diff.gene.pro = DESeq2.result, exp.data.process = DESeq2.normalized_counts, clin.data = clin.data, EnhancedVolcano.plot = F, Box.plot = F, main = "tumor.specific.TFs", meta.signature = T, single.signature = T)

dev.off()

Idents(scRNA.data) <- scRNA.data$cellType_low

tumor.TFs.avgExp <- AverageExpression(scRNA.data, features = tumor.specific.TFs$Name, assays = "RNA", slot = "data")

tumor.TFs.avgExp <- scale(t(tumor.TFs.avgExp$RNA))

pdf("5.Motif/Analysis/tumor.TFs.expression.pdf")

Heatmap(tumor.TFs.avgExp, 

        show_row_dend = F, show_column_dend = F, 

        column_names_gp = gpar(fontsize = 8), row_names_gp = gpar(fontsize = 8), 

        name = "Expression", height = unit(6, "cm"), 

        show_heatmap_legend = T, 

        heatmap_legend_param = list(title = "Expression", direction = "horizontal", title_position = "topcenter"))

dev.off()