<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.sstatic.net/ZTomm.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ZTomm.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2020-07-14T18:52:06.867Z
thank you veryyyyy much..i did look into couple of papers but it didn't occur to me how they did it..i will try this and update you..
kcm
CC BY-SA 4.0
2020-07-14T18:55:05.897Z
"This is the CG, and you can decide how you want to bin your regions:" very naive  question the sequence also contain the "AT " how do you define CG?
kcm
CC BY-SA 4.0
2020-07-14T21:35:05.563Z
sorry I typed too quickly, so in the code above, the result from dinucleotideFrequency() contains the proportion of dinucleotides that are CG, which is what you need, if you are interested in CpG?
StupidWolf
CC BY-SA 4.0
2020-07-14T21:47:20.697Z
yes i need CG for example this "TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC"  "the result from dinucleotideFrequency() contains the proportion of dinucleotides that are CG" this sequence which i have put also contains AT my question is is the function "dinucleotideFrequency" only look for CG only??
kcm
CC BY-SA 4.0
2020-07-14T22:01:47.393Z
of course not.. see when i wrote this ```freq = dinucleotideFrequency(SEQ,as.prob=TRUE)```, freq is a matrix. the next line i called out ```freq[,"CG"]```, did you look at the matrix?
StupidWolf
CC BY-SA 4.0
2020-07-14T22:05:58.587Z
" freq[,"CG"]" im extremely sorry not sure how did i miss that!!!!  now i got it and big thank to make it clear to me...
kcm
CC BY-SA 4.0
2020-07-14T10:11:33.260Z
2020-07-14T10:11:33.260Z
StupidWolf
<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>

CC BY-SA 4.0
2020-07-13T18:23:27.487Z
you can either calculate the GC content? And you can also check overlap with annotated cpg island from UCSC... Not very sure what is the question
StupidWolf
CC BY-SA 4.0
2020-07-13T20:54:22.277Z
can i get the GC content from  atac seq data ? if yes..how that is my question I mean what input do i require
kcm
CC BY-SA 4.0
2020-07-13T21:11:44.930Z
did i get you wrong? It seems like you have annotated peaks. It's a matter of getting the sequences in the annotated peaks and calculating the number of CGs (sorry I typed too quickly previously, it should be CGs instead of GC content)
StupidWolf
CC BY-SA 4.0
2020-07-14T06:37:47.323Z
"It seems like you have annotated peaks." yes correct. "It's a matter of getting the sequences in the annotated peaks and calculating the number of CGs" okay that looks straight forward though can it be done in R ? can you suggest me some library
kcm
<atac-seq>
Chromatin accessibility level on promoters with different CpG ratio calculation
# Chromatin accessibility level on promoters with different CpG ratio calculation

(Asked by: username_0 on 2020-07-11)

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?

**Tags:** atac-seq
> **Comments:**
> **username_1**: you can either calculate the GC content? And you can also check overlap with annotated cpg island from UCSC... Not very sure what is the question  
> **username_0**: can i get the GC content from  atac seq data ? if yes..how that is my question I mean what input do i require  
> **username_1**: did i get you wrong? It seems like you have annotated peaks. It's a matter of getting the sequences in the annotated peaks and calculating the number of CGs (sorry I typed too quickly previously, it should be CGs instead of GC content)  
> **username_0**: "It seems like you have annotated peaks." yes correct. "It's a matter of getting the sequences in the annotated peaks and calculating the number of CGs" okay that looks straight forward though can it be done in R ? can you suggest me some library  

---

## 1 answers

### Answer by username_1 on 2020-07-14. Score: 2 | Accepted answer

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 |
                 |
  [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
                                                
  [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)

> **Comments:**
> **username_0**: thank you veryyyyy much..i did look into couple of papers but it didn't occur to me how they did it..i will try this and update you..  
> **username_0**: "This is the CG, and you can decide how you want to bin your regions:" very naive  question the sequence also contain the "AT " how do you define CG?  
> **username_1**: sorry I typed too quickly, so in the code above, the result from dinucleotideFrequency() contains the proportion of dinucleotides that are CG, which is what you need, if you are interested in CpG?  
> **username_0**: yes i need CG for example this "TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC"  "the result from dinucleotideFrequency() contains the proportion of dinucleotides that are CG" this sequence which i have put also contains AT my question is is the function "dinucleotideFrequency" only look for CG only??  
> **username_1**: of course not.. see when i wrote this ```freq = dinucleotideFrequency(SEQ,as.prob=TRUE)```, freq is a matrix. the next line i called out ```freq[,"CG"]```, did you look at the matrix?  
> **username_0**: " freq[,"CG"]" im extremely sorry not sure how did i miss that!!!!  now i got it and big thank to make it clear to me...  

---

