Re: MEDIPS: how does MEDIPS define a methylated region / cluster?

--_002_3CEBD81408766A4082523CA416461083CA89DEE8MBX08stfnusedus_
Content-Type: text/plain; charset="us-ascii"
Content-Transfer-Encoding: quoted-printable

Dear Lukas,
I went back and tried to find those with Inf but I can't. I don't remember =
where I put the files that contained those values.

I did however find those that contained the NA values and I have attached t=
he file containing some of these.

I was thinking whether I should use the "transf" parameter so as all values=
 are set to between 0 and 1000 for RMS. Perhaps this would not give the NA =
or Inf values

Thanks,
Allen
________________________________________
From: Lukas Chavez [lukas.chavez.mailings-gM/[email protected]]
Sent: Saturday, September 13, 2014 1:12 AM
To: Chong Kim San Allen
Cc: Allen [guest]; [email protected]; MEDIPS Maintainer
Subject: Re: [BioC] MEDIPS: how does MEDIPS define a methylated region / cl=
uster?

Dear Allen,

please see my comments below.

Best,
Lukas

On Fri, Sep 12, 2014 at 3:00 AM, Chong Kim San Allen <[email protected]<ma=
ilto:[email protected]>> wrote:
Dear Lukas,
Firstly, thank you so very much for your kind replies. You have been so pat=
ient and I am very, very grateful for your help.

I have been away and so I'm trying now to gather my thoughts on this again.

1) I did read section 7.7 of the MEDIPS manual. My understanding (correct m=
e if I am wrong) is that MEDIPS identifies differentially methylated region=
s (DMRs) by a sample-against-sample method rather than a group-against-grou=
p method. So, if I have 5 controls (A,B,C,D,E) and 6 (1,2,3,4,5,6) test sam=
ples, then EACH individual control sample is compared against EACH test sam=
ple. So, A vs. 1, A vs. 2, A vs. 3, A vs. 4, A vs. 5, A vs. 6, B vs. 1, B v=
s. 2, ......etc. instead of comparing ALL of the control samples against AL=
L of the test samples. Is my understanding correct?

No that is not correct. MEDIPS identifies differentially enriched regions b=
y comparing two groups of samples.

You are therefore saying that because of this all-against-all method that M=
EDIPS uses, correcting for CpG density is actually not necessary. But that'=
s for identifying differentially methylated regions in MEDIPS. I don't quit=
e understand how this relates to my Principal Component Analysis since I be=
lieve that I need to export my data out from MEDIPS and maybe use a differe=
nt R package to do the PCA. So I am guessing that I need to export the meth=
ylation profile files for each sample and then input them into the PCA anal=
ysis software.

The motivation of CpG normalization is to quantify MeDIP signals at genomic=
 regions with varying CpG density across a genome. The differential enrichm=
ent analysis compares IP-seq derived signals at genomic windows across samp=
les. Therefore, the signal for a specific window will not be affected by di=
fferent CpG densities across samples. The same assumption applies when you =
are comparing different samples with each other employing PCA. Due to a lac=
k of experience I cannot recommend to use CpG normalized methylation estima=
tes or for example variance stabilized counts as done by DESeq for RNA-seq =
data.


2) I have decided to try and do a PCA/hierarchical clustering using only th=
e methylation levels of a specific genomic feature (eg. all gene promoters)=
 to make things simple. I have decided to take your suggestion. I believe I=
 can use the MEDIPS.createROIset() to calculate the methylation levels at a=
ll gene promoters. I can then use the methylation profile for the ROI (in t=
his case, promoter) from each sample to do a a hierarchical clustering or P=
CA.

3) You said: " In case you want to apply a PCA on CpG density methylation e=
stimates anyway, you can either use the rms values ......."...however, I st=
ill don't know if RMS is a linear or logarithmic value. If I have 4 100-bp =
windows where the RMS is 3, 5, 8, 9, can I say that the RMS for the this 40=
0-bp region is 25? If RMS is a log value then obviously I cannot just add t=
hem up.
The rms values are not log transformed. However, the signals in neighbourin=
g windows are not independent. Extended reads spanning more than one window=
 will be counted several times. Therefore, adding up counts or rms values i=
