Chapter 26 Call Peaks and Identify Marker Peaks

library(ArchR)
set.seed(1)

26.1 Description

call peaks using the annotated types and find differential accessible peaks

26.2 Set env and load arrow project

## Section: set default para
##################################################
addArchRThreads(threads = 16) # setting default number of parallel threads


## Section: load object
##################################################
proj <- loadArchRProject(path = "data/ArchR/ArrowProject/Merged/")

26.3 Call peaks

## create pseudo-bulk replicates
proj <- addGroupCoverages(
  ArchRProj = proj,
  groupBy = "predictedGroup_Un",
  minCells = 15,
  maxCells = 700, 
  force = TRUE
)

## call peaks
proj <- addReproduciblePeakSet(
  ArchRProj = proj,
  groupBy = "predictedGroup_Un",
  genomeSize = 2.7e09,
  pathToMacs2 = "/lustre/user/liclab/liuyt/software/anaconda2/bin/macs2",
  force = TRUE
)

## calculate peak matrix
proj <- addPeakMatrix(proj, binarize = TRUE, force = TRUE)

26.4 Identify Marker Peaks

id <- unique(proj$predictedGroup_Un)
id <- id[ id %ni% c('Cycling', 'InPSB')]
markersPeaks <- getMarkerFeatures(
  ArchRProj = proj,
  useMatrix = "PeakMatrix",
  groupBy = "predictedGroup_Un",
  useGroups = id,
  bias = c("TSSEnrichment", "log10(nFrags)"),
  testMethod = "wilcoxon"
)

## ArchR logging to : ArchRLogs/ArchR-getMarkerFeatures-34806b3ca245-Date-2021-07-26_Time-21-02-37.log
## If there is an issue, please report to github with logFile!

## MatrixClass = Sparse.Binary.Matrix

## 2021-07-26 21:02:39 : Matching Known Biases, 0.019 mins elapsed.

## ###########
## 2021-07-26 21:03:12 : Completed Pairwise Tests, 0.566 mins elapsed.
## ###########

## ArchR logging successful to : ArchRLogs/ArchR-getMarkerFeatures-34806b3ca245-Date-2021-07-26_Time-21-02-37.log

markersPeaks

## class: SummarizedExperiment 
## dim: 253044 10 
## metadata(2): MatchInfo Params
## assays(7): Log2FC Mean ... AUC MeanBGD
## rownames(253044): 1 2 ... 253043 253044
## rowData names(4): seqnames idx start end
## colnames(10): End vRG ... IP oRG
## colData names(0):

markerList <-
  getMarkers(markersPeaks, cutOff = "FDR <= 0.01 & Log2FC >= 1")
markerList

## List of length 10
## names(10): End vRG Ex-SP Ex InMGE InCGE Ex-U OPC IP oRG

26.5 Visualize marker features

heatmapPeaks <- markerHeatmap(seMarker = markersPeaks,
                              cutOff = "FDR <= 0.1 & Log2FC >= 0.5",
                              transpose = TRUE)

## ArchR logging to : ArchRLogs/ArchR-plotMarkerHeatmap-2608142ca21-Date-2021-11-12_Time-15-00-04.log
## If there is an issue, please report to github with logFile!

## Identified 3497 markers!

##  [1] "chr1:66562-67062"         "chr1:5049504-5050004"    
##  [3] "chr1:5379884-5380384"     "chr1:5738993-5739493"    
##  [5] "chr1:6114707-6115207"     "chr1:6455298-6455798"    
##  [7] "chr1:6696209-6696709"     "chr1:8039395-8039895"    
##  [9] "chr1:9686053-9686553"     "chr1:11478011-11478511"  
## [11] "chr1:12094155-12094655"   "chr1:12338999-12339499"  
## [13] "chr1:12471345-12471845"   "chr1:12471904-12472404"  
## [15] "chr1:12781828-12782328"   "chr1:18182612-18183112"  
## [17] "chr1:41140409-41140909"   "chr1:55565312-55565812"  
## [19] "chr1:58439231-58439731"   "chr1:59571095-59571595"  
## [21] "chr1:59571722-59572222"   "chr1:59572264-59572764"  
## [23] "chr1:60294996-60295496"   "chr1:60833047-60833547"  
## [25] "chr1:60833605-60834105"   "chr1:62148070-62148570"  
## [27] "chr1:62151940-62152440"   "chr1:62228384-62228884"  
## [29] "chr1:18486708-18487208"   "chr1:28968967-28969467"  
## [31] "chr1:32820427-32820927"   "chr1:33720975-33721475"  
## [33] "chr1:37770442-37770942"   "chr1:38194127-38194627"  
## [35] "chr1:40914397-40914897"   "chr1:59093481-59093981"  
## [37] "chr1:61001656-61002156"   "chr1:62439718-62440218"  
## [39] "chr1:62944073-62944573"   "chr1:66229198-66229698"  
## [41] "chr1:66604925-66605425"   "chr1:66721616-66722116"  
## [43] "chr1:22217012-22217512"   "chr1:23930253-23930753"  
## [45] "chr1:34840533-34841033"   "chr1:40577169-40577669"  
## [47] "chr1:60828755-60829255"   "chr1:66153432-66153932"  
## [49] "chr1:108003704-108004204" "chr1:113490966-113491466"
## [51] "chr1:121328447-121328947" "chr1:134365379-134365879"
## [53] "chr1:150056492-150056992" "chr1:187594701-187595201"
## [55] "chr1:187605442-187605942" "chr1:192613817-192614317"
## [57] "chr1:195832008-195832508" "chr1:19919057-19919557"  
## [59] "chr1:21140521-21141021"   "chr1:32056370-32056870"  
## [61] "chr1:34837701-34838201"   "chr1:35314084-35314584"  
## [63] "chr1:59703058-59703558"   "chr1:65746944-65747444"  
## [65] "chr1:67426762-67427262"   "chr1:70738724-70739224"  
## [67] "chr1:70833615-70834115"   "chr1:70839986-70840486"  
## [69] "chr1:70876229-70876729"   "chr1:94535883-94536383"  
## [71] "chr1:104036944-104037444" "chr1:111968338-111968838"
## [73] "chr1:13367443-13367943"   "chr1:14560867-14561367"  
## [75] "chr1:18510289-18510789"   "chr14:31568379-31568879" 
## [77] "chr1:23073288-23073788"   "chr10:14012106-14012606" 
## [79] "chr10:92085839-92086339"  "chr19:1246202-1246702"   
## [81] "chr2:162974750-162975250" "chr3:35715890-35716390"  
## [83] "chr3:79108027-79108527"   "chr7:67721209-67721709"  
## [85] "chr7:78325753-78326253"

