Using tximport to read in data from sailfish experiment
1
0
Entering edit mode
@linafaller-9082
Last seen 6.9 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

 

tximport sailfish rnaseq • 1.5k views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 2 hours ago
United States

Can you run traceback() after the error and report the output?

ADD COMMENT
0
Entering edit mode

Thanks for the quick reply!

​> traceback()
6: terms.default(object, data = data)
5: terms(object, data = data)
4: stats::model.matrix.default(design, data = as.data.frame(colData(se)))
3: DESeqDataSet(se, design = design, ignoreRank)
2: DESeqDataSetFromMatrix(countData = counts, colData = colData,
       design = design, ...)
1: DESeqDataSetFromTximport(data, colData = metadata, design = model)​
ADD REPLY
0
Entering edit mode

Can you show what happens when you print metadata?

ADD REPLY
0
Entering edit mode
> metadata
               genotype                    timing
Sample_1658940       2B linear feed phase EFT 22h
Sample_1658941       3B linear feed phase EFT 22h
Sample_1658952       2B linear feed phase EFT 27h
Sample_1658953       3B linear feed phase EFT 27h
Sample_1658964       2B linear feed phase EFT 22h
Sample_1658965       3B linear feed phase EFT 22h
Sample_1658976       2B linear feed phase EFT 27h
Sample_1658977       3B linear feed phase EFT 27h
Sample_1658988       2B linear feed phase EFT 22h
Sample_1658989       3B linear feed phase EFT 22h
Sample_1659000       2B linear feed phase EFT 27h
Sample_1659001       3B linear feed phase EFT 27h​
ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

Yes, I think that did the trick! Thanks so much for the help! 

This worked:

# set up the DESeq2 object
model = formula(~genotype)

dds <- DESeqDataSetFromTximport(data, colData = metadata, design = model)
ADD REPLY

Login before adding your answer.

Traffic: 692 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6