s not valid. You can use MEDIPS.createROIset() instead.

4) Also, I was looking at the methylation profile file and I noticed that s=
ome RMS values are given as "Inf" or "NA". I am confused as to how to use o=
r interpret this. I also don't understand why it is "NA"...why not just put=
 a 0 (zero) instead of an "NA"?
What are the CF and count values for windows where rms values are NA or Inf=
? Can you send examples?

5) Thanks for the suggestion to use the MEDIPS.meth() function to remove wi=
ndows with low or zero counts. That's really helpful. I hadn't even conside=
red that.

Have a great weekend.

Warmest regards,
Allen




________________________________________
From: Lukas Chavez [lukas.chavez.mailings-gM/[email protected]<mailto:lukas.chave=
z.mailings-gM/[email protected]>]
Sent: Tuesday, September 09, 2014 1:09 AM
To: Chong Kim San Allen
Cc: Allen [guest]; [email protected]<mailto:bioconductor@r-project=
.org>; MEDIPS Maintainer
Subject: Re: [BioC] MEDIPS: how does MEDIPS define a methylated region / cl=
uster?

Dear Allen,

I was referring to section 7.7 of the MEDIPS Tutorial (or vignette) availab=
le at:
http://www.bioconductor.org/packages/release/bioc/vignettes/MEDIPS/inst/doc=
/MEDIPS.pdf

I recommend to apply the a PCA to the counts at the genomic windows returne=
d by the MEDIPS.meth() function where you probably want to remove windows w=
ith zero or low counts first. Here, the assumption is that CpG density corr=
ection of MeDIP signals might not be necessary as you are comparing differe=
nt samples with each other (please see section 7.7 of the MEDIPS tutorial).=
 In case you want to apply a PCA on CpG density methylation estimates anywa=
y, you can either use the rms values or calculate methylation estimates usi=
ng  BayMeth, a recent tool that appears to be superior for such methylation=
 estimates.

All the best,
Lukas

On Sun, Sep 7, 2014 at 11:10 AM, Chong Kim San Allen <[email protected]<ma=
ilto:[email protected]><mailto:[email protected]<mailto:[email protected]=
g>>> wrote:
Dear Lukas,
Again, I am so grateful for your quick reply.

The only manual for MEDIPS that I know of is this one: http://www.bioconduc=
tor.org/packages/release/bioc/manuals/MEDIPS/man/MEDIPS.pdf
and I have looked and just can't find the section 7.7 that you mentioned.

The reason I asked about how MEDIPS defines a cluster is because I wanted t=
o try and reduce the number of methylated regions that I look at and then u=
se these "consolidated" regions to do a principal component analysis (PCA).=
 I have 12 cell types and I wanted to see if I could either do a PCA or Uns=
upervised Hierarchical Clustering to separate these 12 cell types into subg=
roups, based on their methylation profile.

You said: "When windows with significant differential enrichment between gr=
oups have been identified, the MEDIPS.mergeFrames() function can be used to=
 merge adjacent significant windows into extended regions."
But at this point, I don't have any groupings. Normally, one would have a c=
ontrol and test group and so you can do what you suggested. But I have only=
, for example, a test group and I want to know if this test group can be fu=
rther subdivided into a few subgroups based on their methylation profile. I=
 would think that it would be computationally intensive if I used original =
100bp windows generated by MEDIPS and entered all the rms values for each w=
indow of each cell type into the PCA.

You said: "Instead of merging any values across windows, I recommend to mak=
e use of the MEDIPS.createROIse() function which allows for analyzing speci=
fic regions of interest."
I didn't want my PCA analysis to be biased by a predefined "region of inter=
est", be it whether the ROI is a promoter, CpG, etc. I wanted to just input=
 the raw methylation scores (rms) to see what I could get. So that is why I=
 asked if the rms can be added up (in question 4 of my last email)...also I=
 was curious if it could be done.

Once again., thanks so very much for taking time from your busy schedule to=
 answer my questions.

