hello，大家好，今天分享一个轨迹分析的软件VECTOR，文章在Unsupervised Inference of Developmental Directions for Single Cells Using VECTOR,2020年8月发表于Cell Report，影响因子8.1分，关于轨迹分析的软件太多了，我们这次来看一看示例代码
Unsupervised Inference of Developmental Directions for Single Cells
Step 1. Please prepare a Seurat object with UMAP and 150 PCs.（这个简单，大家用Seurat分析的时候，应该都是这样的结果）。
简单回顾一下Seurat分析
library(Seurat)
# DATA: Expression matrix. Rownames are gene names. Colnames are cell names.
pbmc <- CreateSeuratObject(counts = DATA, project = "pbmc3k", min.cells = 0, min.features = 0)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 5000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc),npcs = 150)
pbmc <- RunUMAP(pbmc, dims = 1:50)
DimPlot(pbmc, reduction = "umap")
saveRDS(pbmc,file='pbmc.RDS')
Step 2. Get UMAP and PCs from Seurat3 object. (pbmc: a Seurat object):(这个也很简单)：
VEC = pbmc@reductions$umap@cell.embeddings
rownames(VEC) = colnames(pbmc)
PCA = pbmc@reductions$pca@cell.embeddings
source('https://raw.githubusercontent.com/jumphone/Vector/master/Vector.R')
# Remove quantile-based colinearity among PCs (new feature in VECTOR 0.0.3):
PCA=vector.rankPCA(PCA)
Step 3. Use VECTOR:
source('https://raw.githubusercontent.com/jumphone/Vector/master/Vector.R')
# Define pixel
OUT=vector.buildGrid(VEC, N=30,SHOW=TRUE)
# Build network
OUT=vector.buildNet(OUT, CUT=1, SHOW=TRUE)
# Calculate Quantile Polarization (QP) score
OUT=vector.getValue(OUT, PCA, SHOW=TRUE)
# Get pixel's QP score
OUT=vector.gridValue(OUT,SHOW=TRUE)
# Find starting point
OUT=vector.autoCenter(OUT,UP=0.9,SHOW=TRUE)
# Infer vector
OUT=vector.drawArrow(OUT,P=0.9,SHOW=TRUE, COL=OUT$COL, SHOW.SUMMIT=TRUE)
# OUT$P.PS : Peseudotime Score (PS) of each cell
图片.png
这个结果跟RNA Velocyto的结果很像。
Additional function 1: Change QP score to a given gene's expression value (e.g. Nes):
NES.EXP = pbmc@assays$RNA@data[which(rownames(pbmc) =='Nes'),]
OUT=vector.buildGrid(VEC, N=30,SHOW=TRUE)
OUT=vector.buildNet(OUT, CUT=1, SHOW=TRUE)
OUT=vector.getValue(OUT, PCA, SHOW=TRUE)
OUT$VALUE=NES.EXP
OUT=vector.showValue(OUT)
OUT=vector.gridValue(OUT, SHOW=TRUE)
OUT=vector.autoCenter(OUT,UP=0.9,SHOW=TRUE)
OUT=vector.drawArrow(OUT,P=0.9,SHOW=TRUE, COL=OUT$COL)
图片.png
看起来很不错。
Additional function 2: Manually select starting point:（人为指定起始位点）
OUT=vector.buildGrid(VEC, N=30,SHOW=TRUE)
OUT=vector.buildNet(OUT, CUT=1, SHOW=TRUE)
OUT=vector.getValue(OUT, PCA, SHOW=TRUE)
OUT=vector.gridValue(OUT,SHOW=TRUE)
OUT=vector.selectCenter(OUT)
OUT=vector.drawArrow(OUT,P=0.9,SHOW=TRUE, COL=OUT$COL)
图片.png
Additional function 3: Manually select region of interest:(选择感兴趣的区域)
OUT=vector.buildGrid(VEC, N=30,SHOW=TRUE)
OUT=vector.buildNet(OUT, CUT=1, SHOW=TRUE)
OUT=vector.getValue(OUT, PCA, SHOW=TRUE)
OUT=vector.gridValue(OUT,SHOW=TRUE)
OUT=vector.autoCenter(OUT,UP=0.9,SHOW=TRUE)
OUT=vector.drawArrow(OUT,P=0.9,SHOW=TRUE, COL=OUT$COL)
#######################
OUT=vector.reDrawArrow(OUT, COL=OUT$COL)
OUT=vector.selectRegion(OUT)
#######################
SELECT_PS=OUT$SELECT_PS #Peseudotime Score (PS) of selected cells
SELECT_INDEX=OUT$SELECT_INDEX #Index of selected cells in the expression matrix
SELECT_COL=OUT$COL[OUT$SELECT_INDEX] #Colors
#######################
# Identify development related genes
EXP=as.matrix(pbmc@assays$RNA@data)[which(rownames(pbmc) %in% VariableFeatures(pbmc)),SELECT_INDEX]
COR=c()
i=1
while(i<=nrow(EXP)){
this_cor=cor(SELECT_PS, EXP[i,],method='spearman')
COR=c(COR,this_cor)
if(i %%100==1){print(i)}
i=i+1}
names(COR)=rownames(EXP)
head(sort(COR),n=10) #Decreasing (top 10)
tail(sort(COR),n=10) #Increasing (top 10)
# Select one gene to draw figure
show_gene=names(head(sort(COR),n=10))[1]
show_gene.exp=EXP[which(rownames(EXP)==show_gene),]
# Smooth expression value along pesudotime order (optional)
show_gene.exp[order(SELECT_PS)]=smooth.spline(show_gene.exp[order(SELECT_PS)], df=5)$y
# Draw figure
plot(jitter(SELECT_PS), show_gene.exp, pch=16,col=SELECT_COL, ylab=show_gene,xlab='PS')
show_gene.fit=lm(show_gene.exp~SELECT_PS)
abline(show_gene.fit,col='black',lwd=1)
图片.png
Other: Get UMAP and PCs from Monocle3. (cds: a Monocle object):
# Get UMAP:
VEC = cds@reducedDims$UMAP
colnames(VEC) = c('UMAP_1','UMAP_2')
# Get 150 PCs
library(Seurat)
DATA=as.matrix(cds@assays$data[[1]])
pbmc <- CreateSeuratObject(counts = DATA, project = "pbmc3k", min.cells = 0, min.features = 0)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 5000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc),npcs = 150)
PCA = pbmc@reductions$pca@cell.embeddings
看一下文献分析出来的图：
图片.png
图片.png
图片.png
大家感觉怎么样，颜值和RNA Velocyto一摸一样啊，不用做RNA速率分析出来这样的图，赶紧试一下吧
A key step in trajectory inference is the determination of starting cells, which is typically done by using manually selected marker genes. In this study, we find that the quantile polarization（分位数极化） of a cell’s principal-component values is strongly associated with their respective states in development hierarchy（这个地方是关键，细胞的主成分值与发育状态相关）, and therefore provides an unsupervised solution for determining the starting cells. Based on this finding, we developed a tool named VECTOR that infers vectors of developmental directions for cells in Uniform Manifold Approximation and Projection (UMAP). In seven datasets of different developmental scenarios, VECTOR correctly identifies the starting cells and successfully infers the vectors of developmental directions. VECTOR is freely available for academic use at https://github.com/jumphone/Vector.。
用起来倒是很方便
生活很好，有你更好
