This is a test version of Biostars. For the public version, visit https://www.biostars.org.
WGBS Data Analysis Related

I am working on WGBS Data and using MethyKIt package. When I used 8 samples in groups, it worked fine. Now on another set of data, where just 1 sample in each group is there. Here I am facing a problem

myDiff=calculateDiffMeth(meth,mc.cores=2) two groups detected: will calculate methylation difference as the difference of treatment (group: 1) - control (group: 0)

NOTE: performing 'fast.fisher' instead of 'F' for two groups testing. Error in mclapply(cntlist, function(x) fast.fisher(matrix(as.numeric(x), : 'mc.cores' > 1 is not supported on Windows

How can I correct this error? So that in later steps which is working on the object myDiff can be obtained.

Regards Shrinka

sequencing

'mc.cores' > 1 is not supported on Windows

then don't use Windows, run you R code in a Linux machine or use Docker

Thanks for your tip. I have done it as I put mc.cores=1, It worked.

after that step I am getting another error. Any clue how to solve this?

> annotateWithGeneParts(as(myDiff25p,"GRanges"),gene.obj)

Error in NSBS(i, x, exact = exact, strict.upper.bound = !allow.append, : subscript contains NAs In addition: Warning messages: 1: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 2: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 3: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 4: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 5: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 6: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.)

Thanks Shrinka

check the contig/chromosome names (for eg. chr1 vs 1 or contig 1 vs 1) in bed and cpg files used.

how did you create gene.obj ? Print few lines from myDiff25p object.

Hello I made gene.obj by using this command, where "refseq.hg18.bed.txt" file is there with MethyKit package

gene.obj<-readTranscriptFeatures("refseq.hg18.bed.txt")

Regards

Shrinka

Here are these. PFA first 10 lines of refseq.hg18.bed.txt file.

methylDiff object with 6 rows

 chr   start     end strand       pvalue       qvalue meth.diff

385 21 9483774 9483774 + 4.462922e-04 0.0068459260 -45.83333 1032 21 9656103 9656103 + 4.325999e-04 0.0066740425 25.85200 1220 21 9674182 9674182 + 4.919951e-06 0.0001579141 38.53448 1298 21 9679300 9679300 + 1.461574e-04 0.0027225173 -48.14815 1597 21 9753103 9753103 + 2.995814e-05 0.0007255432 77.27273

1898 21 9826814 9826814 + 2.030031e-04 0.0035856587 -36.45594

sample.ids: delta21-mod gamma21-mod destranded FALSE assembly: hg18 context: CpG treament: 1 0 resolution: base [1] 11082 7 GRangesList object of length 4: $exons GRanges object with 22359 ranges and 2 metadata columns: seqnames ranges strand | score name <rle> <iranges> <rle> | <numeric> <character> [1] chr20 33506637-33506986 + | 1 NM_007186 [2] chr20 33509539-33509609 + | 2 NM_007186 [3] chr20 33511217-33511339 + | 3 NM_007186 [4] chr20 33513504-33513792 + | 4 NM_007186 [5] chr20 33514814-33514870 + | 5 NM_007186 ... ... ... ... . ... ... [22355] chr22 49561066-49561145 - | 5 NM_001130922 [22356] chr22 49561964-49562043 - | 4 NM_001130922 [22357] chr22 49563246-49563275 - | 3 NM_001130922 [22358] chr22 49567482-49567645 - | 2 NM_001130922 [22359] chr22 49568795-49568953 - | 1 NM_001130922


seqinfo: 3 sequences from an unspecified genome; no seqlengths

... <3 more elements>

I am clueless why the dimension is showing NULL NULL

Regards Shrinka

I see that methyldiff object has chromosome is represented by "21" (ti.e without word "chr") where as bed file (refseq.hg18.bed.txt) seems to have "chr" appended to chromosome name. Remove "chr" from the chromosome names of the bed file and recreate gene.obj and the run the code again.

Formatting of your post is confusing and Please format the data for better understanding.

On GNU-linux machine,

$ sed 's/^chr//g' refseq.hg18.bed.txt > refseq.hg18.modified.bed.txt

Use refseq.hg18.modified.bed.txt, for building gene.obj.

Btw, if you are using refseq.hg18.bed.txt from methylkit package, you must know that it is truncated bed file. Are you using this bed.txt for your work or for reproducing methylkit tutorial?

Hello I am using this bed file for my work, to analyse my data.

Thanks Shrinka

Before I comment further, could you please post the source of the bed.txt file?

I can attach it off course. I mentioned It comes with Methylkit package itself. It is in exdata folder of that package only.

Thanks Shrinka

That bed.txt file has only these chromosomes: chr20, chr21, chr22, chr21_random, chr22_random, chr22_h2_hap1. Full bed file (with all main chromosomes) is supposed to be used in the analysis, instead of truncated bed file (meant for a tutorial only) unless you are looking for DMRs in chr 20,21 and 22 only. Btw, for alignment, did you use bismark? if so, which genome fasta you have used?

No I have not used bismark, as my input file was in matrix form, text file.

Thanks Shrinka

  1. Remove the chr in bed.txt, rebuild the gene.obj with the new bed.txt. Please post here, if it works.
  2. Please inquire how the text (matrix) file is generated. There might be a scientific issue, as I understand esp reference genome part.
  3. The bed file, currently in use, is partial and outdated. Current version is hg38 (hg 38 > hg19 > hg18). Because of these reasons, work will not be scientifically up to date. However, this depends on the reference file used in building the matrix.

Thanks, I will work on that and will keep posting here.

What do you mean by text matrix?

Regards Shrinka

as my input file was in matrix form, text file..this file

1 answer

The mc.cores parameter refers to the use of mclapply, a parallelised version of lapply which does not support forking (so multiple workers) on Windows. Leave it at 1 if you are on Windows or switch to a Unix machine.

Thanks for the tip! But after that step I am getting another error. Any clue how to solve this?

> annotateWithGeneParts(as(myDiff25p,"GRanges"),gene.obj)

Error in NSBS(i, x, exact = exact, strict.upper.bound = !allow.append, : subscript contains NAs In addition: Warning messages: 1: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 2: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 3: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 4: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 5: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.) 6: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.)

Thanks Shrinka

No idea from remote, sorry.

Log in to answer this question.