Question

DESeq2 polysome profiling analysis with 2 interaction terms

1

Entering edit mode

erica ▴ 10

@erica-13193

Last seen 7.6 years ago

Hello,

I am trying to analyse a polysome profiling experiment using DESeq2 and I got stuck deciding on which design would be best to use.

The aims of my experiment are:

1- to investigate the effect of a drug (L-leucine) on translation

2-to look at mRNA translation in patients vs disease samples

I want to test for the ratio of the ratios using a LRT and two interaction terms and these are the groups I want to test:

type: total RNA (RNA), Ribosome-bound RNA (RBR)

treatment: control (D-leucine), treated (L-leucine)

disease: Healthy, 5q syndrome

I worked out a design based on previous posts and DESeq2 vignettes but I am not sure whether it's correct and whether I am extracting the results right.

Here it is:

Design(dds) < - ~type + treatment + disease + type:treatment + disease:treatment

dds <- DESeq(dds, test="LRT", reduced= ~ type + treatment + disease)

To get effect of treatment on ratio between RBR/RNA I do:

Results(dds, contrast=c(0,1,0,0,.5,0))

Then to find the interaction effect of condition:treatment effect across disease:

Results(dds, contrast=c(0 ,0 ,0 ,0 ,1 ,0))

This was based on this previous post https://support.bioconductor.org/p/76966/ but I am wondering:

1- to get the disease effect on the ratio between RBR/RNA should I maybe have this instead?

Design(dds) <- ~condition + treatment + disease + condition:treatment + condition:disease

Perhaps it's simply a different thing than disease:treatment but I am not sure which one to use in this case.

2- to get the disease effect on the ratio between RBR/RNA I used Results(dds, contrast=c(0 ,0 ,0 ,0 ,1 ,0)) but, since I cannot assume that the treatment effect on type is the same across disease, should I have

Results(dds, contrast=c(0 ,1 ,0 ,0 ,0 ,.5)) instead?

It is a quite complex design so any help would be very much appreciated.

All the best,
Erica

deseq2 sequencing polysome profiling • 2.2k views

ADD COMMENT • link updated 7.6 years ago by Michael Love 43k • written 7.6 years ago by erica ▴ 10

0

Entering edit mode

In case it's needed

sessionInfo()

R version 3.2.1 (2015-06-18)
Platform: x86_64-unknown-linux-gnu (64-bit)
Running under: CentOS release 6.9 (Final)

