Question about paired experimental design with correction plus testing interaction effect?

[KristenKrolick]

@haowulab

Hello, I just performed a paired analysis. For my next analysis, I will need the methylation values per individual sample for each of my significant CpG sites in the paired analysis. Since I performed the DMLfit to correct for differences in cell types, I do not want to take the original methylation values (from my bismark cov files) for performing my next analysis. Is there a way to get the DSS corrected methylation values per individual sample back?

Thank you so much for your help and tool!
Kristen

[KristenKrolick]

@haowulab
Nevermind, I just read your paper and FAQ on the DSS manual page about how relative methylation values should not be extracted here.

Instead, Do you know of a way I can use the results from my paired analysis to perform a second statistical test? Basically, I am trying to find a work around here. I found those CpG sites that were significantly different between before and after surgery (after correcting for cell type fractions), and now I want to take those sites and see if they are associated with phenotypes of surgical complications...

Code:
Treatment = factor(rep(c("Preop","Postop"),each = 51))
pair = factor(rep(1:51,2))
design = data.frame(Treatment, pair, covar)

DMLfit = DMLfit.multiFactor(chr2_bismarkBSseq, design, formula = ~Treatment + pair + Mono_3 + NK_3 + B_TOTAL + CD4_TOTAL + CD8_TOTAL + GRANULOCYTES_TOTAL)
dmlTest = DMLtest.multiFactor(DMLfit, term="Treatment")

[haowulab]

Sorry for my tardy reply. You can certainly analyze paired data using DSS multiple factor function. But you need to be careful to setup the design. Please read this for RNA-seq analysis using DESeq2: https://support.bioconductor.org/p/84241/#84256. Basically, the effect you try to test should be an interaction term.

[KristenKrolick]

@haowulab
Thank you so much for responding! I am reading the example you sent, and the experimental design is a very similar example to what I am doing, thank you. However, I am getting a little confused between the DESeq2 example you sent plus and the DSS guide for setting up paired design (see my question below):

Q1) Could you please help with how I setup (the paired design plus test for an interaction of chronic pain response) exactly using DSS multiple factor function for bsseq data?
I have n = 51 preoperation and postoperation patients.
Need to correct for different cell type fractions before surgery and after surgery.
For those CpG sites differently methylated after surgery, need to test their association with chronic pain (so as you said setup as interaction effect somehow). I will call chronic pain the "response". Where would I put the response variable in my equation below? The DESeq2 example was showing their patient.id already setup as a sort of paired design and then the response variable being in the equation like so: ~Patient.ID + treatment * (status + Response), w/ status being what they wanted to correct for.

So far, mine looks like this:

DMLfit = DMLfit.multiFactor(BSobj, design, formula = ~ Treatment + pair + Mono_3 + NK_3 + B_TOTAL + CD4_TOTAL + CD8_TOTAL + GRANULOCYTES_TOTAL)

I had set it up using:

