I was trying to use the bsseq package to analyze methylation. I was able to make a bsseq object but cannot seem to get the methylation values from the object. For example, I create an object (tdmr) with 4 samples:
tdmr
An object of type 'BSseq' with 27027 methylation loci 4 samples has not been smoothed
The object seems to contain the data for the 4 samples, but when I convert this into a GRanges object I can't see the data:
gr <- granges(tdmr)
gr GRanges object with 27027 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chrX [3045200, 3045200] *
[2] chrX [3045201, 3045201] *
[3] chrX [3045249, 3045249] *
[4] chrX [3045250, 3045250] *
[5] chrX [3045279, 3045279] *
... ... ... ...
[27023] chrX [155614920, 155614920] *
[27024] chrX [155615020, 155615020] *
[27025] chrX [155615021, 155615021] *
[27026] chrX [155615169, 155615169] *
[27027] chrX [155615170, 155615170] *
seqinfo: 450 sequences from an unspecified genome; no seqlengths
Have I converted the bsseq object incorrectly? Or do I need to do something else?
thanks!
2 answers
The GRanges object is shared between the samples.
getCoverage(tdmr,type="M")
will give the methylated read count and
getCoverage(tdmr)
will give you overall coverage for each of these samples over these ranges.
sampleNames(tdmr)
gives the sample names.
Thanks, yes I realize I can get it through the getCoverage functions, but should it be in the granges object too (as data)?
...sorry, new to GRanges objects!
Also, when I try to intersect this with genes/exons, I get nothing:
=================
$ library(TxDb.Hsapiens.UCSC.hg38.knownGene)
$ tx <- TxDb.Hsapiens.UCSC.hg38.knownGene
$ genegr <- genes(tx)
$ promUp <- promoters(tx,upstream = 2000)
$ xx <- intersect(granges(tdmr),genegr)
$ xx
GRanges object with 0 ranges and 0 metadata columns: seqnames ranges strand <rle> <iranges> <rle>
seqinfo: 455 sequences (1 circular) from hg38 genome
How should I be doing this to get overlaps between the methylation loci and genes?
Log in to answer this question.
