methylkit from Bismark coverage output

189 views
Skip to first unread message

rthapa

unread,
Nov 2, 2019, 9:24:22 PM11/2/19
to methylkit_discussion
Hi,

I am new to methylkit. I want to estimate the coverage of methylation in the gene-body region. I have coverage output from bismark. I also have bed file of gene region. I am wondering if I could use the bismark coverage output directly or do I need to reformat the data. I would appreciate any suggestion.

Thanks
Screen Shot 2019-11-02 at 8.23.04 PM.png
Screen Shot 2019-11-02 at 8.22.40 PM.png

Alexander Blume

unread,
Nov 4, 2019, 1:57:39 AM11/4/19
to methylkit_discussion
Dear rthapa,

You can use the coverage files directly, please use the methRead() function with "bismarkCoverage" for the pipeline argument. 
Check the help page of methRead for more information. 

Best,
Alex

rthapa

unread,
Nov 4, 2019, 11:56:19 AM11/4/19
to methylkit_discussion
Thank you for your suggestion. I tried to use methRead but got the following error.

myobj=methRead(file.list, sample.id=list("CHG_context_SRR5748817_1"), dbtype = NA, pipeline = "bismarkCoverage")
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘methRead’ for signature ‘"list", "list", "missing"’


I am not sure if that is because of the file format that I used. The bismark coverage output I have has chr, start, end, percentage coverage, number of cytosines (methylated bases) and number of cytosines (unmethylated bases) unlike the one with chr, start, end, number of cytosines (methylated bases) and number of thymines (unmethylated bases) that is used by methRead.

Alexander Blume

unread,
Nov 4, 2019, 2:14:08 PM11/4/19
to methylkit_discussion
Hi,

I guess you missed to fill the assembly argument, sorry that the error is not very informative.

Best, Alex

Ranjita Thapa

unread,
Nov 4, 2019, 2:51:40 PM11/4/19
to methylkit_...@googlegroups.com
I do not have assembly file so I didn't fill the assembly argument. I only have coverage output from bismark and the gene region file where I want to find the coverage. What do you suggest?

myobj=methRead(file.list, sample.id=list("CHG_context_SRR5748817_1"), dbtype = NA, pipeline = "bismarkCoverage")
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘methRead’ for signature ‘"list", "list", "missing"’
--
You received this message because you are subscribed to the Google Groups "methylkit_discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to methylkit_discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/methylkit_discussion/7d5b9d8c-1e99-4bbe-a445-2142efeed8e4%40googlegroups.com.

Alexander Blume

unread,
Nov 4, 2019, 3:14:57 PM11/4/19
to methylkit_...@googlegroups.com
The assembly is only a string, you can put any content you want. 


Best, 
Alex 

Sent from mobile.

You received this message because you are subscribed to a topic in the Google Groups "methylkit_discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/methylkit_discussion/9nPkCbPyKLI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to methylkit_discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/methylkit_discussion/CAPdM%3DRuNJs4AOYkE78gLqNvD7_%3Db6oBVPyqLGDGp0FS9zEeqyQ%40mail.gmail.com.

Ranjita Thapa

unread,
Nov 4, 2019, 5:06:12 PM11/4/19
to methylkit_...@googlegroups.com
I used this script: myobj=methRead(file.list, sample.id=list("CHG_context_SRR5748817_1"), assembly="sb", dbtype = NA, treatment = 1, pipeline = "bismarkCoverage")

and got the output like this;
[[1]]
methylRaw object with 20521740 rows
--------------
    chr start end strand coverage numCs numTs
1 Chr05   419 419      *       15    15     0
2 Chr05   438 438      *       15    15     0
3 Chr05   440 440      *       10    10     0
4 Chr05   441 441      *       15    15     0
5 Chr05   443 443      *       10     9     1
6 Chr05   597 597      *       21    20     1


whereas my Bismark coverage output is different than this;

Chr05 241 241 100 9 0

Chr05 258 258 100 1 0

Chr05 260 260 88.8888888888889 8 1

Chr05 261 261 100 2 0

Chr05 263 263 88.8888888888889 8 1

Chr05 419 419 100 15 0

Chr05 421 421 100 7 0

Chr05 438 438 100 15 0

Chr05 440 440 100 10 0

Chr05 441 441 100 15 0

Chr05 443 443 90 9 1


Also, it looks like some of the rows are escaped during the processing. Is it the correct way to process? 

Alexander Blume

unread,
Nov 5, 2019, 2:42:16 AM11/5/19
to methylkit_...@googlegroups.com
Hi, 

There is a default minimum coverage threshold of 10 bases set with the  mincov argument, so every site with less reads are skipped. 

Best, 
Alex

Sent from mobile.

Reply all
Reply to author
Forward
0 new messages