covar <- read.csv(#path to .csv containing my cell type fractions to correct for next to filename of each patient...

Treatment = factor(rep(c("Preop","Postop"),each = 51))
pair = factor(rep(1:51,2))
design = data.frame(Treatment, pair, covar)

And where my design dataframe ended up looking like this: design.csv

I know I need to add in a "response" column to my dataframe in which 1= pain and 0 = no pain; then somewhere in my equation add in the response. Also not sure but DESeq2 example there were talking about it being important that their paired design was set up 1,1,2,2,3,3,4,4,5,5,6,6 instead of 1,2,3,4,5,6,1,2,3,4,5,6 for some reason? Should I change mine?

Sincerely grateful for your help and time to allow us to perform this analysis,
Kristen

[haowulab]

The paired design described in the DSS manual is a very simple one, where only treatment (surgery in your design) is tested. In your study, I believe you want to test the association between the change of methylation before/after surgery and chronic pain. So you want to study the change (with chronic pain) of change (with surgery), which is an interaction.

Do you have cell type proportions before and after surgery? If so, your formula should be

~ Treatment + pair + Response*( Treatment + Mono_3 + NK_3 + B_TOTAL + CD4_TOTAL + CD8_TOTAL + GRANULOCYTES_TOTAL)

Then you want to test the Response*Treatment term. For that you need to figure out where it is in the design matrix or the name for the term.

I wrote a simple simulation for you to play with. See below. You can study this to understand the concept.

### do a paried test simulation                                                                     
x.bef = rnorm(10)
x.aft = x.bef + 1 + rnorm(10,sd=0.5)
y.bef = rnorm(10, mean=1)
y.aft = y.bef - 1 + rnorm(10,sd=0.5)

boxplot(x.bef,x.aft, y.bef, y.aft)
boxplot(x.aft-x.bef, y.aft-y.bef)

## to test without considering pairing - not significant                                            
t.test(c(x.bef, y.bef), c(x.aft,y.aft))

## to test with considering pairing - significant                                                   
t.test(x.aft-x.bef, y.aft-y.bef)

## do a regression                                                                                  
dat = c(x.bef, y.bef, x.aft, y.aft)

condition = factor(rep(c("Preop","Postop"),each = 20))
pair = factor(rep(1:20,2))
Trt = factor( rep(rep(c("trt1","trt2"), each = 10), 2))
design = data.frame(condition, pair, Trt)

dd = data.frame(dat, design)
fit = lm(dat~pair + condition + Trt + condition:Trt)
summary(fit)
## note, the coefficient for the interaction term is the same as the t-test statistics              

[KristenKrolick]

@haowulab

Thank you Hao! Yes, I have cell type proportions both before and after surgery (please see
design.csv). Went through the toy example you sent to try to understand it a little better, but mainly had to ask our collaborating Statistician for help.

The ~ Treatment + pair + Response*( Treatment + Mono_3 + NK_3 + B_TOTAL + CD4_TOTAL + CD8_TOTAL + GRANULOCYTES_TOTAL) you sent returned the error: Error in solve.default(XTVinv %*% X) : system is computationally singular.... Our collaborating Statistician said that would test pain specific cell type effect (not our question).

She told me to use the following formula instead :
R DMLfit = DMLfit.multiFactor(chr1_bismarkBSseq, design, formula = ~Treatment + pair + response + Mono_3 + NK_3 + B_TOTAL + CD4_TOTAL + CD8_TOTAL + GRANULOCYTES_TOTAL + (Treatment*response)) . I just wanted to confirm with you that this formula will give us the pain-specific pre-post difference correcting for cell types? I am trying to understand which places of before and after tilde correspond to the x and y in "y=Bo + B1x"? Are we sure formula should not be ~pair instead of ~treatment?
If I run this formula, and then run dmlTest = DMLtest.multiFactor(DMLfit, term="Treatment*response"), the following error is returned: Error in makeContrast(DMLfit, term) : Some term(s) to be tested can't be found in the formula.

Thank you for being so helpful to allow us to perform this analysis!
Kristen

[haowulab]

The singularity problem is because the design matrix is not of full rank. Perhaps you have to remove one of the cell type proportions in the model, since they sum up to 1 they are not independent.

For the formula suggested by your statistician: you can use that formula too. But if the cell type proportions are independent with Response, it won't make much difference. However, if the proportions are associated with Response, the interpretation will be different. Even with the interaction, if you test Response:Treatment, it's not testing the "pain specific cell type effect". Such effect can be tested by testing things like Response:Mono_3, etc. You can play with the following codes:

x.bef = rnorm(10)
x.aft = x.bef + 1 + rnorm(10,sd=0.5)
y.bef = rnorm(10, mean=1)
y.aft = y.bef - 1 + rnorm(10,sd=0.5)

boxplot(x.bef,x.aft, y.bef, y.aft)
boxplot(x.aft-x.bef, y.aft-y.bef)

## to test without considering pairing - not significant
t.test(c(x.bef, y.bef), c(x.aft,y.aft))

## to test with considering pairing - significant
t.test(x.aft-x.bef, y.aft-y.bef)

## do a regression
x = c(x.bef, x.aft)
y = c(y.bef, y.aft)
dat = c(x.bef, y.bef, x.aft, y.aft)

pair = factor(rep(1:20,2))
Treatment = factor(rep(c("Preop","Postop"),each = 20))
Response = factor( rep(rep(c("Good","Bad"), each = 10), 2))
## randomly generate proportions, assuming two cell types
ct1 = runif(length(condition))
ct2 = 1-ct1

## now do regression

## without proportion
fit = lm(dat ~ pair + Treatment + Response + Response:Treatment)
summary(fit)
## same as t-test

## with proportion, but without proportion by Response interaction 
fit = lm(dat ~ pair + Treatment + Response + ct1 + ct2 + Response:Treatment)
summary(fit)
## a little different

## with proportion by Response interaction 
dd = data.frame(dat, design)
fit = lm(dat ~ pair + Treatment + Response*(Treatment+ct1+ct2))
summary(fit)
## slightly different 

[haowulab]

By the way, I'm curious how you get the cell type proportions. Did you estimate them from the data?

[KristenKrolick]

@haowulab

Hmm I am in the middle of studying your examples (and multiple regression in general), but I do not understand:
Q1) The 'fit = lm()' function in R seems to be a little different from your 'DMLfit = DMLfit.multifactor function'-- as in the dependent variable (aka the y in y = B0 + B1x) usually goes to the left of the tilde. In your function are we assuming that to the left of the tilde would be the Bsobj, the methylation values?
Q2) I am not understanding how we are testing the difference in methylation between preop and postop (after correcting for blood cell type fractions) from this formula? -- I suppose it's okay if I don't quite understand this one, as long as you and our Statistician confirm this for me.
Q3) Okay so if I use the DMLfit = DMLfit.multiFactor(chr1_bismarkBSseq, design, formula = ~Treatment + pair + response + Mono_3 + NK_3 + B_TOTAL + CD4_TOTAL + CD8_TOTAL + GRANULOCYTES_TOTAL + (Treatment*response)). I get the error, Error in makeContrast(DMLfit, term) : Some term(s) to be tested can't be found in the formula . How would I fix this?

