shenorrlabtrdf / bseqsc Goto Github PK
View Code? Open in Web Editor NEWBulk-Sequence Single-Cell Gene Expression Deconvolution Pipeline
Bulk-Sequence Single-Cell Gene Expression Deconvolution Pipeline
When I tried to use my own data, I got the bug information. I tested using the sample data, it worked. But my data was failed to run.
fit <- bseqsc_proportions(bulk, B, verbose = TRUE)
Rerun with Debug
Error in predict.svm(ret, xhold, decision.values = TRUE) :
Model is empty!
Traceback:
12.
stop("Model is empty!")
11.
predict.svm(ret, xhold, decision.values = TRUE)
10.
predict(ret, xhold, decision.values = TRUE)
9.
na.action(predict(ret, xhold, decision.values = TRUE))
8.
svm.default(X, y, type = "nu-regression", kernel = "linear",
nu = nus, scale = F)
7.
svm(X, y, type = "nu-regression", kernel = "linear", nu = nus,
scale = F) at CIBERSORT.R#58
6.
FUN(X[[i]], ...)
5.
lapply(X, FUN, ...)
4.
mclapply(1:svn_itor, res, mc.cores = 1) at CIBERSORT.R#62
3.
CoreAlg(X, y, absolute, abs_method) at CIBERSORT.R#200
2.
CIBERSORT(xf, yf, ...)
1.
bseqsc_proportions(bulk, B, verbose = TRUE)
In the tutorial for how to do this we see:
eset <- readRDS('https://shenorrlab.github.io/bseqsc/data/GSE50244.rds')
but this file doesn't exist. Can you post it please.
Hi,
It would be great to get access to the bseq-sc code, which was used in A single-cell transcriptomic map of the human and mouse pancreas reveals inter- and intra-cell population structure (PMID:27667365). The paper was published few months ago and the code is not available yet.
Thank you very much for your attention and your great work!
Best,
Javier Perales-Paton
Hi,
The data for the tutorial is missing (ERROR 404):
http://shenorrlab.github.io/
BSeq-sc
/data/[...]
, however this direction does not exist.https://shenorrlab.github.io/
bseqsc/articles/pages
/data.html
), which contains two links for the two data sets. However, there is no /articles
directory in the repository.It seems that you have mixed two github webs (BSeq-sc
and bseqsc
). Or somehow you have forgetten to upload the data files.
The tutorial looks promising! Thank you very much!!
Best,
Javier
Thanks for developing this software!
I have successfully downloaded the test data, but i get an error when i run the pipeline from your tutorial on the sample data sets. It said:
`fit <- bseqsc_proportions(eset, B, verbose = TRUE)
404
There isn't a GitHub Pages site here.
If you're trying to publish one, read the full documentation to learn how to set up GitHub Pages for your repository, organization, or user account.
GitHub Status — @githubstatus
hi developers,
just try to run your pipeline getting following error:
edger_pvals <- mutate(edger_pvals, Regulated = Adjusted <= 0.01 & Adjusted <= Base)
Error in mutate(edger_pvals, Regulated = Adjusted <= 0.01 & Adjusted <= :
could not find function "mutate"
Hi there,
I wish to deconvolve my bulk RNA seq data with single cell RNA seq data from the same tissue. I run in to the following error when running the code line below and I do not understand why. Can any one please help?
csfit <- bseqsc_csdiff(ilc_bulk[genes, ] ~ treatment | fibroblast + macrophage,
verbose = 2, nperms = 100, .rng = 12345)
Then comes the error:
Groups: a_nacl1=4L | a_nacl2=4L | a_nacl3=4L | a_nacl4=4L | a_nacl5=4L | NA0L
Cell type(s): 'fibroblast', 'macrophage' (2 total)
Fitting mode: auto
Data (filtered): 1684 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
Fitting linear interaction model ... OK
Computing FDR using 100 permutations ... 101/100
Alternative 'two.sided' ... OK
Alternative 'greater' ... OK
Alternative 'less' ... OK
OK
Timing:
user system elapsed
2.136 0.153 2.267
Error in array(NA_real_, dim = c(nrow(b), ncol(cc), length(lev)), dimnames = c(rownames(b), :
length of 'dimnames' [1687] must match that of 'dims' [3]
> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=de_DE.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=de_DE.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] edgeR_3.26.6 limma_3.40.6 forcats_0.4.0 stringr_1.4.0
[5] purrr_0.3.2 readr_1.3.1 tidyr_0.8.3 tibble_2.1.3
[9] tidyverse_1.2.1 SoupX_0.3.1 ggplot2_3.2.0 dplyr_0.8.3
[13] DropletUtils_1.4.3 SingleCellExperiment_1.6.0 SummarizedExperiment_1.14.0 DelayedArray_0.10.0
[17] BiocParallel_1.18.0 matrixStats_0.54.0 GenomicRanges_1.36.0 GenomeInfoDb_1.20.0
[21] RColorBrewer_1.1-2 xbioc_0.1.17 AnnotationDbi_1.46.0 IRanges_2.18.1
[25] S4Vectors_0.22.0 BisqueRNA_1.0 Seurat_3.0.2 preprocessCore_1.46.0
[29] e1071_1.7-2 bseqsc_1.0 csSAM_1.4 Rcpp_1.0.2
[33] openxlsx_4.1.0.1 Biobase_2.44.0 BiocGenerics_0.30.0
loaded via a namespace (and not attached):
[1] reticulate_1.13 R.utils_2.9.0 tidyselect_0.2.5 RSQLite_2.1.2 htmlwidgets_1.3
[6] grid_3.6.1 Rtsne_0.15 devtools_2.1.0 munsell_0.5.0 codetools_0.2-16
[11] ica_1.0-2 future_1.14.0 withr_2.1.2 colorspace_1.4-1 rstudioapi_0.10
[16] ROCR_1.0-7 gbRd_0.4-11 listenv_0.7.0 NMF_0.22 Rdpack_0.11-0
[21] labeling_0.3 GenomeInfoDbData_1.2.1 bit64_0.9-7 rhdf5_2.28.0 rprojroot_1.3-2
[26] vctrs_0.2.0 generics_0.0.2 R6_2.4.0 doParallel_1.0.14 rsvd_1.0.2
[31] locfit_1.5-9.1 bitops_1.0-6 assertthat_0.2.1 SDMTools_1.1-221.1 scales_1.0.0
[36] gtable_0.3.0 npsurv_0.4-0 globals_0.12.4 processx_3.4.1 rlang_0.4.0
[41] zeallot_0.1.0 splines_3.6.1 lazyeval_0.2.2 broom_0.5.2 modelr_0.1.4
[46] BiocManager_1.30.4 yaml_2.2.0 reshape2_1.4.3 backports_1.1.4 tools_3.6.1
[51] usethis_1.5.1 gridBase_0.4-7 gplots_3.0.1.1 sessioninfo_1.1.1 ggridges_0.5.1
[56] plyr_1.8.4 zlibbioc_1.30.0 RCurl_1.95-4.12 ps_1.3.0 prettyunits_1.0.2
[61] pbapply_1.4-1 viridis_0.5.1 cowplot_1.0.0 zoo_1.8-6 haven_2.1.1
[66] ggrepel_0.8.1 cluster_2.1.0 fs_1.3.1 magrittr_1.5 data.table_1.12.2
[71] lmtest_0.9-37 RANN_2.6.1 fitdistrplus_1.0-14 pkgload_1.0.2 hms_0.5.0
[76] lsei_1.2-0 xtable_1.8-4 readxl_1.3.1 gridExtra_2.3 testthat_2.2.1
[81] compiler_3.6.1 KernSmooth_2.23-15 crayon_1.3.4 R.oo_1.22.0 htmltools_0.3.6
[86] Formula_1.2-3 lubridate_1.7.4 DBI_1.0.0 MASS_7.3-51.4 Matrix_1.2-17
[91] cli_1.1.0 R.methodsS3_1.7.1 gdata_2.18.0 metap_1.1 igraph_1.2.4.1
[96] pkgconfig_2.0.2 registry_0.5-1 plotly_4.9.0 xml2_1.2.1 foreach_1.4.7
[101] dqrng_0.2.1 rngtools_1.4 pkgmaker_0.28 XVector_0.24.0 rvest_0.3.4
[106] bibtex_0.4.2 callr_3.3.1 digest_0.6.20 sctransform_0.2.0 tsne_0.1-3
[111] cellranger_1.1.0 dendextend_1.12.0 curl_4.0 gtools_3.8.1 nlme_3.1-141
[116] jsonlite_1.6 Rhdf5lib_1.6.0 desc_1.2.0 viridisLite_0.3.0 pillar_1.4.2
[121] lattice_0.20-38 httr_1.4.0 pkgbuild_1.0.3 survival_2.44-1.1 glue_1.3.1
[126] remotes_2.1.0 zip_2.0.3 png_0.1-7 iterators_1.0.12 bit_1.1-14
[131] class_7.3-15 stringi_1.4.3 HDF5Array_1.12.1 blob_1.2.0 caTools_1.17.1.2
[136] memoise_1.1.0 irlba_2.3.3 future.apply_1.3.0 ape_5.3
In order to accomplish downstream analysis, I have to download the CIBERSORT software source code in R.
But I have waited for approval a few days, but download page still showed "Download Link Pending Approval".
Anybody can send me or shared the CIBERSORT software source code in R ?
Hi Thank you for this great package. Is it possible to use DEseq with the package instead of edgeR?
Best wishes
Nurun
Hi. I will be very grateful for any hint on how to overcome the error. I wish to deconvolve my bulk RNA seq data from the lungs of mice using single cell RNA seq data. For practice, I am following your tutorial using my own bulk seq data and the single cell rna seq data you have provided. All seems to go well until when I try to obtain the Cell type-specific differential expression by running the following code. From the beginning I thought the presence of factors(e.g. for th treatment conditions or phenoData) in my data set was the cause but after getting rid of the factors I still get the same error:
csfit <- bseqsc_csdiff(the_best_ilc_bulk_expr_set[genes, ] ~ phenoData| alpha + beta + ductal + acinar,
verbose = 2, nperms = 5000, .rng = 12345)`
Then I get this error:
Groups: bleo=5L | bleo_ko_i=5L | bleo_wt_i=5L | Nacl=5L | NA0L
Cell type(s): 'alpha', 'beta', ..., 'acinar' (4 total)
Fitting mode: auto
Data (filtered): 1202 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
Fitting linear interaction model ... Error in dimnames(covmat.unscaled) <- list(xnames, xnames) :
length of 'dimnames' [1] not equal to array extent
In addition: Warning messages:
1: In lsfit(D, G, intercept = FALSE) : 'X' matrix was collinear
2: In sqrt(((n - p) * stddevmat^2 - resids^2/(1 - hatdiag[good]))/(n - :
NaNs produced
Here are the details of my code:
# inspired from https://github.com/cozygene/bisque/issues/4
# creating expression set from my bulk rna seq data
library(openxlsx)
library(Biobase)
library(bseqsc)
library(tidyverse)
bseqsc_config('CIBERSORT.R')
data(pancreasMarkers)
# read in the raw data
bulk <- read.xlsx("all_samples.xlsx",
sheet = "rawCounts_unrepeated_genes")
head(bulk)
bulk$Gene <- NULL
rownames(bulk) <- toupper(bulk$external_gene_name)
bulk$external_gene_name <- NULL
bulk <- as.matrix(bulk)
head(bulk)
# create expression matrix
featureData <- as.character(rownames(bulk))
featureData <- as(data.frame(featureData,
stringsAsFactors = FALSE), "AnnotatedDataFrame")
rownames(featureData) <- rownames(bulk)
phenoData <- c(rep("bleo", 5), rep("Nacl", 5), rep("bleo_wt_i", 5), rep("bleo_ko_i", 5))
phenoData <- as(as.data.frame(phenoData), "AnnotatedDataFrame")
rownames(phenoData) <- colnames(bulk)
the_best_ilc_bulk_expr_set <- ExpressionSet(assayData = bulk,
phenoData = phenoData, featureData = featureData)
# read in the single cell rna seq data and fit the model
eislet <- readRDS('islet-eset.rds')
B <- bseqsc_basis(eislet, pancreasMarkers,
clusters = 'cellType', samples = 'sampleID', ct.scale = TRUE)
fit <- bseqsc_proportions(the_best_ilc_bulk_expr_set, B, verbose = TRUE)
pData(the_best_ilc_bulk_expr_set) <- cbind(pData(the_best_ilc_bulk_expr_set), t(coef(fit)))
fit_edger<- fitEdgeR(the_best_ilc_bulk_expr_set, ~phenoData,
coef = c("phenoDatableo_ko_i", "phenoDatableo_wt_i", "phenoDataNacl"))
# extended
fit_edger_ext <- fitEdgeR(the_best_ilc_bulk_expr_set, ~ phenoData + beta + ductal + acinar +gamma,
coef = c("phenoDatableo_ko_i", "phenoDatableo_wt_i",
"phenoDataNacl", "beta", "ductal",
"acinar", "gamma"))
fit_edger_ext$Symbol <- rownames(fit_edger_ext)
# gather P-values from both models
df_fit_edger <- as.data.frame(fit_edger, stringsAsFactors = FALSE)
df_fit_edger$Symbol <- rownames(df_fit_edger)
req_df_fit_edger <- df_fit_edger["PValue"]
colnames(req_df_fit_edger) <- "Base"
head(req_df_fit_edger)
req_df_fit_edger$Symbol <- rownames(req_df_fit_edger)
df_fit_edger_ext <- as.data.frame(fit_edger_ext, stringsAsFactors = FALSE)
df_fit_edger_ext$Symbol <- rownames(df_fit_edger_ext)
head(df_fit_edger_ext)
req_df_fit_edger_ext <- df_fit_edger_ext["PValue"]
colnames(req_df_fit_edger_ext) <- "Adjusted"
req_df_fit_edger_ext$Symbol <- rownames(req_df_fit_edger_ext)
edger_pvals <- req_df_fit_edger_ext %>%
inner_join(req_df_fit_edger)
head(edger_pvals)#
rownames(edger_pvals) <- edger_pvals$Symbol
edger_pvals <- mutate(edger_pvals, Regulated = Adjusted <= 0.01 & Adjusted <= Base)
# plot Base vs Adjusted
pvalueScatter(Base ~ Adjusted, edger_pvals, pval.th = 0.01, label.th = 3.5)
# ER genes
genes_ER <- c('HSPA5', 'MAFA', 'HERPUD1', 'DDIT3', 'UCN3', 'NEUROD1')
# Fit on ER stress genes and genes regulated beyond cell type proportion differences
genes <- union(genes_ER, subset(edger_pvals, Regulated)$Symbol)
csfit <- bseqsc_csdiff(the_best_ilc_bulk_expr_set[genes, ] ~ phenoData| alpha + beta + ductal + acinar,
verbose = 2, nperms = 5000, .rng = 12345)
Then I have the error:
Groups: bleo=5L | bleo_ko_i=5L | bleo_wt_i=5L | Nacl=5L | NA0L
Cell type(s): 'alpha', 'beta', ..., 'acinar' (4 total)
Fitting mode: auto
Data (filtered): 1202 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
Fitting linear interaction model ... Error in dimnames(covmat.unscaled) <- list(xnames, xnames) :
length of 'dimnames' [1] not equal to array extent
In addition: Warning messages:
1: In lsfit(D, G, intercept = FALSE) : 'X' matrix was collinear
2: In sqrt(((n - p) * stddevmat^2 - resids^2/(1 - hatdiag[good]))/(n - :
NaNs produced
> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=de_DE.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=de_DE.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] edgeR_3.26.6 limma_3.40.6 preprocessCore_1.46.0 e1071_1.7-2 forcats_0.4.0 stringr_1.4.0
[7] dplyr_0.8.3 purrr_0.3.2 readr_1.3.1 tidyr_0.8.3 tibble_2.1.3 ggplot2_3.2.0
[13] tidyverse_1.2.1 bseqsc_1.0 csSAM_1.4 Rcpp_1.0.2 Biobase_2.44.0 BiocGenerics_0.30.0
[19] openxlsx_4.1.0.1
loaded via a namespace (and not attached):
[1] nlme_3.1-140 lubridate_1.7.4 bit64_0.9-7 doParallel_1.0.14 RColorBrewer_1.1-2 httr_1.4.0
[7] tools_3.6.0 backports_1.1.4 R6_2.4.0 DBI_1.0.0 lazyeval_0.2.2 colorspace_1.4-1
[13] withr_2.1.2 tidyselect_0.2.5 gridExtra_2.3 xbioc_0.1.17 bit_1.1-14 compiler_3.6.0
[19] cli_1.1.0 rvest_0.3.4 xml2_1.2.1 pkgmaker_0.28 scales_1.0.0 NMF_0.22
[25] digest_0.6.20 pkgconfig_2.0.2 bibtex_0.4.2 rlang_0.4.0 readxl_1.3.1 rstudioapi_0.10
[31] RSQLite_2.1.2 generics_0.0.2 jsonlite_1.6 dendextend_1.12.0 zip_2.0.3 magrittr_1.5
[37] Formula_1.2-3 munsell_0.5.0 S4Vectors_0.22.0 viridis_0.5.1 stringi_1.4.3 yaml_2.2.0
[43] plyr_1.8.4 grid_3.6.0 blob_1.2.0 crayon_1.3.4 lattice_0.20-38 haven_2.1.1
[49] splines_3.6.0 hms_0.5.0 locfit_1.5-9.1 zeallot_0.1.0 pillar_1.4.2 rngtools_1.4
[55] reshape2_1.4.3 codetools_0.2-16 stats4_3.6.0 glue_1.3.1 BiocManager_1.30.4 modelr_0.1.4
[61] vctrs_0.2.0 foreach_1.4.7 cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1 gridBase_0.4-7
[67] xtable_1.8-4 broom_0.5.2 class_7.3-15 viridisLite_0.3.0 iterators_1.0.12 AnnotationDbi_1.46.0
[73] registry_0.5-1 memoise_1.1.0 IRanges_2.18.1 cluster_2.1.0
Thanks for developing this software!
I have successfully downloaded and run the pipeline from your tutorial on the sample data sets. When I try to perform the pipeline on my own data I get an error during the 'bseqsc_proportions' functions:
> myfit <- bseqsc_proportions(bulkSet, scB, log=TRUE, verbose = TRUE)
* Data features: 'DPM1', 'SCYL3', ..., 'NEP' (11,561 total)
* Basis features: 'HA', 'RPL39', ..., 'RPL9' (38 total)
* Common features: 'GCLM', 'RPL6', ..., 'M2' (38 total)
* Converting to linear scale
* Writing input files ... OK
* Running CIBERSORT ... Error in if (max(Y) < 50) { : missing value where TRUE/FALSE needed
My 'bulkSet' has 11561 features, 957 samples and contains the corresponding expression matrix of log transformed TPM values. I am assuming there is something wrong with this ExpressionSet that I generated and/or the corresponding expression matrix, but the error message isn't too helpful in helping me identify what that might be. Any ideas? Thanks in advance.
Hi
What kind of pre-processing do you do in the scRNA-seq and RNA-seq? Apart from having consistent gene_ids and no duplicated rows, do you do any normalisation? And if so, which one?
I downloaded pre-processed scRNA-seq data (TPM normalised) folowed the workflow to generate the B matrix. Seemed to work. Then I used VST transformed RNA-seq data and ran Cibersort.. but that didn't work very well.. Maybe I am doing something wrong?
Hi there,
I just started working with your package, but the CIBERSORT dependency has recently been replaced with CIBERSORTx, which has self-reported to have greater capacity for cell-type deconvolution. Will we have to continue to use the classic CIBERSORT package, or is there be an upgrade in the software coming soon?
Hi,
Bseq_Sc is a great tool. I like it!
However, there is a problem with fitEdgeR, if I want to perform analysis using five or more cell types, fitEdgeR return following error;
"Error in glmFit.default(sely, design, offset = seloffset, dispersion = 0.05, :
Design matrix not of full rank. The following coefficients not estimable:
Microglia # (the last of coefficients )"
Nevertheless, its work pretty good with four or less cell types. Is there any bug? Do you have any advice for the analysis with five or more cell types?
Best
Gürkan
We are currently working on putting bseqsc on bioconda, but for that we would need you to create a release. Would that be possible?
Thank you very much in advance!
we are testing our data with bseqsc, while occurs the errror:
> fit <- bseqsc_proportions(skinbulkset, B, verbose = TRUE)
* Data features: 'MAP2K6', 'HK2', ..., 'RP11-412H9.2' (11,262 total)
* Basis features: 'LMX1A', 'CACNA2D1', ..., 'COL4A1' (198 total)
* Common features: 'DSCAML1', 'CACNA1B', ..., 'C12orf55' (133 total)
* Converting to linear scale
* Writing input files ... OK
* Running CIBERSORT ... Error in out[[t]]$coefs : $ operator is invalid for atomic vectors
In addition: Warning message:
In parallel::mclapply(1:svn_itor, res, mc.cores = svn_itor) :
all scheduled cores encountered errors in user code
Any advice provided is appreciated.
hello!How to normalize the gene expression of bulk RNA seq and scRNA seq and calculate the fpkm value?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.