轨迹分析系列：
- 单细胞之轨迹分析-1：RNA velocity
- 单细胞之轨迹分析-2：monocle2 原理解读+实操
- 单细胞之轨迹分析-3：monocle3
- 单细胞之轨迹分析-4：scVelo
- 单细胞之轨迹分析-5：slingshot
- 单细胞之轨迹分析-6：velocyto.R+Seurat
一般要去计算RNA velocity的时候，是已经预先处理过数据了，比如做过了降维，聚类，差异分析等。因此，做RNA velocity的时候，考虑的经常是怎么把之前的结果和RNA velocity的结果合并展示。而不是对同一份数据使用RNA velocity重新做一次降维聚类。
思路：把velocyto生成的loom文件读取之后，和Seurat分析过的数据整合在一起，然后再导出为loom格式，最后用scVelo做velocity分析。
1. Introduction
需要用到的软件：
- scVelo (For RNA Velocity)
- Velocyto or Kallisto Bustools (To produce our initial RNA Velocity Object)
- Anndata (For manipulation of our RNA Velocity object)
- Seurat
- Samtools -- optional (Velocyto will run Samtools sort on unsorted .bam)
2. 生成loom文件
loom文件是从fastq/loom文件中得到的
pip install git+https://github.com/pachterlab/kb_python@devel
kb ref -i index.idx -g t2g.txt -f1 cdna.fa -f2 intron.fa -c1 cdna_t2c.txt -c2 intron_t2c.txt --workflow lamanno -n 4 \
fasta.fa \
gtf.gtf
kb count -i transcriptome.idx -g t2g.txt -x 10xv2 --workflow lamanno --loom -c1 cdna_t2c.txt -c2 intron_t2c.txt read_1.fastq.gz read_2.fastq.gz
#Download dependencies first
conda install numpy scipy cython numba matplotlib scikit-learn h5py click
pip install velocyto
velocyto run -b filtered_barcodes.tsv -o output_path -m repeat_msk_srt.gtf bam_file.bam annotation.gtf
3. 读取Seurat对象和loom文件
需要先转换成h5ad格式，参考Seurat对象、SingleCellExperiment对象和scanpy对象的转化
#数据转换
library(scater)
library(Seurat)
library(SeuratData)
#remotes::install_github("mojaveazure/seurat-disk")
library(SeuratDisk)
library(patchwork)
pbmc <- readRDS("pbmc.rds")
SaveH5Seurat(pbmc, filename = "pbmc.h5Seurat")
Convert("pbmc.h5Seurat", dest = "h5ad")
读取数据Seurat整合对象
import anndata
import scvelo as scv
import pandas as pd
import numpy as np
import matplotlib as plt
import scanpy as sc
%load_ext rpy2.ipython
adata=sc.read_h5ad('pbmc.h5ad')
adata.obs.seurat_clusters=adata.obs.seurat_clusters.astype('category')
读取每个样品的loom文件
data1 = anndata.read_loom("data1.loom")
data2 = anndata.read_loom("data2.loom")
data3 = anndata.read_loom("data3.loom")
4. 根据Seurat对象的细胞ID，修改loom文件细胞ID
barcodes=[bc.split(':')[1] for bc in data1.obs.index.tolist()]
barcodes=[bc[0:len(bc)-1]+ '-1_1' for bc in barcodes]
data1.obs.index=barcodes
data1.var_names_make_unique()
data2和data3的操作相同
5. 整合loom文件
ldata=data1.concatenate([data2,data3])
6. 整合loom文件和metadata
adata=scv.utils.merge(adata,ldata)
画个umap图检查一下
sc.pl.umap(adata, color='celltype', frameon=False, legend_loc='on data', title='', save='_celltypes.pdf')
为不同的细胞类型、样本、细胞类群等设置颜色（可选）
（对应的obs名，然后跟“_colors”）
adata.uns['Group_colors'] = np.array(["#66c2a5", "#8da0cb", "#e78ac3"])
adata.uns['celltype_colors'] = np.array([""#33a02c", "#b2df8a", "#a6cee3", "#fb9a99", "#cab2d6"])
7. scVelo分析
参考scVelo
8. 提取亚群分析
cur_celltypes = ['CD4T', 'CD8T, 'Treg', 'Tnaive']
adata_subset = adata[adata.obs['celltype'].isin(cur_celltypes)]
sc.pl.umap(adata_subset, color=['celltype', 'condition'], frameon=False, title=['', ''])
sc.pp.neighbors(adata_subset, n_neighbors=15, use_rep='X_pca')
# pre-process
scv.pp.filter_and_normalize(adata_subset)
scv.pp.moments(adata_subset)
后续分析同scVelo
参考：
scvelo github网站：https://github.com/theislab/scvelo
scvelo官方文档：https://scvelo.readthedocs.io/index.html
Seurat to RNA-Velocity教程：https://github.com/basilkhuder/Seurat-to-RNA-Velocity#multiple-sample-integration
scvelo实战教程：
https://smorabit.github.io/tutorials/8_velocyto/
RNA velocity：scVelo 应用
