Commit 47279480 authored by Gervaise H. Henry's avatar Gervaise H. Henry 🤠

Add code for Velocyto.R

parent 4387f699
library(methods)
library(optparse)
library(Seurat)
library(readr)
library(readxl)
library(fBasics)
library(pastecs)
library(RColorBrewer)
library(dplyr)
library(viridis)
library(gridExtra)
library(scater)
library(loomR)
library(velocyto.R)
library(SeuratWrappers)
source("./r.scripts/sc-TissueMapper_functions.R")
source("./r.scripts/sc-TissueMapper_process.R")
load("./analysis/sc10x.Rda")
try(
if (as.numeric(substring(sc10x@version,1,1))<3){
sc10x <- UpdateSeuratObject(sc10x)
}
)
Idents(object=sc10x) <- "Merge_Epi.dws_St.go_NE"
sc10x <- RenameIdents(object=sc10x, "Leu" = "Leu")
sc10x <- RenameIdents(object=sc10x, "Endo" = "Endo")
sc10x <- RenameIdents(object=sc10x, "SM" = "SM")
sc10x <- RenameIdents(object=sc10x, "Fib" = "Fib")
sc10x <- RenameIdents(object=sc10x, "NE" = "NE")
sc10x <- RenameIdents(object=sc10x, "OE1" = "Club")
sc10x <- RenameIdents(object=sc10x, "OE2" = "Hillock")
sc10x <- RenameIdents(object=sc10x, "LE" = "LE")
sc10x <- RenameIdents(object=sc10x, "BE" = "BE")
sc10x$id <- Idents(sc10x)
UMAP.aggr <- read.csv("/work/urology/ghenry/RNA-Seq/SingleCell/PIPELINE/RUN/Pd/PUBLISH/Pd_velocyto/analysis/UMAP.aggr.csv", row.names=1)
sc10x[["umap"]] <- CreateDimReducObject(embeddings=as.matrix(UMAP.aggr),key="UMAP_",assay="RNA")
for (i in unique(sc10x$samples)){
if (!dir.exists(paste0("./analysis/DATA/",i,"/"))){
dir.create(paste0("./analysis/DATA/",i,"/"))
}
temp <- colnames(sc10x)[sc10x$samples==i]
temp <- substr(temp,5,nchar(temp))
temp <- substr(temp,1,nchar(temp)-2)
temp <- paste0(temp,"-1")
write.table(temp,file=paste0("./analysis/DATA/",i,"/","barcodes.tsv"),quote=FALSE,sep="\t",row.names=FALSE,col.names=FALSE)
system(paste0("gzip ./analysis/DATA/",i,"/","barcodes.tsv" -f))
file.copy(paste0("./analysis/DATA/",i,"/","barcodes.tsv.gz"),paste0("./analysis/DATA/10x/",i,"/outs/filtered_feature_bc_matrix/","barcodes.tsv.gz"),overwrite=TRUE)
rm(temp)
}
rm(i)
j=1
for (i in unique(sc10x$samples)){
assign(paste0("ldata.",i),ReadVelocity(file=paste0("analysis/DATA/10x/",i,"/velocyto/",i,".loom")))
ldata.temp <- get(paste0("ldata.",i))
temp <- paste0(substr(colnames(get(paste0("ldata.",i))$spliced),1,3),"_",paste0(substr(colnames(get(paste0("ldata.",i))$spliced),14,nchar(colnames(get(paste0("ldata.",i))$spliced))-1),"-",j))
colnames(ldata.temp$spliced) <- temp
colnames(ldata.temp$unspliced) <- temp
colnames(ldata.temp$ambiguous) <- temp
assign(paste0("ldata.",i),ldata.temp)
rm(ldata.temp)
rm(temp)
assign(paste0("bm.",i),as.Seurat(x=get(paste0("ldata.",i))))
assign(paste0("bm.",i),subset(get(paste0("bm.",i)),cells=colnames(sc10x)))
j=j+1
}
rm(i)
rm(j)
bm <- merge(x=get(paste0("bm.",unique(sc10x$samples)[1])),y=mget(ls(pat="bm."))[-1])
sc10x@assays$spliced <- bm@assays$spliced
sc10x@assays$unspliced <- bm@assays$unspliced
sc10x <- SCTransform(sc10x,assay="spliced")
sc10x <- RunVelocity(object=sc10x,deltaT=1,kCells=25,fit.quantile=0.02,reduction="cca.aligned")
Idents(sc10x) <- "id"
ident.colors <- c(rev(brewer.pal(5,"Blues")[2:5]),"grey50",rev(brewer.pal(5,"Reds")[2:5]))
names(x = ident.colors) <- c("BE","LE","Hillock","Club","NE","Fib","SM","Endo","Leu")
cell.colors <- ident.colors[sc10x$id]
names(x = cell.colors) <- colnames(x = sc10x)
gc()
postscript(paste0("./analysis/velocity.eps"))
show.velocity.on.embedding.cor(emb=Embeddings(object=sc10x,reduction="umap"),vel=Tool(object=sc10x,slot="RunVelocity"),n=200,scale="sqrt",cell.colors=ac(x=cell.colors),cex=0.8,arrow.scale=3,show.grid.flow=TRUE,min.grid.cell.mass=0.5,grid.n=40,arrow.lwd=1,do.par=FALSE,n.cores=50)
dev.off()
gc()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment