一般要去计算RNA velocity的时候，是已经有先期处理数据了，比如做过了降维，聚类，差异分析，等。因此，做RNA velocity的时候，考虑的经常是怎么把之前的结果和RNA velocity的结果合并展示。而不是，对同一份数据，之前一套，RNA velocity又重新做一套。
因此，在这里会从数据的整合开始。
讨论的情况是：在R上用Seurat做过了分析，现在要用python的scVelo做RNA velocity分析。
这里选择的方法是把velocyto生成的loom文件读取之后，和Seurat分析过的数据整合在一起，然后再导出为loom格式，最后在jupyter notebook上用scVelo做RNA velocity分析。
需要安装的R包：velocyto.R，loomR，hdf5r
用velocyto.R::read.loom.matrices(file = loomFile, engine = “hdf5r)读取loom文件
以Seurat对象为基础，为loom中的unspliced, spliced, ambiguous三个矩阵建立新的Assay
导出loom
读取loom文件
case 1：单个样本
x <- velocyto.R::read.loom.matrices(file = loomFile, engine = “hdf5r)
object[["spliced"]] <- CreateAssayObject(counts = x[["spliced"]])#object为Seurat object
object[["unspliced"]] <- CreateAssayObject(counts = x[["unspliced"]])
object[["ambiguous"]] <- CreateAssayObject(counts = x[["ambiguous"]])
注意：
1. 如果loom中包括的细胞ID和object中的细胞ID不一致，需要对细胞ID进行调整，选择两者公有的。同理，基因ID。
2. 还有一种情况是细胞ID其实一致，但是加了不同的后缀，同样需要去做调整。
解决方式（以spliced矩阵为例）：
spliced <- x[["spliced"]]
rownames(spliced)，查看spliced矩阵的基因ID
colnames(spliced)，查看spliced矩阵的细胞ID
rownames(object)，查看Seurat对象object的基因ID
colnames(object)，查看Seurat对象object的细胞ID
如果发现spliced矩阵和object的基因或者细胞ID不完全一致，那么就取交集。如果发现一致，但是格式有问题，那么就更改一致。
上述操作完成之后，Seurat对象object中包括了RNA counts ，spliced mRNA counts，unspliced mRNA counts，ambiguous mRNA counts四个矩阵。同时，由于object已经是前期处理过的数据，也会有降维信息等。
如：
An object of class Seurat ***** features across **** samples within 4 assays Active assay: RNA (26334 features, 2000 variable features) 3 other assays present: spliced, unspliced, ambiguous 2 dimensional reductions calculated: pca, umap
case 2：多个样本合并
更普遍的一种情况是，我们有多个样本要一起处理。
用Seurat已经将样本merge合并了，integrate整合了。然后对各个样本分别获得了对应的loom文件。现在需要将不同样本的loom也合并起来然后整理到一个Seurat对象中。
x1 <- velocyto.R::read.loom.matrices(file = loomFile1, engine = “hdf5r)
x2 <- velocyto.R::read.loom.matrices(file = loomFile2, engine = “hdf5r)
spliced <- cbind(x1[["spliced"]], x2[["spliced"]])
object[["spliced"]] <- CreateAssayObject(counts = spliced)
unspliced, ambiguous同理。
注意：
1. Seurat合并样本时，为了让细胞ID不冲突，会给不同样本的细胞ID加后缀或者前缀，因此需要根据自己的数据对细胞ID进行处理。
2. 对于Seurat整合数据，即做了批次校正的数据，Seurat对象object中包含了一个基因只有anchor.feature的Assay，因此，最好是采用merge形式的object，然后把integrate形式的降维数据传递过去。
输出loom文件
因为读取了loom之后，新的object包含了很多Assay，需要在输出loom时把所有的Assay中的矩阵都输出。
代码如下：
x <- object
meta.data <- x[[]]
colnames(x = meta.data) <- gsub(
pattern = '\\.',
replacement = '_',
x = colnames(x = meta.data)
)
meta.data$clusters <- as.integer(x = Idents(object = x)[rownames(x = meta.data)])
meta.data$clusters <- as.character(x = Idents(object = x)[rownames(x = meta.data)])
meta.feature <- x[[assay]][[]]
colnames(x = meta.feature) <- gsub(
pattern = '\\.',
replacement = '_',
x = colnames(x = meta.feature)
)
if (length(x = VariableFeatures(object = x)) > 0) {
meta.feature[VariableFeatures(object = x), 'highly_variable_genes'] <- 1
meta.feature[is.na(x = meta.feature$highly_variable_genes), 'highly_variable_genes'] <- 0
}
DefaultAssay(object = x) <- "RNA"
data <- GetAssayData(object = x, slot = 'counts') # Raw counts matrix
layers = list('spliced' = as.matrix(t(as.data.frame(x@assays$spliced@counts))))
# Make the initial loom object
lfile <- loomR::create(
filename = filename,
data = data,
feature.attrs = as.list(x = meta.feature), # Feature-level metadata
cell.attrs = as.list(x = meta.data), # Cell-level metadata
layers = layers, #spliced data
transpose = TRUE,
calc.count = FALSE
)
lfile$add.layer(layers = list('unspliced' = as.matrix(t(as.data.frame(x@assays$unspliced@counts)))))
lfile$add.layer(layers = list('ambiguous' = as.matrix(t(as.data.frame(x@assays$ambiguous@counts)))))
# Add dimensional reduction information
dim.reducs <- names(x@reductions)
for (dr in dim.reducs) {
message("Adding dimensional reduction information for ", dr)
embeddings <- Embeddings(object = x, reduction = dr)[cell.order, ]
embeddings <- list(embeddings)
names(x = embeddings) <- paste0('X_', dr)
message("Adding cell embedding information for ", dr)
lfile$add.col.attribute(attributes = embeddings)
loadings <- Loadings(
object = x,
reduction = dr,
projected = dim(x = slot(object = x[[dr]], name = 'feature.loadings.projected'))
)
}
# Store assay
hdf5r::h5attr(x = lfile, which = 'assay') <- "RNA"
lfile$close_all()
RNA velocity
参照基础用法
