Hi. I'm very new to coding, but am trying to run the DESeq2 package for an expression analysis and am having trouble that I think should be relatively easily dealt with.

The design is as such: We have two types of tissue - Graves' (G) and normal (N). 10 Graves' samples, 9 normal. For each tissue sample, there are two treatments - ATRA or DMSO.

I would like to see the ATRAvsDMSO comparison for Graves tissue alone, and for normal tissue alone. The way I currently have it set up, my only options are for tissueGvsN, txATRAvsDMSO, or tissueG.txATRA. The way I understand it, this is comparing all the tissue types, regardless of treatment, all treatment regardless of tissue type, or a comparison within a comparison of the two factors.

I'm hoping someone can help me compare treatments within each tissue type separately. One thing I want to keep in mind is pairing - is that considered in the script? the treatments are paired, but the tissue samples are not.

Thanks in advance for any help

OTSCountTable <- read.delim("counts_raw.txt", header = TRUE ,row.names = 1, sep="\t") 
head(OTSCountTable)

#Classify by txs
tx = factor( c( "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA"))
tissue = factor( c("N", "N", "N", "N", "N", "N", "N", "N", "N", "N", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "N", "N", "N", "N", "G", "G", "G", "G", "N", "N", "N", "N", "G", "G", "G", "G"))


#Normalization
library("DESeq")
library("DESeq2")

normalized = newCountDataSet(OTSCountTable, tx)
head(normalized)

normalized = estimateSizeFactors(normalized)
sizeFactors(normalized)
head(counts(normalized, normalized=TRUE))

normalized = estimateDispersions(normalized)


Normalized_coldata= data.frame( 
  row.names = colnames(normalized),
  tissue = c("N", "N", "N", "N", "N", "N", "N", "N", "N", "N", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "N", "N", "N", "N", "G", "G", "G", "G", "N", "N", "N", "N", "G", "G", "G", "G"),
  tx = c("DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA", "DMSO", "ATRA")
  )
Normalized_coldata

write.csv(normalized, file = "Normalized_Counts")


#set DMSO to untreated, reference group
library("magrittr")
Normalized_coldata$tx <-  relevel(Normalized_coldata$tx, "DMSO")
Normalized_coldata$tx

#set N to reference group
Normalized_coldata$tissue <-  relevel(Normalized_coldata$tissue, "N")
Normalized_coldata$tissue

#so now, OTSCountTable is a table with the fragment counts and Normalized_coldata is a table with information about the samples
ddsMat = DESeqDataSetFromMatrix(countData = OTSCountTable,
                               colData = Normalized_coldata,
                              design = ~tissue + tx + tissue:tx)

design(ddsMat) <- ~ tissue + tx + tissue:tx
dds = DESeq(ddsMat)


# find things to compare
resultsNames(dds)

[1] "Intercept"       "tissue_G_vs_N"   "tx_ATRA_vs_DMSO" "tissueG.txATRA"

I would like to see the ATRAvsDMSO comparison for Graves tissue alone, and for normal tissue alone. The way I currently have it set up, my only options are for tissueGvsN, txATRAvsDMSO, or tissueG.txATRA. The way I understand it, this is comparing all the tissue types, regardless of treatment, all treatment regardless of tissue type, or a comparison within a comparison of the two factors.

Actually, no, not with that design. If your design has an interaction term, then tissueGvsN does not mean all G versus all N, it means G versus N in the DMSO alone.

From the vignette

The key point to remember about designs with interaction terms is that, unlike for a design ~genotype + condition, where the condition effect represents the overall effect controlling for differences due to genotype, by adding genotype:condition, the main condition effect only represents the effect of condition for the reference level of genotype (I, or whichever level was defined by the user as the reference level).

You can do it this way for that comparison, but it's (in my opinion) hard to figure out after the fact what you did.

Far simpler is to make a new column that is disease-state concatenated to treatment, and use contrast to compare one grouping to another.

http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions

I strongly recommend you get a table of all the normalized counts, and manually get the different ratios of averages yourself, and see how they compare to what DESeq tells you is the log fold change, (they won't match perfectly, but they should be close) to make sure that you are asking for what you really want.