# How do I select a proper design matrix to use in DRIMSeq?

I am trying to use DRIMseq for DTU with 2 treatments on two different strains of an animal model. How do I select a proper design matrix to use in DRIMSeq ? How I know that one is the correct one to use? And how do I see which transcripts are differentially used among conditions after the pipeline was run? I suspect it is in the plotProportions function.

Also, where do I see the logFoldChange between conditions?

Many thanks in advance as it would help me a lot in my work.

> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

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

other attached packages:
[1] sva_3.30.1          BiocParallel_1.16.6 genefilter_1.64.0   mgcv_1.8-28
[5] nlme_3.1-137        edgeR_3.24.3        limma_3.38.3

loaded via a namespace (and not attached):
[1] Rcpp_1.0.1           compiler_3.5.3       pillar_1.4.2         bitops_1.0-6
[5] tools_3.5.3          digest_0.6.20        zeallot_0.1.0        bit_1.1-14
[9] annotate_1.60.1      RSQLite_2.1.1        memoise_1.1.0        tibble_2.1.3
[13] lattice_0.20-38      pkgconfig_2.0.2      rlang_0.4.0          Matrix_1.2-17
[17] DBI_1.0.0            rstudioapi_0.10      yaml_2.2.0           parallel_3.5.3
[21] IRanges_2.16.0       S4Vectors_0.20.1     vctrs_0.2.0          locfit_1.5-9.1
[25] stats4_3.5.3         bit64_0.9-7          grid_3.5.3           Biobase_2.42.0
[29] AnnotationDbi_1.44.0 survival_2.44-1.1    XML_3.98-1.20        blob_1.2.0
[33] matrixStats_0.54.0   backports_1.1.4      splines_3.5.3        BiocGenerics_0.28.0
[37] xtable_1.8-4         RCurl_1.95-4.12      crayon_1.3.4

• Welcome to the site. The design depends on the question and on the conditions of the experiments. It would be easier for us to help you if you provide a table with the information about each sample sequenced. The strain, treatment, sex, batch and other relevant conditions used of each sample. Also it might (or not) be relevant to know the tissue you are studying (or if you are mixing several tissues on the analysis)
– llrs
Jul 26 '19 at 14:03

From your question I understand that you have two treatments in two mouse strains.
If you are interested in general treatment effects independet of the animal strain you would use the both variables in the design formula (similar to controlling for a batch effect).

meta.data <- data.frame(treatment = c("a", "a", "b", "b"), strain = c("x", "y", "x", "y"))
design_combined <- model.matrix(~ treatment + strain, data = meta.data)


If you want to make the comparisons for each strain separately, the easiest way will be to split your data set by strain and do a set up just for treatment. Alternatively, you can also merge the two columns into a new one in case you also want to also compare between the strains.

meta.data$combined <- paste0(meta.data$treatment, "_", meta.data$strain) # refactor to change the reference level meta.data$combined <- factor(meta.data\$combined, levels = "bx", "by", "ax", "ay")
design_combined <- model.matrix(~ combined, data = meta.data)


To get the fold changes, you can extract the proportions for each transcript using

DRIMSeq::proportions(dataSet)


And from the proportions the log2 fold changes can be calculated with ease.

However, when I look at your plot something seems very off with your analysis right now. In all the examples I have seen and also in my own data the diamonds are representing the mean proportion of a transcript in a condition. See this example from the DRIMSeq manual on Bioconductor:

In your case there are several different means per group and if you look at the line to the very right the two treatments seem to be to be right on top of each other.
So I think that treatment and strain are somehow mixed in your setup and this is causing this kind of merged plot. If you use a design as outlined above you should get much better results. Additionally, you could maybe check that the meta data is correct.

Also in case you have not read it yet, the manual is a great resource:
https://bioconductor.org/packages/release/bioc/vignettes/DRIMSeq/inst/doc/DRIMSeq.pdf

• How do i calculate the LogFC from the proportions? It isn't in the manual.
– user4936
Jul 26 '19 at 10:30
• divide the proportions and take the log2. Log2(proportion A / porportion B) = Log2FC of the proportions with the condition B as the reference level.
– PPK
Jul 26 '19 at 10:49
• Thank you very much for your responses, they are helping me a lot. When i do the command DRIMSeq::proportions(dataSet) , I get proportions for each transcript for all samples. Should I sum the proportions in the samples that belong to category A and the same for samples that belong to category B, and then finally do Log2(proportion A / porportion B) ?
– user4936
Aug 9 '19 at 12:52
• Glad to help. Yes, you can calculate the overall proportion from the single proportions but you should calculate the mean proportion for each experimental group (which is the same as the diamonds in the proportion plot) and then calculate log2(mean(propA)/mean(propB)). Otherwise you run into problems when your sample sizes are not equal.
– PPK
Aug 9 '19 at 13:10
• Very clear! Thank you. I'm now computing the DRIMSeq pipeline for every possible combination of the design matrix I could think of (they are 8) and see any give results like the ones you described above. Despite being computacionally expensive do you think its a good idea, in order to select a proper design matrix ?
– user4936
Aug 9 '19 at 13:57