RNAvelocyto ver.2
Seurat to Velocyto using only R
지난번에 python과 R을 동시에 사용해서 RNA velocity analysis하는 법을 포스팅했었는데, 이번에는 SeuratWrapper libarary 사용해서 python 없이 R만으로 동일한 분석 하는 법을 정리해보려고 한다.
이전 포스팅
Velocyto 사용해서 loom file 제작하는 것까지는 이전 포스팅의 방법과 동일하다. 이후에 Loom file을 R로 불러오고 seurat data의 umap데이터를 합쳐 projection하는 방법에 대해 작성해보고자한다.
(만약에 제 방법에 문제가 있거나 더 좋은 아이디어가 있다면 댓글로 알려주세요!!)
🗂️ Files to prepare before starting analysis
Velocyto돌려서 얻은 loom file- preprocessed(filtering, dimension reduction, clustering, visualization) Seurat data : umap정보 가져와서 그 위에 velocity결과 얹어서 보려고 함.
🔧 Prerequisites to install
- Seurat: install.packages(‘Seurat’)
- velocyto.R: library(devtools); install_github(“velocyto-team/velocyto.R”)
- SeuratWrappers: remotes::install_github(‘satijalab/seurat-wrappers’)
conda로 velocyto.R 설치방법
1
2
3
4
5
git clone https://github.com/velocyto-team/velocyto.R
conda install -c conda-forge boost
conda install -c conda-forge openmp
R CMD build velocyto.R
R CMD INSTALL velocyto.R_0.6.tar.gz
:rotating_light: Error: package or namespace
library(devtools) 했는데 error나면, Error: package or namespace load failed for ‘devtools’ in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]): there is no package called ‘remotes’
1
2
install.packages("glue")
install.packages("devtools")
🚨 ERROR: compilation failed for package ‘hdf5r’
velocyto.R 설치하다 hdf5r과 pcaMethods 때문에 설치가 안된다고 하면 BiocManager::install로 따로 설치해보고, 설치하다 위와 같은 에러가 나면 r-hdf5r을 따로 설치해주고 R package를 다시 깔아보자. 아래는 conda로 해결하는 법.
1
conda install -c conda-forge r-hdf5r
🚨 Error in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]) :
velocyto.R 설치하다
1
2
3
4
5
6
7
8
9
Error in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]) :
there is no package called ‘R.utils’
Calls: <Anonymous> ... loadNamespace -> withRestarts -> withOneRestart -> doWithOneRestart
Execution halted
ERROR: lazy loading failed for package ‘SeuratWrappers’
* removing ‘/home/subin/.conda/envs/Project11/lib/R/library/SeuratWrappers’
Warning message:
In i.p(...) :
installation of package ‘/tmp/RtmpUlGvJj/file3b0d5cab1747/SeuratWrappers_0.3.0.tar.gz’ had non-zero exit status
이런 에러가 난다면, R.utils를 설치한다.
1
install.packages("R.utils")
🚨 Error: Failed to install 'unknown package' from GitHub
에러를 해결하고 seuratwrapper를 설치하려고 했는데, 다음과 같은 에러가 나타난다면..
1
2
3
4
5
6
7
8
9
10
Error: Failed to install 'unknown package' from GitHub:
HTTP error 403.
API rate limit exceeded for 203.255.191.217. (But here's the good news: Authenticated requests get a higher rate limit. Check out the documentation for more details.)
Rate limit remaining: 0/60
Rate limit reset at: 2021-09-14 11:50:10 UTC
To increase your GitHub API rate limit
- Use `usethis::browse_github_pat()` to create a Personal Access Token.
- Use `usethis::edit_r_environ()` and add the token as `GITHUB_PAT`.
reference giuhub api는 시간당 60회의 호출 제한이 걸려 있어서 호출한도에 도달하면 위같은 에러가 나타난다. reference 페이지를 참조하면 api호출 한동 제한을 증가시킬수 있다.
Preprocessing velocity data
만든 loom file을 R로 불러들여 Seurat object로 변환하고 preprocessing과정을 거친다.
1
2
3
library(Seurat)
library(velocyto.R)
library(SeuratWrappers)
1
2
3
4
5
6
7
8
9
10
#Load loom files
sample1_ldat <- ReadVelocity(file = "s1.loom")
sample2_ldat <- ReadVelocity(file = "s2.loom")
..
#Load seurat data
seurat.obj <- readRDS("sample_seurat.rds")
#Get filtered cell ids from seurat data (optional)
sample1_cells <- sapply(colnames(seurat.obj)[grepl("(^(\\w*)-1$)",colnames(seurat.obj))], function(x) paste0("sample1:",substr(x,1,16),"x"))
sample2_cells <- sapply(colnames(seurat.obj)[grepl("(^(\\w*)-2$)",colnames(seurat.obj))], function(x) paste0("sample2:",substr(x,1,16),"x"))
⚠️ 여기서 주의해야할 점은 cell id로 filtering을 진행하기 위해서는 seurat object의 cell id와 loom file로 얻은 데이터의 cell id가 동일해야 한다는 것이다. 내 경우 seurat object의 cell id와 loom file의 cell id의 패턴이 달라서 이를 변환해서 맞는 cell id만 골라 낼 수 있게 한과정을 더 추가해줬다. (굳이 변환해주지 않아도 둘이 동일하다면 skip해도 좋다)
| Seurat cell ids | loom file cell ids |
|---|---|
| AAACCCAGTCCGATCG-1 | sample1:AAACCCAGTCCGATCGx |
| AAACCCAGTGTTGCCG-1 | sample1:AAACCCAGTGTTGCCGx |
| … | … |
| AAACCCATCAGGGTAG-2 | sample2:AAACCCATCAGGGTAGx |
| AAACCCATCCGTAATG-2 | sampel2:AAACCCATCCGTAATGx |
1
2
3
4
5
6
7
8
9
10
11
12
13
#convert loom to seurat and cell filtering
sample1_seurat <- as.Seurat(sample1_ldat)[, sample1_cells]
sample2_seurat <- as.Seurat(sample2_ldat)[, sample2_cells]
#merge the sample files
merged_seurat <- merge(sample1_seurat, sample2_seurat)
#preprocessing
merged_seurat <- SCTransform(object = merged_seurat, assay = "spliced")
merged_seurat <- RunPCA(object = merged_seurat, verbose = FALSE)
merged_seurat <- FindNeighbors(object = merged_seurat, dims = 1:20)
merged_seurat <- FindClusters(object = merged_seurat)
merged_seurat <- RunUMAP(object = merged_seurat, dims = 1:20) #나중에 UMI count로 먼저 계산했던 UMAP으로 덮어씌울 예정
Running RNA velocity
SeuratWrapper로 RNA velocity계산하고 미리 UMI count로 계산했던 UMAP 결과를 velocity seurat object에 추가해 visualization해보자. (이게 맞는 방법인진 확실하지 않음: integrated seurat object와 샘플별 velocity loom file을 단순히 merge해서 계산한 velocity 값을 이런식으로 얹어서 봐도 되는지..)1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# Calculate RNA velocity
## 서버사용하면 ncores옵션 반드시 사용하기 안그러면 속터진다.
merged_seurat <- RunVelocity(object = merged_seurat, deltaT = 1, kCells = 25, fit.quantile = 0.02, ncores = 20)
#add umap coordinate information calculated using UMI counts
numap <- as.matrix(Embeddings(seurat.obj, reduction = "umap"))
rownames(numap) <- rownames(merged_seurat@reductions[["umap"]]@cell.embeddings[,1:2])
merged_seurat@reductions[["numap"]] <- CreateDimReducObject(embeddings = numap, key = "UMAP_", assay = DefaultAssay(merged_seurat))
#cell color
ident.colors <- (scales::hue_pal())(n = length(x = levels(x = merged_seurat)))
names(x = ident.colors) <- levels(x = merged_seurat)
cell.colors <- ident.colors[Idents(object = merged_seurat)]
names(x = cell.colors) <- colnames(x = merged_seurat)
#visualization
#pdf("velocity_analysis.pdf")
show.velocity.on.embedding.cor(emb = Embeddings(object = merged_seurat, reduction = "numap"), vel = Tool(object = merged_seurat,
slot = "RunVelocity"), n = 200, scale = "sqrt", cell.colors = ac(x = cell.colors, alpha = 0.5),
cex = 0.8, arrow.scale = 3, show.grid.flow = TRUE, min.grid.cell.mass = 0.5, grid.n = 40, arrow.lwd = 1,
do.par = FALSE, cell.border.alpha = 0.1, n.cores=20)
#dev.off()