ahmohamed / lipidr

Data Mining and Analysis of Lipidomics datasets in R

https://www.lipidr.org/

Other

27 stars 13 forks source link

Incorporating the multi-level experiments from Limma into Lipidr #42

Closed Ivan-vechetti closed 1 year ago

Ivan-vechetti commented 1 year ago

Hi, I am using lipir to identify differentially expressed lipids in my experiment design, and I have difficulty incorporating some of the analysis used in limma into the lipir. I have an experimental design that goes like that: Group (control, aerobic exercise, and resistance exercise), time-course (pre, 0 min, 30 min), and 8 subjects (all subjects went through all the procedures, paired design). I ran all the initial parts of the package, but the statistical part is giving me trouble.

Design

f <- paste(d_normalized$Group, d_normalized$Time, sep = ".") design <- model.matrix(~0+f, data = colData(d_normalized)) colnames(design)<-levels(f)

Then estimate the correlation between measurements made on the same subject

mat <- assay(d_normalized, "Area") (thanks, Ahmed to show me how to extract the matrix) corfit<-duplicateCorrelation(mat ,design,block=d_normalized$Subject) corfit$consensus

Then this inter-subject correlation is input into the linear model fit:

fit<-lmFit(mat,design,block=d_normalized$Subject,correlation=corfit$consensus).

According to limma package, now I would have to make contrasts and then compute them and moderated t-tests using contrasts.fit and ebayes. How would I perform those steps after I ran the lipir analysis:

de_results= de_design(data=d_normalized, design = design, ALL contrasts, measure="Area")

Thanks in advance

ahmohamed commented 1 year ago

Thanks for creating an issue for your question. The contrasts needed will depend on the biological question you want to ask. Would you mind elaborating more on what you want to investigate, in plain English?

Ivan-vechetti commented 1 year ago

Hi Ahmed, thanks for your answer. I am sorry I wasn't clear before. The contrast, I know how to set up; my question is regarding the contrast.fit and ebayes. Below you will see that I performed all the steps for the analysis, including the addition of the correlation between measurements for paired samples. After that, I made the contrasts, but I need to incorporate the correlation into the results, like in limma section 9.7 Multi-level experiments. This part I am having problems. Could you help me with that?

Complex design with paired samples

condition <- factor(paste(d_normalized$Group, d_normalized$Time, sep = ".")) design <- model.matrix(~0+condition, data = colData(d_normalized)) colnames(design)<-levels(condition)

Then estimate the correlation between measurements made on the same subject

mat <- assay(d_normalized, "Area") corfit<-duplicateCorrelation(mat ,design,block=d_normalized$Subject) corfit$consensus

Then this inter-subject correlation is input into the linear model fit:

fit<-lmFit(mat,design,block=d_normalized$Subject,correlation=corfit$consensus)

design de_results= de_design( data=d_normalized, design = design, PrevsPostForControl= Control.0min-Control.Pre, Prevs30ForControl= Control.30min-Control.Pre, PrevsPostForAerobic= Aerobic.0min-Aerobic.Pre, Prevs30ForAerobic= Aerobic.30min-Aerobic.Pre, PrevsPostForResistance= Resistance.0min-Resistance.Pre, Prevs30ForResistance= Resistance.30min-Resistance.Pre, PostForControlandAerobic=(Aerobic.0min-Aerobic.Pre)-(Control.0min-Control.Pre), Min30ForControlandAerobic=(Aerobic.30min-Aerobic.Pre)-(Control.30min-Control.Pre), PostForControlandResistance=(Resistance.0min-Resistance.Pre)-(Control.0min-Control.Pre), Min30ForControlandResistance=(Resistance.30min-Resistance.Pre)-(Control.30min-Control.Pre), PostForAerobicandResistance=(Resistance.0min-Resistance.Pre)-(Aerobic.0min-Aerobic.Pre), Min30ForAerobicandResistance=(Resistance.30min-Resistance.Pre)-(Aerobic.30min-Aerobic.Pre), measure="Area")

