BSMAP and methyratio.py input but incorrect DM% outputed by methylKit

[Claire Morgan]

Hello,

Im having a problem with methylKit as it is outputting the incorrect % differences between sample comparisons P3.CPG.txt and P7.CpG.txt.  I must be making a mistake somewhere but it's not obvious to me. Any help you can provide me would be greatly appreciated. I have provided the methylkit code below.

I aligned my bisulfite reads using BSMAP and used the methyratio.py program to call the C and CT values.

The methyKit output for position chr13.119601308 between two comparisons is 92% and reads as follows:

"id" "chr" "start" "end" "strand" "pvalue" "qvalue" "meth.diff"
"3256" "chr13.119601308" "chr13" 119601308 119601308 "+" 5.60865078296767e-06 0.00563332884641272 92.8571428571429

However when i extract the information from the input files P3.CpG.txt and P7.CpG.txt respectively these are the scores
[header] chr    pos    strand    context    ratio    total_C    methy_C    CI_lower    CI_upper
[P3.CpG.txt] chr13    119601308    +    CG    0.071    14    1    0.013    0.315
[P7.CpG.txt] chr13    119601308    +    CG    0.240    1247    299    0.217    0.264

So the probem is that the difference in ratios between P3.CpG.txt and P7.CpG.txt is -0.16 not 92.85.

Can anyone tell me if this is an error on my end, have i called the input columns to methylKit wrong or is this perhaps an error in the methylKit code?

Looking forward to your response.

Claire





####################################################################################################
methylKit commands
#####################################################################################################
library(methylKit)
file.list=list( "P3.CpG.txt","P7.CpG.txt")

myobj=read(file.list,sample.id=list("ctrl","test"),assembly="mm10",treatment=c(0,1),context="CpG", header=TRUE, resolution="base",pipeline=list(fraction=TRUE,chr.col=1,start.col=2,end.col=2,coverage.col=6,strand.col=3,freqC.col=5 ))

filtered.myobj=filterByCoverage(myobj,lo.count=5,lo.perc=NULL, hi.count=NULL,hi.perc=99.9)

meth=unite(filtered.myobj,destrand=TRUE)

myDiff=calculateDiffMeth(meth)

myDiff33p=get.methylDiff(myDiff,difference=20,qvalue=0.05)

write.table(getData(myDiff33p),file="myDiff20p.txt")

[Dhruva Chandramohan]

Check the methylBase objects myobj and filtered.myobj to see whether there was a problem during the reading process; you can check the first few rows (just type their names and hit return, or use head), or if you want to examine the specific loci you noticed, use either select or getData.

--
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 post to this group, send email to methylkit_...@googlegroups.com.
Visit this group at http://groups.google.com/group/methylkit_discussion.
For more options, visit https://groups.google.com/d/optout.

[Claire Morgan]

Hello,

Thank you Dhurva. it would appear i did not correctly tab delimit my file following extraction of only CpGs from the methratio output which cause the errors.

Problem solved!

To unsubscribe from this group and stop receiving emails from it, send an email to methylkit_discussion+unsub...@googlegroups.com.

[ABDELJALIL SENHAJI RACHIK]

Hello,

Please can you share with me the commands you used to separate the 3 contexts ?

thank you very much

best, 

Alex