why 150k peaks

[AnjaliC4]

Hi, I have a very simple question: What is the rationale for an upper limit of 150k on peak calling? Can this be related to %genome in the open chromatin?

Thanks

[rcorces]

This is our intuition for a reasonable upper bound on true peaks detectable in a single cell type. This comes from analyzing bulk ATAC-seq data on sorted cells and trying to define a cutoff for what a true peak is.

[AnjaliC4]

Thanks Ryan, From what I read in the literature - the percentage open chromatin can vary between 2-5%. 150k would be about 2.5% so it seems reasonable but could miss some valid OCRs. Also, could you please point me in the direction of the the paper you are referring to?

[rcorces]

No paper reference, just our personal experience after analyzing a lot of bulk ATAC-seq.

Of course some cell types may be different or outliers but we had to draw a cutoff somewhere and this is our intuition for the best default. As with most parameters, we let you change the default if your intuition differs from ours.

[AnjaliC4]

Hi @rcorces ,

Since I have already called peaks (and do not wish to re-call them), I am just trying to get more peaks by changing "maxPeaks" in the addReproduciblePeaks(), but getting the error.

Error:
<simpleError in readRDS(summitFiles[x]): bad 'file' argument>

Code:
temp = list.files(path="/home//scratch/ArchR_outputs/ArchR_BatchTSS_MajorClusters/PeakCalls/ReplicateCalls/", pattern="*.rds",full.names = TRUE)
coverageMetadata <- ArchR:::.getCoverageMetadata(ArchRProj = proj, groupBy = "ClustersMapped", minCells = 0)
outSummits <- unlist(lapply(temp, function(x){readRDS(x)}))
outSummitList <- split(outSummmits, coverageMetadata$Group)
summitNamesList <- split(coverageMetadata$Name, coverageMetadata$Group)

groupPeaks <- .safelapply(seq_along(outSummitList), function(i){
prefix <- sprintf("Rep. Peaks Group (%s of %s) :", i, length(outSummitList))
.logDiffTime(sprintf("%s Creating Reproducible Peaks", prefix), tstart, verbose = FALSE, logFile = logFile)
peaks <- suppressMessages(.identifyReproduciblePeaks(
summitFiles = outSummitList[[i]],
summitNames = summitNamesList[[i]],
reproducibility = reproducibility,
extendSummits = extendSummits,
blacklist = genomeAnnotation$blacklist,
prefix = prefix,
logFile = logFile
))
.logDiffTime(sprintf("%s Annotating and Filtering Peaks", prefix), tstart, verbose = FALSE, logFile = logFile)
peaks <- sort(sortSeqlevels(peaks))
peaks <- subsetByOverlaps(peaks, genomeAnnotation$chromSizes, type = "within")
peaks <- .fastAnnoPeaks(peaks, BSgenome = BSgenome, geneAnnotation = geneAnnotation, promoterRegion = promoterRegion, logFile = logFile)
peaks <- peaks[which(mcols(peaks)$N < 0.001)] #Remove N Containing Peaks
peaks <- peaks[order(peaks$groupScoreQuantile, decreasing = TRUE)]
peaks <- head(peaks, groupSummary[names(outSummitList)[i],"maxPeaks"])
mcols(peaks)$N <- NULL #Remove N Column
print(file.path(outDir, paste0(make.names(names(outSummitList)[i]), "-reproduciblePeaks.gr.rds")))
saveRDS(peaks, file.path(outDir, paste0(make.names(names(outSummitList)[i]), "-reproduciblePeaks.gr.rds")))
return(peaks)
}, threads = threads) %>% GRangesList()
names(groupPeaks) <- names(outSummitList)

I think the outSummitList argument I am passing isnt correct. It looks like this for me (in the manner i created it). If you can help me point what the issue is or show me how outSummitList should like - that would be really great.

$Mic
$Mic$Mic._.4
GRanges object with 250708 ranges and 1 metadata column:
seqnames ranges strand | score
|
[1] chr1 10105 * | 28.27672
[2] chr1 180888 * | 9.88767
[3] chr1 629280 * | 55.96005
[4] chr1 629955 * | 518.61304
[5] chr1 630312 * | 697.62787
... ... ... ... . ...
[250704] chrY 56877765 * | 1.59032
[250705] chrY 56879025 * | 2.90198
[250706] chrY 56880265 * | 1.04491
[250707] chrY 56881099 * | 2.90198
[250708] chrY 56883901 * | 3.64443

seqinfo: 24 sequences from an unspecified genome; no seqlengths

$Mic$Mic._.1
GRanges object with 181543 ranges and 1 metadata column:
seqnames ranges strand | score
|
[1] chr1 10086 * | 10.21162
[2] chr1 268007 * | 8.15063
[3] chr1 629958 * | 506.43799
[4] chr1 634034 * | 3574.29102
[5] chr1 778671 * | 43.86810
... ... ... ... . ...
[181539] chrY 56875066 * | 5.29193
[181540] chrY 56876283 * | 4.41304
[181541] chrY 56878364 * | 2.08017
[181542] chrY 56879849 * | 4.41304
[181543] chrY 56880177 * | 1.44153