Based on Limma

Then compute these contrasts and moderated t-tests:

fit2 <- contrasts.fit(fit, de_results) fit2 <- eBayes(fit2) topTable(fit2, coef="____")

error: Error in contrasts.fit(fit, de_results) : contrasts must be a numeric matrix

ahmohamed commented 1 year ago

I apologize for the late response. It's been a busy couple of months.

Basically de_design wraps the logic for lmFit and ebayes, i.e. it already runs this code. As you might have found, the current implementation does not allow passing duplicate correlations, so you'd need to copy the code from de_deign and add the extra bit about the correlation. You can have a look at the code here https://github.com/ahmohamed/lipidr/blob/master/R/de-analysis.R#L108

Following the logic of the code, your code should look like that:

vfit <- lmFit(assay(d_normalized, "Area"), design, block=d_normalized$Subject,correlation=corfit$consensus)
contr.matrix <- limma::makeContrasts(
  PrevsPostForControl= Control.0min-Control.Pre,
  Prevs30ForControl= Control.30min-Control.Pre,
  PrevsPostForAerobic= Aerobic.0min-Aerobic.Pre,
  Prevs30ForAerobic= Aerobic.30min-Aerobic.Pre,
  PrevsPostForResistance= Resistance.0min-Resistance.Pre,
  Prevs30ForResistance= Resistance.30min-Resistance.Pre,
  PostForControlandAerobic=(Aerobic.0min-Aerobic.Pre)-(Control.0min-Control.Pre),
  Min30ForControlandAerobic=(Aerobic.30min-Aerobic.Pre)-(Control.30min-Control.Pre),
  PostForControlandResistance=(Resistance.0min-Resistance.Pre)-(Control.0min-Control.Pre),
  Min30ForControlandResistance=(Resistance.30min-Resistance.Pre)-(Control.30min-Control.Pre),
  PostForAerobicandResistance=(Resistance.0min-Resistance.Pre)-(Aerobic.0min-Aerobic.Pre),
  Min30ForAerobicandResistance=(Resistance.30min-Resistance.Pre)-(Aerobic.30min-Aerobic.Pre), 
  levels = colnames(design)
)
vfit <- limma::contrasts.fit(vfit, contrasts = contr.matrix)
coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix))

efit <- eBayes(vfit)
dimname_x <- metadata(d_normalized)$dimnames[[1]]
top <- lapply(
  coef, function(x)
    limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x)
) %>%
  bind_rows(.id = "contrast")

top <- lipidr:::to_df(d_normalized, dim = "row") %>%
  dplyr::select(
    one_of("Molecule", "Class", "total_cl", "total_cs", "istd", dimname_x)
  ) %>%
  lipidr:::.left_join_silent(top)
attr(top, 'measure') <- measure

Let me know if this works.

Ivan-vechetti commented 1 year ago

Hi Ahmed, no problem at all.

I know how busy you probably are. Thanks for the code. I tried to run, and faced an error, would you be able to tell me what it is? And then, how should I go about doing the enrichment part?

vfit <- limma::contrasts.fit(fit, contrasts = contr.matrix)

coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit) dimname_x <- metadata(d_normalized)$dimnames[[1]] top <- lapply(

coef, function(x)

limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x)

) %>%

bind_rows(.id = "contrast") Error in bind_rows(., .id = "contrast") : could not find function "bind_rows"

On Mon, Dec 5, 2022 at 6:06 PM Ahmed Mohamed @.***> wrote:

I apologize for the late response. It's been a busy couple of months.

Basically de_design wraps the logic for lmFit and ebayes, i.e. it already runs this code. As you might have found, the current implementation does not allow passing duplicate correlations, so you'd need to copy the code from de_deign and add the extra bit about the correlation. You can have a look at the code here https://github.com/ahmohamed/lipidr/blob/master/R/de-analysis.R#L108

Following the logic of the code, your code should look like that:

vfit <- lmFit(assay(d_normalized, "Area"), design, block=d_normalized$Subject,correlation=corfit$consensus)contr.matrix <- limma::makeContrasts( PrevsPostForControl= Control.0min-Control.Pre, Prevs30ForControl= Control.30min-Control.Pre, PrevsPostForAerobic= Aerobic.0min-Aerobic.Pre, Prevs30ForAerobic= Aerobic.30min-Aerobic.Pre, PrevsPostForResistance= Resistance.0min-Resistance.Pre, Prevs30ForResistance= Resistance.30min-Resistance.Pre, PostForControlandAerobic=(Aerobic.0min-Aerobic.Pre)-(Control.0min-Control.Pre), Min30ForControlandAerobic=(Aerobic.30min-Aerobic.Pre)-(Control.30min-Control.Pre), PostForControlandResistance=(Resistance.0min-Resistance.Pre)-(Control.0min-Control.Pre), Min30ForControlandResistance=(Resistance.30min-Resistance.Pre)-(Control.30min-Control.Pre), PostForAerobicandResistance=(Resistance.0min-Resistance.Pre)-(Aerobic.0min-Aerobic.Pre), Min30ForAerobicandResistance=(Resistance.30min-Resistance.Pre)-(Aerobic.30min-Aerobic.Pre), levels = colnames(design) )vfit <- limma::contrasts.fit(vfit, contrasts = contr.matrix)coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit)dimname_x <- metadata(d_normalized)$dimnames[[1]]top <- lapply( coef, function(x) limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x) ) %>% bind_rows(.id = "contrast") top <- lipidr:::to_df(d_normalized, dim = "row") %>% dplyr::select( one_of("Molecule", "Class", "total_cl", "total_cs", "istd", dimname_x) ) %>% lipidr:::.left_join_silent(top) attr(top, 'measure') <- measure

Let

— Reply to this email directly, view it on GitHub https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1338401575, or unsubscribe https://github.com/notifications/unsubscribe-auth/APD5TZKEW3MWWM4MSDMZIODWLZ7P7ANCNFSM6AAAAAARKDEI6M . You are receiving this because you authored the thread.Message ID: @.***>

ahmohamed commented 1 year ago

Use dplyr::bind_rows and try again. Enrichment takes de_results directly, so it shouldn't be a problem.

On Wed, 7 Dec 2022, 2:03 am Ivan-vechetti, @.***> wrote:

Hi Ahmed, no problem at all.

I know how busy you probably are. Thanks for the code. I tried to run, and faced an error, would you be able to tell me what it is? And then, how should I go about doing the enrichment part?

vfit <- limma::contrasts.fit(fit, contrasts = contr.matrix)

coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit) dimname_x <- metadata(d_normalized)$dimnames[[1]] top <- lapply(

coef, function(x)

limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x)

) %>%

bind_rows(.id = "contrast") Error in bind_rows(., .id = "contrast") : could not find function "bind_rows"

On Mon, Dec 5, 2022 at 6:06 PM Ahmed Mohamed @.***> wrote:

I apologize for the late response. It's been a busy couple of months.

Basically de_design wraps the logic for lmFit and ebayes, i.e. it already runs this code. As you might have found, the current implementation does not allow passing duplicate correlations, so you'd need to copy the code from de_deign and add the extra bit about the correlation. You can have a look at the code here https://github.com/ahmohamed/lipidr/blob/master/R/de-analysis.R#L108

Following the logic of the code, your code should look like that:

