runObservation error

[seynard]

Hi,
The PhD candidate I'm working with, @stacyrousse is trying to use MethylInheritance on RRBS data from 3 consecutive generations but she is facing an error message that we do not seem to be able to correct.
She is loading here 3 methylation files, in methylkit format, putting them into a list but when trying to run runObservation she gets the following error message uniting... Error in unite(filtered.sites, destrand = destrand) : no base were united. try adjusting 'min.per.group'.
We tested the unite function from MethylKit and it works correctly, without error on 1 generation. We also tested to add the option min.per.group to runObservation but it doesn't work.
How can we fix this issue?

Thanks in advance for you help.

Here is the code she wrote

`library(methylKit)
library(data.table)
library(dplyr)
library(methylInheritance)

famfile<-fread("samples_for_methylInheritance_test.tsv", data.table=F)
famfile$line<-factor(famfile$line, levels = c("M", "P"))
levels(famfile$line)<-c(0, 1)
famfile$line<-as.numeric(famfile$line)

list_G0<-famfile%>%filter(pop=='Gen0')%>%select(lib_id, line)%>%mutate(loc=paste0( lib_id, "_intersection.bed"))
list_G1<-famfile%>%filter(pop=='Gen1')%>%select(lib_id, line)%>%mutate(loc=paste0(lib_id, "_intersection.bed"))
list_G2<-famfile%>%filter(pop=='Gen2')%>%select(lib_id, line)%>%mutate(loc=paste0( lib_id, "_intersection.bed"))

gen_0<-methRead(location=as.list(list_G0$loc),
sample.id=as.list(list_G0$lib_id),
assembly="Cjap2.1",
treatment=c(list_G0$line),
context="CpG",
mincov=0)

gen_1<-methRead(location = as.list(list_G1$loc),
sample.id=as.list(list_G1$lib_id),
assembly="Cjap2.1",
treatment=c(list_G1$line),
context="CpG",
mincov=0)

gen_2<-methRead(location = as.list(list_G2$loc),
sample.id=as.list(list_G2$lib_id),
assembly="Cjap2.1",
treatment=c(list_G2$line),
context="CpG",
mincov=0)

final_list <- list(gen_0, gen_1, gen_2)

runObservation(methylKitData = final_list,
type = "sites",
outputDir = runObs1",
nbrCores = 1,
minReads = 10,
minMethDiff = 10,
qvalue = 0.01,
vSeed = 2101)`

[adeschen]

Hi Sonia,

Can you share the entire error output? This error is not generated by the methylIheritance package but by the methylKit package it is using.

How many differentially methylated sites do you find in each generation? Are those sites overlapping between generations?

There is a related issue on the methylKit github site when no Cpg are ovelapping over samples:
al2na/methylKit#63

At this time, I would manually check the site for each generation to see if some are overlaps.

Best,
Astrid

[seynard]

Hi Astrid,

