Differential gene expression: how do I correctly compare three groups to each other?

Question

I'm trying to figure out the right way to do differential gene expression analysis with a discrete variable with three groups. The context is, I want to assess differential expression as a function of whether samples are WT, have heterozygous loss, or homozygous loss for a particular gene. A key part of it is comparing the differences between the WT and het-loss samples with the differences between het-loss and hom-loss samples.

I tried doing just doing differential expression with DESeq2, treating it as a categorical variable, but I got a strange result: differential expression between het and hom samples seemed to usually go in the opposite direction as in wt vs. het (so if a gene was upregulated in het relative to WT, it was downregulated in hom samples relative to het). This was very counterintuitive; why would losing the second copy of a gene consistently have the opposite effect on overall gene expression that losing the first copy has? So to test if something was going awry, I did the same analysis, but randomly assigning samples into three groups, and found the same thing: the log fold change comparing group 3 to 2 was strongly negatively correlated with the log FC comparing group 2 to 1, even though the groups are entirely random. So this must be some artifact of the model being used (note that I observe the same thing when used limma voom instead of DESeq 2). I include an example here from data loaded from TCGA via TCGAbiolinks:

library(TCGAbiolinks)
library(DESeq2)

#Loading KIRC data from TCGA:
query_TCGA = GDCquery(
  project = "TCGA-KICH",
  data.type = "Gene Expression Quantification",
  data.category = "Transcriptome Profiling", # parameter enforced by GDCquery
  experimental.strategy = "RNA-Seq",
  workflow.type = "STAR - Counts")
download <- GDCdownload(query = query_TCGA)
dat <- GDCprepare(query = query_TCGA, save = FALSE)
rna <- as.data.frame(SummarizedExperiment::assay(dat))

#Randomly assigning samples to groups 1, 2, and 3 and running DESeq2:
set.seed(1234)
xvar <- factor(sample(1:3, ncol(rna), replace = TRUE), levels = c(1,2,3))
columnData <- data.frame('Patient'=colnames(rna),'xvar'=xvar)
columnData$xvar <- factor(columnData$xvar, levels = c('1','2','3'))
deseq2Data <- DESeqDataSetFromMatrix(countData=rna, colData=columnData, design= ~     xvar)
deseq2Data <- deseq2Data[rowSums(counts(deseq2Data)) > 5, ]
deseq2Data <- DESeq(deseq2Data)
res.1_v_3 <- as.data.frame(results(deseq2Data, contrast = c('xvar','1','3')))
res.1_v_2 <- as.data.frame(results(deseq2Data, contrast = c('xvar','1','2')))
res.2_v_3 <- as.data.frame(results(deseq2Data, contrast = c('xvar','2','3')))

#Merging comparisons of 1 v 2, 1 v 3, and 2 v 3 into one data frame:
merged <- merge(data.frame('Gene'=rownames(res.1_v_2),     'log2FC_1_v_2'=res.1_v_2$log2FoldChange, 'padj_1_v_2'=res.1_v_2$padj), 
                data.frame('Gene'=rownames(res.2_v_3),     'log2FC_2_v_3'=res.2_v_3$log2FoldChange, 'padj_2_v_3'=res.2_v_3$padj), by='Gene')
merged <- merge(merged, data.frame('Gene'=rownames(res.1_v_3),     'log2FC_1_v_3'=res.1_v_3$log2FoldChange, 
                                   'padj_1_v_3'=res.1_v_3$padj), by='Gene')
#Color coding for illustration:
merged$color <- 'gray'
merged$sig_1_v_2 <- merged$padj_1_v_2 < .05
merged$sig_1_v_3 <- merged$padj_1_v_3 < .05
merged$sig_2_v_3 <- merged$padj_2_v_3 < .05
merged$color[which(merged$padj_1_v_2 < .05 & merged$padj_2_v_3 > .05)] <- 'blue'
merged$color[which(merged$padj_1_v_2 > .05 & merged$padj_2_v_3 < .05)] <- 'red'
merged$color[which(merged$padj_1_v_2 < .05 & merged$padj_2_v_3 < .05)] <- 'purple'
maxim <- max(abs(c(merged$log2FC_1_v_2, merged$log2FC_2_v_3)))
coef <- cor(merged$log2FC_1_v_2, merged$log2FC_2_v_3)
p <- try(round(chisq.test(table(merged[,c('sig_1_v_2','sig_2_v_3')]))$p.value,     digits = 4), silent = TRUE)
if(class(p) == 'try-error'){p <- NA}

plot(merged$log2FC_1_v_2, merged$log2FC_2_v_3, col=merged$color, pch=16,     cex=.5,main='Log2 FC',
     xlim=c(-maxim1,maxim), ylim=c(-maxim,maxim),xlab=c('LogFC:     Het/WT'),ylab=c('LogFC: Biallelic/Het'))
abline(a=0,b=coef1,col='red')
legend('topright',legend = c('Coefficient','Significant Group 2 v. Group 1',
                             'Significant Group 3 vs. group 2','Significant for     both comparisons'), 
       lty=c(1, NA, NA,NA), pch=c(NA, 16, 16,16),     col=c('red','blue','red','purple'))
text(x=-maxim*.8, y=-maxim, labels=paste('p=',p,' (chi-sq. test',sep=''))

Is there something particular I'm doing wrong? And what would be the correct way to do differential expression analysis with an independent variable that is essentially an ordinal variable? Thanks.

Answer 1

I will try to answer this open-ended question with general statements towards the analysis of larger cohorts:

Larger cohorts are usually not homogeneous.

Based on the provided annotations (colData(dat)) there are cancer samples and normal controls in this dataset, and they're not equal in size. It is 66 Primary solid Tumor and 25 Solid Tissue Normal. Running a PCA, for example via plotPCA() after running vst() will demonstrate that this is the main source of variation in this dataset, separating along PC1. That means a random split of the dataset will still have more or fewer cancers or normals in one or the other groups, and at that sample size and based on the expected strong expression differences between cancer and normal this will get DEGs almost certainly.

Here the PCA via plotPCA(vst(dds, blind=TRUE), "definition")

That having said, you need to inspect the data you deal with in depth before throwing any test on it. If not you might test some sort of confounders or factors that bias the analysis without noticing it. Based on the PCA there is also notable heterogeneity in the cancer samples, and the annotations list additional factors that can contribute to this, be it gender, smoking status, tumor stage, age and others. You would need to ensure that your random assignment does not nest any of these factors or the tumor/normal type with the random groups.

Data exploration is key. There was a recent paper that got quite some attention on Twitter claiming that common testing frameworks (incl DESeq2) lead to exaggerated false positive rates. However, the DESeq2 author (analysis here) showed that this paper missed aforementioned proper data exploration, and the dataset infact had unaddressed technical variation and samples forming groups within that dataset. Addressing that (he did it via RUVSeq) then deemphasized the strong claims in that paper.

Summarizing: Large cohorts, especially human ones, often have internal nesting and confounding with other factors that are unwanted. That can be the aforementioned obvious ones like age or gender (sex) but also more subtle ones often not recorded, such as dietary status, unreported drug consumption, amount of physical exercise the person does (if any), genetic predisposition, previous disease or accidents/trauma and so on. That having said, I am quite sure that the "random" DEGs you see are nested unwanted variation rather than anything "wrong" with the testing frameworks.