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.