search

MarioniLab / FurtherMNN2018

Code for further development of the mutual nearest neighbours batch correction method, as implemented in the batchelor package.
22 stars 6 forks source link

Error with using fastMNN #14

Open Zifeng-L opened 4 years ago

Zifeng-L commented 4 years ago

Hi, I tried to use fastMNN to correct batch effects, and the codes were as follows, Cause I used seurat for processing the data, so I merge the data using seurat and transform into singlecellexperiment.

### Merge data sets
Merge.B <-  merge(sobj.B$BL1,c(sobj.B$BL2, sobj.B$BL3, sobj.B$BL4, sobj.B$BNHL1, sobj.B$BNHL2, 
sobj.B$rLN1, sobj.B$rLN2, sobj.B$rLN3),
                     add.cell.ids = c("BL1","BL2","BL3", "BL4", "BNHL1", "BNHL2", "rLN1", "rLN2", "rLN3"))     
# Clean meta data and add information about batch
Merge.B@meta.data <- 
  Merge.B@meta.data %>% 
  rownames_to_column(var = "Barcode") %>%
  select(-contains("res.")) %>% 
  dplyr::mutate(Sample=strsplit(Barcode, split = "_") %>% 
           sapply(., "[[", 1),
         Batch = ifelse(Sample %in% c("BL2", "BL4", "BNHL2",  "rLN1", "rLN2"), 2, 1)) %>%
  column_to_rownames(var = "Barcode") 

