Hi,
I am trying to use DESeq2 to analyze an experiment where I have read counts from individuals assigned to 1 of 2 groups, and each individual is measured at 5 separate time points. I want to determine whether changes over time for read counts of individual bacteria differ between groups (i.e., group-by-time interaction). I have primarily been using the webpages "Analyzing RNA-seq data with DESeq2" and "RNA-seq workflow: gene-level exploratory analysis and differential expression" to figure out the appropriate code.
My code (where diet=study group, studyday=time, and ind distinguishes individuals nested with a group as described on the Analyzing RNA-seq data with DESeq2 website section titled "Group-specific effects, individuals nested within groups):
mat <- Genus.table.reads_NOnames
coldata <- altitude_map_no_cont_9.5
dds <- DESeqDataSetFromMatrix(countData = as.matrix(mat),
colData = DataFrame(coldata),
design = formula(~diet + studyday))
#I have n=8 in group 1 and n=9 in group 2, so use this to remove columns of all zeros.
mm1<-model.matrix(~ diet + studyday + diet:ind + diet:studyday, coldata)
idx<-which(colSums(mm1==0)==nrow(mm1))
mm1<-mm1[,-idx]
mm0<-model.matrix(~ diet + studyday + diet:ind, coldata)
idx<-which(colSums(mm0==0)==nrow(mm0))
mm0<-mm0[,-idx]
ddsOTU <- DESeq(dds, betaPrior=FALSE, full=mm1)
ddsLRT <- DESeq(ddsOTU, test="LRT", betaPrior=FALSE, full=mm1, reduced=mm0)
Everything works fine until the last line of code. Instead of a result, I receive the error message: "$ operator is invalid for atomic vectors."
I have modified the code in multiple different ways, but still receive the same error.
Any help is greatly appreciated.
Thank you,
Phil
Session info -----------------------------------------------------------------------------------------
setting value
version R version 3.0.3 (2014-03-06)
system x86_64, mingw32
ui RStudio (0.98.501)
language (EN)
collate English_United States.1252
tz EDT
Packages ---------------------------------------------------------------------------------------------
package * version date source
annotate * 1.40.1 2014-03-01 Bioconductor
AnnotationDbi * 1.24.0 2013-10-15 Bioconductor
Biobase * 2.22.0 2013-10-15 Bioconductor
BiocGenerics 0.8.0 2013-10-15 Bioconductor
BiocInstaller 1.12.1 2014-04-14 Bioconductor
DBI * 0.3.1 2014-09-24 CRAN (R 3.0.3)
DESeq2 1.2.10 2014-02-05 Bioconductor
devtools * 1.7.0 2015-01-17 CRAN (R 3.0.3)
genefilter * 1.44.0 2013-10-15 Bioconductor
GenomicRanges 1.14.4 2013-12-10 Bioconductor
IRanges 1.20.7 2014-03-06 Bioconductor
lattice * 0.20-27 2014-02-27 CRAN (R 3.0.3)
locfit * 1.5-9.1 2013-04-20 CRAN (R 3.0.3)
RColorBrewer * 1.0-5 2011-06-17 CRAN (R 3.0.0)
Rcpp 0.11.5 2015-03-06 CRAN (R 3.0.3)
RcppArmadillo 0.4.650.1.1 2015-02-26 CRAN (R 3.0.3)
RSQLite * 1.0.0 2014-10-25 CRAN (R 3.0.3)
rstudio * 0.98.501 2015-08-19 local
rstudioapi * 0.3.1 2015-04-07 CRAN (R 3.0.3)
survival * 2.37-7 2014-01-22 CRAN (R 3.0.3)
XML * 3.98-1.1 2013-06-20 CRAN (R 3.0.3)
xtable * 1.7-4 2014-09-12 CRAN (R 3.0.3)
XVector 0.2.0 2013-10-15 Bioconductor
