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, thankyou!
#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.