5_motifAnalysis.Rmd

---
title: "scATAC-seq Rat Metrial Glands"
author: "Ha T. H. Vu"
output: html_document
---
  
```{r setup, include=FALSE}
options(max.print = "75")
knitr::opts_chunk$set(
  echo = TRUE,
  collapse = TRUE,
  comment = "#>",
  fig.path = "Files/",
  fig.width = 15,
  prompt = FALSE,
  tidy = FALSE,
  message = FALSE,
  warning = TRUE
)
knitr::opts_knit$set(width = 75)
```

This is a documentation for analyses of scATAC-seq data, generated from rat metrial gland tissues on gestational day (GD) 15.5 and 19.5. <br>


## File prep:
We first define common peaks to regions enriched in TBC cluster at GD15.5 and GD19.5 that overlap by at least 50\%. <br>
```
module load bedtools2
bedtools intersect -a /work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd15.5-tbc.bed -b /work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-tbc.bed -f 0.5 -r -wa -wb > /work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/commonPeaks.bed
cd /work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/
cut -f1-4 commonPeaks.bed | sort -u > commonPeaks-gd15.5.bed
cut -f5-8 commonPeaks.bed | sort -u > commonPeaks-gd19.5.bed
```

Load libraries
```{r}
library(dplyr)
library(Signac)
library(Seurat)
library(ggplot2)
library(GenomicRanges)
library(patchwork)
set.seed(1234)
```

Associate common peaks to genes:
```{r}
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.rda")
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.ranges.rda")
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-ATACwithRNAlables.rda")
Idents(rats) <- rats@meta.data$predicted.id # change from old identities to those predicted by scRNA-seq
DefaultAssay(rats) <- 'MACSpeaks'
peaks <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/commonPeaks.bed", header=F)
commToGenes <- ClosestFeature(rats, regions = peaks$V8, annotation = rnor6.ranges)
write.table(commToGenes, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/commonPeaksToGenes.txt", quote = F, sep = "\t", row.names = F)
```

