hello刽锤，大家好溃论，今天分享一個(gè)軌跡分析的軟件VECTOR傀缩，文章在Unsupervised Inference of Developmental Directions for Single Cells Using VECTOR,2020年8月發(fā)表于Cell Report驻粟，影響因子8.1分役电，關(guān)于軌跡分析的軟件太多了，我們這次來看一看示例代碼

Unsupervised Inference of Developmental Directions for Single Cells

Step 1. Please prepare a Seurat object with UMAP and 150 PCs.（這個(gè)簡單庶柿，大家用Seurat分析的時(shí)候村怪，應(yīng)該都是這樣的結(jié)果）。

簡單回顧一下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):(這個(gè)也很簡單)：

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

這個(gè)結(jié)果跟RNA Velocyto的結(jié)果很像浮庐。

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

看起來很不錯(cuò)甚负。

Additional function 2: Manually select starting point:（人為指定起始位點(diǎn)）

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:(選擇感興趣的區(qū)域)

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

看一下文獻(xiàn)分析出來的圖：

圖片.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（分位數(shù)極化） of a cell’s principal-component values is strongly associated with their respective states in development hierarchy（這個(gè)地方是關(guān)鍵斑举，細(xì)胞的主成分值與發(fā)育狀態(tài)相關(guān)）, 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.。
用起來倒是很方便
生活很好病涨，有你更好