I would like to investigate DGE and DTE in my samples using the R packages DESeq2 and Fishpond/Swish, respectively. In my subject species, there are three phenotypes and for each of them, I have samples from 13 different tissues for a total of 267 samples. First, I would like to compare the three phenotypes considering all the samples from all the tissues. Then, I would also like to analyze the samples in order to highlight differentially expressed genes and transcripts between the three different phenotypes within each tissue.

In DESeq2, to compare the three phenotypes I set up the experimental design as follows:

design = ~ tissue + phenotype
dds <- DESeqDataSetFromTximport(txi.g, samps, design = design)
and then, after filtering and running DESeq
res1 <- results(dds, contrast = c("phenotype", "pheno1", "pheno2"), alpha = 0.05)
res2 <- results(dds, contrast = c("phenotype", "pheno1", "pheno3"), alpha = 0.05)
res3 <- results(dds, contrast = c("phenotype", "pheno2", "pheno3"), alpha = 0.05)

For the comparisons between the three phenotypes within each tissue I added the interaction term as follows:

design = ~ tissue + phenotype
dds <- DESeqDataSetFromTximport(txi.g, samps, design = design)
dds$group <- factor(paste0(dds$tissue, dds$phenotype))
design(dds) <- ~ group
and then, after filtering and running DESeq (e.g for one of the 13 tissues)
res1 <- results(dds, contrast = c("group", "tissue1pheno1", "tissue1pheno2"), alpha = 0.05)
res2 <- results(dds, contrast = c("group", "tissue1pheno1", "tissue1pheno3"), alpha = 0.05)
res3 <- results(dds, contrast = c("group", "tissue1pheno2", "tissue1pheno3"), alpha = 0.05)

Is this correct as I have done? I am more than satisfied with the results from DESeq2 so far. Now, I would like to perform the same analyses and comparisons in Swish for the DTE analysis. Do you have any recommendations on how I should proceed? I have been following the Swish tutorial, but I have a hard time understanding how to analyze the interactions and contrasts in this package.

DESeq2 analysis looks good.

For Swish, you can treat tissue as the "batch" variable, and it will asses what genes/transcripts are DE consistently across all tissues:

swish(y_1_and_2, x="phenotype", cov="tissue")

You will have to subset the dataset to perform analysis of pairs of phenotypes. DE in Swish is two-group based, although we have other complications like interactions or paired analysis, or continuous covariates. So it's just one more step ahead of swish:

y_1_and_2 <- y[, y$pheno %in% c("1","2") ]
y_1_and_2$pheno <- droplevels(y_1_and_2$pheno)

That will give you something like the first analysis in DESeq2.

Now if you want to perform across phenotype analysis _within_ each tissue (not across all tissues), I would just recommend to perform more subsetting:

y_t1_p1_and_p2 <- y[, y$tissue == "1" & y$pheno %in% c("1","2") ]
y_t1_p1_and_p2$pheno <- droplevels(y_t1_p1_and_p2$pheno)
# then regular swish
y_t1_p1_and_p2 <- swish(y_t1_p1_and_p2, x="pheno")

It cumbersome, but Swish will provide you with strong error control in the presence of large and heterogeneous transcript abundance uncertainty. If it calls a transcript DE, it means it was above and beyond the quantification uncertainty. Also I think you should have enough power to perform these within tissue analyses.