Hi I am running the BiSeq software and have 2 rejected clusters after running testClusters with FDR.cluster=0.1

When I attempt to run trimClusters on these two clusters I get no results unless I use a FDR.loc >= 1.5... I am having some trouble understanding what exactly this FDR.loc number means, in general and one as high as 1.5, and also why I have to use one so high to even see results? Could someone help me with some insight?

Thanks,

Marcus

Hi Katja,

Sorry one more question. Is there a way to write the DMRs object to a csv?

Thanks,

Marcus

Hi Marcus,

Are you sure that you used an FDR.loc > 1? What other values for FDR.loc did you try?

In general, it can happen that testing on CpG-wise level with the same FDR that you used for calling significant regions doesn't result in any significant CpG sites.

Cheers,

Katja

Yes, I tried incrementing the fdr.loc from .1 up to 1.5 by .1 until the results were not null. The way the code is written you are not restricted to any value <1. I tried 10 and got more locations than at 1.5. I imagine using something like 1000000 would just spit out every bp that was in the range of the cluster. This is why I am confused because statistically it seems that 1 should be the highest.

That is true, the FDR.loc should be restricted to be between 0 and 1. I will change this.

To figure out the problem here I would need reproducible data. Would it be possible for you to send me a subset of your data (e.g. only one chromosome) together with your code?

Yes, how should I send it? 

I uploaded the .r script and chr 3 of the data, .bedgraph and .cov files, compressed into a tar.gz to my github @ https://github.com/naymikm/BiSeq_subset

Thank you very much!

Marcus

Hi Katja,

I figured out the issue. There was a problem with 2 of the samples that we decided to ultimately throw out due to sequencing quality issues, without these two samples the analysis runs fine. I do have one trivial question though just to be 100% sure... In the DMR results group 1 and group 2 correspond to the same order as the levels of factor(colData(rrbs)$group)?  In my case factor(colData(rrbs)$group) shows: 

Levels: L O

Hence, L is group 1 and O is group 2?

Thanks again,

Marcus

Hi Katja,

Awesome thank you! Yes that is fine, I will leave them there for a while. To clarify, in the dataset we decided to throw out the samples T01 and  T07 which left 11 samples of L (control) and 9 samples of O (case).

Thanks,

Marcus

Hi Katja,

Another question, I wanted to clarify what the median.p value is. Is this the median p-value for all sites within each DMR from the Wald test when testing for group effect?

Thanks,

Marcus