Hi everyone,

I am new to Chip-seq, I have a control sample and 2 treatment samples which are in .bedGraph. I want to see the binding at enhancer and super-enhancer region of MYCN gene of human in treatment and control sample . How can I achieve that, I am not sure or where to start either. I want to perform this analysis in R.

Following is what I have done so far. I have gotten a reference genome of human hg38 MYCN gene coordinates tried to plot how my sample tracks look compared to the ref genome for which I get an error which I have pasted at the end.

So, I have 3 questions.

1. How to over come the error when making plottracks?
2. I have .bedGraph file and want to see if they have enhancer and super enhancer regions annotated?
3. I want to the binding for my samples at the enhancer and super enhancer region of MYCN gene?

Any help would greatly be appreciated, thank you!

#loading the libraries
library(Gviz)
library(GenomicRanges)
library(rtracklayer)
library(IRanges)
library(biomaRt)

#laoding the bed.graph files for analysis 
S1C = import.bedGraph('/home/bedGraph_bigWig_files_for_genomeBrowser/1C.coverage.bedGraph')
S2D1 = import.bedGraph('/home/bedGraph_bigWig_files_for_genomeBrowser/2D1.coverage.bedGraph')
S3D2 = import.bedGraph('/home/bedGraph_bigWig_files_for_genomeBrowser/3D2.coverage.bedGraph')
S4C = import.bedGraph('/home/bedGraph_bigWig_files_for_genomeBrowser/4C.coverage.bedGraph')
S5D1 = import.bedGraph('/home/bedGraph_bigWig_files_for_genomeBrowser/5D1.coverage.bedGraph')
S6D2 = import.bedGraph('/home/bedGraph_bigWig_files_for_genomeBrowser/6D2.coverage.bedGraph')

seqnames(S6D2)
ranges(S6D2)
mcols(S6D2)

plot( 200000:201000, S6D2$score[200000:201000], 
      xlab="chr2", ylab="counts per bin" )

#getting the human pre-annotated model set
library(BSgenome.Hsapiens.UCSC.hg38)
genome = BSgenome.Hsapiens.UCSC.hg38
si = seqinfo(genome)
si = si[ paste0('chr', c(1:19, 'X', 'Y'))]


library(biomaRt)
# Set up the connection to the Ensembl BioMart for human (archived version)
mart = useMart(biomart = "ENSEMBL_MART_ENSEMBL", 
               dataset = "hsapiens_gene_ensembl", 
               host = "www.ensembl.org")

# Retrieve the feature map using Gviz
fm = Gviz:::.getBMFeatureMap()
# Modify the feature map for gene symbols
fm["symbol"] = "ensembl_gene_id"

# List all available attributes for the current BioMart dataset
attributes = listAttributes(mart)
head(attributes)

bm = BiomartGeneRegionTrack(chromosome='chr2', genome="hg38", 
                            start=15940566, end = 15947007,
                            biomart=mart, 
                            size=4, name="RefSeq", utr5="red3", utr3="red3", 
                            protein_coding="black", col.line=NULL, cex=7,
                            collapseTranscripts="longest",
                            featureMap=fm)

gen <- "hg38"
chr <- "chr2"
atrack <- bm

gtrack <- GenomeAxisTrack()
plotTracks(list(gtrack, atrack))
itrack <- IdeogramTrack(genome = gen, chromosome = chr)
plotTracks(list(itrack, gtrack, atrack))

plotTracks(c(itrack, bm),
           start=15940550, end = 15947004,
           transcriptAnnotation="symbol", window="auto", 
           cex.title=1, fontsize=10 )


input.S1C = DataTrack(S1C, 
                      genome="hg38", name='Peaks Rep. 1',
                      chromosome='chr2',
                      shape='box',fill='blue3',size=2)

input.S2D1 = DataTrack(S2D1, 
                       strand="*", genome="hg38", col.histogram='steelblue',
                       fill.histogram='black',chromosome = "chr2", name="S2D1", col.axis='steelblue',
                       cex.axis=0.4, ylim=c(0,150))

input.S3D2 = DataTrack(S3D2, 
                       strand="*", genome="hg38", col.histogram='sienna2',
                       fill.histogram='black',chromosome = "chr2", name="S3D2", col.axis='sienna2',
                       cex.axis=0.4, ylim=c(0,150))


input.S4C = DataTrack(S4C, 
                      strand="*", genome="hg38", col.histogram='maroon3',
                      fill.histogram='black',chromosome = "chr2", name="S4C", col.axis='maroon3',
                      cex.axis=0.4, ylim=c(0,150))


input.S5D1 = DataTrack(S5D1, 
                       strand="*", genome="hg38", col.histogram='yellowgreen',
                       fill.histogram='black',chromosome = "chr2", name="S5D1", col.axis='yellowgreen',
                       cex.axis=0.4, ylim=c(0,150))


input.S6D2 = DataTrack(S6D2, 
                       strand="*", genome="hg38", col.histogram='darksalmon',
                       fill.histogram='black',chromosome = "chr2", name="S6D2", col.axis='darksalmon',
                       cex.axis=0.4, ylim=c(0,150))

plotTracks(list(bm, input.S1C),
           start=15940550, end = 15941004,
           transcriptAnnotation="symbol", window="auto", 
           type="histogram", cex.title=0.7, fontsize=10 )


Error in .local(GdObject, ...) :    Too many stacks to draw. Either
increase the device size or limit the drawing to a smaller region.