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.

&gt; sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows &gt;= 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
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 &lt;- data.frame(treatment = c("a", "a", "b", "b"), strain = c("x", "y", "x", "y"))
design_combined &lt;- 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 &lt;- paste0(meta.data$treatment, "_", meta.data$strain)
# refactor to change the reference level
meta.data$combined &lt;- factor(meta.data$combined, levels = "bx", "by", "ax", "ay")
design_combined &lt;- 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