drisso / zinbwave

Clone of the Bioconductor repository for the zinbwave package, see https://bioconductor.org/packages/zinbwave
43 stars 10 forks source link

Error if add gene-level covariates #34

Open rtian30 opened 6 years ago

rtian30 commented 6 years ago

Hi, I added V matrix (gene-level covariates) in ZinbFit, but got an error message: Error in validObject(.Object) : invalid class “ZinbModel” object: beta_mu must have J columns!

This is my code: data<-read.csv(file ="./new_Expression_matrix.csv", header=T, row.names = 1, sep = "") data<-as.data.frame(t(data)) #row names are genes

construct SummartizedExperiment class

ExperimentalDesign<-read.csv(file="./Experimental_Design.csv",header=T,row.names=1,sep=",") UMI<- as.matrix(data[-1, ]) #small subset, change later colData <- DataFrame(ExperimentalDesign,row.names=colnames(UMI)) #experimental design scRNA<-SummarizedExperiment(assays=list(counts=UMI), colData=colData) #the function uses assay named "counts"

identify the 2000 most variable genes, and keep in scRNA, you could also increase the number of variable genes kept

assay(scRNA) %>% log1p %>% rowVars -> vars names(vars) <- rownames(scRNA) vars <- sort(vars, decreasing = TRUE) head(vars) scRNA <- scRNA[names(vars)[1:1000],]

create gene-level covariates

mart <- useMart("ensembl") mart <- useDataset("hsapiens_gene_ensembl", mart = mart) #get access to ensembl,hsapiens dataset bm <- getBM(attributes=c('hgnc_symbol', 'start_position', 'end_position', 'percentage_gene_gc_content'), filters = 'hgnc_symbol', values = rownames(scRNA), mart = mart)

bm$length <- bm$end_position - bm$start_position len <- tapply(bm$length, bm$hgnc_symbol, mean) len <- len[rownames(scRNA)] gcc <- tapply(bm$percentage_gene_gc_content, bm$hgnc_symbol, mean) gcc <- gcc[rownames(scRNA)] rowData(scRNA) <- data.frame(gccontent = gcc, length = len)

Zinbfit, fit a zinb regression model

fit<-zinbFit(Y=scRNA, X="~Run", V="~gccontent+ log(length)", verbose=T, K = 2, epsilon=1000)

It is very similar as vignette.

And I also run the code from vignette, the same error message happened.

Can you check if there is a problem? Thanks!

drisso commented 6 years ago

Hi @rtian30 ,

what is the exact error that you're getting? What version of the package are you using? Copy/pasting the result of sessionInfo() as well as the error message would be useful.

erbader commented 6 years ago

UPDATE: My issue was with my biomaRt query. Some of the HUGO gene symbols were not in the ensembl database, meaning the matrix of my input data would have a larger J than the matrix V.

Has this been resolved? I'm having the same issue, but I ran the code in the vignette and it worked.

Create model: Error in validObject(.Object) : invalid class “ZinbModel” object: beta_mu must have J columns!

R version 3.4.4 (2018-03-15) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Sierra 10.12.6

Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

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

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

other attached packages: [1] Rtsne_0.13 magrittr_1.5
[3] scRNAseq_1.4.0 biomaRt_2.34.2
[5] dplyr_0.7.6 scone_1.2.0
[7] zinbwave_1.0.0 SingleCellExperiment_1.0.0 [9] SummarizedExperiment_1.8.1 DelayedArray_0.4.1
[11] matrixStats_0.54.0 Biobase_2.38.0
[13] GenomicRanges_1.30.3 GenomeInfoDb_1.14.0
[15] IRanges_2.12.0 S4Vectors_0.16.0
[17] BiocGenerics_0.24.0 Seurat_2.3.4
[19] Matrix_1.2-14 cowplot_0.9.3
[21] ggplot2_3.0.0

