SymSim

Overview

SymSim is a method simulating single cell data. I find it comprehensive, generative and mimic the experimental procedure.

github, tutorial

Considers three variations: extrinsic variation(cell type, cell state), intrinsic variation(promoter on/off, mRNA synthesis rate, degradation), and technical variation(library preparation and sequencing).

Geneate counts: 1. generate transcript counts; 2. generate observed expressions.

The true transcript counts are generated using classical promoter kinetic model, with three parameters: promotoer on rate(\(k_{on}\)), promoter off rate(\(k_{off}\)) and RNA synthesis rate(\(s\)). The value of these parameters are determined by the product of gene effects and cell specific extrinsic variability EVF(indicate the cell state, low dimensional manifold). EVF value is determined by a tree stucture(dictates cell-cell similarity map). For homogenous population(a single location in a tree, like one cell type at a certain state/time), EVFs are drawn iid from a distribution.

Generate transcript counts

Two-state kinetic model: promoter switches between an on and an off states with certain probabilities. Let \(s\) be the transcript rate and \(d\) be the degradation rate fixed at 1. The stationary distribution for each gene analytically is a Beta-Poisson mixture.

Suppose there are \(m\) genes and \(n\) cells. Then there are three \(m\times n\) matrices: \(K_{on}\), \(K_{off}, S\), each is a product of two low rank matrix. When the gene is on it is transcribed with transcription rate \(s\). Given the three matrices, transcript counts are drawn from a beta-Poisson model.

Draw transcript counts: for each gene in a cell, 1. generate \(p\) from \(Beta(k_{on},k_{on})\); 2, generate transcript counts from \(Poisson(p*s)\)

Generate observed expressions

cell efficiency, amplification bias, fragmentation, sequencing.

library(SymSim)
cal_amp_bias <- function(lenslope, nbins, gene_len, amp_bias_limit){
  
  ngenes <- length(gene_len)
  len_bias_bin <- (-c(1:nbins))*lenslope
  len_bias_bin <- len_bias_bin-median(len_bias_bin)
  if (max(len_bias_bin) > amp_bias_limit[2]) {
    stop("The lenslope parameter is too large.")
  }
  max_rand_bias <- amp_bias_limit[2] - max(len_bias_bin)
  
  rand_bias <- rnorm(ngenes, mean=0, sd=max_rand_bias)
  rand_bias[rand_bias > max_rand_bias] <- max_rand_bias
  rand_bias[rand_bias < -max_rand_bias] <- -max_rand_bias
  #rand_bias <- runif(ngenes, -max_rand_bias,  max_rand_bias)
  
  binsize <- floor(ngenes/nbins)
  genes_in_bins <- vector("list", nbins)
  bin4genes <- numeric(ngenes)
  for (ibin in 1:(nbins-1)){
    genes_in_bins[[ibin]] <- order(gene_len)[((ibin-1)*binsize+1) : (ibin*binsize)]
    bin4genes[genes_in_bins[[ibin]]] <- ibin
  }
  genes_in_bins[[nbins]] <- order(gene_len)[((nbins-1)*binsize+1) : ngenes]
  bin4genes[genes_in_bins[[nbins]]] <- nbins
  
  len_bias <- numeric(ngenes); len_bias <- len_bias_bin[bin4genes]
  amp_bias <- rand_bias+len_bias
  return(amp_bias)
}

expand2binary <- function(true_counts_1cell){
  expanded_vec <- rep(1, sum(true_counts_1cell))
  trans_idx <- sapply(which(true_counts_1cell>0), 
                      function(igene){return(rep(igene, true_counts_1cell[igene]))})
  trans_idx <- unlist(trans_idx)
  return(list(expanded_vec, trans_idx))
}


data("gene_len_pool")
true_counts_res <- SimulateTrueCounts(ncells_total=100, ngenes=100, evf_type="one.population", Sigma=0.4, randseed=0)
true_counts_1cell = true_counts_res$counts[,1]
protocol = 'nonUMI' 
rate_2cap= 0.1
gene_len = sample(gene_len_pool,100)
amp_bias = cal_amp_bias(0.02, 20, gene_len, c(-0.2, 0.2))
rate_2PCR=0.8 
nPCR1=16 
nPCR2=10 
LinearAmp=FALSE 
N_molecules_SEQ = 1e4

