Commit 8de3b69e authored by Gervaise Henry's avatar Gervaise Henry 🤠
Browse files

Correct r.scripts

parent ce300d01
library(Seurat)
library(biomaRt)
options(future.globals.maxSize= 1000000000000)
setwd("../")
sc10x <- readRDS("./analysis/huPr_PdPb_id_all.rds")
sc10x.fmst <- readRDS("./analysis/huPr_PdPb_id_fmst.rds")
DefaultAssay(sc10x.fmst) <- "SCT"
Idents(sc10x.fmst) <- "resolution_0.2"
sc10x.fmst <- RenameIdents(sc10x.fmst,"6"="Prostate SM")
sc10x.fmst <- RenameIdents(sc10x.fmst,"2"="Prostate SM")
sc10x.fmst <- RenameIdents(sc10x.fmst,"1"="Pericytes")
sc10x.fmst <- RenameIdents(sc10x.fmst,"3"="vSM")
sc10x.fmst <- RenameIdents(sc10x.fmst,"5"="Interstitial Fibroblasts")
sc10x.fmst <- RenameIdents(sc10x.fmst,"0"="Interstitial Fibroblasts")
sc10x.fmst <- RenameIdents(sc10x.fmst,"4"="Peri-epithelial Fibroblasts")
sc10x.fmst$Populations <- Idents(sc10x.fmst)
Idents(sc10x.fmst) <- "Populations"
sc10x$pop[names(sc10x.fmst$Populations)] <- as.character(sc10x.fmst$Populations)
rm(sc10x.fmst)
as_matrix <- function(mat){
tmp <- matrix(data=0L, nrow = mat@Dim[1], ncol = mat@Dim[2])
row_pos <- mat@i+1
col_pos <- findInterval(seq(mat@x)-1,mat@p[-1])+1
val <- mat@x
for (i in seq_along(val)){
tmp[row_pos[i],col_pos[i]] <- val[i]
}
row.names(tmp) <- mat@Dimnames[[1]]
colnames(tmp) <- mat@Dimnames[[2]]
return(tmp)
}
count_norm <- as.data.frame(as_matrix(GetAssayData(sc10x,assay="SCT",slot="data")))
count_norm$Symbol <- rownames(count_norm)
ensembl <- useMart("ensembl",dataset="hsapiens_gene_ensembl",host="https://sep2019.archive.ensembl.org")
symbol2ID <- getBM(mart = ensembl,attributes = c("external_gene_name","ensembl_gene_id"),filters="external_gene_name",values=count_norm$Symbol,uniqueRows=TRUE)
colnames(symbol2ID)[colnames(symbol2ID)=="ensembl_gene_id"] <- "Genes"
count_norm <- merge(x=count_norm,y=symbol2ID,by.x="Symbol",by.y="external_gene_name",all.x=TRUE)
count_norm$Symbol <- NULL
count_norm <- count_norm[,c(which(colnames(count_norm)=="Genes"),which(colnames(count_norm)!="Genes"))]
write.table(count_norm,'./analysis/cellphonedb_count.txt',sep='\t',quote=FALSE,row.names=FALSE,col.names=TRUE)
meta_data <- data.frame(Cell=colnames(sc10x),cell_type=as.character(sc10x$pop))
write.table(meta_data,'./analysis/cellphonedb_meta.txt',sep='\t',quote=FALSE,row.names=FALSE,col.names=TRUE)
This diff is collapsed.
This diff is collapsed.
gc()
library(Seurat)
library(RColorBrewer)
library(viridis)
library(gridExtra)
library(ggplot2)
library(qusage)
library(VennDiagram)
library(corrplot)
setwd("../")
source("./r.scripts/sc-TissueMapper_functions.R")
source("./r.scripts/sc-TissueMapper_process.R")
sc10x.fmst <- readRDS("./analysis/muPrUr_id_fmst.rds")
DefaultAssay(sc10x.fmst) <- "SCT"
sc10x.fmst$Region <- sc10x.fmst$samples
pr <- unique(sc10x.fmst$Region)[grepl("PrF",unique(sc10x.fmst$Region))]
ur <- unique(sc10x.fmst$Region)[grepl("UrF",unique(sc10x.fmst$Region))]
sc10x.fmst$Region[sc10x.fmst$Region %in% pr] <- "Prostate"
sc10x.fmst$Region[sc10x.fmst$Region %in% ur] <- "Urethra"
Idents(sc10x.fmst) <- "resolution_0.1"
sc10x.fmst <- RenameIdents(sc10x.fmst,"4"="SM")
sc10x.fmst <- RenameIdents(sc10x.fmst,"2"="SM")
sc10x.fmst <- RenameIdents(sc10x.fmst,"3"="Ductal Fibroblasts")
sc10x.fmst <- RenameIdents(sc10x.fmst,"0"="Urethral Fibroblasts")
sc10x.fmst <- RenameIdents(sc10x.fmst,"1"="Prostate Fibroblasts")
sc10x.fmst$Populations <- Idents(sc10x.fmst)
Idents(sc10x.fmst) <- "Populations"
sc10x.fib <- subset(sc10x.fmst,idents=c("Prostate Fibroblasts","Urethral Fibroblasts","Ductal Fibroblasts"))
postscript("./analysis/UMAP.muPrUr_fmst.pops.eps")
DimPlot(sc10x.fmst,group.by="Populations",reduction="umap",label=TRUE,repel=TRUE,label.size=10,pt.size=2)+theme(legend.position="none")+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
postscript("./analysis/UMAP.muPrUr_fmst.region.eps")
DimPlot(sc10x.fmst,group.by="Region",reduction="umap",label=FALSE,repel=TRUE,label.size=10,pt.size=2)+theme(legend.position="right",legend.text=element_text(size=25))+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
postscript("./analysis/DotPlot.muPrUr_fmst.eps")
DotPlot(sc10x.fmst,features=c("Dcn","C3","Lgr5","Wnt2","Acta2"),cols=rev(heat.colors(2)))+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25,angle=45,hjust=1),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
postscript("./analysis/DotPlot.muPrUr_fib.eps")
DotPlot(sc10x.fib,features=c("C3","Igf1","Sult1e1","Gpx3","Lgr5","Apoe","Osr1","Sfrp2","Wnt2","Rorb","Srd5a2","Wif1"),cols=rev(heat.colors(2)))+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25,angle=45,hjust=1),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
postscript("./analysis/VlnPlot.muPrUr_fib_Tgfb1.eps")
VlnPlot(sc10x.fib,features="Tgfb1",pt.size=2)+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
postscript("./analysis/VlnPlot.muPrUr_fib_Tgfb2.eps")
VlnPlot(sc10x.fib,features="Tgfb2",pt.size=2)+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
postscript("./analysis/VlnPlot.muPrUr_fib_Tgfb3.eps")
VlnPlot(sc10x.fib,features="Tgfb3",pt.size=2)+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
genes.wnt <- c("SFRP1","SFRP2","SFRP4","SFRP5",'WNT10A',
"WNT11","WNT16","WNT2","WNT2B","WNT3A",
"WNT4","WNT5A","WNT5B","WNT6",
"WNT7B","WNT9A","WIF1","AXIN2",
"FZD1","FZD10","FZD2","FZD3", "FZD4",
"FZD5","FZD6","FZD7","FZD8","FZD9",
"APC","CTNNB1","DVL1","LEF1","PORCN") #WNT5A-AS1,WNT7A,DKK1
genes.wnt <- paste0(toupper(substr(genes.wnt,1,1)),substr(tolower(genes.wnt),2,nchar(genes.wnt)))
postscript("./analysis/DotPlot.muPrUr_fib_wnt.eps")
DotPlot(sc10x.fib,features=genes.wnt,cols=rev(heat.colors(2)),assay="SCT")+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25,angle=45,hjust=1),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
postscript("./analysis/Heatmap.muPrUr_fib_wnt.eps")
DoHeatmap(sc10x.fib,features=genes.wnt,assay = "SCT")
dev.off()
genes.tgfb <- c("BMP2","BMP4","BMP7","INHBA",
"TGFB1","TGFB2","TGFB3","TGFBR1","TGFBR2",
"TGFBR3")#,"MIS"
genes.tgfb <- paste0(toupper(substr(genes.tgfb,1,1)),substr(tolower(genes.tgfb),2,nchar(genes.tgfb)))
postscript("./analysis/DotPlot.muPrUr_fib_tgfb.eps")
DotPlot(sc10x.fib,features=genes.tgfb,cols=rev(heat.colors(2)),assay="SCT")+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25,angle=45,hjust=1),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
postscript("./analysis/Heatmap.muPrUr_fib_tgfb.eps")
DoHeatmap(sc10x.fib,features=genes.tgfb,assay = "SCT")
dev.off()
genes.wnts <- c(rownames(sc10x.fib)[startsWith(rownames(sc10x.fib),"Wnt")],rownames(sc10x.fib)[startsWith(rownames(sc10x.fib),"Rspo")])
sc10x.fib <- AddModuleScore(sc10x.fib,features=genes.wnts,name="Wnt_Score",assay="SCT",nbin=15)
postscript("./analysis/DotPlot.muPrUr_fib_wnt.score.eps")
DotPlot(sc10x.fib,group.by="Populations",features="Wnt_Score1",cols=rev(heat.colors(2)),assay="SCT")+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25,angle=45,hjust=1),axis.text.y=element_text(size=25,angle=45),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
postscript("./analysis/VlnPlot.muPrUr_fib_wnt.score.eps")
VlnPlot(sc10x.fib,group.by="Populations",features="Wnt_Score1",pt.size=0.1)+theme(plot.title=element_text(size=50),axis.text.x=element_text(size=25),axis.text.y=element_text(size=25),axis.title.x=element_text(size=40),axis.title.y=element_text(size=40))
dev.off()
deg.fmst <- FindAllMarkers(sc10x.fmst,assay="SCT",slot="data",logfc.threshold=0,min.pct=0.25,test.use="MAST",only.pos=TRUE)
write.table(deg.fmst,file=paste0("./analysis/DEG.muPrUr_fmst.csv"),row.names=FALSE,col.names=TRUE,append=FALSE,sep=",")
deg.fib <- FindAllMarkers(sc10x.fib,assay="SCT",slot="data",logfc.threshold=0,min.pct=0.25,test.use="MAST",only.pos=TRUE)
write.table(deg.fib,file=paste0("./analysis/DEG.muPrUr_fib.csv"),row.names=FALSE,col.names=TRUE,append=FALSE,sep=",")
sc10x.fib <- RenameIdents(sc10x.fib,"Ductal Fibroblasts"="dFibroblasts")
sc10x.fib <- RenameIdents(sc10x.fib,"Urethral Fibroblasts"="uFibroblasts")
sc10x.fib <- RenameIdents(sc10x.fib,"Prostate Fibroblasts"="pFibroblasts")
sc10x.fib$Pops <- Idents(sc10x.fib)
Idents(sc10x.fib) <- "Pops"
gene.set.c5.bp <- read.gmt("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/genesets_Mouse_GOBP_AllPathways_no_GO_iea_March_01_2021_symbol.gmt")
neg <- names(gene.set.c5.bp)[grepl("NEGATIVE",names(gene.set.c5.bp))]
neg <- neg[! neg %in% names(gene.set.c5.bp)[grepl("GRAM NEGATIVE",names(gene.set.c5.bp))]]
gene.set.c5.bp <- gene.set.c5.bp[! names(gene.set.c5.bp) %in% neg]
rm(neg)
results <- scQuSAGE(sc10x.fib,gs=gene.set.c5.bp,save=FALSE,type="lg",id="Pops",ds=0,nm="fib",print="0")
results.cor.c5.bp <- results[[2]]
results.clust.c5.bp <- results[[3]]
rm(results)
write.table(results.cor.c5.bp,file=paste0("./analysis/CorResults.muPrUr_fib.go_bp.fib.csv"),row.names=FALSE,col.names=TRUE,append=FALSE,sep=",")
ortholog <- read.delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/genesets_Ensemble.mus-hum_GRCm38.p6.txt")
ortholog <- ortholog[order(ortholog[,"Human.gene.name"]),]
ortholog <- ortholog[!duplicated(ortholog$Human.gene.name),]
sc10x.human.fmst <- readRDS("../huPr_PdPb/analysis/huPr_PdPb_id_fmst.rds")
DefaultAssay(sc10x.human.fmst) <- "SCT"
bph <- unique(sc10x.human.fmst$samples)[startsWith(unique(sc10x.human.fmst$samples),"BPH")]
donor <- unique(sc10x.human.fmst$samples)[startsWith(unique(sc10x.human.fmst$samples),"D")]
sc10x.human.fmst$Phenotype <- sc10x.human.fmst$samples
Idents(sc10x.human.fmst) <- "Phenotype"
sc10x.human.fmst$Phenotype[sc10x.human.fmst$Phenotype %in% bph] <- "BPH"
sc10x.human.fmst$Phenotype[sc10x.human.fmst$Phenotype %in% donor] <- "Donor"
Idents(sc10x.human.fmst) <- "Phenotype"
sc10x.human.fmst.donor <- subset(sc10x.human.fmst,idents="Donor")
genes.mouse <- rownames(sc10x.fmst)
genes.human <- rownames(sc10x.human.fmst)
genes.human.donor <- rownames(sc10x.human.fmst.donor)
genes.mouse.filtered <- genes.mouse[genes.mouse %in% ortholog$Gene.name]
genes.human.filtered <- genes.human[genes.human %in% ortholog$Human.gene.name]
genes.human.donor.filtered <- genes.human.donor[genes.human.donor %in% ortholog$Human.gene.name]
genes.mouse.filtered.humanified <- ortholog$Human.gene.name[match(genes.mouse.filtered,ortholog$Gene.name)]
genes.mouse.filtered.humanified <- genes.mouse.filtered.humanified[genes.mouse.filtered.humanified %in% genes.human.filtered]
genes.human.filtered <- genes.human.filtered[genes.human.filtered %in% genes.mouse.filtered.humanified]
genes.human.donor.filtered <- genes.human.donor.filtered[genes.human.donor.filtered %in% genes.mouse.filtered.humanified]
genes.mouse.filtered <- ortholog$Gene.name[match(genes.mouse.filtered.humanified,ortholog$Human.gene.name)]
ortholog <- ortholog[ortholog$Gene.name %in% genes.mouse.filtered,]
ortholog <- ortholog[ortholog$Human.gene.name %in% genes.human.filtered,]
data.mouse.fmst <- as.data.frame(AverageExpression(sc10x.fmst,assays="SCT",slot="scale.data")$SCT)
data.mouse.fmst <- data.mouse.fmst[rownames(data.mouse.fmst) %in% genes.mouse.filtered,]
rownames(data.mouse.fmst) <- ortholog$Human.gene.name[match(rownames(data.mouse.fmst),ortholog$Gene.name)]
data.mouse.fmst <- data.mouse.fmst[ortholog$Human.gene.name,]
Idents(sc10x.human.fmst) <- "resolution_0.2"
sc10x.human.fmst <- RenameIdents(sc10x.human.fmst,"6"="pSM")
sc10x.human.fmst <- RenameIdents(sc10x.human.fmst,"2"="pSM")
sc10x.human.fmst <- RenameIdents(sc10x.human.fmst,"1"="Pericytes")
sc10x.human.fmst <- RenameIdents(sc10x.human.fmst,"3"="vSM")
sc10x.human.fmst <- RenameIdents(sc10x.human.fmst,"5"="iFib")
sc10x.human.fmst <- RenameIdents(sc10x.human.fmst,"0"="iFib")
sc10x.human.fmst <- RenameIdents(sc10x.human.fmst,"4"="peFib")
sc10x.human.fmst$Populations <- Idents(sc10x.human.fmst)
Idents(sc10x.human.fmst) <- "Populations"
data.human.fmst <- as.data.frame(AverageExpression(sc10x.human.fmst,assays="SCT",slot="scale.data")$SCT)
data.human.fmst <- data.human.fmst[rownames(data.human.fmst) %in% genes.human.filtered,]
data.human.fmst <- data.human.fmst[ortholog$Human.gene.name,]
Idents(sc10x.human.fmst.donor) <- "resolution_0.2"
sc10x.human.fmst.donor <- RenameIdents(sc10x.human.fmst.donor,"6"="pSM")
sc10x.human.fmst.donor <- RenameIdents(sc10x.human.fmst.donor,"2"="pSM")
sc10x.human.fmst.donor <- RenameIdents(sc10x.human.fmst.donor,"1"="Pericytes")
sc10x.human.fmst.donor <- RenameIdents(sc10x.human.fmst.donor,"3"="vSM")
sc10x.human.fmst.donor <- RenameIdents(sc10x.human.fmst.donor,"5"="iFib")
sc10x.human.fmst.donor <- RenameIdents(sc10x.human.fmst.donor,"0"="iFib")
sc10x.human.fmst.donor <- RenameIdents(sc10x.human.fmst.donor,"4"="peFib")
sc10x.human.fmst.donor$Populations <- Idents(sc10x.human.fmst.donor)
Idents(sc10x.human.fmst.donor) <- "Populations"
data.human.fmst.donor <- as.data.frame(AverageExpression(sc10x.human.fmst.donor,assays="SCT",slot="scale.data")$SCT)
data.human.fmst.donor <- data.human.fmst.donor[rownames(data.human.fmst.donor) %in% genes.human.filtered,]
data.human.fmst.donor <- data.human.fmst.donor[ortholog$Human.gene.name,]
correlation.fmst <- cor(data.mouse.fmst,data.human.fmst,use="complete.obs")
correlation.fmst.scale <- ((correlation.fmst-min(correlation.fmst))/(max(correlation.fmst)-min(correlation.fmst)))*(1-(-1))-1
postscript("./analysis/Corr.muhu_fmst.eps")
corrplot(t(correlation.fmst.scale),method="color",addCoef.col="white",tl.srt=45,tl.col="black",title="Human:Mouse FMSt Correlations",mar=c(0,0,2,0),is.corr=FALSE,col=colorRampPalette(c("blue","black","red"))(100))
dev.off()
correlation.fmst.donor <- cor(data.mouse.fmst,data.human.fmst.donor,use="complete.obs")
correlation.fmst.donor.scale <- ((correlation.fmst.donor-min(correlation.fmst.donor))/(max(correlation.fmst.donor)-min(correlation.fmst.donor)))*(1-(-1))-1
postscript("./analysis/Corr.muhu_fmst.donor.eps")
corrplot(t(correlation.fmst.donor.scale),method="color",addCoef.col="white",tl.srt=45,tl.col="black",title="Human:Mouse FMSt Correlations",mar=c(0,0,2,0),is.corr=FALSE,col=colorRampPalette(c("blue","black","red"))(100))
dev.off()
postscript("./analysis/Corr.muhu_fmst.donor_unscaled.eps")
corrplot(t(correlation.fmst.donor),method="color",addCoef.col="white",tl.srt=45,tl.col="black",title="Human:Mouse FMSt Correlations",mar=c(0,0,2,0),is.corr=FALSE,col=colorRampPalette(c("blue","black","red"))(100))
dev.off()
load("./analysis/sc10x.raw.RData")
filter <- data.frame(raw=unlist(counts.cell.raw),umi=unlist(counts.cell.filtered.1umi),mito=unlist(counts.cell.filtered.2mito),ribo=unlist(counts.cell.filtered.3ribo),gene=unlist(counts.cell.filtered.4gene),stress=unlist(counts.cell.destress))
write.table(filter,'./analysis/cells.txt',sep='\t',quote=FALSE,row.names=TRUE,col.names=NA)
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