Best regards,
Allen

________________________________________
From: Lukas Chavez [lukas.chavez.mailings-gM/[email protected]<mailto:lukas.chave=
z.mailings-gM/[email protected]><mailto:lukas.chavez.mailings-gM/[email protected]<mail=
to:lukas.chavez.mailings-gM/[email protected]>>]
Sent: Saturday, September 06, 2014 11:09 AM
To: Allen [guest]
Cc: [email protected]<mailto:[email protected]><mailto:bi=
[email protected]<mailto:[email protected]>>; Chong Kim San=
 Allen; MEDIPS Maintainer
Subject: Re: [BioC] MEDIPS: how does MEDIPS define a methylated region / cl=
uster?

Dear Allen,

please see my comments below. Moreover, please consider reading section 7.7=
 of the MEDIPS manual.

Best,
Lukas


On Fri, Sep 5, 2014 at 1:14 AM, Allen [guest] <[email protected]<mailt=
o:[email protected]><mailto:[email protected]<mailto:guest@biocon=
ductor.org>><mailto:[email protected]<mailto:[email protected]><m=
ailto:[email protected]<mailto:[email protected]>>>> wrote:
Hi,
I have aligned my sequencing data for one sample (filename: Ca2_MAPQ20.bam =
--- because I filtered out everything with MAPQ20 or less) and I've put it =
through MEDIPS. At the moment, I am not making a comparison between samples=
. Therefore, the result file shows the methylation profile for only one sam=
ple. In this results file, I noticed that the data is organized as such:
chr     start   stop    CF      Ca2_MAPQ20.bam.counts   Ca2_MAPQ20.bam.rpkm=
     Ca2_MAPQ20.bam.rms      Ca2_MAPQ20.bam.prob     MSets1.counts.mean    =
  MSets1.rpkm.mean        MSets1.rms.mean MSets1.prob.mean

1) I tried reading the Down et al. (2008) paper that explains the concept o=
f coupling factor and I think it is, simply put, a measure of local CpG den=
sity. I am not sure if my understanding of 'coupling factor' is correct?

That is correct.

2) This lead to question how or what MEDIPS defines as a "region or cluster=
"? That is to say, on Chromosome 1, I have reads aligning from position "10=
02501" to "1003300" but then there is a region of 1400 bp (from "1003301" t=
o "1004700") where there are no reads aligned (rpkm =3D 0). Then, again, fr=
om "1004701" to "1005400", there reads aligning to this region. So my quest=
ion is does MEDIPS consider this to be a case of 2 methylated regions, that=
 is, "1002501-1003300" and "1004701-1005400" or does MEDIPS try to consolid=
ate these 2 methylated regions into one methylated region/cluster (that is,=
 from "1002501-1005400") since they are fairly close to one another.

Each window is considered as one region; CpG density, counts, rpkm, rms, an=
d differential coverage will be calculated for each window separately. When=
 windows with significant differential enrichment between groups have been =
identified, the MEDIPS.mergeFrames() function can be used to merge adjacent=
 significant windows into extended regions.

3) I used "uniq=3DTRUE" to get rid of stacked/clonal reads and so for each =
100 bp bin, I am usually getting rpkm=3D1 within each bin but occasionally,=
 it may be as many rpkm=3D3. This lead me to wonder how the relative methyl=
ation score ("Ca2_MAPQ20.bam.rms") is calculated? For example, 3 bins have =
a value of 1 rpkm each, but one has an rms value of "1660.409928", another =
an rms of "7122.10254" and the third is only "679.1814523". How is that pos=
sible that 3 bins that have the same rpkm can have such varying rms values?

This is due to the different CpG densities in these windows.

4)Can the rms be added up for a region so as to represent the cumulative me=
thylation level for that region. I am asking because I do not know if the r=
ms value is a log value or what? So, in my above question (3), if the 3 bin=
s are adjacent to one another and I decide to cluster them together and con=
sider them as one methylated region, then would the rms value for this new =
300bp bin that I have created be 1660.409928 + 7122.10254 + 679.1814523 =3D=
 9461.68.

