sandhya212 / biscuit_singlecell_imm_icml_2016 Goto Github PK

Hi Sandhya,

Please see below error:

for(r_all in 1:runs.all){

• print(paste("Beginning of batch ",r_all));
  

• results.per.MCMC <- foreach(r=((1+(num_gene_sub_batches*(r_all-1))):(num_gene_sub_batches*r_all)),.packages=c("Rtsne","lattice","MASS","bayesm","robustbase","chron","mnormt","MCMCpack", "coda","Matrix", "mvtnorm")) %dopar% {
  

•     values <- list()
  

•     values <- main.fun(r)
  

•     return(list(unlist(values[[1]]),unlist(values[[2]]),unlist(values[[3]]),values[[4]],values[[5]],values[[6]])) ##z, alpha, beta, mu and Sigma
  

• }
  

• results.all.MCMC[[r_all]] <- results.per.MCMC
• print(paste("End of batch ",r_all));
• }
  [1] "Beginning of batch 1"
  Error in { : task 10 failed - "missing value where TRUE/FALSE needed"

I am not sure how to resolve this error. Thanks.

Dear Developer,

I attempt to run BISCUIT on the scRNA-seq data of infiltrating immune cells to solid tumor obtained from two patients. 1303 and 637 cells are obtained from patient #1 and #2, respectively.

Cells are divided into 8 clusters by BISCUIT. One cluster is composed of cells from only patient #2. I think that batch effect remains in spite of normalization by BISCUIT.

I set parameter alpha = 1.0 and select variable 600 genes for analysis. Could you teach me how the batch effect is canceled if there is any information?

Hi, sandhya,
I wonder if it would be appropriate to plug in Biscuit normalized and imputed data to Seurat for downstream analysis?
Thanks!

Two questions.

1.What should be the correct format of input file if "z_true_labels_avl" set be TRUE?
First 4 lines of your example file 'expression_mRNA_17-Aug-2014.txt' is the cell annotation, is that will be recognized as the 'true labels of cells' mentioned in readme file?

2. I know the gene id through "gene_selected.csv" one of output files,but what about the cell id selected by program? The "post_imputed_tSNE_coord.txt" only provide the tSNE coordinates of cells but I need corresponding cell-id,how can I find it ?

Thanks.

Hello,
I am getting the following error messages when I try to load my data set:

Error in read.table(file = file, header = header, sep = sep, quote = quote, :
duplicate 'row.names' are not allowed

and;

Error in plot.window(...) : need finite 'xlim' values
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf

I have attached an exert of the data set so that you can see the format. Can you let me know the correct format?
Also when I run Biscuit do I need to put these conditions or is it optional?

is_format_genes_cells <- FALSE

choose_cells <- 300
choose_genes <- 150
gene_batch <- 50;
num_iter <- 20;
z_true_labels_avl <- TRUE
num_cells_batch <- 100;
alpha <- 1;
Thank you!

Insiya
final_matrix.txt

i am running this on a cluster with 8 cores and 4 GB RAM per core. it's about 3K cells and about 10K genes

after 3 days, Rprof.out is 4.3 GB and every line is like this:

$ tail -3 Rprof.out
:1641588:41371563:101158008:0:"solve.default" "solve" "standardGeneric" "solve" "IMM.MCMC" "main.fun" "eval" "eval" "doTryCatch" "tryCatchOne" "tryCatchList" "tryCatch" "fun" "do.call" "docall" "doTryCatch" "tryCatchOne" "tryCatchList" "tryCatch" "doTryCatch" "tryCatchOne" "tryCatchList" "tryCatch" "slaveLoop" "eval" "eval" "eval" "eval" "eval.parent" "local"
:1641588:41371563:101158008:0:"solve.default" "solve" "standardGeneric" "solve" "IMM.MCMC" "main.fun" "eval" "eval" "doTryCatch" "tryCatchOne" "tryCatchList" "tryCatch" "fun" "do.call" "docall" "doTryCatch" "tryCatchOne" "tryCatchList" "tryCatch" "doTryCatch" "tryCatchOne" "tryCatchList" "tryCatch" "slaveLoop" "eval" "eval" "eval" "eval" "eval.parent" "local"
:1641588:41371563:101158008:0:"solve.default" "solve" "standardGeneric" "solve" "IMM.MCMC" "main.fun" "eval" "eval" "doTryCatch" "tryCatchOne" "tryCatch

what else can i check? none of the other files have been recently modified since the first day.

I'm getting an error running the code on the sample data at line 133 in BISCUIT_IMM_Gibbs_MCMC_parallel.R:

path <- working_path;

working_path is not defined. I replaced it with

path=getwd()

and it seems to run. Just wanted to let you know and make sure that fix doesn't break anything else.

One other thing is that in the final line of BISCUIT_main.R, f1 is not defined.

Error in socketSelect(socklist) : not a list of sockets???

Hi,
Just wondering what's the typical memory usage for biscuit on a dataset with around 4000 genes and 3000 cells? I've been running into memory issues all the time.
Thanks!

Every time I run the code, when the program reaches here "Beginning of Batch 1", and after a few hours it will end up with this "Error in unserialize(socklist[[n]]) : 'connection' must be a connection"

but at the same time, the parallel processes are still running,

so I don't know if the error actually put an end to the running program.
By the way, the cores to use paramater in "start_file.R" seems make no difference. I set it to 18 this time but it still consumes nearly 100 cores while parallel processing.

Hi,
Is there any way that the 10x data can be directly integrate to this pipeline or format the 10x data for biscuit.
Thanks,
Shams

Hello,

I am attempting to run Biscuit analysis on the Zeisel data before moving on to our other dataset. I was wondering if I have the settings correct for the Start File. From your paper you mention setting choose_cells <- 3005 and choose_genes <- 558 and in the Readme you say to set "alpha = 0.005 or less" Are these the only changes that need to be made?

I've attempted to run the data with these settings using 16 cores with 121G of Memory. The program runs to completion but I receive a warning messages. I was wondering if the messages (below) are to be expected and if not how I can correct them? Also is the time to be expected for this set up?

Thank you,

Jared

[1] 3005
Time difference of 2.497088 hours
There were 50 or more warnings (use warnings() to see the first 50)

warnings()
Warning messages:
1: In write.table(df_indicator, file = paste0(getwd(), "/", ... :
appending column names to file
2: In write.table("Cells selected are ", file = paste0(getwd(), ... :
appending column names to file
3: In write.table(((num_cells_batch * divsr + 1):(num_cells_batch * ... :
appending column names to file
4: In write.table(f, file = paste0(getwd(), "/debug_CM.txt"), ... :
appending column names to file
5: In write.table(f, file = paste0(getwd(), "/debug_CM.txt"), ... :
appending column names to file
6: In write.table(f, file = paste0(getwd(), "/debug_CM.txt"), ... :
appending column names to file
7: In write.table(f, file = paste0(getwd(), "/debug_CM.txt"), ... :
appending column names to file
8: In write.table(f, file = paste0(getwd(), "/debug_CM.txt"), ... :
appending column names to file
9: Quick-TRANSfer stage steps exceeded maximum (= 150250)
10: In write.table(df_indicator, file = paste0(getwd(), "/", ... :
appending column names to file
11: In write.table("Cells selected are ", file = paste0(getwd(), ... :
appending column names to file
12: In write.table(((num_cells_batch * divsr + 1):(num_cells_batch * ... :
appending column names to file
13: In sqrt(varcov) : NaNs produced
14: In sqrt(varcov) : NaNs produced
15: In sqrt(varcov) : NaNs produced

Hi dear:
I have encountered some problems when I run Biscuit. I need your help sincerely.
I set the parameter as following picture.

When I running Biscuit,I got the following output.

Finally, it have the following error.

The R version is 3.5.1. When I runing small dataset(3000cell,150gene),it did not have such error and run successfully.

How can I deal with such problem? Is there someting wrong with the parameter?

Thank you very much. Hope your answer.

Parry

I cant' seem to be able to run start_file.R using the data provided in the file expression_mRNA_17-Aug-2014.txt. Here is the error:

Error in read.table(file = file, header = header, sep = sep, quote = quote, : duplicate 'row.names' are not allowed

Hello,

How do you compute the F-measure shown in Figure 6 of your paper? In the supplementary text you talk about creating "a binary vector h encoding the ’true’ z of each cell". How do you do this binarization?

When I performed BISCUIT in the step as following, the code gave rise to an error "Error in as.vector(data) : no method for coercing this S4 class to a vector". Would you help me with this problem?
Thank you!

"""
......
[1] "Monitor debug_CM.txt for parallel process diagnostics and log_CM.txt for MAP progression"
[1] "Time for creating the MAP estimates"
Time difference of 6.888685 secs
[1] "Resample moments per cluster"
[1] "cluster is 1"
"""

Dear developer,

Actually I have several questions.

1. I got most of the output files and they look normal. Four clusters can be seen in the tSNE plot. Unfortunately the "cluster_posterior_probabilities.csv" shows zero for all four clusters for all cells.
   I wonder what caused this issue?
2. I really need to know the inferred cluster id for each cell. I have the "Imputed_Y_logspace.txt" and "post_imputed_tSNE_coord.txt". Is there a way to generate the cluster ids from these files? I cannot afford to rerun it, for the reason I will mention in 3.
3. It took me two weeks to run this on a dataset with 17K genes and 20K cells on a typical Linux node. Is it normal? Or any way to speed it up? By the way I saw some socket connection warnings.

Thank you very much!

Hi,

I got the following error while trying to run BISCUIT_main.R on my dataset:

Error in serialize(data, node$con) : error writing to connection

It looks like it specifically occurs during the BISCUIT_IMM_Gibbs_MCMC_parallel.R step at lines 593-604:

for(r_all in 1:runs.all){

print(paste("Beginning of batch ",r_all));

results.per.MCMC <- foreach(r=((1+(num_gene_sub_batches*(r_all-1))):(num_gene_sub_batches*r_all)),.packages=c("Rtsne","lattice","MASS","bayesm","robustbase","chron","mnormt","MCMCpack", "coda","Matrix", "mvtnorm")) %dopar% {
    values <- list()
    values <- main.fun(r)
    return(list(unlist(values[[1]]),unlist(values[[2]]),unlist(values[[3]]),values[[4]],values[[5]],values[[6]])) ##z, alpha, beta, mu and Sigma
}
results.all.MCMC[[r_all]] <- results.per.MCMC
print(paste("End of batch ",r_all));

}

It looks like it may be related to a cluster failing because of memory? It happens at gene batch 15. I used the numCores <-detectCores()-1 feature. And is there any way to solve this?

Thank you

Hello,
I am writing to ask for help or suggestion concerning the error as in topic. I am testing BISCUIT on the cell cycle related dataset. I am processing my single cell RNAseq dataset (not normalized, genes in rows, cells in columns) with parameters: choose_cells <- 200, choose_genes <- 300, gene_batch <- 100, num_iter <- 30, num_cells_batch <- 100, alpha <- 1. I am getting the error: kmeans: more cluster centers than distinct data points. I tried to vary three first parameters but each time getting the same error. Could you suggest me if the issue is due to the parameters or maybe to the character of dataset (Im using only cell cycle related genes).
Best,
Monika Krzak

Greetings,

What is the appropriate parallel environment (MPI, etc.) to use for this program?
Thanks!

Hello,
I have a question concerning MCMC convergence. I run the example data and I wish to plot the iteration versus number of clusters. With default parameters there should be three MCMC running in parallel so I assume I have to 3 plots of convergences. Is there any output that I can access the number of clusters for each iteration? I see something it in the algorithm (nr of clusters are used in the loop to produce tsne plots for each iteration) however Im not sure if I can access those values..I would be grateful for help!
Best,
Monika

Dear Sandhya,

I have question concerning the run of BISCUIT on Zeisel dataset (558 genes, 3005 cells) that you present in the article. How long it took ? Did you parallelize the run with number of gene batches , if yes how many you chose ?

Best,
Monika

Sorry

I am interested in BLOOD samples from your paper I therefore downloaded GSE114725_rna_raw.csv.gz and separated cells in BLOOD part

I named the matrix of resulted counts loaded in dropbox

https://www.dropbox.com/s/a1ef6c5e6wxo1h1/expression_mRNA_17-Aug-2014.txt?dl=0

After about an hour running, I got this error

> source("start_file.R")
[1] "Loading Data"
[1] "Calculating the Fiedler vector of the data"
 Error in full.data.1[, order(f.vec, decreasing = TRUE)] : 
  subscript out of bounds 

Thank you for any help

Hello!
I am working on a project comparing different clustering methods for single-cell, and I am using the Zeisel et al. 2015 data set. The results presented in Figure 8 look really good, and I am looking to reproduce the BISCUIT results found in Figure 8 Panel 2. Unfortunately, I can't seem to get the results using 3005 cells, 558 genes, and an alpha = 0.005. What were the parameters used to get this? When I put 558 genes in the parameters, does it select the same 558 genes as the ones used to generate Figure 8's results?
Thanks for such a well-documented repo, too! It was so easy to get up and running!

Hi,
When I run start_file.R on around 3000 cells, the above error reported. How could I continue? Thank you.
Best,
Elaine

Hello,

I am attempting a run of the BISCUIT algorithm with a gene X cell matrix, obtained from Seurat, of dimensions ~40,000 X ~60,000. I am running this through the BISCUIT docker installation on Ubuntu 18.04 with a 32 core instance containing 256 GB of RAM. I am receiving the error below after the load data step runs for ~24 hours. It seems as if I need more RAM. I am able to run the test data and have no issues. Should I just increase the RAM of the instance or do you have any other recommendations? Can you provide any suggestions to resolve this error?

Thanks,
Michael

Attaching package: 'snow'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
    parCapply, parLapply, parRapply, parSapply, splitIndices,
    stopCluster


Attaching package: 'bayesm'

The following object is masked from 'package:gtools':

    rdirichlet

The following object is masked from 'package:MCMCpack':

    rdirichlet


Attaching package: 'chron'

The following object is masked from 'package:foreach':

    times

[1] "Loading Data"
[1] "Calculating the Fiedler vector of the data"
Error: cannot allocate vector of size 7.2 Gb

UPDATE: I re-ran this matrix with the following settings (mostly default) and was able to get past this initial error:

input_data_tab_delimited <- TRUE; #set to TRUE if the input data is tab-delimited

is_format_genes_cells <-  TRUE; #set to TRUE if input data has rows as genes and columns as cells

choose_cells <- 3000; #comment if you want all the cells to be considered

choose_genes <- 150; #comment if you want all the genes to be considered

gene_batch <- 50; #number of genes per batch, therefore num_batches = choose_genes (or numgenes)/gene_batch. Max value is 150

num_iter <- 20; #number of iterations, choose based on data size.

num_cores <- detectCores() - 4; #number of cores for parallel processing. Ensure that detectCores() > 1 for parallel processing to work, else set num_cores to 1.

z_true_labels_avl <- FALSE; #set this to TRUE if the true labels of cells are available, else set it to FALSE. If TRUE, ensure to populate 'z_true' with the true labels in 'BISCUIT_process_data.R'

num_cells_batch <- 1000; #set this to 1000 if input number of cells is in the 1000s, else set it to 100.

alpha <- 1; #DPMM dispersion parameter. A higher value spins more clusters whereas a lower value spins lesser clusters.

output_folder_name <- "output"; #give a name for your output folder.

I then received a different error:

1] "Loading Data"
[1] "Calculating the Fiedler vector of the data"
[1] "Ensuring entire data is numeric and then log transforming it"
[1] "numcells is 3000"
[1] "numgenes is 150"
[1] "Number of gene batches is 3"
[1] "Number of gene subbatches is 3"
[1] "Ensuring user-specified data is numeric"
[1] "Computing t-sne projection of the data"
[1] "Monitor log.txt and outputs/plots/ folder for outputs"
[1] "floor(num_gene_batches/num_gene_sub_batches): 1"
[1] "MCMC begins"
[1] "Begin parallel processing of gene splits"
[1] "Beginning of batch  1"
Error in serialize(data, node$con) : error writing to connection

Do you have any suggestions to resolve this issue? Do you have any recommendations for settings so that I can include all cells and genes in the matrix (~40,000 X ~60,000) for the analysis?

Hello!

I am trying BISCUIT on a data, however got the following error:

[1] "Loading Data"
[1] "Calculating the Fiedler vector of the data"
[1] "Ensuring entire data is numeric and then log transforming it"
[1] "numcells is 74"
[1] "numgenes is 9377"
[1] "Number of gene batches is 62"
[1] "Number of gene subbatches is 2"
[1] "Ensuring user-specified data is numeric"
[1] "Computing t-sne projection of the data"
Error in .check_tsne_params(nrow(X), dims = dims, perplexity = perplexity,  : 
  perplexity is too large for the number of samples
In addition: There were 19 warnings (use warnings() to see them)

Data comes from Grun et al 2014.

The parameter setting is as following:

## 21st Dec 2016
## BISCUIT R implementation
## Start_file with user inputs
## 
## Code author SP


###
###

############## packages required ##############

library(MCMCpack)
library(mvtnorm)
library(ellipse)
library(coda)
library(Matrix)
library(Rtsne)
library(gtools)
library(foreach)
library(doParallel)
library(doSNOW)
library(snow)
library(lattice)
library(MASS)
library(bayesm)
library(robustbase)
library(chron)
library(mnormt)
library(schoolmath)
library(RColorBrewer)

#############################################


#input_file_name <- "expression_mRNA_17-Aug-2014.txt";

input_data_tab_delimited <- TRUE; #set to TRUE if the input data is tab-delimited

is_format_genes_cells <-  TRUE; #set to TRUE if input data has rows as genes and columns as cells

#choose_cells <- 3000; #comment if you want all the cells to be considered

#choose_genes <- 150; #comment if you want all the genes to be considered

gene_batch <- 150; #number of genes per batch, therefore num_batches = choose_genes (or numgenes)/gene_batch. Max value is 150

num_iter <- 5; #number of iterations, choose based on data size.

num_cores <- detectCores() - 4; #number of cores for parallel processing. Ensure that detectCores() > 1 for parallel processing to work, else set num_cores to 1.

z_true_labels_avl <- FALSE; #set this to TRUE if the true labels of cells are available, else set it to FALSE. If TRUE, ensure to populate 'z_true' with the true labels in 'BISCUIT_process_data.R'

num_cells_batch <- 1000; #set this to 1000 if input number of cells is in the 1000s, else set it to 100.

alpha <- 0.1; #DPMM dispersion parameter. A higher value spins more clusters whereas a lower value spins lesser clusters.

#output_folder_name <- "output"; #give a name for your output folder.

## call BISCUIT
source("BISCUIT_main.R")

The data I used can be found here: https://github.com/WT215/Raw_data (Grun_2i.txt)

Thank you for your help!

Best wishes,
Wenhao

I am trying to replicate the experiments done in section 6.1 of your paper.

1. I have a question about the identifiability constraints you put on the alpha and beta scaling parameters: how do you decide what values of maximum and minimum to impose on those parameters for synthetic data? I see that in your code, for real gene expression data, you use the library size in order to make this decision. But what about in synthetic data generated according to your model?

2. What was the run time for fitting BISCUIT and HDPMM to the generated data with 100 samples, 50 dimensions and 3 components? Did you enforce some kind of regularization on the covariance matrices to reduce complexity, like making them diagonal? Inverting high-dimensional matrices takes quite a while...

3. In your code I don't see how you are using the Hierarchical Conditionally Conjugate DPMM model with cell-specific scalings. Where are the hyperparameter updates of psi, mu_prime, Sigma_prime, H_prime and sigma_prime, as in the plate model in Figure 4??

Hello!
I keep encountering this error when I run Biscuit. I'm relatively new to coding so it maybe something small that I'm missing or haven't done.
I set the parameters as the following:

And when I ran Biscuit I got the following output and error:

I have R version 4.0.2 but ran the program in RStudio.

Please let me know how I can fix this issue or if I'm doing something incorrectly.

Thank you!
Rhea

Hi,

I am trying to re-analyze a subset of public scRNA-Seq data from Chung et al breast cancer paper. Attached is my (expected) count matrix file sce_immune_subset.txt for 65 mostly immune cells that I am interested in.

From other clustering methods, e.g. SC3 or kmeans, I am getting about 4-5 clusters: macrophages, T cells, B cells, maybe some pDCs and some unclassified. This is also roughly what the authors report in their paper.

When I try BISCUIT on this data set I am getting 60+ clusters depending on various options. I tried choose_genes to include either all genes or just top 1024. I tried alpha=0.01 as well. Shown below is my latest output, BISCUIT also aborts at kmeans step. Do I do everything correctly here?

[1] "Loading Data"
[1] "Calculating the Fiedler vector of the data"
[1] "Ensuring entire data is numeric and then log transforming it"
[1] "numcells is 65"
[1] "numgenes is 1024"
[1] "Number of gene batches is 16"
[1] "Number of gene subbatches is 8"
[1] "Ensuring user-specified data is numeric"
[1] "Computing t-sne projection of the data"
[1] "Monitor log.txt and outputs/plots/ folder for outputs"
[1] "floor(num_gene_batches/num_gene_sub_batches): 2"
[1] "MCMC begins"
[1] "Begin parallel processing of gene splits"
[1] "Beginning of batch 1"
[1] "End of batch 1"
[1] "Beginning of batch 2"
[1] "End of batch 2"
[1] "End of parallel runs"
Time difference of 8.947547 mins
[1] "Merging gene splits"
[1] "Computing the global confusion matrix"
[1] "Monitor log_CM.txt in outputs folder and debug_CM.txt"
[1] "Time for creating the upper triangular global confusion matrix is "
Time difference of 7.58055 secs
[1] "Time for creating the overall upper triangular global confusion matrix is "
Time difference of 7.581766 secs
[1] "Number of clusters over gene splits: 57" "Number of clusters over gene splits: 63"
[3] "Number of clusters over gene splits: 61" "Number of clusters over gene splits: 62"
[5] "Number of clusters over gene splits: 62" "Number of clusters over gene splits: 59"
[7] "Number of clusters over gene splits: 60" "Number of clusters over gene splits: 56"
[9] "Number of clusters over gene splits: 57" "Number of clusters over gene splits: 60"
[11] "Number of clusters over gene splits: 46" "Number of clusters over gene splits: 56"
[13] "Number of clusters over gene splits: 51" "Number of clusters over gene splits: 51"
[15] "Number of clusters over gene splits: 54" "Number of clusters over gene splits: 55"
[1] "Mean number of clusters over gene splits: 63"
[1] "Cluster the global CM by kmeans"

Error in kmeans(ave_CM, mean_num_clusters, nstart = 1) :
more cluster centers than distinct data points.