For my blood fraction in silico estimation:

1. I used the Reference provided by Salas et al. 2022 (https://www.nature.com/articles/s41467-021-27864-7), actual Ref available:(https://github.com/immunomethylomics/FlowSorted.BloodExtended.EPIC/tree/main).
2. I converted this reference from the Illumina Array sites into UCSC chromosomal coordinates using the manifest provided by Illumina.
3. I used Epidish to perform the deconvolution.
   All of this is outlined and the code I used here: https://github.com/KristenKrolick/Chidambaran_Lab_Manuscript/blob/main/n)%20Epidish
   p.s. We had samples that we had flowcyto data for as well as whole-genome meth seq, so I was able to perform a Pearson correlation, and (at least for those cells whose same cell markers were used in our flow data vs. salas ref cell markers he used to make his Ref), they correlated! So, I really liked his Reference.

[KristenKrolick]

Oh! And Q4) Your DSS manual outlines reasoning for why one would not need to subset their whole-genome bisulfite seq data by minimum coverage threshold: does the same apply for the general hypothesis testing design that I am using??

...I would think I would want to subset my data by at least mincov 4 or so, just to qualify that I am getting an accurate methylation call since my seq. data is from whole blood (mixed cells)...

Thank you!

[haowulab]

See my comments below.

On Nov 12, 2023, at 12:24 AM, KristenKrolick ***@***.***> wrote: @haowulab <https://github.com/haowulab> Hmm I am in the middle of studying your examples (and multiple regression in general), but I do not understand: Q1) The 'fit = lm()' function in R seems to be a little different from your 'DMLfit = DMLfit.multifactor function'-- as in the dependent variable (aka the y in y = B0 + B1x) usually goes to the left of the tilde. In your function are we assuming that to the left of the tilde would be the Bsobj, the methylation values?

