WGBS Data Analysis Related
1
0
Entering edit mode
4.5 years ago

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 • 2.3k views
ADD COMMENT
0
Entering edit mode

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

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

Sorry, not understood, where should I check this?

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

Post few lines from mydiff25p object and refseq.hg18.bed.txt @ shrinka.genetics

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

ok, how will I do that? by what command?

Regards Shrinka

ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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

Thanks Shrinka

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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

Thanks Shrinka

ADD REPLY
0
Entering edit mode
  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.
ADD REPLY
0
Entering edit mode

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

What do you mean by text matrix?

Regards Shrinka

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

Yes it worked

Thanks Shrinka

ADD REPLY

Login before adding your answer.

Traffic: 2594 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