Instead of merging any values across windows, I recommend to make use of th=
e MEDIPS.createROIse() function which allows for analyzing specific regions=
 of interest.

5)Finally, what does "Ca2_MAPQ20.bam.prob" represent? I figured that is som=
e sort of a probability score but I am not sure of what.

As mentioned in the manual, please disregard the prob values which will be =
removed in a future version.

Thanks in advance for any assistance in answering my questions.

Best regards,
Allen

 -- output of sessionInfo():

error

--
Sent via the guest posting facility at bioconductor.org<http://bioconductor=
.org><http://bioconductor.org><http://bioconductor.org>.

_______________________________________________
Bioconductor mailing list
[email protected]<mailto:[email protected]><mailto:Biocon=
[email protected]<mailto:[email protected]>><mailto:Bioconducto=
[email protected]<mailto:[email protected]><mailto:Bioconductor@r-pr=
oject.org<mailto:[email protected]>>>
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatic=
s.conductor

Important: This email is confidential and may be privileged. If you are not=
 the intended recipient, please delete it and notify us immediately; you sh=
ould not copy or use it for any purpose, nor disclose its contents to any o=
ther person. Thank you.

Important: This email is confidential and may be privileged. If you are not=
 the intended recipient, please delete it and notify us immediately; you sh=
ould not copy or use it for any purpose, nor disclose its contents to any o=
ther person. Thank you.


Important: This email is confidential and may be privileged. If you are not=
 the intended recipient, please delete it and notify us immediately; you sh=
ould not copy or use it for any purpose, nor disclose its contents to any o=
ther person. Thank you.

--_002_3CEBD81408766A4082523CA416461083CA89DEE8MBX08stfnusedus_
Content-Type: text/plain; name="To Lukas.txt"
Content-Description: To Lukas.txt
Content-Disposition: attachment; filename="To Lukas.txt"; size=505;
	creation-date="Mon, 15 Sep 2014 10:20:41 GMT";
	modification-date="Mon, 15 Sep 2014 10:20:41 GMT"
Content-Transfer-Encoding: base64

DQpjaHIJc3RhcnQJc3RvcAlDRglTUlI1MTc3MjguYmFtLmNvdW50cwlTUlI1MTc3MjguYmFtLnJw
a20JU1JSNTE3NzI4LmJhbS5ybXMJU1JSNTE3NzI4LmJhbS5wcm9iCU1TZXRzMS5jb3VudHMubWVh
bglNU2V0czEucnBrbS5tZWFuCU1TZXRzMS5ybXMubWVhbglNU2V0czEucHJvYi5tZWFuDQoNCmNo
cjIyCTE2Mjg4NjAxCTE2Mjg4NzAwCTMxCTE1CTcuMjEzMzk0ODcwMzI0MDgJTkEJMS41NDQwODcw
NTA0NDI0ZS0wNwkxNQk3LjIxMzM5NDg3MDMyNDA4CU5BCTEuNTQ0MDg3MDUwNDQyNGUtMDcNCmNo
cjIyCTE5NzQ4NTAxCTE5NzQ4NjAwCTI2CTAJMAlOQQkwCTAJMAlOQQkwDQpjaHIyMgk0MjY2Njcw
MQk0MjY2NjgwMAkyNQkwCTAJTkEJMS4xMTAyMjMwMjQ2MjUxNmUtMTYJMAkwCU5BCTEuMTEwMjIz
MDI0NjI1MTZlLTE2DQpjaHIyMgk0NDcyNzcwMQk0NDcyNzgwMAkyNgkwCTAJTkEJMAkwCTAJTkEJ
MA0KY2hyMjIJNDg4ODUyMDEJNDg4ODUzMDAJMjkJMAkwCU5BCTAJMAkwCU5BCTANCg==

--_002_3CEBD81408766A4082523CA416461083CA89DEE8MBX08stfnusedus_
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

_______________________________________________
Bioconductor mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
--_002_3CEBD81408766A4082523CA416461083CA89DEE8MBX08stfnusedus_--