vfit <- lmFit(assay(d_normalized, "Area"), design, block=d_normalized$Subject,correlation=corfit$consensus)contr.matrix <- limma::makeContrasts( PrevsPostForControl= Control.0min-Control.Pre, Prevs30ForControl= Control.30min-Control.Pre, PrevsPostForAerobic= Aerobic.0min-Aerobic.Pre, Prevs30ForAerobic= Aerobic.30min-Aerobic.Pre, PrevsPostForResistance= Resistance.0min-Resistance.Pre, Prevs30ForResistance= Resistance.30min-Resistance.Pre,

PostForControlandAerobic=(Aerobic.0min-Aerobic.Pre)-(Control.0min-Control.Pre),

Min30ForControlandAerobic=(Aerobic.30min-Aerobic.Pre)-(Control.30min-Control.Pre),

PostForControlandResistance=(Resistance.0min-Resistance.Pre)-(Control.0min-Control.Pre),

Min30ForControlandResistance=(Resistance.30min-Resistance.Pre)-(Control.30min-Control.Pre),

PostForAerobicandResistance=(Resistance.0min-Resistance.Pre)-(Aerobic.0min-Aerobic.Pre),

Min30ForAerobicandResistance=(Resistance.30min-Resistance.Pre)-(Aerobic.30min-Aerobic.Pre), levels = colnames(design) )vfit <- limma::contrasts.fit(vfit, contrasts = contr.matrix)coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit)dimname_x <- metadata(d_normalized)$dimnames[[1]]top <- lapply( coef, function(x) limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x) ) %>% bind_rows(.id = "contrast") top <- lipidr:::to_df(d_normalized, dim = "row") %>% dplyr::select( one_of("Molecule", "Class", "total_cl", "total_cs", "istd", dimname_x) ) %>% lipidr:::.left_join_silent(top) attr(top, 'measure') <- measure

Let

— Reply to this email directly, view it on GitHub https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1338401575, or unsubscribe < https://github.com/notifications/unsubscribe-auth/APD5TZKEW3MWWM4MSDMZIODWLZ7P7ANCNFSM6AAAAAARKDEI6M

. You are receiving this because you authored the thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339518533, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABQWVVFN2XTM4CSHC3AZ36DWL5IVTANCNFSM6AAAAAARKDEI6M . You are receiving this because you commented.Message ID: @.***>

Ivan-vechetti commented 1 year ago

Thank you so much. It worked, with the exception of:

attr(top, 'measure') <- measure Error: object 'measure' not found

On Tue, Dec 6, 2022 at 9:13 AM Ahmed Mohamed @.***> wrote:

Use dplyr::bind_rows and try again. Enrichment takes de_results directly, so it shouldn't be a problem.

On Wed, 7 Dec 2022, 2:03 am Ivan-vechetti, @.***> wrote:

Hi Ahmed, no problem at all.

I know how busy you probably are. Thanks for the code. I tried to run, and faced an error, would you be able to tell me what it is? And then, how should I go about doing the enrichment part?

vfit <- limma::contrasts.fit(fit, contrasts = contr.matrix)

coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit) dimname_x <- metadata(d_normalized)$dimnames[[1]] top <- lapply(

coef, function(x)

limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x)

) %>%

bind_rows(.id = "contrast") Error in bind_rows(., .id = "contrast") : could not find function "bind_rows"

On Mon, Dec 5, 2022 at 6:06 PM Ahmed Mohamed @.***> wrote:

I apologize for the late response. It's been a busy couple of months.

Basically de_design wraps the logic for lmFit and ebayes, i.e. it already runs this code. As you might have found, the current implementation does not allow passing duplicate correlations, so you'd need to copy the code from de_deign and add the extra bit about the correlation. You can have a look at the code here https://github.com/ahmohamed/lipidr/blob/master/R/de-analysis.R#L108

Following the logic of the code, your code should look like that:

