velocyto預(yù)處理

###安裝velocyto###
#安裝依賴(lài)包
conda install numpy scipy cython numba matplotlib scikit-learn h5py click
#安裝velocyto包
pip install velocyto

Run velocyto

nohup velocyto run10x /.../P60-loupe/scVelo/loupe/ /.../refdata-gex-GRCh38-2020-A/genes/genes.gtf &

從Seurat對(duì)象導(dǎo)出細(xì)胞、基因、表達(dá)信息，供 scVelo 使用

#安裝R包
library(Seurat)
library(SpatialCPie)
library(Spaniel)
library(SingleR)
library(infercnv)
library(clustree)
library(clusterProfiler)
library(org.Hs.eg.db)
library(fgsea)
library(tidyverse)
library(ggplot2)
library(ggpubr)
library(devtools)
library(Matrix)
library(cowplot)
library(SeuratData)
library(patchwork)
library(dplyr)
library(hdf5r)

#讀空轉(zhuǎn)數(shù)據(jù)
 sRCC_P1<- Load10X_Spatial(
     data.dir="/.../P1-loupe/outs/",
     filename = "filtered_feature_bc_matrix.h5",
     assay = "Spatial",
     slice = "slice1",
     filter.matrix = TRUE,
     to.upper = FALSE,
)
Idents(sRCC_P1)<-"orig.ident"
sRCC_P1 <- RenameIdents(sRCC_P1, 'SeuratProject' = "sRCC_P1")
sRCC_P1$orig.ident<-Idents(sRCC_P1)

#標(biāo)準(zhǔn)化
sRCC_P1 <- SCTransform(sRCC_P1, assay = "Spatial", verbose = FALSE)

#降維和聚類(lèi)
# umap
sRCC_P1 <- RunPCA(sRCC_P1, assay = "SCT", verbose = FALSE, dims = 1:30)
sRCC_P1 <- FindNeighbors(sRCC_P1, reduction = "pca", dims = 1:30)
sRCC_P1 <- FindClusters(sRCC_P1, verbose = FALSE, resolution = 0.6)
sRCC_P1 <- RunUMAP(sRCC_P1, reduction = "pca", dims = 1:30)

#構(gòu)造函數(shù)Seurat to scVelo
Seurat2scVelo=function(sce, outputRoot){
  message("output to:", outputRoot)
  # save metadata table:
  sce$barcode <- colnames(sce)
  sce$UMAP_1 <- sce@reductions$umap@cell.embeddings[,1]
  sce$UMAP_2 <- sce@reductions$umap@cell.embeddings[,2]
  #write.csv(data.frame(sce@meta.data)[,c("barcode","seurat_clusters", "UMAP_1","UMAP_2")], 
  write.csv(sce@meta.data, 
            file=paste0(outputRoot, 'metadata.csv'), quote=F, row.names=F)
  
  # write expression counts matrix
  library(Matrix)
  counts_matrix <- GetAssayData(sce, assay='Spatial', slot='counts')
  writeMM(counts_matrix, file=paste0(outputRoot, 'counts.mtx'))
  
  # write dimesnionality reduction matrix, in this example case pca matrix
  write.csv(sce@reductions$pca@cell.embeddings, 
            file=paste0(outputRoot, 'pca.csv'), quote=F, row.names=F)
  
  # write gene names
  write.table(
    data.frame('gene'=rownames(counts_matrix)),
    file=paste0(outputRoot, 'gene_names.csv'),
    quote=F,row.names=F,col.names=F
  )
  return(invisible(NULL))
}

# 調(diào)用:
Seurat2scVelo(sRCC_P1, "/.../scVelo/1_sRCC-P1-loupe/velocyto/Seurat2scVelo/")

使用 Python 預(yù)處理數(shù)據(jù)

###讀取 Seurat 導(dǎo)出的數(shù)據(jù)###
import os
os.chdir("/.../scVelo/1_sRCC-P1-loupe/velocyto/Seurat2scVelo/")
os.getcwd()

import scanpy as sc
import anndata 
from scipy import io
from scipy.sparse import coo_matrix, csr_matrix
import numpy as np
import os
import pandas as pd

#(1) load sparse matrix:
X = io.mmread(f"./counts.mtx")
# create anndata object
adata = anndata.AnnData(
    X=X.transpose().tocsr()
)
#(2) load cell metadata:
cell_meta = pd.read_csv(f"./metadata.csv")
# load gene names:
with open(f"./gene_names.csv", 'r') as f:
    gene_names = f.read().splitlines()
# set anndata observations and index obs by barcodes, var by gene names
adata.obs = cell_meta
adata.obs.index = adata.obs['barcode']