ngenes <- length(gene_len)
if (protocol=="nonUMI"){data(len2nfrag)} 
inds <- vector("list",2)
expanded_res <- expand2binary(c(true_counts_1cell,1))
expanded_vec <- expanded_res[[1]]; trans_idx <- expanded_res[[2]]
inds[[1]] <- which(expanded_vec > 0); expanded_vec <- expanded_vec[inds[[1]]]
trans_idx <- trans_idx[inds[[1]]]

captured_vec <- expanded_vec; captured_vec[runif(length(captured_vec)) > rate_2cap] <- 0
captured_vec[length(captured_vec)] <- 1
inds[[2]] <- which(captured_vec > 0); captured_vec <- captured_vec[inds[[2]]]
trans_idx <- trans_idx[inds[[2]]]
amp_rate <- c((rate_2PCR+amp_bias[trans_idx[1:(length(trans_idx)-1)]]),1)

## what does this step do?

temp <- runif(length(captured_vec)) < amp_rate
    temp <- temp*2+captured_vec-temp
    for (iPCR in 2:nPCR1){
      eff <- runif(length(temp))*amp_rate
      v1 <- temp*(1-eff)
      round_down <- (v1-floor(v1)) < runif(length(v1))
      v1[round_down] <- floor(v1[round_down]); v1[!round_down] <- ceiling(v1[!round_down])
      temp <- v1 + 2*(temp-v1)
    }
    PCRed_vec <- temp

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] SymSim_0.0.0.9000           SummarizedExperiment_1.12.0
 [3] DelayedArray_0.8.0          BiocParallel_1.16.0        
 [5] matrixStats_0.54.0          GenomicRanges_1.34.0       
 [7] GenomeInfoDb_1.18.1         IRanges_2.16.0             
 [9] S4Vectors_0.20.1            repr_0.17                  
[11] phytools_0.6-99             maps_3.3.0                 
[13] roxygen2_6.1.1              stringi_1.2.4              
[15] MASS_7.3-51.1               ape_5.2                    
[17] Biobase_2.42.0              BiocGenerics_0.28.0        
[19] RColorBrewer_1.1-2          reshape_0.8.8              
[21] Rtsne_0.15                  ggplot2_3.1.1              
[23] plyr_1.8.4                 

loaded via a namespace (and not attached):
 [1] gtools_3.8.1            assertthat_0.2.0       
 [3] expm_0.999-3            animation_2.5          
 [5] GenomeInfoDbData_1.2.0  yaml_2.2.0             
 [7] numDeriv_2016.8-1       pillar_1.3.1           
 [9] backports_1.1.2         lattice_0.20-38        
[11] glue_1.3.0              quadprog_1.5-5         
[13] phangorn_2.5.5          digest_0.6.18          
[15] promises_1.0.1          XVector_0.22.0         
[17] colorspace_1.3-2        htmltools_0.3.6        
[19] httpuv_1.4.5            Matrix_1.2-15          
[21] pkgconfig_2.0.2         zlibbioc_1.28.0        
[23] purrr_0.3.2             scales_1.0.0           
[25] whisker_0.3-2           later_0.7.5            
[27] git2r_0.26.1            tibble_2.1.1           
[29] combinat_0.0-8          withr_2.1.2            
[31] lazyeval_0.2.1          mnormt_1.5-5           
[33] magrittr_1.5            crayon_1.3.4           
[35] evaluate_0.12           fs_1.3.1               
[37] nlme_3.1-137            xml2_1.2.0             
[39] tools_3.5.1             stringr_1.3.1          
[41] munsell_0.5.0           plotrix_3.7-4          
[43] compiler_3.5.1          clusterGeneration_1.3.4
[45] rlang_0.4.0             RCurl_1.95-4.11        
[47] igraph_1.2.2            bitops_1.0-6           
[49] base64enc_0.1-3         rmarkdown_1.10         
[51] gtable_0.2.0            R6_2.3.0               
[53] knitr_1.20              dplyr_0.8.0.1          
[55] fastmatch_1.1-0         commonmark_1.6         
[57] workflowr_1.6.2         rprojroot_1.3-2        
[59] Rcpp_1.0.4.6            scatterplot3d_0.3-41   
[61] tidyselect_0.2.5        coda_0.19-2