Liger batch effect elimination
Ehsan
05/05/2022
library(rliger)## Loading required package: cowplot## Loading required package: Matrix## Loading required package: patchwork##
## Attaching package: 'patchwork'## The following object is masked from 'package:cowplot':
##
## align_plotslibrary(Seurat)## Attaching SeuratObjectlibrary(SeuratData)## ── Installed datasets ───────────────────────────────────── SeuratData v0.2.2 ──## ✓ ifnb 3.1.0 ✓ pbmcsca 3.0.0
## ✓ panc8 3.0.2## ────────────────────────────────────── Key ─────────────────────────────────────## ✓ Dataset loaded successfully
## > Dataset built with a newer version of Seurat than installed
## ❓ Unknown version of Seurat installedlibrary(SeuratWrappers)ifnb <- readRDS("/mnt/nectar_volume/home/eraz0001/new/IFNB.rds")head(ifnb)## orig.ident nCount_RNA nFeature_RNA stim seurat_annotations
## AAACATACATTTCC.1 IMMUNE_CTRL 3017 877 CTRL CD14 Mono
## AAACATACCAGAAA.1 IMMUNE_CTRL 2481 713 CTRL CD14 Mono
## AAACATACCTCGCT.1 IMMUNE_CTRL 3420 850 CTRL CD14 Mono
## AAACATACCTGGTA.1 IMMUNE_CTRL 3156 1109 CTRL pDC
## AAACATACGATGAA.1 IMMUNE_CTRL 1868 634 CTRL CD4 Memory T
## AAACATACGGCATT.1 IMMUNE_CTRL 1581 557 CTRL CD14 Mono
## AAACATACTGCGTA.1 IMMUNE_CTRL 2747 980 CTRL T activated
## AAACATACTGCTGA.1 IMMUNE_CTRL 1341 581 CTRL CD4 Naive T
## AAACATTGAGTGTC.1 IMMUNE_CTRL 2155 880 CTRL CD8 T
## AAACATTGCTTCGC.1 IMMUNE_CTRL 2536 669 CTRL CD14 MonoWe can use the following functions in LIGER pckg or even in seurat. we used the seurat’s ones.
pbmcsca <- NormalizeData(pbmcsca)
pbmcsca <- FindVariableFeatures(pbmcsca)
pbmcsca <- ScaleData(pbmcsca, split.by = "Method", do.center = FALSE)## Scaling data matrix## Scaling data from split 10x Chromium (v2)## Scaling data from split 10x Chromium (v2) A## Scaling data from split 10x Chromium (v2) B## Scaling data from split 10x Chromium (v3)## Scaling data from split CEL-Seq2## Scaling data from split Drop-seq## Scaling data from split inDrops## Scaling data from split Seq-Well## Scaling data from split Smart-seq2head(pbmcsca)## orig.ident nCount_RNA nFeature_RNA nGene nUMI
## pbmc1_SM2_Cell_108 pbmc1 437125 2200 2200 437125
## pbmc1_SM2_Cell_115 pbmc1 335596 2438 2438 335596
## pbmc1_SM2_Cell_133 pbmc1 302204 1874 1874 302204
## pbmc1_SM2_Cell_142 pbmc1 377420 2480 2480 377420
## pbmc1_SM2_Cell_143 pbmc1 385514 2196 2196 385514
## pbmc1_SM2_Cell_144 pbmc1 304994 2216 2216 304994
## pbmc1_SM2_Cell_146 pbmc1 410522 2331 2331 410522
## pbmc1_SM2_Cell_148 pbmc1 335187 2075 2075 335187
## pbmc1_SM2_Cell_152 pbmc1 485530 2830 2830 485530
## pbmc1_SM2_Cell_153 pbmc1 280213 2263 2263 280213
## percent.mito Cluster CellType Experiment
## pbmc1_SM2_Cell_108 0.0297434465355702 0 Cytotoxic T cell pbmc1
## pbmc1_SM2_Cell_115 0.0311521658159055 0 Cytotoxic T cell pbmc1
## pbmc1_SM2_Cell_133 0.0431128105727693 0 Cytotoxic T cell pbmc1
## pbmc1_SM2_Cell_142 0.0260323569927476 0 Cytotoxic T cell pbmc1
## pbmc1_SM2_Cell_143 0.0404759383962183 0 Cytotoxic T cell pbmc1
## pbmc1_SM2_Cell_144 0.023409951391094 0 Cytotoxic T cell pbmc1
## pbmc1_SM2_Cell_146 0.0266065031418581 0 Cytotoxic T cell pbmc1
## pbmc1_SM2_Cell_148 0.0446061115049967 0 Cytotoxic T cell pbmc1
## pbmc1_SM2_Cell_152 0.0234517593805909 5 Cytotoxic T cell pbmc1
## pbmc1_SM2_Cell_153 0.0346036948333838 0 Cytotoxic T cell pbmc1
## Method
## pbmc1_SM2_Cell_108 Smart-seq2
## pbmc1_SM2_Cell_115 Smart-seq2
## pbmc1_SM2_Cell_133 Smart-seq2
## pbmc1_SM2_Cell_142 Smart-seq2
## pbmc1_SM2_Cell_143 Smart-seq2
## pbmc1_SM2_Cell_144 Smart-seq2
## pbmc1_SM2_Cell_146 Smart-seq2
## pbmc1_SM2_Cell_148 Smart-seq2
## pbmc1_SM2_Cell_152 Smart-seq2
## pbmc1_SM2_Cell_153 Smart-seq2tail(pbmcsca)## orig.ident nCount_RNA nFeature_RNA
## pbmc2_inDrops_1_GTGATAAA.TTATGCGA.CTTCTTCG pbmc2 1082 704
## pbmc2_inDrops_1_GTGGTGCT.TTATGCGA.ACCTTGCC pbmc2 832 522
## pbmc2_inDrops_1_GTTACTAG.GAGCCTTA.GCACCTCT pbmc2 1463 807
## pbmc2_inDrops_1_TAAACCGA.TTATGCGA.ATGCGGAT pbmc2 917 561
## pbmc2_inDrops_1_TAGATCAA.TTATGCGA.CTTCGCAC pbmc2 807 556
## pbmc2_inDrops_1_TAGTCTCT.GAGCCTTA.ATCCGCTA pbmc2 717 453
## pbmc2_inDrops_1_TCCAGAAG.TTATGCGA.TAAGACGG pbmc2 938 592
## pbmc2_inDrops_1_TGAATCCT.GAGCCTTA.CCCAAGCA pbmc2 662 406
## pbmc2_inDrops_1_TGAATCCT.TTATGCGA.CATCTCCC pbmc2 2066 1001
## pbmc2_inDrops_1_TGAGCACA.GAGCCTTA.CGAGTCTG pbmc2 253 184
## nGene nUMI percent.mito
## pbmc2_inDrops_1_GTGATAAA.TTATGCGA.CTTCTTCG 704 1082 0.0286506469500924
## pbmc2_inDrops_1_GTGGTGCT.TTATGCGA.ACCTTGCC 522 832 0.0240384615384615
## pbmc2_inDrops_1_GTTACTAG.GAGCCTTA.GCACCTCT 807 1463 0.050580997949419
## pbmc2_inDrops_1_TAAACCGA.TTATGCGA.ATGCGGAT 561 917 0.0588876772082879
## pbmc2_inDrops_1_TAGATCAA.TTATGCGA.CTTCGCAC 556 807 0.0247831474597274
## pbmc2_inDrops_1_TAGTCTCT.GAGCCTTA.ATCCGCTA 453 717 0.097629009762901
## pbmc2_inDrops_1_TCCAGAAG.TTATGCGA.TAAGACGG 592 938 0.035181236673774
## pbmc2_inDrops_1_TGAATCCT.GAGCCTTA.CCCAAGCA 406 662 0.138972809667674
## pbmc2_inDrops_1_TGAATCCT.TTATGCGA.CATCTCCC 1001 2066 0.0556631171345595
## pbmc2_inDrops_1_TGAGCACA.GAGCCTTA.CGAGTCTG 184 253 0.146245059288538
## Cluster CellType
## pbmc2_inDrops_1_GTGATAAA.TTATGCGA.CTTCTTCG 11 Plasmacytoid dendritic cell
## pbmc2_inDrops_1_GTGGTGCT.TTATGCGA.ACCTTGCC 11 Plasmacytoid dendritic cell
## pbmc2_inDrops_1_GTTACTAG.GAGCCTTA.GCACCTCT 11 Plasmacytoid dendritic cell
## pbmc2_inDrops_1_TAAACCGA.TTATGCGA.ATGCGGAT 11 Plasmacytoid dendritic cell
## pbmc2_inDrops_1_TAGATCAA.TTATGCGA.CTTCGCAC 7 Dendritic cell
## pbmc2_inDrops_1_TAGTCTCT.GAGCCTTA.ATCCGCTA 11 Plasmacytoid dendritic cell
## pbmc2_inDrops_1_TCCAGAAG.TTATGCGA.TAAGACGG 11 Plasmacytoid dendritic cell
## pbmc2_inDrops_1_TGAATCCT.GAGCCTTA.CCCAAGCA 11 Plasmacytoid dendritic cell
## pbmc2_inDrops_1_TGAATCCT.TTATGCGA.CATCTCCC 11 Plasmacytoid dendritic cell
## pbmc2_inDrops_1_TGAGCACA.GAGCCTTA.CGAGTCTG 11 Plasmacytoid dendritic cell
## Experiment Method
## pbmc2_inDrops_1_GTGATAAA.TTATGCGA.CTTCTTCG pbmc2 inDrops
## pbmc2_inDrops_1_GTGGTGCT.TTATGCGA.ACCTTGCC pbmc2 inDrops
## pbmc2_inDrops_1_GTTACTAG.GAGCCTTA.GCACCTCT pbmc2 inDrops
## pbmc2_inDrops_1_TAAACCGA.TTATGCGA.ATGCGGAT pbmc2 inDrops
## pbmc2_inDrops_1_TAGATCAA.TTATGCGA.CTTCGCAC pbmc2 inDrops
## pbmc2_inDrops_1_TAGTCTCT.GAGCCTTA.ATCCGCTA pbmc2 inDrops
## pbmc2_inDrops_1_TCCAGAAG.TTATGCGA.TAAGACGG pbmc2 inDrops
## pbmc2_inDrops_1_TGAATCCT.GAGCCTTA.CCCAAGCA pbmc2 inDrops
## pbmc2_inDrops_1_TGAATCCT.TTATGCGA.CATCTCCC pbmc2 inDrops
## pbmc2_inDrops_1_TGAGCACA.GAGCCTTA.CGAGTCTG pbmc2 inDropspbmcsca <- RunOptimizeALS(pbmcsca, k = 20, lambda = 5, split.by = "Method")##
|
| | 0%
|
|== | 3%
|
|===== | 7%
|
|======= | 10%
|
|========= | 13%
|
|============ | 17%
|
|============== | 20%
|
|================ | 23%
|
|=================== | 27%
|
|===================== | 30%
|
|======================= | 33%
|
|========================== | 37%
|
|============================ | 40%
|
|============================== | 43%
|
|================================= | 47%
|
|=================================== | 50%
|
|===================================== | 53%
|
|======================================== | 57%
|
|========================================== | 60%
|
|============================================ | 63%
|
|=============================================== | 67%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|====================================================== | 77%
|
|======================================================== | 80%
|
|========================================================== | 83%
|
|============================================================= | 87%
|
|=============================================================== | 90%
|
|================================================================= | 93%
|
|==================================================================== | 97%
|
|======================================================================| 100%
## Finished in 16.9668 mins, 30 iterations.
## Max iterations set: 30.
## Final objective delta: 1.718198e-05.
## Best results with seed 1.## Warning: No columnames present in cell embeddings, setting to 'riNMF_1:20'pbmcsca <- RunQuantileNorm(pbmcsca, split.by = "Method")pbmcsca <- FindNeighbors(pbmcsca, reduction = "iNMF", dims = 1:20)## Computing nearest neighbor graph## Computing SNNpbmcsca <- FindClusters(pbmcsca, resolution = 0.3)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 31021
## Number of edges: 1037536
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9258
## Number of communities: 12
## Elapsed time: 8 secondspbmcsca <- RunUMAP(pbmcsca, dims = 1:ncol(pbmcsca[["iNMF"]]), reduction = "iNMF")## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 13:28:55 UMAP embedding parameters a = 0.9922 b = 1.112## 13:28:55 Read 31021 rows and found 20 numeric columns## 13:28:55 Using Annoy for neighbor search, n_neighbors = 30## 13:28:55 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 13:29:00 Writing NN index file to temp file /tmp/RtmpVKEUNG/file2da55b1402a26b
## 13:29:00 Searching Annoy index using 1 thread, search_k = 3000
## 13:29:17 Annoy recall = 100%
## 13:29:17 Commencing smooth kNN distance calibration using 1 thread
## 13:29:19 Initializing from normalized Laplacian + noise
## 13:29:21 Commencing optimization for 200 epochs, with 1339728 positive edges
## 13:29:37 Optimization finishedDimPlot(pbmcsca, group.by = c("Method", "ident", "CellType"), ncol = 3)