Differential expression analysis in RNA-seq

DEG_RNA-seq.R

# Source: Additional file 1 from Łabaj, Paweł P., and David P. Kreil. “Sensitivity, Specificity, and Reproducibility of RNA-Seq Differential Expression Calls.” Biology Direct 11, no. 1 (December 20, 2016): 66. doi:10.1186/s13062-016-0169-7.
# https://static-content.springer.com/esm/art%3A10.1186%2Fs13062-016-0169-7/MediaObjects/13062_2016_169_MOESM1_ESM.pdf

DEfun <- function(counts, design) {
    DE <- list()
    ## limma
    gene.dge <- DGEList(counts = counts, group = factor(rep(1:2, each = 4)))
    gene.dge.norm <- calcNormFactors(gene.dge)
    gene.dge.norm2 <- gene.dge.norm
    gene.dge.norm2$counts <- gene.dge.norm2$counts - 0.5
    rpm <- voom(gene.dge.norm2, design)
    cm <- makeContrasts(AvsC = As - Cs, levels = design)
    fit <- lmFit(rpm, design)
    cf <- contrasts.fit(fit, cm)
    fit.eB <- eBayes(cf, proportion = 0.3, robust = TRUE)
    adjust.meth <- "BH"
    tt.limma <- topTable(fit.eB, adjust.method = adjust.meth, number = Inf)
    FC <- tt.limma$logFC
    AE <- tt.limma$AveExpr
    qval <- tt.limma$adj.P.Val
    names(FC) <- names(AE) <- names(qval) <- rownames(tt.limma)
    DE.list <- DElist(FC, AE, qval, 0.05)
    DE[["limma"]] <- list(DE.list = DE.list)
    
    ## edgeR
    y <- estimateCommonDisp(gene.dge.norm)
    y <- estimateTagwiseDisp(y)
    et <- exactTest(y)
    tt.edgeR <- topTags(et, n = Inf, adjust.method = adjust.meth)
    FC <- -tt.edgeR$table$logFC
    AE <- tt.edgeR$table$logCPM
    qval <- tt.edgeR$table$FDR
    names(FC) <- names(AE) <- names(qval) <- rownames(tt.edgeR$table)
    DE.list <- DElist(FC, AE, qval, 0.05)
    ## DE[['edgeR']] <- list(fit=et, tt=tt.edgeR$table, DE.list=DE.list)
    DE[["edgeR"]] <- list(DE.list = DE.list)
    
    ## DESeq2
    colData <- data.frame(condition = rep(c("A", "C"), each = 4), type = "paired-end")
    rownames(colData) <- colnames(counts)
    
    dds <- DESeqDataSetFromMatrix(countData = round(counts), colData = colData, 
        design = ~condition)
    
    
    dds <- DESeq(dds)
    tt.DESeq2 <- DESeq2::results(dds, contrast = c("condition", "A", "C"), 
        pAdjustMethod = "BH", )
    FC <- tt.DESeq2$log2FoldChange
    FC[is.na(FC)] <- 0
    AE <- log2(tt.DESeq2$baseMean + 0.5)
    qval <- tt.DESeq2$padj
    qval[is.na(qval)] <- 1
    names(FC) <- names(AE) <- names(qval) <- rownames(tt.DESeq2)
    DE.list <- DElist(FC, AE, qval, 0.05)
    DE[["DESeq2"]] <- list(DE.list = DE.list)
    
    return(DE)
}