Converging multiple methylraw objects with 5hmc data into one

[hchetia]

Hi, I have created a bunch of adjusted 5hmc values using corresponding BS and oxBS values (derived using bismark coverage files).

file.list=as.list(list.files(pattern` = ".*.cov"))

OXBS_21=methRead("5268_ctrl.cov",sample.id="5268_ctrl",assembly="hg38",header=F,skip=0,sep="\t",context="CpG", pipeline="bismarkCoverage", mincov = 1)
BS_21=methRead("5269_ctrl.cov",sample.id="5269_ctrl",assembly="hg38",header=F,skip=0,sep="\t",context="CpG", pipeline= "bismarkCoverage", mincov = 1)

ADJUSTED_5HMC_21=adjustMethylC(BS_21, OXBS_21)

OXBS_22=methRead("5270_test.cov",sample.id="5270_test",assembly="hg38",header=F,skip=0,sep="\t",context="CpG",pipeline="bismarkCoverage",mincov = 1)
BS_22=methRead("5271_test.cov",sample.id="5271_test.cov",assembly="hg38",header=F,skip=0,sep="\t",context="CpG",pipeline="bismarkCoverage",mincov = 1)

ADJUSTED_5HMC_22=adjustMethylC(BS_22, OXBS_22)

How do I converge ADJUSTED_5HMC_21 (being control) and ADJUSTED_5HMC_22 (test) into a single methylkit object so that I can carry out further differential analysis?

Regards, Hasna

[alexg9010]

Hi @hcheatia,

You should be able to create the methylrawlist object for further analysis by using the methylRawList() constructor function. This takes a list of methylraw objects and a treatment vector.

Best, Alex

[hchetia]

Hi @alexg9010 I get the foll. error while doing methylRawList("ADJUSTED_5HMC_21, ADJUSTED_5HMC_22", treatment=c(1,0))

Error in methylRawList("ADJUSTED_5HMC_21, ADJUSTED_5HMC_22", treatment = c(1, : all elements in '...' must be methylRaw objects

[alexg9010]

H @hchetia,

As the error message suggests, you should pass the methylRaw objects. Thus, it should work if you leave the extra quotes out.

methylRawList(ADJUSTED_5HMC_21, ADJUSTED_5HMC_22, treatment=c(1,0))

Best Alex

[hchetia]

Hi @alexg9010 I passed the methylraw objects as suggested.

meth_5hmC_cov5=methylRawList(ADJUSTED_5HMC_21, ADJUSTED_5HMC_22, treatment=c(1,0))

Then I ran the unite function meth_5hmC_cov5 = unite(myobj_5hmC_cov5, destrand=FALSE, mc.cores = 4)

The object looks the way it should look-

However, on running clustering or PCA, only the BS samples are showing up on the plots. Now, I am not sure if somehow the object only ended up with the BS data (which should be impossible given the previous steps I hav described in my earlier issues) or is it more of a labelling issue?

Best, Hasna

[hchetia]

@al2na @alexg9010 looking for a response here.

[alexg9010]

Hi @hchetia, I think the output looks like expected.

Given your code, you loaded the BS (5mC+5hmC) and oxBS (only 5mC) data for each of your samples and estimated 5hmC methylation by adjusting BS with oxBS:

ADJUSTED_5HMC_21=adjustMethylC(BS_21, OXBS_21)

Here, the features (sample.id, assembly, ... ) from BS_21 are kept and counts are updated, thus only the BS samples remain but with 5hmC counts.

What were you trying to achieve instead?

[hchetia]

Hi @alexg9010 Okay that clarifies my doubt. Is there a way to change the features for the plot?

[alexg9010]

Hi @hchetia

what exactly do you want to change?

If you only want to update the sample names, you may use something like this: getSampleID(meth_5hmC_cov5) <- c("new_first","new_second","new_third", .... )

If you want to inspect the pca object, you may do so with: pca_res = PCASamples(meth_5hmC_cov5, obj.return = TRUE)

[hchetia]

This worked. Thanks