Commit da218911 authored by Gervaise Henry's avatar Gervaise Henry 🤠
Browse files

Finish shiny output for raw and id including deg

parent a68d0d81
......@@ -8,6 +8,7 @@ library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(monocle)
library(dplyr)
library(viridis)
library(gridExtra)
......@@ -208,6 +209,12 @@ for (i in c("epi","fmst","st","leu")){
rm(list=ls(pattern="^sc10x.se"))
scShinyOutput(sc10x,anal="id")
scShinyOutput(sc10x.epi,anal="id.epi")
scShinyOutput(sc10x.fmst,anal="id.fmst")
scShinyOutput(sc10x.st,anal="id.st")
scShinyOutput(sc10x.leu,anal="id.leu")
save(list=ls(pattern="^singler"),file='./analysis/singler_objects.RData')
save(list=ls(pattern="^sc10x"),file='./analysis/sc10x.id.rda')
#save.image(file="./analysis/sc10x.id.RData")
......@@ -15,6 +15,7 @@ library(gridExtra)
library(SingleR)
library(sctransform)
library(autothresholdr)
library(ggplot2)
options(bitmapType="cairo")
options(future.globals.maxSize= 1000000000000)
......@@ -140,7 +141,7 @@ sc10x <- scCluster(sc10x,res=res,red="pca",dim=pc.use.poststress,print="2",folde
DefaultAssay(object=sc10x) <- "SCT"
scShinyOutput()
scShinyOutput(sc10x,anal="raw")
save(sc10x,file=paste0("./analysis/sc10x.raw.rda"))
save.image(file="./analysis/sc10x.raw.RData")
......@@ -32,38 +32,43 @@ scFolders <- function(){
if (!dir.exists("./analysis/shiny")){
dir.create("./analysis/shiny")
}
if (!dir.exists("./analysis/shiny/raw")){
dir.create("./analysis/shiny/raw")
}
if (!dir.exists("./analysis/shiny/raw/outs")){
dir.create("./analysis/shiny/raw/outs")
}
if (!dir.exists("./analysis/shiny/raw/outs/analysis")){
dir.create("./analysis/shiny/raw/outs/analysis")
}
if (!dir.exists("./analysis/shiny/raw/outs/analysis/clustering")){
dir.create("./analysis/shiny/raw/outs/analysis/clustering")
}
if (!dir.exists("./analysis/shiny/raw/outs/analysis/diffexp")){
dir.create("./analysis/shiny/raw/outs/analysis/diffexp")
}
if (!dir.exists("./analysis/shiny/raw/outs/analysis/pca")){
dir.create("./analysis/shiny/raw/outs/analysis/pca")
}
if (!dir.exists("./analysis/shiny/raw/outs/analysis/pca/10_components")){
dir.create("./analysis/shiny/raw/outs/analysis/pca/10_components")
}
if (!dir.exists("./analysis/shiny/raw/outs/analysis/tsne")){
dir.create("./analysis/shiny/raw/outs/analysis/tsne")
}
if (!dir.exists("./analysis/shiny/raw/outs/analysis/tsne/2_components")){
dir.create("./analysis/shiny/raw/outs/analysis/tsne/2_components")
}
if (!dir.exists("./analysis/shiny/raw/outs/analysis/umap")){
dir.create("./analysis/shiny/raw/outs/analysis/umap")
}
if (!dir.exists("./analysis/shiny/raw/outs/analysis/umap/2_components")){
dir.create("./analysis/shiny/raw/outs/analysis/umap/2_components")
for (i in c("raw","id","id.epi","id.fmst","id.st","id.leu")){
if (!dir.exists(paste0("./analysis/shiny/",i))){
dir.create(paste0("./analysis/shiny/",i))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs"))){
dir.create(paste0("./analysis/shiny/",i,"/outs"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/filtered_feature_bc_matrix"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/filtered_feature_bc_matrix"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/analysis"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/analysis"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/analysis/clustering"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/analysis/clustering"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/analysis/diffexp"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/analysis/diffexp"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/analysis/pca"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/analysis/pca"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/analysis/pca/10_components"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/analysis/pca/10_components"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/analysis/tsne"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/analysis/tsne"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/analysis/tsne/2_components"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/analysis/tsne/2_components"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/analysis/umap"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/analysis/umap"))
}
if (!dir.exists(paste0("./analysis/shiny/",i,"/outs/analysis/umap/2_components"))){
dir.create(paste0("./analysis/shiny/",i,"/outs/analysis/umap/2_components"))
}
}
}
......@@ -971,7 +976,21 @@ scQuSAGE <- function(sc10x,gs,save=FALSE,type,id,ds=0,nm="pops",print="tsne"){
return(results)
}
scShinyOutput <- function(){
scShinyOutput <- function(sc10x,anal="raw"){
write_delim(as.data.frame(colnames(sc10x)),path=paste0("./analysis/shiny/",anal,"/outs/filtered_feature_bc_matrix/barcodes.tsv.gz"),delim="\t",col_names=FALSE)
features <- rownames(sc10x)
features <- c(features,c("nFeature","nCount","percent.mito","percent.ribo","Stress.score"))
features <- data.frame(ENSG=features,Feature=features,Label="feature")
write_delim(features,path=paste0("./analysis/shiny/",anal,"/outs/filtered_feature_bc_matrix/features.tsv.gz"),delim="\t",col_names=FALSE)
exp <- GetAssayData(sc10x,slot="scale.data")
exp.extra <- matrix(nrow=5,ncol=ncol(sc10x))
exp.extra[1,] <- as.numeric(sc10x$nFeature_RNA)
exp.extra[2,] <- as.numeric(sc10x$nCount_RNA)
exp.extra[3,] <- as.numeric(sc10x$percent.mito)
exp.extra[4,] <- as.numeric(sc10x$percent.ribo)
exp.extra[5,] <- as.numeric(sc10x$Stress1)
exp <- rbind(exp,exp.extra)
Matrix::writeMM(as(exp,"dgCMatrix"),file=paste0("./analysis/shiny/",anal,"/outs/filtered_feature_bc_matrix/matrix.mtx.gz"))
for (i in c("pca","tsne","umap")){
dr <- Embeddings(sc10x,i)
if (i != "pca"){
......@@ -984,23 +1003,38 @@ scShinyOutput <- function(){
dr <- dr[,c(3,1,2)]
dr <- as.data.frame(dr,row.names=FALSE)
if (i != "pca"){
write_csv(dr,paste0("./analysis/shiny/raw/outs/analysis/",i,"/2_components/projection.csv"),col_names=TRUE)
write_csv(dr,paste0("./analysis/shiny/",anal,"/outs/analysis/",i,"/2_components/projection.csv"),col_names=TRUE)
} else {
write_csv(dr,paste0("./analysis/shiny/raw/outs/analysis/",i,"/10_components/projection.csv"),col_names=TRUE)
write_csv(dr,paste0("./analysis/shiny/",anal,"/outs/analysis/",i,"/10_components/projection.csv"),col_names=TRUE)
}
}
for (i in c("samples",paste0("integrated_snn_res.",res))){
if (!dir.exists(paste0("./analysis/shiny/raw/outs/analysis/clustering/",gsub("integrated_snn_res.","res_",i)))){
dir.create(paste0("./analysis/shiny/raw/outs/analysis/clustering/",gsub("integrated_snn_res.","res_",i)))
sc10x <- NormalizeData(sc10x,assay="RNA")
clusters <- c("samples",paste0("integrated_snn_res.",res),"lin","pops","leu","scDWSpr")
clusters <- intersect(clusters,names(sc10x@meta.data))
for (i in clusters){
if (!dir.exists(paste0("./analysis/shiny/",anal,"/outs/analysis/clustering/",gsub("integrated_snn_res.","res_",i)))){
dir.create(paste0("./analysis/shiny/",anal,"/outs/analysis/clustering/",gsub("integrated_snn_res.","res_",i)))
}
clust <- as.matrix(sc10x[[i]])
colnames(clust) <- "Cluster"
clust <- cbind(clust,Barcode=rownames(clust))
clust <- clust[,c(2,1)]
clust <- as.data.frame(clust,row.names=FALSE)
write_csv(clust,paste0("./analysis/shiny/raw/outs/analysis/clustering/",gsub("integrated_snn_res.","res_",i),"/clusters.csv"),col_names=TRUE)
clust[,2] <- paste0("Cluster ",clust[,2])
write_csv(clust,paste0("./analysis/shiny/",anal,"/outs/analysis/clustering/",gsub("integrated_snn_res.","res_",i),"/clusters.csv"),col_names=TRUE)
if (!dir.exists(paste0("./analysis/shiny/",anal,"/outs/analysis/diffexp/",gsub("integrated_snn_res.","res_",i)))){
dir.create(paste0("./analysis/shiny/",anal,"/outs/analysis/diffexp/",gsub("integrated_snn_res.","res_",i)))
}
Idents(sc10x) <- i
deg <- FindAllMarkers(sc10x,assay="RNA",slot="data",logfc.threshold=0,test.use="MAST",min.pct=0.25,min.diff.pct=0.25,max.cells.per.ident=500)
dexp <- data.frame("Feature ID"=unique(deg$gene),"Feature Name"=unique(deg$gene))
for (cluster in unique(deg$cluster)){
dexp[,paste0("Cluster.",cluster,".Mean.Counts")] <- deg$pct.1[deg$cluster==cluster][match(dexp$Feature.ID,deg$gene[deg$cluster==cluster])]
dexp[,paste0("Cluster.",cluster,".Log2.fold.change")] <- deg$avg_logFC[deg$cluster==cluster][match(dexp$Feature.ID,deg$gene[deg$cluster==cluster])]
dexp[,paste0("Cluster.",cluster,".Adjusted.p.value")] <- deg$p_val_adj[deg$cluster==cluster][match(dexp$Feature.ID,deg$gene[deg$cluster==cluster])]
}
colnames(dexp) <- gsub("\\."," ",colnames(dexp))
write_csv(dexp,paste0("./analysis/shiny/",anal,"/outs/analysis/diffexp/",gsub("integrated_snn_res.","res_",i),"/differential_expression.csv"),col_names=TRUE)
}
#Idents(sc10x) <- "integrated_snn_res.0.1"
#deg <- FindAllMarkers(sc10x,assay="SCT",slot="scale.data",logfc.threshold=0,test.use="MAST")
}
}
\ No newline at end of file
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