# Chromatin accessibility level on promoters with different CpG ratio calculation

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?

• 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 Jul 13 '20 at 18:23
• 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
Jul 13 '20 at 20:54
• 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) Jul 13 '20 at 21:11
• "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
Jul 14 '20 at 6:37

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)

extraCols <- c(signalValue = "numeric", pValue = "numeric",
qValue = "numeric", peak = "integer")
gr <- import(peaks_narrowbed, format = "BED",extraCols = extraCols)

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>


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)

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)


• 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
Jul 14 '20 at 18:52
• "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
Jul 14 '20 at 18:55
• 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? Jul 14 '20 at 21:35
• 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? Jul 14 '20 at 22:01
• " 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
Jul 14 '20 at 22:05