text: user285: I want to categorize promoters based upon high, intermediate and low-CpG content (high-CpG-density promoters (HCPs), intermediate-CpG-density promoters(ICPs), and low-CpG-density promoters (LCPs)).
  So the data I have is for promoter which i have annotated 1000+/- around TSS and taken them as something less than 1kb is promoter region and beyond distal.
  In term of tool i have used chipseeker to annotate and do the above step.
  So now if i have to find the CpG density as i have mentioned above how do i get that?
  
  user5077: Below I import the peaks call in narrowbed format from macs2 into a GRanges object. If you have a bed file, you can just use rtracklayer for it:
  library(rtracklayer)
  library(BSgenome.Hsapiens.UCSC.hg38)
  
  peaks_narrowbed = "https://www.encodeproject.org/files/ENCFF846JAO/@@download/ENCFF846JAO.bed.gz"
  
  extraCols <- c(signalValue = "numeric", pValue = "numeric",
                            qValue = "numeric", peak = "integer")
  gr <- import(peaks_narrowbed, format = "BED",extraCols = extraCols)
  
  head(gr)
  
  GRanges object with 6 ranges and 6 metadata columns:
        seqnames        ranges strand |
           <Rle>     <IRanges>  <Rle> |
    [1]     chr1   10049-10326      * |
    [2]     chr1 180686-181040      * |
    [3]     chr1 186648-186903      * |
    [4]     chr1 191293-191604      * |
    [5]     chr1 267872-268146      * |
    [6]     chr1 586110-586324      * |
                                                     name     score signalValue
                                              <character> <numeric>   <numeric>
    [1] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_1       135     7.10907
    [2] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_2       309    11.65888
    [3] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_3        81     5.40289
    [4] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_4        56     4.54981
    [5] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_5       145     7.39343
    [6] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_6       116     6.54035
  
  Then we extract the sequences that are within these regions, you can of course redefine the region by extending it from the peak:
  library(Biostrings)
  SEQ = getSeq(Hsapiens,gr)
  head(SEQ)
  
  A DNAStringSet instance of length 6
      width seq
  [1]   278 TAACCCTAACCCTAACCCTAACCCTAACCCTAAC...CCCCAACCCCAACCCCAACCCCAACCCCAACCC
  [2]   355 CTAACCCCTAACCCCTAACCCTAACCCTACCCTA...CTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGA
  [3]   256 GCAACCACCTGAGCGCGGGCATCCTGTGTGCAGA...CCCCTTCTTTCCATTGGTTTAATTAGGAACGGG
  [4]   312 CAATCAGCAGGGACCGTGCACTCTCTTGGAGCCA...CCCTGCCCCTGTCTCTTCCGTGCAGGAGGAGCA
  [5]   275 TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC...TTGTTTAGTTTGCGTTGTGTTTCTCCAACTTTG
  [6]   215 GAGGCTGAAGGAGACTGATGTGGTTTCTCCTCAG...CCAACTTTGTGCCTCATCAGGAAAAGCTTTGGA
  
  This is the CG, and you can decide how you want to bin your regions:
  freq = dinucleotideFrequency(SEQ,as.prob=TRUE)
  hist(freq[,"CG"],br=100)
  
  
  
original_text: 285: <p>I want to categorize promoters based upon high, intermediate and low-CpG content (high-CpG-density promoters (HCPs), intermediate-CpG-density promoters(ICPs), and low-CpG-density promoters (LCPs)).</p>
  <p>So the data I have is for promoter which i have annotated 1000+/- around TSS and taken them as something less than 1kb is promoter region and beyond distal.</p>
  <p>In term of tool i have used chipseeker to annotate and do the above step.</p>
  <p>So now if i have to find the CpG density as i have mentioned above how do i get that?</p>
   5077: <p>Below I import the peaks call in narrowbed format from macs2 into a GRanges object. If you have a bed file, you can just use rtracklayer for it:</p>
  <pre><code>library(rtracklayer)
  library(BSgenome.Hsapiens.UCSC.hg38)
  
  peaks_narrowbed = &quot;https://www.encodeproject.org/files/ENCFF846JAO/@@download/ENCFF846JAO.bed.gz&quot;
  
  extraCols &lt;- c(signalValue = &quot;numeric&quot;, pValue = &quot;numeric&quot;,
                            qValue = &quot;numeric&quot;, peak = &quot;integer&quot;)
  gr &lt;- import(peaks_narrowbed, format = &quot;BED&quot;,extraCols = extraCols)
  
  head(gr)
  
  GRanges object with 6 ranges and 6 metadata columns:
        seqnames        ranges strand |
           &lt;Rle&gt;     &lt;IRanges&gt;  &lt;Rle&gt; |
    [1]     chr1   10049-10326      * |
    [2]     chr1 180686-181040      * |
    [3]     chr1 186648-186903      * |
    [4]     chr1 191293-191604      * |
    [5]     chr1 267872-268146      * |
    [6]     chr1 586110-586324      * |
                                                     name     score signalValue
                                              &lt;character&gt; &lt;numeric&gt;   &lt;numeric&gt;
    [1] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_1       135     7.10907
    [2] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_2       309    11.65888
    [3] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_3        81     5.40289
    [4] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_4        56     4.54981
    [5] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_5       145     7.39343
    [6] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_6       116     6.54035
  </code></pre>
  <p>Then we extract the sequences that are within these regions, you can of course redefine the region by extending it from the peak:</p>
  <pre><code>library(Biostrings)
  SEQ = getSeq(Hsapiens,gr)
  head(SEQ)
  
  A DNAStringSet instance of length 6
      width seq
  [1]   278 TAACCCTAACCCTAACCCTAACCCTAACCCTAAC...CCCCAACCCCAACCCCAACCCCAACCCCAACCC
  [2]   355 CTAACCCCTAACCCCTAACCCTAACCCTACCCTA...CTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGA
  [3]   256 GCAACCACCTGAGCGCGGGCATCCTGTGTGCAGA...CCCCTTCTTTCCATTGGTTTAATTAGGAACGGG
  [4]   312 CAATCAGCAGGGACCGTGCACTCTCTTGGAGCCA...CCCTGCCCCTGTCTCTTCCGTGCAGGAGGAGCA
  [5]   275 TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC...TTGTTTAGTTTGCGTTGTGTTTCTCCAACTTTG
  [6]   215 GAGGCTGAAGGAGACTGATGTGGTTTCTCCTCAG...CCAACTTTGTGCCTCATCAGGAAAAGCTTTGGA
  </code></pre>
  <p>This is the CG, and you can decide how you want to bin your regions:</p>
  <pre><code>freq = dinucleotideFrequency(SEQ,as.prob=TRUE)
  hist(freq[,&quot;CG&quot;],br=100)
  </code></pre>
  <p><a href="https://i.stack.imgur.com/ZTomm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ZTomm.png" alt="enter image description here" /></a></p>
  