Entering edit mode
lina.faller
▴
10
@linafaller-9082
Last seen 7.7 years ago
Hi all,
I am just getting started with RNAseq DE analyses. I processed my data using sailfish and would like to use DESeq2 to further process the results. Here is a snippet of my code
library('tximport') library('DESeq2') sf_dirs <- dir() files = paste(sf_dirs, 'quant.sf', sep = '/') tx2gene = read.table('ecoli_transcript_to_gene_mapping.txt', header = FALSE, stringsAsFactors = FALSE, sep = '\t') data = tximport(files, type = "sailfish", tx2gene = tx2gene) metadata = read.table('metadata.txt', header = TRUE, stringsAsFactors = FALSE, sep = '\t', row.names = 1) colnames(metadata) = c('genotype', 'timing') # set up the DESeq2 object model <- "~ genotype" dds <- DESeqDataSetFromTximport(data, colData = metadata, design = model)
Unfortunately, when I run the last line, I get the following error message:
Error: $ operator is invalid for atomic vectors My tximport object looks like this: > str(data) List of 4 $ abundance : num [1:4325, 1:12] 0 1028.1 623.6 1086.6 62.5 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:4325] "b0001" "b0002" "b0003" "b0004" ... .. ..$ : NULL $ counts : num [1:4325, 1:12] 0 331.01 65.01 168.01 1.01 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:4325] "b0001" "b0002" "b0003" "b0004" ... .. ..$ : NULL $ length : num [1:4325, 1:12] 24.9 2262.6 732.6 1086.6 113.6 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:4325] "b0001" "b0002" "b0003" "b0004" ... .. ..$ : NULL $ countsFromAbundance: chr "no" > sessionInfo() R version 3.3.1 (2016-06-21) Platform: x86_64-apple-darwin13.4.0 (64-bit) Running under: OS X 10.11.6 (El Capitan) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] parallel stats4 stats graphics grDevices utils datasets methods base other attached packages: [1] DESeq2_1.14.0 SummarizedExperiment_1.4.0 Biobase_2.34.0 GenomicRanges_1.26.1 [5] GenomeInfoDb_1.10.0 IRanges_2.8.0 S4Vectors_0.12.0 BiocGenerics_0.20.0 [9] tximport_1.2.0 BiocInstaller_1.24.0 shiny_0.14.2 wasabi_0.1 [13] sleuth_0.28.1 dplyr_0.5.0 ggplot2_2.1.0 loaded via a namespace (and not attached): [1] locfit_1.5-9.1 Rcpp_0.12.7 lattice_0.20-33 tidyr_0.6.0 assertthat_0.1 [6] digest_0.6.10 mime_0.5 R6_2.2.0 plyr_1.8.4 chron_2.3-47 [11] acepack_1.4.1 RSQLite_1.0.0 zlibbioc_1.20.0 lazyeval_0.2.0 data.table_1.9.6 [16] annotate_1.52.0 rpart_4.1-10 Matrix_1.2-6 labeling_0.3 splines_3.3.1 [21] BiocParallel_1.8.1 geneplotter_1.52.0 stringr_1.1.0 foreign_0.8-66 RCurl_1.95-4.8 [26] munsell_0.4.3 httpuv_1.3.3 htmltools_0.3.5 nnet_7.3-12 tibble_1.2 [31] gridExtra_2.2.1 htmlTable_1.7 Hmisc_4.0-0 XML_3.98-1.4 bitops_1.0-6 [36] grid_3.3.1 jsonlite_1.1 xtable_1.8-2 gtable_0.2.0 DBI_0.5-1 [41] magrittr_1.5 scales_0.4.0 stringi_1.1.2 XVector_0.14.0 reshape2_1.4.2 [46] genefilter_1.56.0 latticeExtra_0.6-28 Formula_1.2-1 rjson_0.2.15 RColorBrewer_1.1-2 [51] tools_3.3.1 survival_2.40-1 AnnotationDbi_1.36.0 colorspace_1.2-7 rhdf5_2.18.0 [56] cluster_2.0.4 knitr_1.14
If anyone has any advice as to how I can fix my code, please let me know.
Thanks in advance!
~Lina
Thanks for the quick reply!
Can you show what happens when you print
metadata
?I think I've got it. It's passing a character string instead of a formula as the design. Can you try just using a formula for design?
Yes, I think that did the trick! Thanks so much for the help!
This worked: