Issues with interpretation of DESeq2 time series analysis results

Question

[Maksim Nesterenko]

Dear colleagues, My research dedicated to transcriptomic analysis of regeneration process in two different sites (Head and Tail) of one invertebrate animal. Experiment design includes RNA-seq data from different time point: 0, 4, 12, 24, 48 and 96 hours. DESeq2 library was used to differential expression analysis:

samps <- read.csv2("Samples_and_conditions.csv")
samps$Site <- factor(samps$Site) 
samps$TPA <- factor(samps$TPA) 
samps$Replicates <- factor(samps$Replicates) 
files <- file.path(samps$X, "quant.sf")
txi.tx <- tximport(files, type="salmon", txOut = TRUE)
ddsTxi <- DESeqDataSetFromTximport(txi.tx, colData=samps, design = ~ Site + TPA + Site:TPA)
dds <- ddsTxi[rowSums(counts(ddsTxi)) > 10, ] 
dds$TPA <- relevel(dds$TPA, ref="0") 
dds <- DESeq(dds, test="LRT", reduced=~Site + TPA)
res <- results(dds)

And I have several questions: 1) Do I understand correctly that function “relevel” sets as a reference point only 1 time point (0 hours) in one site (in my case, Head)? Or will each site have its own reference (0 hours in Head and 0 hours in Tail)?

2) How I should intepretate Sitetail.TPA*(4,12,24,48,96) in output table and heatmaps? If I understand analysis right, these values are log2 fold changes of expression in Tail at *(4, 12,24,48,96) time point regarding to *(4,12,24,48, 96) time point in Head. Please, correct me if I wrong.

3) If I want to find genes which demonstrate differential expression in Tail at any time point regarding to 0 hour in Tail, should I use, for example, next command: results(dds, contrast=list(c(“TPA_4_vs_0”, “Sitetail.TPA4”))) ?

I will be grateful for any help!

[ATpoint]

Is there a specific reason you use txOut = TRUE instead of aggregating the transcripts to gene level? DESeq2 is intended for gene rather than tx-level analysis.