华为云的linux系统ubuntu R中跑velocyto.r

华为云的linux系统ubuntu 跑velocyto.r

实战时候，只需要改动subset_data 和path

getwd()
path="/root/yll/code_test/silicosis/yll/0419"
dir.create(path,recursive = TRUE)
setwd(path)

library(Seurat)
load("macrophage_clean.rds")
as.matrix(table(Idents(subset_data), subset_data$stim))
levels(Idents(subset_data))
table(subset_data$stim)
unique(subset_data$stim)
table(Idents(subset_data))
subset_data$cluster = Idents(subset_data)
Idents(subset_data) = subset_data$stim

##	1
library(velocyto.R)
library(sccore)
dir.create(paste(path, "RNAVelocity", sep = "/"))
setwd(paste(path, "RNAVelocity", sep = "/"))
stims = data.frame(stim = unique(subset_data$stim), 
                   loom = c("/root/yll/SiO2_7d.loom", "/root/yll/SiO2_56d_2.loom",
                            "/root/yll/NS_7d.loom", "/root/yll/NS_56d_2.loom"), s1 = c(9, 12, 7, 10))
stims
for(i in 2:nrow(stims)){
  print(stims[i,"stim"])
  subset_data_stim = subset(subset_data, ident = stims[i,"stim"], invert = FALSE)
  Idents(subset_data_stim) = subset_data_stim$cluster
  ldat <- read.loom.matrices(stims[i,"loom"])##	list: spliced, unspliced, ambiguous
  emat.temp <- ldat$spliced
  colnames(emat.temp) = paste(substr(colnames(emat.temp), stims[i, "s1"], (stims[i, "s1"]+15)), i, sep = ".")
  cells = intersect(colnames(emat.temp), colnames(subset_data_stim))
  emat.temp <- emat.temp[,cells]
  nmat.temp <- ldat$unspliced
  colnames(nmat.temp) = paste(substr(colnames(nmat.temp), stims[i, "s1"], (stims[i, "s1"]+15)), i, sep = ".")
  cells = intersect(colnames(nmat.temp), colnames(subset_data_stim))
  nmat.temp <- nmat.temp[,cells]
  emat = emat.temp
  nmat = nmat.temp
  cluster.label <- Idents(subset_data_stim)
  cell.colors <- sccore::fac2col(cluster.label)
  emat <- filter.genes.by.cluster.expression(emat,cluster.label,min.max.cluster.average = 0.5)
  nmat <- filter.genes.by.cluster.expression(nmat,cluster.label,min.max.cluster.average = 0.05)
  emb <- subset_data_stim$tsne@cell.embeddings
  cell.dist <- as.dist(1-armaCor(t(emb)))
  fit.quantile <- 0.02
  rvel.cd <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells=20,cell.dist=cell.dist,fit.quantile=fit.quantile)
  pdf(paste("RNA_Velocity",stims[i, "stim"], "tsne.pdf", sep = "_"))
  p = show.velocity.on.embedding.cor(emb,rvel.cd,n=300,scale='sqrt',cell.colors=ac(cell.colors,alpha=0.5),cex=0.8,arrow.scale=5,show.grid.flow=TRUE,min.grid.cell.mass=0.5,grid.n=40,arrow.lwd=1,do.par=F,cell.border.alpha = 0.1)
  print(p)
  dev.off()
  emb <- subset_data_stim$umap@cell.embeddings
  cell.dist <- as.dist(1-armaCor(t(emb)))
  fit.quantile <- 0.02
  rvel.cd <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells=20,cell.dist=cell.dist,fit.quantile=fit.quantile)
  pdf(paste("RNA_Velocity",stims[i, "stim"], "umap.pdf", sep = "_"))
  p = show.velocity.on.embedding.cor(emb,rvel.cd,n=300,scale='sqrt',cell.colors=ac(cell.colors,alpha=0.5),cex=0.8,arrow.scale=5,show.grid.flow=TRUE,min.grid.cell.mass=0.5,grid.n=40,arrow.lwd=1,do.par=F,cell.border.alpha = 0.1)
  print(p)
  dev.off()
}