Merge.B <- NormalizeData(Merge.B, normalization.method = "LogNormalize")
Merge.B <- FindVariableFeatures(Merge.B, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
Merge.B <- ScaleData(Merge.B, verbose = FALSE)
# Set batch to main ident
Merge.B <- SetIdent(Merge.B, cells = WhichCells(Merge.B, expression = Batch %in% c('1')), 1)
Merge.B <- SetIdent(Merge.B, cells = WhichCells(Merge.B, expression = Batch %in% c('2')), 2)
#### Prepare data
Subset_Merge.B <- Merge.B@meta.data %>% rownames_to_column(var="Barcode") %>% 
  group_by(Sample) %>% sample_n(600) %>% pull(Barcode)
Merge.B <- subset(Merge.B, cells = Subset_Merge.B)
# Convert to single cell experiment
sce_B <- as.SingleCellExperiment(Merge.B)
#### Compute SumFactors and Normalize
clusters <- quickCluster(sce_B, min.size=200, get.spikes=FALSE, method="igraph", use.ranks=T)
sce_B <- computeSumFactors(sce_B)
sce_B <- normalize(sce_B)
#### Variance modelling
alt.fit2 <- trendVar(sce_B, use.spikes=FALSE, block=sce_B$ident)
alt.decomp2 <- decomposeVar(sce_B, alt.fit2, block=sce_B$ident)
top.hvgs <- order(alt.decomp2$bio, decreasing=TRUE)
#### Find correlated genes
null.dist <- correlateNull(ncol(sce_B))
cor.pairs <- correlatePairs(sce_B, subset.row=top.hvgs[1:1000], null.dist=null.dist)
null.dist2 <- correlateNull(block=sce_B$ident, iter=1e5)
cor.pairs2 <- correlatePairs(sce_B, subset.row=top.hvgs[1:1000], 
    null.dist=null.dist2, block=sce_B$ident)
cor.genes <- correlatePairs(sce_B, subset.row=top.hvgs[1:1000], 
    null.dist=null.dist, per.gene=TRUE)
#### Separate batches
batch1 <- Merge.B@meta.data %>%
  rownames_to_column(var = "Barcode") %>%
  dplyr::filter(Batch==1) %>% 
  pull(Barcode)
batch2 <- Merge.B@meta.data %>%
  rownames_to_column(var = "Barcode") %>%
  dplyr::filter(Batch==2) %>% 
  pull(Barcode)
sce_B1 <- sce_B[, colnames(sce_B) %in% batch1]
sce_B2 <- sce_B[, colnames(sce_B) %in% batch2]
#### Run MNN
out <- fastMNN(sce_B1, sce_B2)

And the error appeared.

'fastMNN' is deprecated.
Use 'batchelor::fastMNN' instead.
See help("Deprecated")'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available 
Error in env[[as.character(i)]] <- value : 
  wrong args for environment subassignment
sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

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

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

other attached packages:
 [1] cowplot_1.0.0               scran_1.12.1                scater_1.12.2               SingleCellExperiment_1.6.0  SummarizedExperiment_1.14.1
 [6] DelayedArray_0.10.0         BiocParallel_1.18.1         Biobase_2.44.0              GenomicRanges_1.36.1        GenomeInfoDb_1.20.0        
[11] IRanges_2.18.3              S4Vectors_0.22.1            BiocGenerics_0.30.0         ggbeeswarm_0.6.0            scales_1.1.1               
[16] ggpmisc_0.3.5               ggrepel_0.8.2               lemon_0.4.5                 RColorBrewer_1.1-2          viridis_0.5.1              
[21] viridisLite_0.3.0           matrixStats_0.56.0          reshape2_1.4.4              Matrix_1.2-17               readxl_1.3.1               
[26] forcats_0.5.0               stringr_1.4.0               dplyr_1.0.0                 purrr_0.3.4                 readr_1.3.1                
[31] tidyr_1.1.0                 tibble_3.0.1                ggplot2_3.3.2               tidyverse_1.3.0             Seurat_3.2.0               

loaded via a namespace (and not attached):
  [1] backports_1.1.8          plyr_1.8.6               igraph_1.2.4.2           lazyeval_0.2.2           splines_3.6.0           
  [6] listenv_0.8.0            digest_0.6.25            htmltools_0.5.0          fansi_0.4.1              magrittr_1.5            
 [11] tensor_1.5               cluster_2.0.8            ROCR_1.0-11              limma_3.40.6             globals_0.12.5          
 [16] modelr_0.1.8             colorspace_1.4-1         blob_1.2.1               rvest_0.3.6              rappdirs_0.3.1          
 [21] haven_2.3.1              xfun_0.15                crayon_1.3.4             RCurl_1.95-4.12          jsonlite_1.7.0          
 [26] spatstat_1.64-1          spatstat.data_1.4-3      survival_3.2-3           zoo_1.8-8                ape_5.4                 
 [31] glue_1.4.1               polyclip_1.10-0          gtable_0.3.0             zlibbioc_1.30.0          XVector_0.24.0          
 [36] leiden_0.3.3             BiocSingular_1.0.0       future.apply_1.5.0       abind_1.4-5              edgeR_3.26.8            
 [41] DBI_1.1.0                miniUI_0.1.1.1           Rcpp_1.0.4.6             xtable_1.8-4             dqrng_0.2.1             
 [46] reticulate_1.16          rsvd_1.0.3               htmlwidgets_1.5.1        httr_1.4.1               ellipsis_0.3.1          
 [51] ica_1.0-2                pkgconfig_2.0.3          uwot_0.1.8               dbplyr_1.4.4             deldir_0.1-28           
 [56] locfit_1.5-9.4           dynamicTreeCut_1.63-1    tidyselect_1.1.0         rlang_0.4.6              later_1.1.0.1           
 [61] munsell_0.5.0            cellranger_1.1.0         tools_3.6.0              cli_2.0.2                generics_0.0.2          
 [66] broom_0.7.0              ggridges_0.5.2           fastmap_1.0.1            goftest_1.2-2            knitr_1.29              
 [71] fs_1.5.0                 fitdistrplus_1.1-1       RANN_2.6.1               pbapply_1.4-2            future_1.17.0           
 [76] nlme_3.1-139             mime_0.9                 xml2_1.3.2               compiler_3.6.0           rstudioapi_0.11         
 [81] beeswarm_0.2.3           plotly_4.9.2.1           png_0.1-7                spatstat.utils_1.17-0    reprex_0.3.0            
 [86] statmod_1.4.34           stringi_1.4.6            lattice_0.20-38          vctrs_0.3.1              pillar_1.4.4            
 [91] lifecycle_0.2.0          lmtest_0.9-37            BiocNeighbors_1.2.0      RcppAnnoy_0.0.16         data.table_1.12.8       
 [96] bitops_1.0-6             irlba_2.3.3              httpuv_1.5.4             patchwork_1.0.1          R6_2.4.1                
[101] promises_1.1.1           KernSmooth_2.23-15       gridExtra_2.3            vipor_0.4.5              codetools_0.2-16        
[106] MASS_7.3-51.6            assertthat_0.2.1         withr_2.2.0              sctransform_0.2.1        GenomeInfoDbData_1.2.1  
[111] mgcv_1.8-28              hms_0.5.3                rpart_4.1-15             DelayedMatrixStats_1.6.1 Rtsne_0.15              
[116] shiny_1.5.0              lubridate_1.7.9

Can you help me? thanks a lot!

LTLA commented 4 years ago

Several points.

  1. This is the wrong repository. The correct repository is https://github.com/LTLA/batchelor/.
  2. You are using a version of Bioconductor that is at least 1 year out of date. Even if this was a bug in the code, I would only be able to apply bugfixes to the latest version of Bioconductor; you should consider updating.
  3. I am reasonably confident in saying that env[[as.character(i)]] <- value does not occur anywhere in my code. If I had to guess, I would say that this is a problem with the BiocParallel package based on the warnings.

My best suggestion is to update to R 4.0 and Bioconductor 3.11 and try again.