## Motif analysis:
We carry out motif analysis with TF motifs from rat, mouse and human database. Firstly, run with rat TFs:
```{r, eval=F}
library(JASPAR2020)
library(TFBSTools)
library(BSgenome.Rnorvegicus.Ensembl.rn6)

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-ATACwithRNAlables.rda")
Idents(rats) <- rats@meta.data$predicted.id # change from old identities to those predicted by scRNA-seq
DefaultAssay(rats) <- 'MACSpeaks'
peaks <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/commonPeaks.bed", header=F)
top.peaks <- unique(peaks$V8)
open.peaks <- AccessiblePeaks(rats, assay = "MACSpeaks", min.cells = 1)
rats <- RegionStats(rats[['MACSpeaks']], genome=BSgenome.Rnorvegicus.Ensembl.rn6)
meta.feature <- GetAssayData(rats, assay = "MACSpeaks", slot = "meta.features")
peaks.matched <- MatchRegionStats(meta.feature = meta.feature[open.peaks, ], query.feature = meta.feature[top.peaks, ], n = 50000, features.match = c("GC.percent", "sequence.length"))
#rats
pfm <- getMatrixSet(x = JASPAR2020, opts = list(species = "Rattus norvegicus", all_versions = FALSE))
motif.matrix <- CreateMotifMatrix(features = granges(rats), pwm = pfm, genome = BSgenome.Rnorvegicus.Ensembl.rn6, p.cutoff = 0.00005)
motif.positions <- motifmatchr::matchMotifs(pwms = pfm, subject = granges(rats), out = 'positions', genome = BSgenome.Rnorvegicus.Ensembl.rn6)
motif <- CreateMotifObject(data = motif.matrix, positions = motif.positions, pwm = pfm)
rats <- SetAssayData(object = rats, slot = 'motifs', new.data = motif)
enriched.motifs <- FindMotifs(object=rats, features=top.peaks, background=peaks.matched)
#write.table(enriched.motifs, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/common19.5-ratsMotifs.txt", quote=F)
#save(rats, file = "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/gd19.5_ATACwithRatPFM.rda")
```
Mouse motifs:
```{r, eval=F}
pfm <- getMatrixSet(x = JASPAR2020, opts = list(species = "Mus musculus", all_versions = FALSE))
motif.matrix <- CreateMotifMatrix(features = granges(rats), pwm = pfm, genome = BSgenome.Rnorvegicus.Ensembl.rn6, p.cutoff = 0.00005)
motif.positions <- motifmatchr::matchMotifs(pwms = pfm, subject = granges(rats), out = 'positions', genome = BSgenome.Rnorvegicus.Ensembl.rn6)
motif <- CreateMotifObject(data = motif.matrix, positions = motif.positions, pwm = pfm)
rats <- SetAssayData(object = rats, slot = 'motifs', new.data = motif)
enriched.motifs <- FindMotifs(object=rats, features=top.peaks, background=peaks.matched)
#write.table(enriched.motifs, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/common19.5-mouseMotifs.txt", quote=F)
#save(rats, file = "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/gd19.5_ATACwithMousePFM.rda")
```
Human motifs:
```{r, eval=F}
pfm <- getMatrixSet(x = JASPAR2020, opts = list(species = "Homo sapiens", all_versions = FALSE))
motif.matrix <- CreateMotifMatrix(features = granges(rats), pwm = pfm, genome = BSgenome.Rnorvegicus.Ensembl.rn6, p.cutoff = 0.00005)
motif.positions <- motifmatchr::matchMotifs(pwms = pfm, subject = granges(rats), out = 'positions', genome = BSgenome.Rnorvegicus.Ensembl.rn6)
motif <- CreateMotifObject(data = motif.matrix, positions = motif.positions, pwm = pfm)
rats <- SetAssayData(object = rats, slot = 'motifs', new.data = motif)
enriched.motifs <- FindMotifs(object=rats, features=top.peaks, background=peaks.matched)
#write.table(enriched.motifs, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/common19.5-humanMotifs.txt", quote=F)
#save(rats, file = "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/gd19.5_ATACwithHumanPFM.rda")
```

We compile the enriched motifs to one table.
```{r}
rat19.5 <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/common19.5-ratsMotifs.txt", header = T)
rat19.5$adj_pval <- p.adjust(rat19.5$pvalue, "BH")
rat19.5 <- rat19.5[rat19.5$adj_pval <= 0.05 & rat19.5$fold.enrichment >= 1.5,]
rat19.5$type <- "rat"

mouse19.5 <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/common19.5-mouseMotifs.txt", header = T)
mouse19.5$adj_pval <- p.adjust(mouse19.5$pvalue, "BH")
mouse19.5 <- mouse19.5[mouse19.5$adj_pval <= 0.05 & mouse19.5$fold.enrichment >= 1.5,]
mouse19.5$type <- "mouse"

human19.5 <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/common19.5-humanMotifs.txt", header = T)
human19.5$adj_pval <- p.adjust(human19.5$pvalue, "BH")
human19.5 <- human19.5[human19.5$adj_pval <= 0.05 & human19.5$fold.enrichment >= 1.5,]
human19.5$type <- "human"

motifs19.5 <- rbind(rat19.5, mouse19.5, human19.5)
```

Number of motifs
```{r}
length(motifs19.5$motif)
```

Next, filter and keep only transcription factors (TFs) with enriched motifs and have expression levels > 0.5 at both GD15.5 and GD19.5.
```{r}
mouseOrthologs <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/rat-mouseOrthologs.txt", header = T, sep = "\t")
mouseOrthologs <- as.data.frame(apply(mouseOrthologs, 2, toupper))
humanOrthologs <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/rat-humanOrthologs.txt", header = T, sep = "\t")
humanOrthologs <- as.data.frame(apply(humanOrthologs, 2, toupper))

#transforming functions
subVar2 <- function(x) {
  ifelse(grepl("var.2", x, ignore.case = TRUE), gsub("\\(var.2\\)", "", x, ignore.case = TRUE), x)
}
subVar3 <- function(x) {
  ifelse(grepl("var.3", x, ignore.case = TRUE), gsub("\\(var.3\\)", "", x, ignore.case = TRUE), x)
}
subColon <- function(x) {
  ifelse(grepl("::", x, ignore.case = TRUE), strsplit(x, "::"), x)
}

enriched <- motifs19.5[, c("motif", "motif.name", "type")]
enriched2 <- data.frame(motif=NA, motif.name=NA, type=NA, TFgenes=NA)
for (i in 1:nrow(enriched)) {
  if (enriched[i, "type"] == "rat") {
    temp <- enriched[i,]
    temp$TFgenes <- enriched[i, "motif.name"]
    enriched2 <- rbind(enriched2, temp)
  } else if (enriched[i, "type"] == "mouse") {
    temp <- enriched[i,]
    names <- sapply(enriched[i, "motif.name"], subVar2)
    names <- sapply(names, subVar3)
    names <- sapply(names, subColon)
    names <- unlist(names)
    df <- data.frame(tf=toupper(names))
    df <- inner_join(df, mouseOrthologs[, c("Gene.name", "Gene.name.1")], by = c("tf" = "Gene.name.1"))
    temp2 <- temp[rep(seq_len(nrow(temp)), each = length(df$Gene.name)), ]
    temp2$TFgenes <- df$Gene.name
    enriched2 <- rbind(enriched2, temp2)
  } else {
    temp <- enriched[i,]
    names <- sapply(enriched[i, "motif.name"], subVar2)
    names <- sapply(names, subVar3)
    names <- sapply(names, subColon)
    names <- unlist(names)
    df <- data.frame(tf=toupper(names))
    df <- inner_join(df, mouseOrthologs[, c("Gene.name", "Gene.name.1")], by = c("tf" = "Gene.name.1"))
    temp2 <- temp[rep(seq_len(nrow(temp)), each = length(df$Gene.name)), ]
    temp2$TFgenes <- df$Gene.name
    enriched2 <- rbind(enriched2, temp2)
  }
}
enriched2 <- enriched2[!is.na(enriched2$motif),]

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd19.5-4-5-6-7_res0.8.rda")
exp19.5 <- AverageExpression(data2)
exp19.5 <- as.data.frame(exp19.5$RNA)
rownames(exp19.5) <- toupper(rownames(exp19.5))
exp19.5$gene <- rownames(exp19.5)
g2 <- rownames(exp19.5[rownames(exp19.5) %in% enriched2$TFgenes & exp19.5$`6` > 0.5,])
enriched2 <- left_join(enriched2, exp19.5[,c("gene", "6")], by = c("TFgenes" = "gene"))
enriched2[is.na(enriched2$`6`), "6"] <- 0

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/scRNA-seqFiles/gd15.5_res0.8.rda")
exp15.5 <- AverageExpression(data)
exp15.5 <- as.data.frame(exp15.5$RNA)
rownames(exp15.5) <- toupper(rownames(exp15.5))
exp15.5$gene <- rownames(exp15.5)
enriched2 <- left_join(enriched2, exp15.5[,c("gene", "23")], by = c("TFgenes" = "gene"))
enriched2[is.na(enriched2$`23`), "23"] <- 0

temp <- left_join(motifs19.5, enriched2)
#write.table(temp, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/tableS3_1.txt",  sep = "\t", quote = F, row.names = F)

exp15.5 <- exp15.5[g2,]
exp15.5 <- exp15.5[exp15.5$`23` > 0.5,]

enriched2 <- enriched2[enriched2$TFgenes %in% rownames(exp15.5),]

motifs19.5 <- motifs19.5[motifs19.5$motif.name %in% enriched2$motif.name,]
#write.table(motifs19.5, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/finalMotifs-commonPeaksGD19.5coor.txt", sep = "\t", quote = F)
```

Number of motifs retained
```{r}
length(motifs19.5$motif)
```

Find peaks that contain the enriched motifs of the expressed TFs. Since the types of the retained motifs were mouse and human, we will load the related objects only.
```{r}
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.ranges.rda")
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-ATACwithRNAlables.rda")
Idents(rats) <- rats@meta.data$predicted.id # change from old identities to those predicted by scRNA-seq
DefaultAssay(rats) <- 'MACSpeaks'
Annotation(rats) <- rnor6.ranges
peaks <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/commonPeaks-gd19.5.bed", header=F)
closest_genes <- ClosestFeature(rats, regions = peaks$V4)

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/gd19.5_ATACwithMousePFM.rda")

temp <- motifs19.5[motifs19.5$type == "mouse",]
motifPeaks <- data.frame(motif=NA, peak=NA)
for (i in 1:nrow(temp)) {
  p <- intersect(closest_genes$query_region, names(which(GetMotifData(rats)[, temp$motif[i]] == 1)))
  temp1 <- data.frame(motif=rep(temp$motif[i], length(p)), 
                      peak=p)
  motifPeaks <- rbind(motifPeaks, temp1)
}

load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/GD19.5_motifAnalysis/gd19.5_ATACwithHumanPFM.rda")

temp <- motifs19.5[motifs19.5$type == "human",]
for (i in 1:nrow(temp)) {
  p <- intersect(closest_genes$query_region, names(which(GetMotifData(rats)[, temp$motif[i]] == 1)))
  temp1 <- data.frame(motif=rep(temp$motif[i], length(p)), 
                      peak=p)
  motifPeaks <- rbind(motifPeaks, temp1)
}

motifPeaks <- motifPeaks[!is.na(motifPeaks$motif),]
#motifPeaks$peak <- paste0("chr", motifPeaks$peak)
#motifPeaks$peak <- sub("-", ":", motifPeaks$peak)
#write.table(motifPeaks, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/commonPeaksGD19.5withFinalMotifs.txt", sep = "\t", row.names = F, quote = F)
```

Number of peaks that contained the retained motifs
```{r}
length(unique(motifPeaks$peak))
```

## Compare with TFAP2C data from mouse stem cells
```
cd /work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks
ml bedtools2

# common peaks with motifs and overlapping with ChIP
sed 's/^/chr/g' commonPeaks-gd19.5-hg38Coor.bed | bedtools intersect -a - -b extraData/GSM3019344_dTSC_Tfap2c.peak-mm10toHg38.txt -wa -u | cut -f 4 | grep -wf - commonPeaksGD19.5withFinalMotifs.txt | grep -E "MA0524.2|MA0524.2|MA0814.2" | cut -f 2 | sort -u | wc -l
#208

# common peaks with motifs
grep -E "MA0524.2|MA0524.2|MA0814.2" commonPeaksGD19.5withFinalMotifs.txt | cut -f 2 | grep -wf - commonPeaks-gd19.5-hg38Coor.bed | sed 's/^/chr/g' | wc -l
#439

# common peaks with ChIP
sed 's/^/chr/g' commonPeaks-gd19.5-hg38Coor.bed | bedtools intersect -a - -b extraData/GSM3019344_dTSC_Tfap2c.peak-mm10toHg38.txt -wa -u | wc -l
#394

# common peaks
wc -l commonPeaks-gd19.5-hg38Coor.bed
#911
```

```{r}
dat <- data.frame(`Rat_peaks_overlapping_with_TFAP2C_ChIP-seq`=c(208,394-208), `Rat_peaks_not_overlapping_with_TFAP2C_ChIP-seq`=c(439-208,911-439-(394-208)))
rownames(dat) <- c("Rat peaks with TFAP2C motifs", "Rat peaks without TFAP2C motifs")

fisher.test(dat, alternative = "greater")
```

