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.logmarkersPeaks## 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 oRG26.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.logdraw(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.logenrichMotifs## 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)