Hi,

I am analyzing time course RNA seq data.

My design looks like this:

To sum up: 3 Genotypes, 4 time points.

I mapped the bam files with star aligner, and extracted raw counts with genomicalignments

I run deseq2 on it with mixed success. My design formula looks like this:

ddstc <- DESeqDataSet(se,~Genotype + SA.treatment + Genotype:SA.treatment)
ddstc <- DESeq(ddstc_universal, test="LRT", reduced = ~Genotype + SA.treatment)

If I look at the pvalue histograms of the individual comparisons, I see, that things are not going all right, here an example pvalue histogramm:

(I switched the filters off, cause there is a heave expression induction upon treatment, and the filter filter all of these)

res <- results(ddstc, name='Genotype_dsr_vs_jc66', test='Wald', independentFiltering = F, cooksCutoff = F) 

hist(res$pvalue, col = "lavender", main = n, xlab = "p-values")

I found that this could be cause by the Wald test and I can run a false discovery rate correction, like explained here:

http://www-huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#inspection-and-correction-of-pvalues

and it indeed improves my pvalue histogram and i can add a corrected adjusted pvalue (not perfect, but already a lot better):

res <- res[ !is.na(res$padj), ]
res <- res[ !is.na(res$pvalue), ]
res <- res[, -which(names(res) == "padj")]
FDR.res <- fdrtool(res$stat, statistic= "normal", plot = T)
FDR.res$param[1, "sd"]
res[,"padj"]  <- p.adjust(FDR.res$pval, method = "BH")

hist(FDR.res$pval, col = "royalblue4",main = paste0(n,' - corrected'), xlab = "CORRECTED p-values")

This is all fine and good, but can I apply the same technique for the whole time course experiment?

I took the time course workflow from the RNAseq workflow:

https://www.bioconductor.org/help/workflows/rnaseqGene/#time-course-experiments

resTC <- results(ddsTC)
​betas <- coef(ddstc) %>% na.omit()
colnames(betas)
topGenes <- head(order(resTC$padj),40)
mat <- betas[topGenes, -c(1)]
thr <- 3
mat[mat < -thr] <- -thr
mat[mat > thr] <- thr
pheatmap(mat, breaks=seq(from=-thr, to=thr, length=101),
         cluster_col=FALSE)

If I understand it correctly, it takes the results of the last resultsName, and identifies the most significant genes, by ordering them with their adjusted pvalue (which is the same in all results, as it's globally Hochberg corrected). Then extracting the coefficients for these genes and plotting it into a heatmap.

But can I run the FDR tool on this result and will it correct the pvalues for all result contrasts? My assumption is that the FDR will correct the pvalue only for the last called results (resTC) and adjust the pvalue based on that and will not take any other contrast into account. Which then in turn would mean, that I have different adjusted pvalues in different contrasts, which then in turn means, that the topgenes I want to extract, are not the real topgenes.

So, can I do a FDR correction on the whole experiment ? Does it make sense to do it?

And I have a second question about the interacting terms:

If I understand correctly the results of our example here:

results(ddsTC , name='Genotypedsr.SA.treatment6h') 

this shows the results for my Genotype dsr at 6h compared to the Genotype dsr at 0h, minus the effect on the control at the same time points ?

Is there a way to call the results without reducing it by the control levels ?

at least that's my interpretation of:

# the condition effect for genotype II
# this is, by definition, the main effect *plus* the interaction term
# (the extra condition effect in genotype II compared to genotype I).
results(dds, list( c("condition_B_vs_A","genotypeII.conditionB") ))

Cheer

Jakob

hi,

So my answer to this would be, to not use fdrtool without filtering, because this downstream method requires certain properties of the distribution of p-values.

The very small count genes do not produce smooth, uniform p-values (this was noted in the DESeq2 paper, Supp Fig 26).

You should filter first, if you want to achieve a smooth p-value distribution to then operate on. You can do this like so:

dds <- estimateSizeFactors(dds)
keep <- ...
dds[keep,]
dds <- DESeq(dds, ...)
res <- results(dds)