Yes this is correct.

Q2) I am not understanding how we are testing the difference in methylation between preop and postop (after correcting for blood cell type fractions) from this formula? -- I suppose it's okay if I don't quite understand this one, as long as you and our Statistician confirm this for me.

Yes. consult your statistician.

Q3) Okay so if I use the DMLfit = DMLfit.multiFactor(chr1_bismarkBSseq, design, formula = ~Treatment + pair + response + Mono_3 + NK_3 + B_TOTAL + CD4_TOTAL + CD8_TOTAL + GRANULOCYTES_TOTAL + (Treatment*response)). I get the error, Error in makeContrast(DMLfit, term) : Some term(s) to be tested can't be found in the formula . How would I fix this?

hmmm, this error should not appear in DMLfit.multiFactor function. It should be in DMLtest.multiFactor function, because you are trying to test something. What did you use for that? To choose the proper term to test is a bit complicated, especially in your design. Please read the manual, and consult your statistician.

For my blood fraction in silico estimation: I used the Reference provided by Salas et al. 2022 (https://www.nature.com/articles/s41467-021-27864-7), actual Ref available:(https://github.com/immunomethylomics/FlowSorted.BloodExtended.EPIC/tree/main). I converted this reference from the Illumina Array sites into UCSC chromosomal coordinates using the manifest provided by Illumina. I used Epidish to perform the deconvolution. All of this is outlined and the code I used here: https://github.com/KristenKrolick/Chidambaran_Lab_Manuscript/blob/main/n)%20Epidish p.s. We had samples that we had flowcyto data for as well as whole-genome meth seq, so I was able to perform a Pearson correlation, and (at least for those cells whose same cell markers were used in our flow data vs. salas ref cell markers he used to make his Ref), they correlated! So, I really liked his Reference.

Okay thanks. I’m not sure if using epidish will be good for BS-seq. What’s the correlation?

…

— Reply to this email directly, view it on GitHub <#41 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AD7H4N7OSVYKSJ2N5H3XTZLYD6RDJAVCNFSM6AAAAAA6472JXKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMBWHA2TONRQHE>. You are receiving this because you were mentioned.

[haowulab]

You can certainly filter the sites by coverage, but it’s not necessary. DSS will consider the coverage in the model. If the coverage is very shallow, it’s very unlikely that the site will be called as DML. This said, it doesn’t hurt to filter.

…

On Nov 12, 2023, at 12:36 AM, KristenKrolick ***@***.***> wrote: Oh! And Q4) Your DSS manual outlines reasoning for why one would not need to subset their whole-genome bisulfite seq data by minimum coverage threshold: does the same apply for the general hypothesis testing design that I am using?? ...I would think I would want to subset my data by at least mincov 4 or so, just to qualify that I am getting an accurate methylation call since my seq. data is from whole blood (mixed cells)... Thank you! — Reply to this email directly, view it on GitHub <#41 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AD7H4NYJPWGOJBPVKEJFIL3YD6SPRAVCNFSM6AAAAAA6472JXKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMBWHA3DAMZWGQ>. You are receiving this because you were mentioned.

[KristenKrolick]

@haowulab
In response to Q3-- sorry it was the same question from my comment 4 days ago, and I did not write it out fully again:

So the design matrix is the same as above (4 days ago comment).

Here was my code:

DMLfit = DMLfit.multiFactor(chr1_bismarkBSseq, design, formula = ~Treatment + pair + response + Mono_3 + NK_3 + B_TOTAL + CD4_TOTAL + CD8_TOTAL + GRANULOCYTES_TOTAL + (Treatment*response))
dmlTest = DMLtest.multiFactor(DMLfit, term="(Treatment*response)")
#Error in makeContrast(DMLfit, term) : Some term(s) to be tested can't be found in the formula.

treatment*response was clearly in the formula? Should I send all the beginning parts of my code again in me making the design dataframe?