vfit <- lmFit(assay(d_normalized, "Area"), design, block=d_normalized$Subject,correlation=corfit$consensus)contr.matrix <- limma::makeContrasts( PrevsPostForControl= Control.0min-Control.Pre, Prevs30ForControl= Control.30min-Control.Pre, PrevsPostForAerobic= Aerobic.0min-Aerobic.Pre, Prevs30ForAerobic= Aerobic.30min-Aerobic.Pre, PrevsPostForResistance= Resistance.0min-Resistance.Pre, Prevs30ForResistance= Resistance.30min-Resistance.Pre,

PostForControlandAerobic=(Aerobic.0min-Aerobic.Pre)-(Control.0min-Control.Pre),

Min30ForControlandAerobic=(Aerobic.30min-Aerobic.Pre)-(Control.30min-Control.Pre),

PostForControlandResistance=(Resistance.0min-Resistance.Pre)-(Control.0min-Control.Pre),

Min30ForControlandResistance=(Resistance.30min-Resistance.Pre)-(Control.30min-Control.Pre),

PostForAerobicandResistance=(Resistance.0min-Resistance.Pre)-(Aerobic.0min-Aerobic.Pre),

Min30ForAerobicandResistance=(Resistance.30min-Resistance.Pre)-(Aerobic.30min-Aerobic.Pre),

levels = colnames(design) )vfit <- limma::contrasts.fit(vfit, contrasts = contr.matrix)coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit)dimname_x <- metadata(d_normalized)$dimnames[[1]]top <- lapply( coef, function(x) limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x) ) %>% bind_rows(.id = "contrast") top <- lipidr:::to_df(d_normalized, dim = "row") %>% dplyr::select( one_of("Molecule", "Class", "total_cl", "total_cs", "istd", dimname_x) ) %>% lipidr:::.left_join_silent(top) attr(top, 'measure') <- measure

Let

— Reply to this email directly, view it on GitHub <https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1338401575 , or unsubscribe <

https://github.com/notifications/unsubscribe-auth/APD5TZKEW3MWWM4MSDMZIODWLZ7P7ANCNFSM6AAAAAARKDEI6M

. You are receiving this because you authored the thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339518533, or unsubscribe < https://github.com/notifications/unsubscribe-auth/ABQWVVFN2XTM4CSHC3AZ36DWL5IVTANCNFSM6AAAAAARKDEI6M

. You are receiving this because you commented.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339533403, or unsubscribe https://github.com/notifications/unsubscribe-auth/APD5TZMKXK4RDBQG33LSSMDWL5J2PANCNFSM6AAAAAARKDEI6M . You are receiving this because you authored the thread.Message ID: @.***>

ahmohamed commented 1 year ago

attr(top, 'measure') <- "Area"

On Wed, 7 Dec 2022, 2:41 am Ivan-vechetti, @.***> wrote:

Thank you so much. It worked, with the exception of:

attr(top, 'measure') <- measure Error: object 'measure' not found

On Tue, Dec 6, 2022 at 9:13 AM Ahmed Mohamed @.***> wrote:

Use dplyr::bind_rows and try again. Enrichment takes de_results directly, so it shouldn't be a problem.

On Wed, 7 Dec 2022, 2:03 am Ivan-vechetti, @.***> wrote:

Hi Ahmed, no problem at all.

I know how busy you probably are. Thanks for the code. I tried to run, and faced an error, would you be able to tell me what it is? And then, how should I go about doing the enrichment part?

vfit <- limma::contrasts.fit(fit, contrasts = contr.matrix)

coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit) dimname_x <- metadata(d_normalized)$dimnames[[1]] top <- lapply(

coef, function(x)

limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x)

) %>%

bind_rows(.id = "contrast") Error in bind_rows(., .id = "contrast") : could not find function "bind_rows"

On Mon, Dec 5, 2022 at 6:06 PM Ahmed Mohamed @.***> wrote:

I apologize for the late response. It's been a busy couple of months.

Basically de_design wraps the logic for lmFit and ebayes, i.e. it already runs this code. As you might have found, the current implementation does not allow passing duplicate correlations, so you'd need to copy the code from de_deign and add the extra bit about the correlation. You can have a look at the code here https://github.com/ahmohamed/lipidr/blob/master/R/de-analysis.R#L108

Following the logic of the code, your code should look like that:

vfit <- lmFit(assay(d_normalized, "Area"), design, block=d_normalized$Subject,correlation=corfit$consensus)contr.matrix <- limma::makeContrasts( PrevsPostForControl= Control.0min-Control.Pre, Prevs30ForControl= Control.30min-Control.Pre, PrevsPostForAerobic= Aerobic.0min-Aerobic.Pre, Prevs30ForAerobic= Aerobic.30min-Aerobic.Pre, PrevsPostForResistance= Resistance.0min-Resistance.Pre, Prevs30ForResistance= Resistance.30min-Resistance.Pre,

PostForControlandAerobic=(Aerobic.0min-Aerobic.Pre)-(Control.0min-Control.Pre),

Min30ForControlandAerobic=(Aerobic.30min-Aerobic.Pre)-(Control.30min-Control.Pre),

PostForControlandResistance=(Resistance.0min-Resistance.Pre)-(Control.0min-Control.Pre),

Min30ForControlandResistance=(Resistance.30min-Resistance.Pre)-(Control.30min-Control.Pre),

PostForAerobicandResistance=(Resistance.0min-Resistance.Pre)-(Aerobic.0min-Aerobic.Pre),

Min30ForAerobicandResistance=(Resistance.30min-Resistance.Pre)-(Aerobic.30min-Aerobic.Pre),

levels = colnames(design) )vfit <- limma::contrasts.fit(vfit, contrasts = contr.matrix)coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit)dimname_x <- metadata(d_normalized)$dimnames[[1]]top <- lapply( coef, function(x) limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x) ) %>% bind_rows(.id = "contrast") top <- lipidr:::to_df(d_normalized, dim = "row") %>% dplyr::select( one_of("Molecule", "Class", "total_cl", "total_cs", "istd", dimname_x) ) %>% lipidr:::.left_join_silent(top) attr(top, 'measure') <- measure

Let

— Reply to this email directly, view it on GitHub < https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1338401575 , or unsubscribe <

https://github.com/notifications/unsubscribe-auth/APD5TZKEW3MWWM4MSDMZIODWLZ7P7ANCNFSM6AAAAAARKDEI6M

. You are receiving this because you authored the thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub <https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339518533 , or unsubscribe <

https://github.com/notifications/unsubscribe-auth/ABQWVVFN2XTM4CSHC3AZ36DWL5IVTANCNFSM6AAAAAARKDEI6M

. You are receiving this because you commented.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339533403, or unsubscribe < https://github.com/notifications/unsubscribe-auth/APD5TZMKXK4RDBQG33LSSMDWL5J2PANCNFSM6AAAAAARKDEI6M

. You are receiving this because you authored the thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339570922, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABQWVVDEECRJ6Z2LR6BIUJDWL5NBFANCNFSM6AAAAAARKDEI6M . You are receiving this because you commented.Message ID: @.***>

Ivan-vechetti commented 1 year ago

Thanks again, Ahmed.

The last thing, I promise. Since the source code you sent has no de_results, I am a little confused about how to do the enrichment or generate the de_results.

I appreciate all your help and help you have a great week.

On Tue, Dec 6, 2022 at 10:56 AM Ahmed Mohamed @.***> wrote:

attr(top, 'measure') <- "Area"

On Wed, 7 Dec 2022, 2:41 am Ivan-vechetti, @.***> wrote:

Thank you so much. It worked, with the exception of:

attr(top, 'measure') <- measure Error: object 'measure' not found

On Tue, Dec 6, 2022 at 9:13 AM Ahmed Mohamed @.***> wrote:

Use dplyr::bind_rows and try again. Enrichment takes de_results directly, so it shouldn't be a problem.

On Wed, 7 Dec 2022, 2:03 am Ivan-vechetti, @.***> wrote:

Hi Ahmed, no problem at all.

I know how busy you probably are. Thanks for the code. I tried to run, and faced an error, would you be able to tell me what it is? And then, how should I go about doing the enrichment part?

vfit <- limma::contrasts.fit(fit, contrasts = contr.matrix)

coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit) dimname_x <- metadata(d_normalized)$dimnames[[1]] top <- lapply(

coef, function(x)

limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x)

) %>%

bind_rows(.id = "contrast") Error in bind_rows(., .id = "contrast") : could not find function "bind_rows"

On Mon, Dec 5, 2022 at 6:06 PM Ahmed Mohamed @.***> wrote:

I apologize for the late response. It's been a busy couple of months.

Basically de_design wraps the logic for lmFit and ebayes, i.e. it already runs this code. As you might have found, the current implementation does not allow passing duplicate correlations, so you'd need to copy the code from de_deign and add the extra bit about the correlation. You can have a look at the code here

https://github.com/ahmohamed/lipidr/blob/master/R/de-analysis.R#L108

Following the logic of the code, your code should look like that:

vfit <- lmFit(assay(d_normalized, "Area"), design, block=d_normalized$Subject,correlation=corfit$consensus)contr.matrix <- limma::makeContrasts( PrevsPostForControl= Control.0min-Control.Pre, Prevs30ForControl= Control.30min-Control.Pre, PrevsPostForAerobic= Aerobic.0min-Aerobic.Pre, Prevs30ForAerobic= Aerobic.30min-Aerobic.Pre, PrevsPostForResistance= Resistance.0min-Resistance.Pre, Prevs30ForResistance= Resistance.30min-Resistance.Pre,

PostForControlandAerobic=(Aerobic.0min-Aerobic.Pre)-(Control.0min-Control.Pre),

Min30ForControlandAerobic=(Aerobic.30min-Aerobic.Pre)-(Control.30min-Control.Pre),

PostForControlandResistance=(Resistance.0min-Resistance.Pre)-(Control.0min-Control.Pre),

Min30ForControlandResistance=(Resistance.30min-Resistance.Pre)-(Control.30min-Control.Pre),

PostForAerobicandResistance=(Resistance.0min-Resistance.Pre)-(Aerobic.0min-Aerobic.Pre),

Min30ForAerobicandResistance=(Resistance.30min-Resistance.Pre)-(Aerobic.30min-Aerobic.Pre),

levels = colnames(design) )vfit <- limma::contrasts.fit(vfit, contrasts = contr.matrix)coef <- setNames(seq_len(ncol(contr.matrix)), colnames(contr.matrix)) efit <- eBayes(vfit)dimname_x <- metadata(d_normalized)$dimnames[[1]]top <- lapply( coef, function(x) limma::topTable(efit, number = Inf, coef = x) %>% tibble::rownames_to_column(dimname_x) ) %>% bind_rows(.id = "contrast") top <- lipidr:::to_df(d_normalized, dim = "row") %>% dplyr::select( one_of("Molecule", "Class", "total_cl", "total_cs", "istd", dimname_x) ) %>% lipidr:::.left_join_silent(top) attr(top, 'measure') <- measure

Let

— Reply to this email directly, view it on GitHub < https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1338401575 , or unsubscribe <

https://github.com/notifications/unsubscribe-auth/APD5TZKEW3MWWM4MSDMZIODWLZ7P7ANCNFSM6AAAAAARKDEI6M

. You are receiving this because you authored the thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub < https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339518533 , or unsubscribe <

https://github.com/notifications/unsubscribe-auth/ABQWVVFN2XTM4CSHC3AZ36DWL5IVTANCNFSM6AAAAAARKDEI6M

. You are receiving this because you commented.Message ID: @.***>

— Reply to this email directly, view it on GitHub <https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339533403 , or unsubscribe <

https://github.com/notifications/unsubscribe-auth/APD5TZMKXK4RDBQG33LSSMDWL5J2PANCNFSM6AAAAAARKDEI6M

. You are receiving this because you authored the thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339570922, or unsubscribe < https://github.com/notifications/unsubscribe-auth/ABQWVVDEECRJ6Z2LR6BIUJDWL5NBFANCNFSM6AAAAAARKDEI6M

. You are receiving this because you commented.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/ahmohamed/lipidr/issues/42#issuecomment-1339676283, or unsubscribe https://github.com/notifications/unsubscribe-auth/APD5TZOAP6DPJSLI2TY7X4DWL5V3LANCNFSM6AAAAAARKDEI6M . You are receiving this because you authored the thread.Message ID: @.***>

ahmohamed commented 1 year ago

de_results <- top