Question

Setting up contrasts with 'limma', patient data, small number of repeats

0

Entering edit mode

uridavid.akavia • 0

@uridavidakavia-22012

Last seen 5.1 years ago

Hello,

I'm trying to run limma-vooom on a Nanostring experiment. The experiment design is multiple patients, each assessed before and after treatment. So the targets structure look something like this (only more lines):

  SampleID patient Sex           Status
1    C15-1     C15   F Before Treatment
2    C15-2     C15   F    After Relapse
3    C30-1     C30   M Before Treatment
4    C30-2     C30   M    After Relapse
5    C40-1     C40   M Before Treatment
6    C40-2     C40   M    After Relapse

I've set up my design matrix the following manner

mm <- model.matrix(~0 + patient + Status, sampleDetails)

Then I've processed the data using

y <- voom(digital, design = mm, plot=T)

And I plan to run the rest of the limma pipeline (lmFit, topTable) to find differentially expressed genes by the treatment.

Three questions 0) Are there any mistakes I've made so far? 1) Do I need to set up a contrast matrix to focus on treatment? I don't think so, since it is already in my model 2) I am interested in seeing if there are treatment effects specific to Sex, but I'm not sure how to tease them out besides doing a manual comparison. I think I can focus on it using a contrast matrix, but I'm not sure how to do so. How should I go about doing that?

Thank you very much.

Yours,

Uri David Akavia

=========

sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] stringr_1.4.0 readxl_1.3.1  readr_1.3.1   dplyr_0.8.3   edgeR_3.24.3  limma_3.38.3  ggplot2_3.2.1

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.2         pillar_1.4.2       compiler_3.5.1     cellranger_1.1.0   BiocManager_1.30.4
 [6] tools_3.5.1        zeallot_0.1.0      tibble_2.1.3       gtable_0.3.0       lattice_0.20-38   
[11] pkgconfig_2.0.3    rlang_0.4.0        cli_1.1.0          rstudioapi_0.10    yaml_2.2.0        
[16] xfun_0.9           withr_2.1.2        knitr_1.25         vctrs_0.2.0        hms_0.5.1         
[21] locfit_1.5-9.1     grid_3.5.1         tidyselect_0.2.5   glue_1.3.1         R6_2.4.0          
[26] fansi_0.4.0        purrr_0.3.2        magrittr_1.5       scales_1.0.0       backports_1.1.4   
[31] assertthat_0.2.1   colorspace_1.4-1   utf8_1.1.4         stringi_1.4.3      lazyeval_0.2.2    
[36] munsell_0.5.0      crayon_1.3.4

limma nanostring contrasts • 1.2k views

ADD COMMENT • link updated 5.1 years ago by Gordon Smyth 51k • written 5.1 years ago by uridavid.akavia • 0

score 0 · Answer 1 · 2019-09-27

0

Entering edit mode

Gordon Smyth 51k

@gordon-smyth

Last seen 4 hours ago

WEHI, Melbourne, Australia

If you are happy for the treatment effect to be the same for F and M, then the design matrix you have is fine and you do not need to form a contrast.

If you want separate treatment effects for F and M, I would do it like this:

design.patients <- model.matrix(~patient)
F.Status <- Sex=="F" & Status=="After relapse"
M.Status <- Sex=="M" & Status=="After relapse"
design <- cbind(design.patients, F.Status=F.Status, M.Status=M.Status)

ADD COMMENT • link 5.1 years ago Gordon Smyth 51k

0

Entering edit mode

Thank you Gordon!

What should I do if I don't know if there are any gender specific treatment effects? Should I run the analysis twice, once with my original design and once with your design, and then compare the coefficients for F.Status and M.Status? Can you suggest code on how to do this?

I think that if I try to combine all the design matrices, using something like

design <- cbind(mm, F.Status=F.Status, M.Status=M.Status)

it might work, but I'm not sure. My thinking is that F.Status and M.Status will be zeroed out UNLESS there are gender specific effects. Have I got this right?

Thank you!

Uri David

ADD REPLY • link 5.1 years ago uridavid.akavia • 0

0

Entering edit mode

I already gave you the correct design matrix. To test whether the treatments effects are the same for M and F, simply use contrasts.fit to compare the two coefficients.

Your combination design matrix is over-parametrized and will give uninterpretable results.

ADD REPLY • link 5.1 years ago Gordon Smyth 51k

0

Entering edit mode

Thank you for your help. I'm using makeContrasts, and my code looks like this

design.patients <- model.matrix(~patient, sampleDetails)
F.Status <- sampleDetails$Sex=="F" & sampleDetails$Status=="After_Relapse"
M.Status <- sampleDetails$Sex=="M" & sampleDetails$Status=="After_Relapse"
designMatrix <- cbind(design.patients, F.Status=F.Status, M.Status=M.Status)
y2 <- voom(digital, designMatrix, plot=T)
fit2 <- lmFit(y2, designMatrix)
contrastMatrix <- makeContrasts("(M.Status + F.Status)/2", "M.Status - F.Status", levels=colnames(coef(fit2)))
fitWithContrasts <- contrasts.fit(fit2, contrastMatrix)
fitWithContrasts <- eBayes(fitWithContrasts)

Have I got this right?

Second, I also have other factors, such as cancer type and so on, so my sample Matrix looks like head(sampleDetails)

  SampleID patient nanostring_run Sex GCB tDLBCL ASCT      Age prior_lines Immune_stimulation           Status
1    C15-1     C15              1   F   1      0    0 73.99589           1              FALSE Before_Treatment
2    C15-2     C15              1   F   1      0    0 73.99589           1              FALSE    After_Relapse
3    C30-1     C30              1   M   1      1    1 29.17180           6              FALSE Before_Treatment
4    C30-2     C30              1   M   1      1    1 29.17180           6              FALSE    After_Relapse
5    C40-1     C40              1   M   0      1    0 71.58658           1              FALSE Before_Treatment
6    C40-2     C40              1   M   0      1    0 71.58658           1              FALSE    After_Relapse

I'm running chi-Sq test between the various properties to see that they're not correlated. Assuming I want to check if there are any effects specific to GCB, my understanding is that I should add it to the design matrix (if it is not correlated to the others) in this way

design.patients <- model.matrix(~patient, sampleDetails)
F.Status <- sampleDetails$Sex=="F" & sampleDetails$Status=="After_Relapse"
M.Status <- sampleDetails$Sex=="M" & sampleDetails$Status=="After_Relapse"
GCB.Status <- sampleDetails$GCB== 1 & sampleDetails$Status=="After_Relapse"
designMatrix <- cbind(design.patients, F.Status=F.Status, M.Status=M.Status, GCB.Status)

And then to the contrast matrix using

contrastMatrix <- makeContrasts("(M.Status + F.Status)/2", "M.Status - F.Status", "GCB.Status", levels=colnames(coef(fit2)))

Is that correct?

Thank you very much, and thank you for writing limma and related packages in the first place!

Yours,

Uri David

ADD REPLY • link 5.1 years ago uridavid.akavia • 0