Riborex equivalence to using an interaction term

suppressMessages(library(riborex))
suppressMessages(library(DESeq2))

Is riborex equivalent to using an interaction term?

data(riborexdata)
RNACntTable <- riborexdata$rna
RiboCntTable <- riborexdata$ribo
rnaCond <- c("control", "control", "treated", "treated")
riboCond <- c("control", "control", "treated", "treated")
res.riborex <- as.data.frame(riborex(RNACntTable, RiboCntTable, rnaCond, riboCond))

DESeq2 mode selected
combining design matrix
applying DESeq2 to modified design matrix
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

res.riborex <- res.riborex[order(res.riborex$padj),]
res.riborex.sig <- subset(res.riborex, padj<0.1)
res.riborex

summary(res.riborex)

    baseMean       log2FoldChange         lfcSE             stat               pvalue             padj       
 Min.   :      1   Min.   :-7.56081   Min.   :0.1776   Min.   :-21.14617   Min.   :0.00000   Min.   :0.0000  
 1st Qu.:     99   1st Qu.:-0.32806   1st Qu.:0.2462   1st Qu.: -0.76849   1st Qu.:0.09693   1st Qu.:0.3482  
 Median :    453   Median : 0.03521   Median :0.3106   Median :  0.09082   Median :0.40558   Median :0.7877  
 Mean   :   3607   Mean   :-0.01939   Mean   :0.5390   Mean   :  0.05475   Mean   :0.42014   Mean   :0.6333  
 3rd Qu.:   1420   3rd Qu.: 0.32434   3rd Qu.:0.5423   3rd Qu.:  0.88267   3rd Qu.:0.70605   3rd Qu.:0.9325  
 Max.   :5058124   Max.   : 6.37661   Max.   :4.3364   Max.   : 18.53963   Max.   :0.99994   Max.   :0.9999  
                                                                                             NA's   :540     

combinedCounts <- cbind(RNACntTable, RiboCntTable)
assay <- factor( c( rep('RNA', ncol(RNACntTable) ), rep( 'Ribo', ncol(RiboCntTable) ) ), levels = c("RNA", "Ribo") )
condition <- factor(c(rnaCond, riboCond), levels = c("control", "treated"))
colData <- data.frame('assay'=assay, 'condition'=condition)
dds <- DESeqDataSetFromMatrix(combinedCounts, colData,  design= ~ assay + condition + assay:condition)
dds <- DESeq(dds, test="LRT", reduced= ~ assay + condition)

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

res.interaction <- as.data.frame(results(dds))
res.interaction <- res.interaction[order(res.interaction$padj),]
res.interaction.sig <- subset(res.interaction, padj<0.1)
res.interaction

summary(res.interaction)

    baseMean       log2FoldChange         lfcSE             stat              pvalue             padj       
 Min.   :      1   Min.   :-7.56081   Min.   :0.1776   Min.   : -0.0002   Min.   :0.00000   Min.   :0.0000  
 1st Qu.:     99   1st Qu.:-0.32806   1st Qu.:0.2462   1st Qu.:  0.1423   1st Qu.:0.09654   1st Qu.:0.3449  
 Median :    453   Median : 0.03521   Median :0.3106   Median :  0.6935   Median :0.40498   Median :0.7875  
 Mean   :   3607   Mean   :-0.01939   Mean   :0.5390   Mean   :  7.2746   Mean   :0.41970   Mean   :0.6326  
 3rd Qu.:   1420   3rd Qu.: 0.32434   3rd Qu.:0.5423   3rd Qu.:  2.7618   3rd Qu.:0.70597   3rd Qu.:0.9324  
 Max.   :5058124   Max.   : 6.37661   Max.   :4.3364   Max.   :424.3174   Max.   :1.00000   Max.   :0.9999  
                                                                                            NA's   :540     

Do we see similar number of DE genes at 0.1 FDR?

Riborex

dim(res.riborex.sig)

[1] 2324    6

Using interaction

dim(res.interaction.sig)

[1] 2336    6

So the p-values differ, probably because of the way size factor estimation is handled.

Do the logFold changes look similar?

plot(res.interaction$log2FoldChange, res.riborex$log2FoldChange, cex=0.2, xlab='Using interaction', ylab='Using riborex')

Correlation

correlation.pearson <- cor(res.interaction$log2FoldChange, res.riborex$log2FoldChange)
correlation.spearman <- cor(res.interaction$log2FoldChange, res.riborex$log2FoldChange, method = 'spearman')

correlation.pearson

[1] 0.2510842

correlation.spearman

[1] 0.1476345