We recommend to use DESeq2 on sequencing time course data with the following settings: Firstly, we recommend the usage of a natural cubic spline basis with less degrees of freedom than the number of sampled time points. Secondly, we recommend less stringent dispersion outlier calling to avoid missing differentially expressed genes with observations that are zero. Thirdly, we recommend the usage of condition-wise batch correction in case-control analysis. Below, we supply example code for this adapted DESeq2 model on simulated data in the case-only and the case-control scenario both with and without batch correction:
Without batch correction:
suppressPackageStartupMessages(library(ImpulseDE2))
suppressPackageStartupMessages(library(splines))
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(BiocParallel))
## simulate data
lsSim <- simulateDataSetImpulseDE2(
vecTimePointsA=rep(1:6, 3),
vecTimePointsB=NULL,
vecBatchesA=rep("B_NULL", 6*3),
vecBatchesB=NULL,
scaNConst=100,
scaNImp=50,
scaNLin=50,
scaNSig=50,
scaNRand=0,
scaMumax=1000,
scaSeedInit = 1)
## run DESeq2
dfAnnotationDESeq2 <- lsSim$dfAnnotation
# create a natural cubic spline basis with 4 degrees of freedom
matTimeSplineBasis <- ns(
dfAnnotationDESeq2$Time, df=4)
colnames(matTimeSplineBasis) <-
paste0("spline", seq(1, dim(matTimeSplineBasis)[2]))
# add spline basis into annotation data frame
dfAnnotationDESeq2 <- cbind(dfAnnotationDESeq2, matTimeSplineBasis)
# create DESeq2 object with spline basis vectors as individual linear predictors
dds <- DESeqDataSetFromMatrix(
countData = lsSim$matObservedCounts,
colData = dfAnnotationDESeq2,
design = ~spline1 + spline2 + spline3 + spline4)## converting counts to integer mode# estimate dispersion with less stringent dispersion outlier handling
# the outlier calling is controlled via outlierSD of estimateDispersionsMAP()
dds <- estimateSizeFactors(dds)
dds <- estimateDispersionsGeneEst(dds)
dds <- estimateDispersionsFit(dds)## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.dds <- estimateDispersionsMAP(dds, outlierSD = 10)
# perform a log-likelihood ratio test
dds <- nbinomLRT(dds, full = ~spline1 + spline2 + spline3 + spline4,
reduced = ~1)
# extract results
print(head( results(dds) ))## log2 fold change (MLE): spline4
## LRT p-value: '~ spline1 + spline2 + spline3 + spline4' vs '~ 1'
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## gene1 300.1373 -0.07164473 0.3442949 1.0974737 0.894672827 0.95177960
## gene2 366.7376 -0.43020095 0.2489286 16.1016613 0.002885739 0.01145135
## gene3 588.3016 -0.25562817 0.4340152 0.6982777 0.951541147 0.98097820
## gene4 952.4941 0.34247246 0.3869531 2.9938291 0.558858622 0.73923098
## gene5 209.4934 -0.20025409 0.3405095 3.4749894 0.481691432 0.68035513
## gene6 977.1628 0.05911703 0.2799341 11.5164745 0.021333532 0.06666729With batch correction:
suppressPackageStartupMessages(library(ImpulseDE2))
suppressPackageStartupMessages(library(splines))
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(BiocParallel))
## simulate data
lsSim <- simulateDataSetImpulseDE2(
vecTimePointsA=rep(1:6, 3),
vecTimePointsB=NULL,
vecBatchesA=do.call(c, lapply(
seq(1,3), function(i) rep(paste0("B",i), 6) )),
vecBatchesB=NULL,
scaMuBatchEffect=1,
scaSDBatchEffect=0.2,
scaNConst=100,
scaNImp=50,
scaNLin=50,
scaNSig=50,
scaNRand=0,
scaMumax=1000,
scaSeedInit = 1)
## run DESeq2
dfAnnotationDESeq2 <- lsSim$dfAnnotation
# create a natural cubic spline basis with 4 degrees of freedom
matTimeSplineBasis <- ns(
dfAnnotationDESeq2$Time, df=4)
colnames(matTimeSplineBasis) <-
paste0("spline", seq(1, dim(matTimeSplineBasis)[2]))
# add spline basis into annotation data frame
dfAnnotationDESeq2 <- cbind(dfAnnotationDESeq2, matTimeSplineBasis)
# create DESeq2 object with spline basis vectors as individual linear predictors
dds <- DESeqDataSetFromMatrix(
countData = lsSim$matObservedCounts,
colData = dfAnnotationDESeq2,
design = ~Batch + spline1 + spline2 + spline3 + spline4)## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors# estimate dispersion with less stringent dispersion outlier handling
# the outlier calling is controlled via outlierSD of estimateDispersionsMAP()
dds <- estimateSizeFactors(dds)
dds <- estimateDispersionsGeneEst(dds)
dds <- estimateDispersionsFit(dds)## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.dds <- estimateDispersionsMAP(dds, outlierSD = 10)
# perform a log-likelihood ratio test
dds <- nbinomLRT(dds, full = ~Batch + spline1 + spline2 + spline3 + spline4,
reduced = ~Batch)
# extract results
print(head( results(dds) ))## log2 fold change (MLE): spline4
## LRT p-value: '~ Batch + spline1 + spline2 + spline3 + spline4' vs '~ Batch'
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## gene1 332.2238 0.0879040309 0.2924181 0.9434907 0.91824014 0.96535590
## gene2 422.1556 0.9170551773 0.3212048 13.1827163 0.01041667 0.03364355
## gene3 780.6804 0.2937970080 0.3761408 0.9062017 0.92366941 0.96535590
## gene4 1211.7727 -0.4093646633 0.3156516 5.8675786 0.20925825 0.35368283
## gene5 239.4681 -0.1719123830 0.3563992 1.8938956 0.75526629 0.89486528
## gene6 1102.1868 0.0006391691 0.2693686 12.9481841 0.01153177 0.03649294Without batch correction:
suppressPackageStartupMessages(library(ImpulseDE2))
suppressPackageStartupMessages(library(splines))
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(BiocParallel))
## simulate data
lsSim <- simulateDataSetImpulseDE2(
vecTimePointsA=rep(1:6, 3),
vecTimePointsB=rep(1:6, 3),
vecBatchesA=rep("B_NULL", 6*3),
vecBatchesB=rep("B_NULL", 6*3),
scaNConst=100,
scaNImp=50,
scaNLin=50,
scaNSig=50,
scaNRand=0,
scaMumax=1000,
scaSeedInit = 1)
## run DESeq2
dfAnnotationDESeq2 <- lsSim$dfAnnotation
# create a natural cubic spline basis with 4 degrees of freedom
matTimeSplineBasis <- ns(
dfAnnotationDESeq2$Time, df=4)
colnames(matTimeSplineBasis) <-
paste0("spline", seq(1, dim(matTimeSplineBasis)[2]))
# add spline basis into annotation data frame
dfAnnotationDESeq2 <- cbind(dfAnnotationDESeq2, matTimeSplineBasis)
# create DESeq2 object with spline basis vectors as individual linear predictors
dds <- DESeqDataSetFromMatrix(
countData = lsSim$matObservedCounts,
colData = dfAnnotationDESeq2,
design = ~Condition + Condition:spline1 + Condition:spline2 +
Condition:spline3 + Condition:spline4)## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors# estimate dispersion with less stringent dispersion outlier handling
# the outlier calling is controlled via outlierSD of estimateDispersionsMAP()
dds <- estimateSizeFactors(dds)
dds <- estimateDispersionsGeneEst(dds)
dds <- estimateDispersionsFit(dds)## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.dds <- estimateDispersionsMAP(dds, outlierSD = 10)
# perform a log-likelihood ratio test
dds <- nbinomLRT(dds, full = ~Condition + Condition:spline1 + Condition:spline2 +
Condition:spline3 + Condition:spline4,
reduced = ~spline1 + spline2 + spline3 + spline4)
# extract results
print(head( results(dds) ))## log2 fold change (MLE): Conditioncontrol.spline4
## LRT p-value: '~ Condition + Condition:spline1 + Condition:spline2 + Condition:spline3 + Condition:spline4' vs '~ spline1 + spline2 + spline3 + spline4'
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## gene1 253.1908 -0.1318652 0.3924393 7.732633 0.17159937 0.40121224
## gene2 388.7622 -0.3678561 0.3906130 8.315500 0.13968363 0.34920908
## gene3 589.7440 0.6272186 0.3207137 13.386425 0.02001451 0.07942264
## gene4 915.7521 0.5030112 0.3572772 6.974332 0.22255627 0.46365890
## gene5 192.9654 0.4259325 0.3920925 6.767761 0.23849309 0.48251801
## gene6 908.7567 -0.3362782 0.3145794 5.713672 0.33508413 0.58180898With batch correction:
suppressPackageStartupMessages(library(ImpulseDE2))
suppressPackageStartupMessages(library(splines))
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(BiocParallel))
## simulate data
lsSim <- simulateDataSetImpulseDE2(
vecTimePointsA=rep(1:6, 3),
vecTimePointsB=rep(1:6, 3),
vecBatchesA=do.call(c, lapply(
seq(1,3), function(i) rep(paste0("B",i), 6) )),
vecBatchesB=do.call(c, lapply(
seq(1,3), function(i) rep(paste0("B",i), 6) )),
scaMuBatchEffect=1,
scaSDBatchEffect=0.2,
scaNConst=100,
scaNImp=50,
scaNLin=50,
scaNSig=50,
scaNRand=0,
scaMumax=1000,
scaSeedInit = 1)
## run DESeq2
dfAnnotationDESeq2 <- lsSim$dfAnnotation
# create a natural cubic spline basis with 4 degrees of freedom
matTimeSplineBasis <- ns(
dfAnnotationDESeq2$Time, df=4)
colnames(matTimeSplineBasis) <-
paste0("spline", seq(1, dim(matTimeSplineBasis)[2]))
# add spline basis into annotation data frame
dfAnnotationDESeq2 <- cbind(dfAnnotationDESeq2, matTimeSplineBasis)
# create DESeq2 object with spline basis vectors as individual linear predictors
dds <- DESeqDataSetFromMatrix(
countData = lsSim$matObservedCounts,
colData = dfAnnotationDESeq2,
design = ~Condition + Condition:spline1 + Condition:spline2 +
Condition:spline3 + Condition:spline4 + Condition:Batch)## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors# estimate dispersion with less stringent dispersion outlier handling
# the outlier calling is controlled via outlierSD of estimateDispersionsMAP()
dds <- estimateSizeFactors(dds)
dds <- estimateDispersionsGeneEst(dds)
dds <- estimateDispersionsFit(dds)
dds <- estimateDispersionsMAP(dds, outlierSD = 10)
# perform a log-likelihood ratio test
dds <- nbinomLRT(dds, full = ~Condition + Condition:spline1 + Condition:spline2 +
Condition:spline3 + Condition:spline4 + Condition:Batch,
reduced = ~spline1 + spline2 + spline3 + spline4 + Batch)
# extract results
print(head( results(dds) ))## log2 fold change (MLE): Conditioncontrol.BatchB3
## LRT p-value: '~ Condition + Condition:spline1 + Condition:spline2 + Condition:spline3 + Condition:spline4 + Condition:Batch' vs '~ spline1 + spline2 + spline3 + spline4 + Batch'
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## gene1 255.5417 0.07621837 0.2898356 11.563099 0.11587328 0.3656402
## gene2 389.9618 0.03258513 0.2737443 12.353518 0.08951636 0.3108207
## gene3 595.1599 -0.06389875 0.2317662 10.649002 0.15467904 0.4210152
## gene4 925.8095 -0.02921463 0.2457753 6.665662 0.46450065 0.7826958
## gene5 200.5167 -0.33910048 0.2928535 13.963302 0.05183793 0.2124505
## gene6 911.3917 0.25502562 0.2275118 7.573716 0.37167789 0.6659283