[AnjaliC4]

UPDATE: I think I figured this one out. If you could please review the code (just to make sure it makes sense with ArchR format) - I will very much appreciate it. I have pasted it here - may also help anyone else who wishes to do this analysis.
Basically, I think this must be the format of outSummitList:
$Mic
[1] "Mic..1-summits.rds" "Mic..2-summits.rds" "Mic..3-summits.rds"
[4] "Mic..4-summits.rds"
$Oli
[1] "Oli..1-summits.rds" "Oli..2-summits.rds" "Oli..3-summits.rds"
[4] "Oli..4-summits.rds"

`library(ArchR)
proj <- loadArchRProject("/home/scratch/ArchR_outputs/ArchR_Batch")
genomeAnnotation = getGenomeAnnotation(proj)
geneAnnotation = getGeneAnnotation(proj)
ArchR:::.requirePackage(genomeAnnotation$genome)
ArchR:::.requirePackage("Biostrings",source="bioc")
BSgenome <- eval(parse(text = genomeAnnotation$genome))
BSgenome <- validBSgenome(BSgenome)
promoterRegion = c(2000, 100)

outDir ="~/scratch/"
#######

outSummit <- list.files(path="/home/scratch/ArchR_outputs/ArchR_Batch/PeakCalls/ReplicateCalls/", pattern="*.rds")
outSummit <- do.call("rbind",strsplit(outSummitList, "-"))
coverageMetadata <- ArchR:::.getCoverageMetadata(ArchRProj = proj, groupBy = "ClustersMapped", minCells = 0)
coverageMetadata <- coverageMetadata[match(outSummit[,1], coverageMetadata$Name),]
summitNamesList <- split(coverageMetadata$Name, coverageMetadata$Group)
outSummitList <- list("Ast" = c(outSummitList[1],outSummitList[2],outSummitList[3],outSummitList[4]), "End" =c(outSummitList[5],outSummitList[6],outSummitList[7],outSummitList[8]), "ExN" = c(outSummitList[9],outSummitList[10],outSummitList[11],outSummitList[12]),"InN" = c(outSummitList[13],outSummitList[14],outSummitList[15],outSummitList[16]),"Mic" = c(outSummitList[17],outSummitList[18],outSummitList[19],outSummitList[20]),"Oli" = c(outSummitList[21],outSummitList[22],outSummitList[23],outSummitList[24]), "OPC" =c(outSummitList[25],outSummitList[26],outSummitList[27],outSummitList[28]))
logFile = createLogFile("addReproduciblePeakSet")
logFile="~/projects/ArchRLogs/ArchR-addReproduciblePeakSet-c6059aa3cf5-Date-2022-03-16_Time-14-07-19.log"
groupPeaks <- ArchR:::.safelapply(seq_along(outSummitList), function(i){
prefix <- sprintf("Rep. Peaks Group (%s of %s) :", i, length(outSummitList))
ArchR:::.logDiffTime(sprintf("%s Creating Reproducible Peaks", prefix), tstart, verbose = FALSE, logFile = ".")
peaks <- ArchR:::.identifyReproduciblePeaks(
summitFiles = outSummitList[[i]],
summitNames = summitNamesList[[i]],
reproducibility = "2",
extendSummits = 250,
blacklist = genomeAnnotation$blacklist,
prefix = prefix,
logFile = logFile
)
#print(peaks)
ArchR:::.logDiffTime(sprintf("%s Annotating and Filtering Peaks", prefix), tstart, verbose = FALSE, logFile = logFile)
peaks <- sort(sortSeqlevels(peaks))
peaks <- subsetByOverlaps(peaks, genomeAnnotation$chromSizes, type = "within")
peaks <- ArchR:::.fastAnnoPeaks(peaks, BSgenome = BSgenome, geneAnnotation = geneAnnotation, promoterRegion = promoterRegion, logFile = logFile)
peaks <- peaks[which(mcols(peaks)$N < 0.001)] #Remove N Containing Peaks
peaks <- peaks[order(peaks$groupScoreQuantile, decreasing = TRUE)]
#peaks <- head(peaks, groupSummary[names(outSummitList)[i],"maxPeaks"])
peaks <- head(peaks, max(300000, length(peaks)))
mcols(peaks)$N <- NULL #Remove N Column
print(peaks)
print(file.path(outDir, paste0(make.names(names(outSummitList)[i]), "-reproduciblePeaks.gr.rds")))
#saveRDS(peaks, file.path(outDir, paste0(make.names(names(outSummitList)[i]), "-reproduciblePeaks.gr.rds")))
return(peaks)
}, threads = 10) %>% GRangesList()

`