Thanks so much for your answer.
As a baseline, we have 10 tested CpG positions that are shared by all the individuals for the 3 generations. These positions have been identified as interesting for us using a different tool than MethylKit.
We tried running a differential analysis with MethylKit and obtain the following output (I reformatted everything to have a table for further analysis.
chr start end strand pval qval diff gen 1 chr21 1042814 1042814 + 2.854197e-01 4.451395e-01 3.3776198 g0 2 chr21 1042825 1042825 + 1.164322e-01 2.269841e-01 -4.8379353 g0 3 chr21 1042835 1042835 + 8.332745e-01 7.219842e-01 0.7085092 g0 4 chr21 1042838 1042838 + 6.654676e-01 6.486629e-01 -1.3171970 g0 5 chr21 1042844 1042844 + 5.026929e-01 5.599984e-01 2.1033868 g0 6 chr21 1042849 1042849 + 5.458710e-02 1.418897e-01 -4.8419794 g0 7 chr21 1042859 1042859 + 4.833831e-01 5.599984e-01 -2.4458883 g0 8 chr21 1042864 1042864 + 9.668909e-01 7.539796e-01 -0.1422475 g0 9 chr27 3171990 3171990 + 1.077950e-02 8.405831e-02 -6.9504386 g0 10 chr27 3172028 3172028 + 3.466959e-02 1.351764e-01 -4.8895659 g0 11 chr21 1042814 1042814 + 3.452313e-03 -2.043392e-16 5.5308851 g1 12 chr21 1042825 1042825 + 3.013722e-04 -2.043392e-16 7.0048130 g1 13 chr21 1042835 1042835 + 1.512248e-03 -2.043392e-16 6.9434532 g1 14 chr21 1042838 1042838 + 2.338911e-02 -2.043392e-16 4.3134051 g1 15 chr21 1042844 1042844 + 5.493896e-02 -2.043392e-16 3.9067533 g1 16 chr21 1042849 1042849 + 1.371103e-02 -2.043392e-16 4.2644416 g1 17 chr21 1042859 1042859 + 1.798440e-06 -2.043392e-16 10.6029336 g1 18 chr21 1042864 1042864 + 1.041140e-04 -2.043392e-16 9.1067426 g1 19 chr27 3171990 3171990 + 3.171208e-03 -2.043392e-16 6.3597283 g1 20 chr27 3172028 3172028 + 5.842354e-04 -2.043392e-16 6.5322526 g1 21 chr21 1042814 1042814 + 3.452313e-03 -2.043392e-16 5.5308851 g2 22 chr21 1042825 1042825 + 3.013722e-04 -2.043392e-16 7.0048130 g2 23 chr21 1042835 1042835 + 1.512248e-03 -2.043392e-16 6.9434532 g2 24 chr21 1042838 1042838 + 2.338911e-02 -2.043392e-16 4.3134051 g2 25 chr21 1042844 1042844 + 5.493896e-02 -2.043392e-16 3.9067533 g2 26 chr21 1042849 1042849 + 1.371103e-02 -2.043392e-16 4.2644416 g2 27 chr21 1042859 1042859 + 1.798440e-06 -2.043392e-16 10.6029336 g2 28 chr21 1042864 1042864 + 1.041140e-04 -2.043392e-16 9.1067426 g2 29 chr27 3171990 3171990 + 3.171208e-03 -2.043392e-16 6.3597283 g2 30 chr27 3172028 3172028 + 5.842354e-04 -2.043392e-16 6.5322526 g2

Does your tool need a minimum differential or a pvalue threshold to be able to estimate the transmissibility across generations? Because in our case generation 1 and 2 are really similar in differential between or test and treatment but generation 0 is somehow different.

To answer the part on full error message, here it is as a screen shot :

Thanks in advance for your help. I really hope that we will find how to use the tool you developed, we are really interested in having it for our analysis.

Best regards, Sonia

[adeschen]

Hi Sonia,

In the vignette, there is a demo dataset that is used as input for the runObservation() method:

library(methylInheritence)
## Load dataset containing information over three generations
data(demoForTransgenerationalAnalysis)


## All samples related to one generation are contained in a methylRawList 
## object.
## The methylRawList object contains two Slots:
## 1- treatment: a numeric vector denoting controls and cases.
## 2- .Data: a list of methylRaw objects. Each object stores the raw 
##           mehylation data of one sample.

Have you validated that your input is in the same format?

In your associated paper, we had fewer than 50 shared sites between the three generations and managed to run the simulations successfully. However, it is unnecessary to run the simulation steps if you have zero shared sites.

Best,
Astrid

[seynard]

Thanks Astrid for the feedback. Yes our data set looks exactly like the test data.

However, we tested only 10 sites but they are all shared by all the individuals and all the generations.
I'm not sure that I understand what you're suggesting to check?

Thanks again for your help.
Cheers, Sonia

[adeschen]

Hi Sonia,

I am unsure that only looking at 10 sites makes sense, as the method includes a normalization by coverage.

You should look at the internal method methylInheritance:::runOnePermutationOnAllGenerations and run it step by step. The function only runs methylKit to identify differentially methylated sites on each generation separately.

    for (i in seq_len(nbrGenerations)) {

        allSamplesForOneGeneration <- methylRawForAllGenerations[[i]]

        ## SITES
        if (doSites && !readySites) {

            ## Filter sites by coverage
            filtered.sites <- filterByCoverage(allSamplesForOneGeneration,
                                                lo.count = minReads,
                                                lo.perc = NULL,
                                                hi.count = NULL,
                                                hi.perc = maxPercReads)

            ## Normalize coverage
            filtered.sites <- normalizeCoverage(filtered.sites, "median")

            ## Merge all samples to one table
            meth.sites <- unite(filtered.sites, destrand = destrand)

            if (length(meth.sites@.Data[[1]]) == 0) {
                stop("meth.sites IS EMPTY")
            }

            ## Get differentially methylated sites
            allSites <- suppressMessages(suppressWarnings(
                calculateDiffMeth(meth.sites, mc.cores = nbrCoresDiffMeth)))

            permutationList[["SITES"]][[i]] <- suppressWarnings(
                getMethylDiff(allSites, difference = minMethDiff,
                                qvalue = qvalue))
        }

Best,
Astrid

[seynard]

Hi Astrid,

Thanks for the suggestions.
I tried the code you sent me and after some adjustments on the values provided it runs on our dataset. However, if I try to run 'runPermutation' with the same parameters I again reach the same issue with the error message
uniting... Error in unite(filtered.sites, destrand = destrand) : no base were united. try adjusting 'min.per.group'.
I really do not understand the issue ...

I was wondering if you could schedule a video call to be able to show you live what we are doing and to have your imputes directly? (we can contact each other outside of github for that if that's ok with you).

Thanks in advance for your help.
Best regards. Sonia

[adeschen]

Hi Sonia,

If you need to modify your dataset to prevent methylKit from crashing, it indicates that the dataset is not fully compatible with methylKit. In the previous link to the associated issue, the problem was solved by aligning the reads "directionally".

Unfortunately, I don't feel a Zoom meeting will be helpful at this step. Access to the real dataset with time devoted to debugging is the next logical step. If you are interested in establishing a collaboration with CSHL, I will send you the contact info for the first two authors of the package.

Best,
Astrid

[seynard]

Ok Astrid,

Thanks for your answer.
I don't really agree with your response as our data set worked perfectly fine in methylkit to run differential analysis on 1 generation at a time. Also, our reads are already aligned directionally.
From our trials I feel that there is an issue with a parameter integrated in the runObservation and runPermutation function that is unclear to us (if there is a constrain on the minimum number of sites to use for the estimation that should be probably mentioned somewhere).

The idea of a Zoom was so we could show you live what we were running, our data and the error messages obtained. If you believe that another contact could help me out understand what we miss to run your package successfully then I would be really grateful if you could provide me their names.

Cheers, Sonia

[adeschen]

Hi Sonia,

The method has been developed to utilize all available sites, as it performs permutations on all chromosomes (see the associated paper).

I still don't see how running the code would be more helpful than seeing the command line you use.

Best,
Astrid

[seynard]

Ok, well noted. Thanks for the support.
Best, Sonia