## Adding Annotations..

## Preparing Main Heatmap..

## 'magick' package is suggested to install to give better rasterization.
## 
## Set `ht_opt$message = FALSE` to turn off this message.

## ArchR logging successful to : ArchRLogs/ArchR-plotMarkerHeatmap-2608142ca21-Date-2021-11-12_Time-15-00-04.log

draw(heatmapPeaks,
     heatmap_legend_side = "bot",
     annotation_legend_side = "bot")

26.6 Motif enrichment for marker features

# add motif presence
proj <- addMotifAnnotations(ArchRProj = proj,
                            motifSet = "cisbp",
                            name = "Motif",
                            species = "homo sapiens",
                            force = TRUE)

# motif enrichment in marker peaks
enrichMotifs <- peakAnnoEnrichment(
  seMarker = markersPeaks,
  ArchRProj = proj,
  peakAnnotation = "Motif",
  cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
)

## ArchR logging to : ArchRLogs/ArchR-peakAnnoEnrichment-3480110840a9-Date-2021-07-26_Time-21-03-34.log
## If there is an issue, please report to github with logFile!

## 2021-07-26 21:03:43 : Computing Enrichments 1 of 10, 0.152 mins elapsed.

## 2021-07-26 21:03:45 : Computing Enrichments 2 of 10, 0.184 mins elapsed.

## 2021-07-26 21:03:47 : Computing Enrichments 3 of 10, 0.208 mins elapsed.

## 2021-07-26 21:03:48 : Computing Enrichments 4 of 10, 0.239 mins elapsed.

## 2021-07-26 21:04:14 : Computing Enrichments 5 of 10, 0.672 mins elapsed.

## 2021-07-26 21:04:27 : Computing Enrichments 6 of 10, 0.877 mins elapsed.

## 2021-07-26 21:04:32 : Computing Enrichments 7 of 10, 0.971 mins elapsed.

## 2021-07-26 21:04:46 : Computing Enrichments 8 of 10, 1.197 mins elapsed.

## 2021-07-26 21:05:19 : Computing Enrichments 9 of 10, 1.755 mins elapsed.

## 2021-07-26 21:06:07 : Computing Enrichments 10 of 10, 2.556 mins elapsed.

## ArchR logging successful to : ArchRLogs/ArchR-peakAnnoEnrichment-3480110840a9-Date-2021-07-26_Time-21-03-34.log

enrichMotifs

## class: SummarizedExperiment 
## dim: 870 10 
## metadata(0):
## assays(10): mlog10Padj mlog10p ... CompareFrequency feature
## rownames(870): TFAP2B_1 TFAP2D_2 ... TBX18_869 TBX22_870
## rowData names(0):
## colnames(10): End vRG ... IP oRG
## colData names(0):

26.7 Visualize motif enrichment for marker features

# visu
heatmapEM <-
  plotEnrichHeatmap(enrichMotifs, n = 7, transpose = TRUE)

## ArchR logging to : ArchRLogs/ArchR-plotEnrichHeatmap-260819bc8eef-Date-2021-11-12_Time-15-00-15.log
## If there is an issue, please report to github with logFile!

## Adding Annotations..

## Preparing Main Heatmap..

## 'magick' package is suggested to install to give better rasterization.
## 
## Set `ht_opt$message = FALSE` to turn off this message.

ComplexHeatmap::draw(heatmapEM,
                     heatmap_legend_side = "bot",
                     annotation_legend_side = "bot")

26.8 Save arrow project

saveArchRProject(ArchRProj = proj)