Tweet
I want to generate an outfile containing or C site (CG, CHG or CHH) that are differentially methylated (DM) using the R package methylKit
The problem is that I got a file that is empty.
I understand that outfiles listing hypermethylated or Hypomethylated can be empty simply because there are not DM sites between the compared samples. However I'm sure this is not the case and please let me explain before go to my question.
I have 4 samples ( control vs treatment with 2 replicates each). If I compare tiles of 100 bp with a step.size of 50 I get a bunch of DMRs. However if I try to found difference with single base resolution I got nothing. This do not have sense to me and I suspect I'm doing something wrong during the analysis. I willl appreciate if someone can read the commands i used and help me to identify my error.
#### Methyl seq sample detail
SAMPLE
C4_ctrl BC1 REP1 control
C7_ctrl BC5 REP2 control
C4_PEG BC3 REP1 treatment
C7_PEG BC7 REP2 treatment
### In this example I'm comparing control=ctrl samples against the
## treatment=PEG in any context (CG, CHG, CHH) from; chr1
> library(methylKit);
> library(graphics);
> comparisonName="C_ctrl_vs_C_PEG_anyContext_site_by_site";ou
> chrom="chr1";
> C4_ctrl="BC1"; C7_ctrl="BC5";
> C4_PEG="BC3"; C7_PEG="BC7";
> ctrl1=paste(C4_ctrl, chrom, "MethylKit.inputdata.anyContext.txt", sep=".");
> ctrl2=paste(C7_ctrl, chrom, "MethylKit.inputdata.anyContext.txt", sep=".");
>PEG1=paste(C4_PEG, chrom, "MethylKit.inputdata.anyContext.txt", sep=".");
>PEG2=paste(C7_PEG, chrom, "MethylKit.inputdata.anyContext.txt", sep=".");
> file.list <- list(ctrl1,ctrl2,PEG1,PEG2);
> myobj <- read(file.list, sample.id = list(ctrl1,ctrl2,PEG1,PEG2), assembly="TAIR10", treatment=c(1,1,0,0), context="anyContext");
filtered.myobj=filterByCoverage(myobj,lo.count=10,lo.perc=NULL,hi.count=NULL,hi.perc=99.9)
>meth <- unite(myobj, destrand = FALSE);
>myDiff <- calculateDiffMeth(meth, num.cores=4);
#filtering for significant differential methylation
> myDiff.hyper <- get.methylDiff(myDiff, difference = 10, qvalue = 0.01, type = "hyper");
> myDiff.hypo <- get.methylDiff(myDiff, difference = 10, qvalue = 0.01, type = "hypo");
#output data for hypermethylated, hypomethylated, and all tiles
> outfile1 = paste(comparisonName, chrom, "hyper", "txt", sep=".");
write.table(as.data.frame(myDiff.hyper), file=outfile1, sep="\t", col.names=TRUE, row.names=FALSE);
> outfile2 = paste(comparisonName, chrom, "hypo", "txt", sep=".");
write.table(as.data.frame(myDiff.hypo), file=outfile2, sep="\t", col.names=TRUE, row.names=FALSE);
> outfile3 = paste(comparisonName, chrom, "united", "txt", sep=".");
write.table(as.data.frame(meth), file=outfile3, sep="\t", col.names=TRUE, row.names=FALSE);
###the result is that outfile 1 and outfile2 are empty.. What I did wrong?
###I will appreciate any help
The problem is that I got a file that is empty.
I understand that outfiles listing hypermethylated or Hypomethylated can be empty simply because there are not DM sites between the compared samples. However I'm sure this is not the case and please let me explain before go to my question.
I have 4 samples ( control vs treatment with 2 replicates each). If I compare tiles of 100 bp with a step.size of 50 I get a bunch of DMRs. However if I try to found difference with single base resolution I got nothing. This do not have sense to me and I suspect I'm doing something wrong during the analysis. I willl appreciate if someone can read the commands i used and help me to identify my error.
#### Methyl seq sample detail
SAMPLE
C4_ctrl BC1 REP1 control
C7_ctrl BC5 REP2 control
C4_PEG BC3 REP1 treatment
C7_PEG BC7 REP2 treatment
### In this example I'm comparing control=ctrl samples against the
## treatment=PEG in any context (CG, CHG, CHH) from; chr1
> library(methylKit);
> library(graphics);
> comparisonName="C_ctrl_vs_C_PEG_anyContext_site_by_site";ou
> chrom="chr1";
> C4_ctrl="BC1"; C7_ctrl="BC5";
> C4_PEG="BC3"; C7_PEG="BC7";
> ctrl1=paste(C4_ctrl, chrom, "MethylKit.inputdata.anyContext.txt", sep=".");
> ctrl2=paste(C7_ctrl, chrom, "MethylKit.inputdata.anyContext.txt", sep=".");
>PEG1=paste(C4_PEG, chrom, "MethylKit.inputdata.anyContext.txt", sep=".");
>PEG2=paste(C7_PEG, chrom, "MethylKit.inputdata.anyContext.txt", sep=".");
> file.list <- list(ctrl1,ctrl2,PEG1,PEG2);
> myobj <- read(file.list, sample.id = list(ctrl1,ctrl2,PEG1,PEG2), assembly="TAIR10", treatment=c(1,1,0,0), context="anyContext");
filtered.myobj=filterByCoverage(myobj,lo.count=10,lo.perc=NULL,hi.count=NULL,hi.perc=99.9)
>meth <- unite(myobj, destrand = FALSE);
>myDiff <- calculateDiffMeth(meth, num.cores=4);
#filtering for significant differential methylation
> myDiff.hyper <- get.methylDiff(myDiff, difference = 10, qvalue = 0.01, type = "hyper");
> myDiff.hypo <- get.methylDiff(myDiff, difference = 10, qvalue = 0.01, type = "hypo");
#output data for hypermethylated, hypomethylated, and all tiles
> outfile1 = paste(comparisonName, chrom, "hyper", "txt", sep=".");
write.table(as.data.frame(myDiff.hyper), file=outfile1, sep="\t", col.names=TRUE, row.names=FALSE);
> outfile2 = paste(comparisonName, chrom, "hypo", "txt", sep=".");
write.table(as.data.frame(myDiff.hypo), file=outfile2, sep="\t", col.names=TRUE, row.names=FALSE);
> outfile3 = paste(comparisonName, chrom, "united", "txt", sep=".");
write.table(as.data.frame(meth), file=outfile3, sep="\t", col.names=TRUE, row.names=FALSE);
###the result is that outfile 1 and outfile2 are empty.. What I did wrong?
###I will appreciate any help