gse52202 script.R

gse52202 script.R

fc.counts <- read.csv("fccounts.csv",row.names=1)
pheno <- read.csv("pheno.csv",row.names=1)
pheno$sam <- paste0(pheno$run,"Aligned.out.sam")

library(DESeq2)
dds <- DESeqDataSetFromMatrix(fc.counts, 
        DataFrame(pheno[match(colnames(fc.counts),pheno$sam),]),
        ~ condition)
dds$condition <- relevel(dds$condition,"Normal control")
levels(dds$condition)

# i think these replicates look like duplicates...
dds <- dds[,1:8 * 2 - 1]

mcols(dds)$ENTREZID <- rownames(dds)
library(org.Hs.eg.db)
map <- select(org.Hs.eg.db,mcols(dds)$ENTREZID,"SYMBOL","ENTREZID")
symbols <- map$SYMBOL[match(mcols(dds)$ENTREZID,map$ENTREZID)]
symbols[duplicated(symbols)] <- NA
symbols[is.na(symbols)] <- paste0("NA",seq_len(sum(is.na(symbols))))
rownames(dds) <- symbols

dds <- DESeq(dds)
res <- results(dds)

# CBLN4 is a "dispersion outlier", has high dispersion
# and then the log fold change is moderated more toward 0
idx <- c("DPP6","SLITRK2",paste0("CBLN",c(1,2,4)))
res[idx,]
round(counts(dds,norm=TRUE)[idx,])
mcols(dds,use.names=TRUE)[idx,"dispOutlier"]

plotMA(res, alpha=.5, ylim=c(-2,2))
points(as.matrix(res[idx,c("baseMean","log2FoldChange")]),
       cex=2,lwd=4,col="green")

# there is a parameter which tunes how high a dispersion
# must be, so that the gene does not use the posterior estimate (MAP)
dds2 <- dds
dds2 <- estimateDispersionsMAP(dds2,outlierSD=4)
dds2 <- nbinomWaldTest(dds2)
res2 <- results(dds2)

# after changing this param, i get more similar p-value to the paper
# log fold change is still moderated but less
mcols(dds2,use.names=TRUE)[idx,"dispOutlier"]
res2[idx,]

plotMA(res2, alpha=.5, ylim=c(-3,3))
points(as.matrix(res2[idx,c("baseMean","log2FoldChange")]),
       cex=2,lwd=4,col="green")