Output differential expression and supporting data between clusters C12 and C13/C14 from Bach et al.

Stevenson_DE_clusters.r

dataList <- readRDS("ExpressionList_QC_norm_clustered_clean.rds")

m <- dataList[[1]]
pD <- dataList[[2]]
fD <- dataList[[3]]

library(dplyr)
library(edgeR)

comps <- list(C12v14=c("C12","C14"), C13v14=c("C13","C14"))

for (cname in c("C12v14","C13v14")) {
    comp <- comps[[cname]]
    pD.sub <- pD[pD$SubClusterNumbers %in% comp,]
    m.sub <- m[,as.character(pD.sub$barcode)]
    keep <- rowMeans(m.sub) > 0.1
    m.sub <- m.sub[keep,]
    fD.sub <- fD[keep,]
    rownames(m.sub) <- fD.sub$symbol

    # DGEList object
    nf <- log(pD.sub$sf/pD.sub$UmiSums)
    pD.sub$nf <- exp(nf-mean(nf))
    y <- DGEList(counts=m.sub,
         samples=pD.sub,
         genes=fD.sub,
         norm.factors=pD.sub$nf)

    # DE
    choice <- comp[1] 
    cluster <- factor(as.numeric(pD.sub$SubClusterNumbers==choice))
    de.design <- model.matrix(~cluster)
    y <- estimateDisp(y, de.design, prior.df=0, trend="none")
    fit <- glmFit(y, de.design)
    res <- glmTreat(fit, lfc=1)
    resTab <- topTags(res, n=Inf, sort.by="PValue")
    topTab <- resTab$table

    # Write DE table
    topTab <- select(topTab, id, symbol, logFC, unshrunk.logFC, logCPM, PValue, FDR)

    fname_de <- sprintf("%s.csv",cname)
    write.csv(topTab,file=fname_de, row.names=FALSE, quote=F)

    # Write sample table for plots
    fname_samples <- sprintf("cpm.%s.csv",cname)
    write.csv(y$samples, file=fname_samples, quote=F)

    # Write cpm table for plots
    fname_cpm <- sprintf("samples.%s.csv",cname)
    write.csv(cpm(y$counts), file=fname_cpm, quote=F)
}