Hi All,

I want to compare Deseq2 with sSeq and I am getting similar log2foldchange values but the direction is opposite. I am confident with the DEseq2 because I use the >relevel to make sure the control is the first level, but I am not sure how does sSeq does it.

Does anyone has an idea of what I am writing wrong or why this is happening?

DEseq2:

countdata=as.matrix(read.table("2014_04_06_6h_Co2.csv",header=T, sep=",", row.names=1))

coldata= read.table ("targets.csv", header = T, sep=",",row.names=1)

dds <- DESeqDataSetFromMatrix(
  countData = countdata,
  colData = coldata,
  design = ~ Subject + Treatment)
dds

dds$Treatment <- factor(dds$Treatment,
                        levels=c("PBS", "LPS"))

dds$Treatment <- relevel (dds$Treatment,"PBS")                             

dds$Subject <- factor(dds$Subject)

colData (dds)

ddsLRT <- DESeq (dds, test ="LRT", reduced= ~Subject)

resLRT<- results(ddsLRT)

sSeq:

countdata=as.matrix(read.table("2014_04_06_6h_Co2.csv",header=T, sep=",", row.names=1))

conds.data=c("PBS","PBS","PBS","PBS","LPS","LPS","LPS","LPS")

coLevels=data.frame(subjects=c(1,2,4,5,1,2,4,5))

coLevels

## Filtering

nSamples <- ncol(countdata)
minSamples <- floor(nSamples / 2)
keep <- countdata > 4
keep <- rowSums(keep) > minSamples
countdata <- countdata[keep,]
res2 = nbTestSH(countdata, conds.data,"PBS", "LPS",
                coLevels= coLevels,  pairedDesign=TRUE, pairedDesign.dispMethod="pooled")

sSeq vs DESeq2 first two columns (after Row.names) are from sSeq and last three columns are from DESeq2

Thanks for your help,

Catalina

> sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DESeq2_1.4.5              RcppArmadillo_0.4.450.1.0 Rcpp_0.11.3               GenomicRanges_1.16.4      GenomeInfoDb_1.0.2       
 [6] IRanges_1.22.10           BiocGenerics_0.10.0       sSeq_1.2.0                RColorBrewer_1.0-5        caTools_1.17.1           

loaded via a namespace (and not attached):
 [1] annotate_1.42.1      AnnotationDbi_1.26.1 Biobase_2.24.0       bitops_1.0-6         DBI_0.3.1            genefilter_1.46.1   
 [7] geneplotter_1.42.0   grid_3.1.0           lattice_0.20-29      locfit_1.5-9.1       RSQLite_0.11.4       splines_3.1.0       
[13] stats4_3.1.0         survival_2.37-7      tools_3.1.0          XML_3.98-1.1         xtable_1.7-4         XVector_0.4.0      

In the help page for ?nbTestSH in sSeq you have the answer. It's good to check the function help pages as they typically explain their output.

Usage

nbTestSH(countsTable, conds, condA = "A", condB = "B",
numPart = 1, SHonly = FALSE, propForSigma = c(0, 1),
shrinkTarget = NULL, shrinkQuantile = NULL, plotASD = FALSE,
coLevels = NULL, contrast = NULL, keepLevelsConsistant = FALSE,
useMMdisp = FALSE, addRawData = FALSE, shrinkVariance = FALSE,
pairedDesign = FALSE, pairedDesign.dispMethod = "per-pair",
useFisher = FALSE, Dispersions = NULL, eSlope = 0.05, lwd_ASD = 4.5,
cex_ASD = 1.2)

...

log2FoldChange
The per-gene fold Changes between condition A and B in the log2 scale