loaded via a namespace (and not attached): [1] reticulate_1.10 R.utils_2.7.0
[3] tidyselect_0.2.4 RSQLite_2.1.1
[5] AnnotationDbi_1.40.0 htmlwidgets_1.2
[7] grid_3.4.4 trimcluster_0.1-2.1
[9] BiocParallel_1.12.0 DESeq_1.30.0
[11] munsell_0.5.0 codetools_0.2-15
[13] ica_1.0-2 withr_2.1.2
[15] colorspace_1.3-2 energy_1.7-5
[17] knitr_1.20 rstudioapi_0.7
[19] pspline_1.0-18 ROCR_1.0-7
[21] bayesm_3.1-0.1 robustbase_0.93-1.1
[23] dtw_1.20-1 gbRd_0.4-11
[25] Rdpack_0.9-0 labeling_0.3
[27] lars_1.2 GenomeInfoDbData_1.0.0
[29] hwriter_1.3.2 bit64_0.9-7
[31] rhdf5_2.22.0 EDASeq_2.12.0
[33] diptest_0.75-7 R6_2.2.2
[35] locfit_1.5-9.1 hdf5r_1.0.0
[37] flexmix_2.3-14 bitops_1.0-6
[39] assertthat_0.2.0 SDMTools_1.1-221
[41] scales_1.0.0 nnet_7.3-12
[43] gtable_0.2.0 npsurv_0.4-0
[45] RUVSeq_1.12.0 rlang_0.2.2
[47] genefilter_1.60.0 splines_3.4.4
[49] rtracklayer_1.38.3 lazyeval_0.2.1
[51] acepack_1.4.1 hexbin_1.27.2
[53] checkmate_1.8.5 reshape2_1.4.3
[55] GenomicFeatures_1.30.3 backports_1.1.2
[57] Hmisc_4.1-1 RMySQL_0.10.15
[59] tensorA_0.36.1 tools_3.4.4
[61] gplots_3.0.1 RColorBrewer_1.1-2
[63] proxy_0.4-22 stabledist_0.7-1
[65] ggridges_0.5.1 Rcpp_0.12.18
[67] plyr_1.8.4 base64enc_0.1-3
[69] progress_1.2.0 zlibbioc_1.24.0
[71] purrr_0.2.5 RCurl_1.95-4.11
[73] prettyunits_1.0.2 rpart_4.1-13
[75] pbapply_1.3-4 zoo_1.8-4
[77] cluster_2.0.7-1 RSpectra_0.13-1
[79] data.table_1.11.6 lmtest_0.9-36
[81] RANN_2.6 mvtnorm_1.0-8
[83] fitdistrplus_1.0-11 gsl_1.9-10.3
[85] aroma.light_3.8.0 hms_0.4.2
[87] lsei_1.2-0 xtable_1.8-3
[89] XML_3.98-1.16 mclust_5.4.1
[91] gridExtra_2.3 compiler_3.4.4
[93] tibble_1.4.2 KernSmooth_2.23-15
[95] crayon_1.3.4 R.oo_1.22.0
[97] htmltools_0.3.6 segmented_0.5-3.0
[99] pcaPP_1.9-73 Formula_1.2-3
[101] snow_0.4-3 tidyr_0.8.1
[103] geneplotter_1.56.0 DBI_1.0.0
[105] MASS_7.3-50 fpc_2.1-11.1
[107] boot_1.3-20 compositions_1.40-2
[109] ShortRead_1.36.1 R.methodsS3_1.7.1
[111] gdata_2.18.0 metap_1.0
[113] bindr_0.1.1 igraph_1.2.2
[115] pkgconfig_2.0.2 GenomicAlignments_1.14.2 [117] numDeriv_2016.8-1 foreign_0.8-71
[119] foreach_1.4.4 rARPACK_0.11-0
[121] annotate_1.56.2 XVector_0.18.0
[123] bibtex_0.4.2 stringr_1.3.1
[125] digest_0.6.17 tsne_0.1-3
[127] copula_0.999-18 ADGofTest_0.3
[129] softImpute_1.4 Biostrings_2.46.0
[131] htmlTable_1.12 edgeR_3.20.9
[133] curl_3.2 kernlab_0.9-27
[135] Rsamtools_1.30.0 gtools_3.8.1
[137] modeltools_0.2-22 nlme_3.1-137
[139] jsonlite_1.5 bindrcpp_0.2.2
[141] limma_3.34.9 pillar_1.3.0
[143] lattice_0.20-35 httr_1.3.1
[145] DEoptimR_1.0-8 survival_2.42-6
[147] glue_1.3.0 png_0.1-7
[149] prabclus_2.2-6 iterators_1.0.10
[151] glmnet_2.0-16 bit_1.1-14
[153] class_7.3-14 stringi_1.2.4
[155] mixtools_1.1.0 blob_1.1.1
[157] doSNOW_1.0.16 latticeExtra_0.6-28
[159] caTools_1.17.1.1 memoise_1.1.0
[161] irlba_2.3.2 ape_5.2

yeroslaviz commented 4 years ago

Has someone solved this problem?

I have the same error with my data.