## Shared genes
First, merge families that shared any motif members:
```{r}
library(dplyr)
library(ggplot2)
library(ComplexHeatmap)


###merge family
toGenes <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/commonPeaks-gd19.5-toGenes.txt", header = T, sep = "\t")

motifPeaks <- inner_join(motifPeaks, motifs19.5[,c("motif", "motif.name")], by = c("motif" = "motif"))

jaspar <- read.table("/work/LAS/geetu-lab/hhvu/JASPAR-library.txt", header = T, sep = "\t")
motifPeaks <- inner_join(motifPeaks, jaspar, by = c("motif" = "ID"))


subColon <- function(x) {
  ifelse(grepl("::", x, ignore.case = TRUE), strsplit(x, "::"), x)
}
keep <- motifPeaks$Family
keep <- subColon(keep)
for (i in 1:length(keep)) {
  keep[i] <- list(unique(keep[i][[1]]))
  for (j in setdiff(1:length(keep), which(keep %in% keep[i]))) {
    if (sum(keep[j][[1]] %in% keep[i][[1]]) > 0) {
      keep[j] <- keep[i]
    }
  }
}
keep2 <- unlist(lapply(keep, toString))
motifPeaks$Family2 <- keep2

motifPeaks <- inner_join(motifPeaks, toGenes, by = c("peak" = "query_region"))
#write.table(motifPeaks, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/5_integrating2timepoints/annotation_scRNAseq/commonPeaks/tableS3_2.txt",  sep = "\t", quote = F, row.names = F)
```
Build a data frame for testing
```{r}
df = data.frame(from = rep(unique(motifPeaks$Family2), times = 11),
                to = rep(unique(motifPeaks$Family2), each = 11),
                q = NA,
                stringsAsFactors = FALSE)

#test significance
#number of samples
k <- c()
for (i in df$from) {
  sub <- motifPeaks[motifPeaks$Family2 == i,]
  k <- c(k, length(unique(sub$gene_name)))
}
df$k <- k #number of genes regulated by family A

#number of successes
m <- c()
for (i in df$to) {
  sub <- motifPeaks[motifPeaks$Family2 == i,]
  m <- c(m, length(unique(sub$gene_name)))
}
df$m <- m #number of genes regulated by family B

#number of failures
df$n <- length(unique(toGenes$gene_name)) - df$m #number of peaks without binding sites of family B

#number of genes regulated by both families
q <- c()
for (i in unique(df$to)) {
  for (j in unique(df$from)) {
    sub1 <- motifPeaks[motifPeaks$Family2 == i,]
    sub2 <- motifPeaks[motifPeaks$Family2 == j,]
    q <- c(q, length(intersect(sub1$gene_name, sub2$gene_name)))
  }
}
df$q <- q #number of genes regulated by both families

df <- df[df$from != df$to,]
sig <- phyper(df$q, df$m, df$n, df$k, lower.tail = FALSE, log.p = FALSE)
df$upper.adj.pval <- p.adjust(sig, "BH")
```

Build data frame for plotting
```{r}
test <- data.frame(matrix(ncol=11, nrow=11))
rownames(test) <- sort(unique(df$from))
colnames(test) <- sort(unique(df$to))
for (i in rownames(test)) {
  for (j in colnames(test)) {
    if (i != j) {
      test[i, j] <- df[df$from == i & df$to == j, "upper.adj.pval"] 
    } else {
      test[i, j] <- NA
    }
  }
}
test <- apply(test, 1, as.numeric)
rownames(test) <- sort(unique(df$from))
keepForPlot <- test
```

## Shared locations
```{r}
dat <- data.frame(matrix(ncol = length(unique(motifPeaks$Family2)), nrow = 0))

for (i in unique(motifPeaks$Family2)) {
  sub <- motifPeaks[motifPeaks$Family2 == i,]
  sub2 <- motifPeaks[motifPeaks$peak %in% unique(sub$peak),]
  t <- rep(0, length(unique(motifPeaks$Family2)))
  names(t) <- unique(motifPeaks$Family2)
  for (j in unique(motifPeaks$Family2)) {
    sub3 <- sub2[sub2$Family2 == j,]
    t[j] <- length(unique(sub3$peak))
  }
  
  dat <- rbind(dat, t)
}
rownames(dat) <- unique(motifPeaks$Family2)
colnames(dat) <- unique(motifPeaks$Family2)

dat <- as.matrix(dat)

df = data.frame(from = rep(rownames(dat), times = ncol(dat)),
                to = rep(colnames(dat), each = nrow(dat)),
                value = as.vector(dat),
                stringsAsFactors = FALSE)
df2 <- df[df$from != df$to,]

#test significance
#number of samples
k <- c()
for (i in df2$from) {
  sub <- motifPeaks[motifPeaks$Family2 == i,]
  k <- c(k, length(unique(sub$peak)))
}
df2$k <- k #number of peaks with binding sites of family A

#number of successes
m <- c()
for (i in df2$to) {
  sub <- motifPeaks[motifPeaks$Family2 == i,]
  m <- c(m, length(unique(sub$peak)))
}
df2$m <- m #number of peaks with binding sites of family B

#number of failures
df2$n <- 1242 - df2$m #number of peaks without binding sites of family B

#df2$value = number of peaks with binding sites of both families

sig <- phyper(df2$value, df2$m, df2$n, df2$k, lower.tail = FALSE, log.p = FALSE)
df2$upper.adj.pval <- p.adjust(sig, "BH")

test <- data.frame(matrix(ncol=11, nrow=11))
rownames(test) <- sort(unique(df2$from))
colnames(test) <- sort(unique(df2$to))
for (i in rownames(test)) {
  for (j in colnames(test)) {
    if (i != j) {
      test[i, j] <- df2[df2$from == i & df2$to == j, "upper.adj.pval"] 
    } else {
      test[i, j] <- NA
    }
  }
}
test <- apply(test, 1, as.numeric)
rownames(test) <- sort(unique(df2$from))
```

Combine two analyses in one plot
```{r}
test[upper.tri(test)] <- keepForPlot[upper.tri(keepForPlot)]
col1 = circlize::colorRamp2(c(1, 0), c("#fcae91", "#a50f15"))
col2 = circlize::colorRamp2(c(1, 0), c("#eff3ff", "#084594"))

ht = Heatmap(test, rect_gp = gpar(type = "none"), show_heatmap_legend = FALSE,
             cluster_rows = FALSE, cluster_columns = FALSE,
             layer_fun = function(j, i, x, y, w, h, fill) {
               l = i > j
               grid.rect(x[l], y[l], w[l], h[l], 
                         gp = gpar(fill = col1(pindex(test, i[l], j[l])), col = NA))
               
               l = i < j
               grid.rect(x[l], y[l], w[l], h[l], 
                         gp = gpar(fill = col2(pindex(test, i[l], j[l])), col = NA))
               
               l = i == j
               grid.rect(x[l], y[l], w[l], h[l], 
                         gp = gpar(fill = "grey", col = NA))
               
               v = pindex(test, i, j)
               grid.text(sprintf("%.3f", v), x, y, gp = gpar(fontsize = 15))
             })

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/fig3b-heatmap-combined.pdf", width=12, height=10)
draw(ht, heatmap_legend_list = list(
  Legend(title = "Shared binding location representation", col_fun = col1),
  Legend(title = "Shared gene representation", col_fun = col2)
))
#dev.off()
```

Distribution of number of TF fams per peak
```{r}
test <- motifPeaks[,c("peak", "Family2")]
test <- distinct(test)
t <- as.data.frame(table(test$peak))

#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figure3c-famPerPeak.pdf", width = 5, height = 6)
ggplot(t, aes(x = Freq)) + geom_histogram(bins = 8, fill = "#c6dbef", color="grey") +
  scale_x_continuous(breaks=c(0:8)) +
  stat_bin(aes(y=..count.., label=ifelse(..count..==0,"",..count..)), geom="text", vjust=-.5) +
  theme_bw() +
  ggtitle("Distribution of family number per peak")
#dev.off()
```