Commit 8713aaa9 authored by Gervaise H. Henry's avatar Gervaise H. Henry 🤠

Merge branch 'develop' into 'master'

Merge ds_D17 from develop into master

See merge request sc-TissueMapper_Pr!4
parents 25690b09 68fdeb80
#!/bin/bash
#SBATCH --job-name R_FullAnalysis
#SBATCH -p 256GB,256GBv1
#SBATCH -N 1
#SBATCH -t 7-0:0:0
#SBATCH -o job_%j.out
#SBATCH -e job_%j.out
#SBATCH --mail-type ALL
#SBATCH --mail-user gervaise.henry@utsouthwestern.edu
module load R/3.4.1-gccmkl
Rscript ../r.scripts/sc-TissueMapper_RUN.DS_D17.R
......@@ -103,14 +103,17 @@ scCellCycle <- function(sc10x,sub=FALSE){
sc10x <- ScaleData(object=sc10x,display.progress=FALSE,do.par=TRUE,num.cores=45)
sc10x <- CellCycleScoring(object=sc10x,s.genes=genes.s,g2m.genes=genes.g2m,set.ident=TRUE)
postscript(paste0(folder,"Ridge_cc.Raw.eps"))
plot <- RidgePlot(object=sc10x,features.plot=c("PCNA","TOP2A","MCM6","MKI67"),y.log=TRUE,nCol=2,do.return=TRUE)
plot(plot)
dev.off()
postscript(paste0(folder,"Violin_cc.Raw.eps"))
plot <- VlnPlot(object=sc10x,features.plot=c("PCNA","TOP2A","MCM6","MKI67"),nCol=2,point.size.use=1,size.title.use=20,x.lab.rot=TRUE)
plot(plot)
dev.off()
tryCatch({
postscript(paste0(folder,"Ridge_cc.Raw.eps"))
plot <- RidgePlot(object=sc10x,features.plot=c("PCNA","TOP2A","MCM6","MKI67"),y.log=TRUE,nCol=2,do.return=TRUE)
plot(plot)
dev.off()
postscript(paste0(folder,"Violin_cc.Raw.eps"))
plot <- VlnPlot(object=sc10x,features.plot=c("PCNA","TOP2A","MCM6","MKI67"),nCol=2,point.size.use=1,size.title.use=20,x.lab.rot=TRUE)
plot(plot)
dev.off()
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
sc10x <- RunPCA(object=sc10x,pc.genes=c(genes.s,genes.g2m),do.print=FALSE,pcs.store=2)
postscript(paste0(folder,"PCA_cc.Raw.eps"))
plot <- PCAPlot(object=sc10x,do.return=TRUE)
......@@ -401,7 +404,7 @@ scStress <- function(sc10x,stg="go",res.use=1,pc.use=10,cut=0.95){
dev.off()
#Subsample all cells (+Stress) to better visualize their clustering
if (ncol(sc10x@data)<2500){
if (ncol(sc10x@data)>2500){
rnd <- sample(1:ncol(sc10x@data),2500)
} else {
rnd <- 1:ncol(sc10x@data)
......@@ -898,7 +901,7 @@ scNE <- function(sc10x,neg="EurUro",cut=0.95){
dev.off()
#Subsample all cells (+NE) to better visualize their clustering
if (ncol(sc10x@data)<2500){
if (ncol(sc10x@data)>2500){
rnd <- sample(1:ncol(sc10x@data),2500)
} else {
rnd <- 1:ncol(sc10x@data)
......
This diff is collapsed.
gc()
library(methods)
library(optparse)
library(Seurat)
library(readr)
library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(monocle)
library(dplyr)
library(viridis)
library(reshape2)
library(NMI)
source("../r.scripts/sc-TissueMapper.R")
setwd("../")
load("./analysis/sc10x.Rda")
sc10x.All <-sc10x
rm(sc10x)
downsample <- c("All","350","300","250","200","150","100","075","050","037","025","012","007","005","002")
for (i in downsample[-1]){
load(paste0("../../",i,"/sc-TissueMapper_Pr/analysis/sc10x.Rda"))
assign(paste0("sc10x.",i),sc10x)
rm(sc10x)
}
all.cells <- NULL
for (i in downsample){
all.cells <- c(all.cells,get(paste0("sc10x.",i))@data@Dimnames[[2]])
}
all.cells <- unique(all.cells)
shared.cells <- all.cells
shared.cells.no002<- all.cells
for (i in downsample){
shared.cells <- intersect(shared.cells,get(paste0("sc10x.",i))@data@Dimnames[[2]])
if (i != "002"){shared.cells <- intersect(shared.cells.no002,get(paste0("sc10x.",i))@data@Dimnames[[2]])}
}
for (i in downsample){
assign(paste0("sc10x.",i),SetAllIdent(get(paste0("sc10x.",i)),id="Merge_Epi.dws.sc_St.dws.sc"))
assign(paste0("cluster.",i),data.frame(Barcodes=names(get(paste0("sc10x.",i))@ident),Cluster=get(paste0("sc10x.",i))@ident))
assign(paste0("cluster.",i,".filter"),get(paste0("cluster.",i))[get(paste0("cluster.",i))$Barcodes %in% sc10x.All@data@Dimnames[[2]],])
}
nmi <- data.frame(Sample=character(),value=double())
for (i in downsample[-1]){
nmi <- rbind(nmi,data.frame(Sample=i,value=NMI(cluster.All.filter,get(paste0("cluster.",i,".filter")))))
}
nmi$Sample <- as.numeric(levels(nmi$Sample))
png(paste0("./analysis/NMI.png"),width=1000,height=500,type="cairo")
plot.nmi <- ggplot(nmi,aes(x=Sample,y=value))+geom_point()+geom_smooth(method='loess',formula=y~log(x))+labs(x="Sample (Million Reads)",y="NMI")
model.nmi <- loess(value~log(Sample),data=nmi)
fit.nmi.y <- 0.9
fit.nmi.x <- approx(x=predict(model.nmi),y=nmi$Sample,xout=fit.nmi.y)$y
plot.nmi <- plot.nmi+geom_vline(xintercept=fit.nmi.x)+geom_hline(yintercept=fit.nmi.y)
plot(plot.nmi)
dev.off()
for (i in downsample[-1]){
assign(paste0("rpc.",i),read_csv(paste0("../../../../count/",i,"M_D17PrTzF_Via/outs/metrics_summary.csv"))[,2])
}
rpc <- data.frame(Sample=character(),value=double())
for (i in downsample[-1]){
rpc <- rbind(rpc,data.frame(Sample=i,value=get(paste0("rpc.",i))))
}
colnames(rpc)[2] <- "value"
rpc$Sample <- as.numeric(levels(rpc$Sample))
png(paste0("./analysis/RPC.png"),width=1000,height=500,type="cairo")
plot.rpc <- ggplot(rpc,aes(x=Sample,y=value))+geom_point()+geom_smooth(method='lm',formula=y~x)+labs(x="Sample (Million Reads)",y="Mean Reads Per Cell")
model.rpc <- lm(value~Sample,data=rpc)
fit.rpc.y <- approx(x=rpc$Sample,y=predict(model.rpc),xout=fit.nmi.x)$y
plot.rpc <- plot.rpc+geom_vline(xintercept=fit.nmi.x)+geom_hline(yintercept=fit.rpc.y)
plot(plot.rpc)
dev.off()
comb <- cbind(nmi,rpc[,2])
colnames(comb) <- c("Sample","NMI","RPC")
nmi.rpc <- merge(nmi,rpc,by="Sample")
nmi.rpc <- nmi.rpc[,-1]
colnames(nmi.rpc) <- c("NMI","RPC")
nmi.rpc$NMI <- round(as.numeric(nmi.rpc$NMI),2)
postscript("./analysis/RPC+NMI.eps")
plot.comb <- ggplot(nmi.rpc,aes(x=RPC,y=NMI))+geom_point(colour="blue",size=4)
plot.comb <- plot.comb+geom_smooth(method='loess',formula=y~log(x),size=2)
model <- loess(NMI~RPC,data=nmi.rpc)
fit.y <- 0.9
fit.x <- approx(y=nmi.rpc$RPC,x=predict(model),xout=fit.y)$y
plot.comb <- plot.comb+geom_vline(xintercept=fit.x,linetype=2,size=1.5)+geom_hline(yintercept=fit.y,linetype=2,size=1.5)
plot.comb <- plot.comb+labs(x="Mean Reads Per Cell",y="NMI")
plot.comb <- plot.comb+scale_x_continuous(expand=c(0,0),limits=c(0,80000),breaks=c(seq(0,100000,25000),round(fit.x,0)))+scale_y_continuous(expand=c(0,0),limits=c(0,1),breaks=c(seq(0,1,0.2),fit.y))
plot(plot.comb)
dev.off()
save.image(file="./analysis/NMI.RData")
......@@ -312,7 +312,7 @@ names(gene.set1) <- "Leu"
gene.set <- c(gene.set,gene.set1)
rm(gene.set1)
gc()
min.st <- min(table(sc10x.Epi@meta.data[,paste0("res",opt$res.poststress)]))
min.st <- min(table(sc10x.St@meta.data[,paste0("res",opt$res.poststress)]))
results <- scQuSAGE(sc10x.St,gs=gene.set,res.use=0.2,ds=min.st,nm="St.go",folder="st")
sc10x.St <- results[[1]]
results.cor.St.go <- results[[2]]
......
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