Call peaks

Call peaks from a bedGraph or bigWig file using multiple methods. By default, it automatically infers a reasonable cutoff threshold as well.
Note : MACS3/MACSr is not currently compatible with Windows (see here for details).

call_peaks(
  bedgraph_path,
  call_peaks_method = "MACSr",
  cutoff = NULL,
  minlen = 200L,
  maxgap = 30L,
  call_summits = TRUE,
  trackline = TRUE,
  log = TRUE,
  outdir = tempdir(),
  outputfile = "MACSr.peaks.bed",
  return_path = FALSE,
  verbose = TRUE
)

Arguments

bedgraph_path

Path to bedGraph file. Can instead provide a bigWig file, but this will first be converted to bedGraph format, which can take some time if trying to convert data from across the entire genome.

call_peaks_method

Method to call peaks with:

• "MACSr" : Uses MACS3 via bdgpeakcall.

cutoff

Cutoff depends on which method you used for score track. If the file contains pvalue scores from MACS3, score 5 means pvalue 1e-5. If NULL, a reasonable cutoff value will be inferred through a cutoff_analysis.

minlen

minimum length of peak, better to set it as d value. DEFAULT: 200", default = 200.

maxgap

maximum gap between significant points in a peak, better to set it as tag size. DEFAULT: 30", default = 30.

call_summits

If set, MACS will use a more sophisticated approach to find all summits in each enriched peak region DEFAULT: False",default=False.

trackline

Tells MACS not to include trackline with bedGraph files. The trackline is required by UCSC.

log

Whether to capture logs.

outdir

Directory to store cutoff_analysis report and peak file in.

outputfile

Name of the peak output file (stored in BED format).

return_path

Whether to return the path to the saved peak file, or the peak data itself as a GRanges object.

verbose

Print messages.

Value

GRanges or path to save peaks file.

Examples

#### Get bedGraph subset ####
## Normally, you'd call peaks on the entire chromosome, 
## or even the whole genome. But for demo purposes we'll just use one locus.
gsm <- "GSM4703766" 
links <- PeakyFinders:::get_geo_links(gsm = gsm) 
#> Determining available file types.
#> Found file link(s) for 1 category.
#> bedgraph : 
#> >>> https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4703nnn/GSM4703766/suppl//GSM4703766_GAWSU510-A549-H3K4me3-ChIP-rep2_AGTCAA_GAWSU510_CKDL200157349-1a-13_HC3LFBBXX_L2_PE_R1-2.fastq.sam_noDup_RPBM.bedgraph.gz
query_granges <- GenomicRanges::GRanges("chr6:165169213-167169213")
gr <- rtracklayer::import(con = links$bedgraph, 
                          which = query_granges)
tmp <- tempfile(fileext = ".bedgraph")
rtracklayer::export.bedGraph(object = gr, con = tmp)

#### Call peaks #### 
if(.Platform$OS.type!="windows"){
peaks <- PeakyFinders::call_peaks(bedgraph_path = tmp)
}
#> File already in bedGraph format.
#> Analyzing cutoff thresholds.
#> INFO:root:Read and build bedGraph...
#> INFO:root:Analyze cutoff vs number of peaks/total length of peaks/average length of peak
#> INFO:root:Write report...
#> INFO:root:Done
#> 
#> Calling peaks.
#> INFO:root:Read and build bedGraph...
#> INFO:root:Call peaks from bedGraph...
#> INFO:root:Write peaks...
#> INFO:root:Done
#> 
#> 11 peaks called.