bsseq - access for methylation data
2
0
Entering edit mode
7.8 years ago
bsmith030465 ▴ 240

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!

bsseq methylation • 2.1k views
ADD COMMENT
0
Entering edit mode
7.8 years ago
hoseintoossi ▴ 40

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.

ADD COMMENT
0
Entering edit mode
7.8 years ago
bsmith030465 ▴ 240

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?

ADD COMMENT

Login before adding your answer.

Traffic: 2485 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6