adata.var.index = gene_names
#(3) load dimensional reduction:
pca = pd.read_csv(f"./pca.csv")
pca.index = adata.obs.index
# set pca and umap
adata.obsm['X_pca'] = pca.to_numpy()
adata.obsm['X_umap'] = np.vstack((adata.obs['UMAP_1'].to_numpy(), adata.obs['UMAP_2'].to_numpy())).T
# plot a UMAP colored by sampleID to test:
sc.pl.umap(adata, color='seurat_clusters', frameon=False, save=True)

###讀取 loom 數(shù)據(jù)(unsplied/spliced data)###
import scvelo as scv
import scanpy as sc
#import cellrank as cr
import numpy as np
import pandas as pd
import anndata as ad

# load loom files for spliced/unspliced matrices for each sample:
sample_loom_1="/.../scVelo/1_sRCC-P1-loupe/velocyto/sRCC.loom"
ldata1 = scv.read( sample_loom_1, cache=True)
ldata1
#加上cell ID
ldata1.obs.index[0:2]
barcodes = [bc.split(':')[1] for bc in ldata1.obs.index.tolist()]
barcodes = [bc[0:len(bc)-1] + '-1_1' for bc in barcodes]
ldata1.obs.index = barcodes
ldata1.obs.index[0:5]
ldata1.var.head()
ldata1.var_names_make_unique()
### 若有多個(gè)樣本就merge一下
#ldata = ldata1.concatenate([ldata2])
#ldata

#和Seurat數(shù)據(jù)合并
# merge matrices into the original adata object
adata = scv.utils.merge(adata, ldata1)
adata
adata.write_h5ad(f'/.../scVelo/1_sRCC-P1-loupe/velocyto/adata_ldata.h5ad')

scVelo計(jì)算RNA速率

#scVelo對(duì)數(shù)據(jù)預(yù)處理
import scvelo as scv
import scanpy as sc
import numpy as np
import pandas as pd
import anndata as ad
import os
os.chdir("/.../scVelo/1_sRCC-P1-loupe/velocyto/")
os.getcwd()


adata = sc.read(f'/.../scVelo/1_sRCC-P1-loupe/velocyto/adata_ldata.h5ad')
adata
###引入過(guò)Seurat信息后可以不做過(guò)濾等操作
scv.tl.velocity(adata, mode='stochastic')
scv.tl.velocity_graph(adata)
#######直接跳到可視化步驟#######
#scv.pp.filter_and_normalize(adata, min_shared_counts=5, min_shared_cells=3, log=True)
scv.pp.filter_and_normalize(adata)

#可選操作
## clean some genes##
#import re
#flag = [not bool(re.match('^Rp[ls]', i)) for i in adata.var_names]
#adata = adata[:,flag]
#adata = adata[:,adata.var_names != "Malat1"]
#adata

scv.set_figure_params()
scv.pp.moments(adata)

# this step will take a long while
import gc
gc.collect()
#
temp_pre= f"sRCC_nue.in_process2"
if False==os.path.exists(f"{temp_pre}.velo.gz.h5ad"):
    scv.tl.recover_dynamics(adata, var_names='all', n_jobs=64)
    scv.tl.velocity(adata, mode='dynamical')
    adata.write(f"{temp_pre}.velo.gz.h5ad", compression='gzip')
    print(">>Write to file")
else:
    adata = sc.read(f"{temp_pre}.velo.gz.h5ad", compression='gzip', ext="h5ad")
    print(">>read from file")

可視化

embedding

scv.pl.velocity_embedding(adata, basis = 'umap', title="sRCC",
                          save=f"./embedding.pdf",
                          color="seurat_clusters")

grid

scv.settings.set_figure_params('scvelo', dpi=300, dpi_save=300)
fig, ax = plt.subplots()
ax.set_aspect(1)

scv.pl.velocity_embedding_grid(adata, basis='umap',color='seurat_clusters', title='sRCC',
                               arrow_size=1, arrow_length=2, arrow_color="#D2691E",
                               alpha=0.1,
                               #density=0.9,
                               legend_loc='right margin',legend_fontsize=5,
                               show=True,
                               save=f"./grid.pdf", #figsize=(10,10),
                               xlim=[-10,10],ylim=[-10,10], ax=ax)

stream

import matplotlib.pyplot as plt
scv.settings.set_figure_params('scvelo', dpi=300, dpi_save=300)
#fig, ax = plt.subplots()
#ax.set_aspect(1)
scv.pl.velocity_embedding_stream(adata, basis='umap',color='seurat_clusters', title='sRCC',
                               #arrow_size=1, ##arrow_length=2, 
                               #arrow_color="#D2691E",
                               #alpha=0.01, density=0.9,
                               legend_loc='right margin',legend_fontsize=5,
                               show=True,
                               save=f"./stream.pdf")
                                #, #figsize=(10,10),
                               #xlim=[-10,10],ylim=[-10,10], 
                               #  ax=ax)