locale:
[1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C LC_TIME=en_GB.UTF-8
[4] LC_COLLATE=en_GB.UTF-8 LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
[7] LC_PAPER=en_GB.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C

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

other attached packages:
[1] DESeq2_1.8.2 RcppArmadillo_0.6.500.4.0 Rcpp_0.12.10
[4] GenomicRanges_1.20.8 GenomeInfoDb_1.4.3 IRanges_2.2.9
[7] S4Vectors_0.6.6 BiocGenerics_0.14.0 tidyr_0.6.3
[10] pheatmap_1.0.8 dplyr_0.5.0

loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-2 futile.logger_1.4.1 plyr_1.8.4
[4] XVector_0.8.0 futile.options_1.0.0 tools_3.2.1
[7] rpart_4.1-10 RSQLite_1.0.0 annotate_1.46.1
[10] tibble_1.3.1 gtable_0.2.0 lattice_0.20-33
[13] rlang_0.1.1 DBI_0.6-1 gridExtra_2.2.1
[16] genefilter_1.50.0 cluster_2.0.3 locfit_1.5-9.1
[19] nnet_7.3-12 grid_3.2.1 Biobase_2.28.0
[22] R6_2.2.1 AnnotationDbi_1.30.1 XML_3.98-1.3
[25] survival_2.38-3 BiocParallel_1.2.22 foreign_0.8-66
[28] latticeExtra_0.6-26 Formula_1.2-1 geneplotter_1.46.0
[31] ggplot2_2.2.1 lambda.r_1.1.7 magrittr_1.5
[34] scales_0.4.1 Hmisc_3.17-1 splines_3.2.1
[37] assertthat_0.2.0 xtable_1.8-2 colorspace_1.3-2
[40] acepack_1.3-3.3 lazyeval_0.2.0 munsell_0.4.3

ADD REPLY • link 7.6 years ago erica ▴ 10

score 0 · Answer 1 · 2017-06-07

0

Entering edit mode

Michael Love 43k

@mikelove

Last seen 5 days ago

United States

If you have two treatments, can you avoid having to report some 'average' effect of disease and just report the two disease effects specific for each treatment? There may be no such thing as an average disease effect after all for some genes, if the treatments do induce differences.

ADD COMMENT • link 7.6 years ago Michael Love 43k

0

Entering edit mode

Hello Michael,

thank you for your reply.

The reason I am interested in the disease effect is because the syndrome is characterised by a defect in mRNA translation (from previous evidence) and therefore I expect the see an "average" effect of disease across treatments. It might also be that the treatment increases translation in the patients to a level that is comparable to translation in untreated controls but I don't know that.

What would be the best way of investigating this in your opinion? Is any of the designs I proposed above appropriate?

Many thanks for your help.

All the best,
Erica

ADD REPLY • link 7.6 years ago erica ▴ 10

0

Entering edit mode

If 'type' gives the ribosomal vs total effect, and you want to see differences in the disease effect across treatment, you will need an interaction between all three. I think you want ~treatment*disease*type. Can you reply back with the resultsNames(dds) you get. The last interaction term will be the test of a difference across treatment for the disease effect on Ribo/total ratio. As with other interaction models, you add main effects and interactions to achieve the comparisons for specific groups, or use 0.5 to get an 'average' effect across e.g. treatment.

ADD REPLY • link 7.6 years ago Michael Love 43k

0

Entering edit mode

If I use

Design(dds) <- ~type + treatment + disease + type:treatment + disease:treatment

dds <- DESeq(dds, test="LRT", reduced= ~ type + treatment + disease)

This gives me 6 contrasts:

> resultsNames(dds)

[1] "Intercept"

[2] "type_RBR_vs_RNA"

[3] "treatment_L_vs_D"

[4] "disease_5q..syndrome_vs_Healthy.Control"

[5] "typeRBR.treatmentL"

[6] "treatmentL.disease5q..syndrome"

Is this correct or should I have Design(dds) <- ~type + treatment + disease + type:treatment + type:disease ? If I know that the disease effect is the same across treatment, should I simply test for average type effect across disease (with 0.5 for type:disease)?

Thank you very much Michael this is really helpful for me!

Best,
Erica

ADD REPLY • link 7.6 years ago erica ▴ 10

0

Entering edit mode

You need the interaction between all three factors, try the code I suggested above.

ADD REPLY • link 7.6 years ago Michael Love 43k

0

Entering edit mode

I am sorry Michael, I misinterpreted your reply above.

Should the design be:

dds <-DESeqDataSetFromMatrix(countData = countdata_all, colData=colData, design = ~type*treatment*disease)?

Can I still do a LRT in this case? I am not entirely sure what the reduced formula would be.

I tried it and the dds <- DESeq(dds) gives me this error

Error in designAndArgChecker(object, betaPrior) :
interactions higher than 1st order with log fold change shrinkage
has not been implemented. Either combine factors into individual groups using
factor(paste(...)) and use a design of ~group, or set betaPrior=FALSE

I am not sure that's what you meant. I also tried:

dds <-DESeqDataSetFromMatrix(countData = countdata_all, colData=colData, design = ~type + treatment + disease+ type*treatment*disease)
dds <- DESeq(dds, "LRT", reduced = (~type + treatment + disease))

and resultsNames(dds) gives me:

[1] "Intercept" "type_RBR_vs_RNA"
[3] "treatment_L_vs_D" "disease_5q..syndrome_vs_Healthy.Control"
[5] "typeRBR.treatmentL" "typeRBR.disease5q..syndrome"
[7] "treatmentL.disease5q..syndrome" "typeRBR.treatmentL.disease5q..syndrome"

Was this what you were asking? Which one is the right contrast to call results?

Many thanks for your patience. I am very new to this and this is a quite complicated design.

Best,
Erica

ADD REPLY • link 7.6 years ago erica ▴ 10

0

Entering edit mode

You should update to the latest version of R and Bioconductor. There have been 4 releases of DESeq2 since April 2015. The code would work with the current release, and that way you'll benefit from the past two years of bug fixes and software improvements.

ADD REPLY • link 7.6 years ago Michael Love 43k

0

Entering edit mode

Hello Michael,

I have updated and it now runs fine. I did

dds <-DESeqDataSetFromMatrix(countData = countdata_all, colData=colData, design = ~type*treatment*disease)

dds <- DESeq(dds)

resultsNames(dds)

> resultsNames(dds)
[1] "Intercept" "type_RBR_vs_RNA"
[3] "treatment_L_vs_D" "disease_5q..syndrome_vs_Healthy.Control"
[5] "typeRBR.treatmentL" "typeRBR.disease5q..syndrome"
[7] "treatmentL.disease5q..syndrome" "typeRBR.treatmentL.disease5q..syndrome"

As I said, is it possible at all to run a LRT with this design? If so, what would be the reduced formula?

With these resultsNames what is the right contrast for the two questions I am asking?

Thanks

Best,
Erica

ADD REPLY • link 7.6 years ago erica ▴ 10

1

Entering edit mode

You can remove whichever term you want to test for a LRT. I tend to use the Wald tests in DESeq2 as they give similar results and do not require fitting the reduced formula. You can obtain the test of differences across treatment by specifying the interaction term with all 3 factors. For the average disease effect you would put a 1 for typeRBR.disease5q..syndrome and 0.5 for the interaction term between all 3 factors and so on.