Add D17 downsample aggrggrigate code

r.scripts/sc-TissueMapper_RUN.DS_D17.aggr.R

+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% shared.cells,])
+  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 <- rbind(nmi,data.frame(Sample=i,value=NMI(cluster.All.filter,get(paste0("cluster.",i)))))
+}
+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+geom_text(aes(40000,0.5,label=paste0("Mean Read/Cell = ",round(fit.x,0)),vjust=3))
+#plot.comb <- plot.comb+geom_text(aes(40000,0.5,label=paste0("NMI = ",fit.y),vjust=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()
\ No newline at end of file