# How to create a DESeqDataSet and define experiment design before varianceStabilizingTransformation?

I have an RNAseq count matrix consisting of 2 groups (high, low) with 6 timepoints per group (T1,T2,...,T6) and 3 replicates per timepoint (rep1, rep2, rep3). So a 2-factor design with 36 samples in total (2x6x3=36).

I want to perform a varianceStabilizingTransformation (with blind=FALSE) of the matrix to later perform rhythm analysis (over the timepoints, but separately for each group). Now I wonder how to create the DESeqDataSet that the varianceStabilizingTransformation command needs as input and how to appropriately define the design formula associated with the DESeqDataSet.

Here a scheme of my count table (the numbers are random), which I would save in a .txt or .csv (named "matrix") and import with matrix <- read.table(FILE.txt, header=TRUE, row.names=1):

ID            high_T1_rep1   high_T1_rep2   high_T1_rep3   high_T2_rep1 ...
transcript1   1.23           1.45           1.67           1.89
transcript2   5.32           5.54           4.76           5.98
transcript3   3.22           3.44           3.66           3.88
transcript4   7.33           7.55           7.77           7.99
...


For the design I would now create a table (also txt. or .csv) as follows and import it with designdata <- read.table(FILE.txt, header=TRUE, row.names=1).

sample         group   timepoint
high_T1_rep1   high    T1
high_T1_rep2   high    T1
high_T1_rep3   high    T1
high_T2_rep1   high    T2
...            ...     ...
low_T5_rep3    low     T5
low_T6_rep1    low     T6
low_T6_rep2    low     T6
low_T6_rep3    low     T6


I went through several DESeq2 manuals, but still do not understand how to find the right definition for the design=... term when creating the DESeqDataSet. The varianceStabilizingTransformation seems to use the design information (if blind=FALSE), and thus I want to make sure, I provide it correctly.

From what I found out so far, my command for the DESeqDataSet would be something like this:

dds <- DESeqDataSetFromMatrix(countData=matrix, colData=designdata, design=???)


As options for the design I thought of e.g. "~ group + timepoint" or "~ group + timepoint + group:timepoint", but I have no idea how one or the other option would affect the following varianceStabilizingTransformation.

So while I have at least some rough idea what to do, my big concern is handling the data incorrectly and thereby distorting it without even noticing. Any help you could give me is greatly appreciated. Thanks a lot in advance!

• Are your real counts...counts? Are they all whole numbers? You realize that the vst values won't be used at all in DE calculations? And what design you want depends on exactly what questions you want answered. – swbarnes2 Jun 6 at 23:38
• Thanks for your thoughts. These are TMM-counts (trimmed mean of M-values) I got from the salmon mapper and they have decimals. I also have raw counts (whole numbers), but those do not take differences in sequencing depth between samples into account, so I'm hesitant to use them. I want use the vst values for a different analysis (R-package RAIN) and performing this transformation beforehand was recommended to me by a specialist on this method. In the analysis, rhythmicity is determined from the changes over the time points (T1, T2,...). This is done separately for the 2 groups (high, low). – colormashed Jun 10 at 16:03
• The DESeq documentation is quite clear. It wants raw counts. Library size correction can be done on those in DESeq. – swbarnes2 Jun 10 at 16:51
• Thanks, then I will use the raw counts. Does DESeq automatically determine/consider sequencing depth differences between samples, or does info need to be provided? (at which step?) – colormashed 20 hours ago
• Have you read the vignette? Of course it will do a library size correction, or you can give it numbers to use, if you think your method is better. – swbarnes2 12 hours ago