Commit 7fbfe85b authored by Gervaise H. Henry's avatar Gervaise H. Henry 🤠

Add prelim HTO exploration

parent 6e7d9f40
library(Matrix)
library(ggplot2)
exp <- readMM("matrix.mtx.gz")
features <- read.table("features.tsv.gz", quote="\"")
ab <- as.data.frame(as.matrix(t(exp[features$V4 == "Capture",])))
colnames(ab) <- features[features$V4 == "Capture",2]
ab$total <- rowSums(ab)
ab <- transform(ab,AP.ratio=(Anterior_Prostate/total))
ab <- transform(ab,VP.ratio=(Ventral_Prostate/total))
ab <- transform(ab,DL.ratio=(Dorsolateral_Prostate/total))
ab <- transform(ab,cutoff=total*0.50)
ab$AP.ident <- (ab$Anterior_Prostate >= ab$cutoff)+0
ab$VP.ident <- (ab$Ventral_Prostate >= ab$cutoff)+0
ab$DL.ident <- (ab$Dorsolateral_Prostate >= ab$cutoff)+0
ab <- transform(ab,Ab.ident=Ab1.ident+Ab2.ident)
ab$Ab.ident[ab$ident==0] <- "Ambiguous"
ab$Ab.ident[ab$ident==2] <- "Ambiguous"
ab$Ab.ident[ab$ident==3] <- "Ambiguous"
ab$Ab.ident[ab$ident==1] <-
ab$Ab.ident[ab$Ab.ident==1] <- "Ab1"
ab$Ab.ident[ab$Ab.ident==-1] <- "Ab2"
ggplot(ab,aes(x=Ab1,y=Ab2,col=total))+geom_point()+scale_x_continuous(trans='log2')+scale_y_continuous(trans='log2')+scale_color_gradient(low="blue",high="red")
ggplot(ab,aes(x=Ab1,y=Ab2,col=Ab.ident))+geom_point()+scale_x_continuous(trans='log2')+scale_y_continuous(trans='log2')
table(ab$Ab.ident)
gc()
library(methods)
library(optparse)
library(Seurat)
library(readr)
library(readxl)
library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(monocle)
library(dplyr)
library(viridis)
library(gridExtra)
library(SingleR)
library(sctransform)
library(autothresholdr)
options(bitmapType="cairo")
setwd("../")
source("./r.scripts/sc-TissueMapper_functions.R")
source("./r.scripts/sc-TissueMapper_process.R")
project.name="musAdPrF"
scFolders()
# sc10x <- scLoad(p=project.name,cellranger=3,ref="mm10")
sc10x.data <- Read10X(data.dir=paste0("./analysis/DATA//10x/filtered_feature_bc_matrix/"))
#sc10x.data$`Antibody Capture`@Dimnames[[1]] <- c("musAd004","musAd005")
sc10x <- CreateSeuratObject(counts=sc10x.data$`Gene Expression`,project="musAd004n5_UrF")
sc10x[["HTO"]] <- CreateAssayObject(counts=sc10x.data$`Antibody Capture`)
results <- scQC(list(all=sc10x),sp="mu")
sc10x <- results[[1]]
counts.cell.raw <- results[[2]]
counts.gene.raw <- results[[3]]
counts.cell.filtered <- results[[4]]
counts.gene.filtered <- results[[5]]
rm(results)
sc10x <- sc10x$all
sc10x <- NormalizeData(object=sc10x,assay="HTO",normalization.method="CLR")
sc10x <- HTODemux(sc10x,assay="HTO",positive.quantile=.99)
# FeatureScatter(sc10x,"hto_musAd004","hto_musAd005",slot="data",group.by = "hash.ID")
# table(sc10x$hash.ID)
HTOHeatmap(sc10x, assay = "HTO")
plot1 <- FeatureScatter(sc10x,feature1 = "Anterior-Prostate",feature2 = "Ventral-Prostate")
plot2 <- FeatureScatter(sc10x,feature1 = "Anterior-Prostate",feature2 = "Dorsolateral-Prostate")
plot3 <- FeatureScatter(sc10x,feature1 = "Dorsolateral-Prostate",feature2 = "Ventral-Prostate")
grid.arrange(plot1,plot2,plot3,ncol=1)
Idents(sc10x) <- "hash.ID"
sc10x <- subset(sc10x,idents="Doublet",invert=TRUE)
sc10x <- RenameIdents(sc10x,"Negative"="musAd004+5_UrF")
sc10x <- RenameIdents(sc10x,"musAd004"="musAd004_UrF")
sc10x <- RenameIdents(sc10x,"musAd005"="musAd005_UrF")
sc10x$samples.dmltpx <- Idents(sc10x)
# results <- scCellCycle(sc10x,sub=FALSE,sp="mu")
# sc10x <- results[[1]]
# genes.s <- results[[2]]
# genes.g2m <- results[[3]]
# rm(results)
sc10x.l <- list()
sc10x.l[["mu1"]] <- subset(sc10x,subset= samples=="musAd001_PrF")
sc10x.l[["mu2"]] <- subset(sc10x,subset= samples=="musAd002_PrF")
sc10x.l[["mu3"]] <- subset(sc10x,subset= samples=="musAd003_PrF_St")
sc10x <- scCCA(sc10x.l)
rm(sc10x.l)
sc10x@project.name <- project.name
# gc()
# sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nFeature_RNA","percent.mito","S.Score","G2M.Score"),do.par=TRUE,num.cores=45,verbose=FALSE)
sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nFeature_RNA","percent.mito"),do.par=TRUE,num.cores=45,verbose=FALSE)
# gc()
results <- scPC(sc10x,pc=50,hpc=0.9,file="pre.stress",print="2",cca=TRUE)
sc10x <- results[[1]]
pc.use.prestress <- results[[2]]
rm(results)
# sc10x <- scCluster(sc10x,res=0.5,red="pca",dim=pc.use.prestress,print="2",folder="pre.stress")
# genes.stress <- read_delim("./genesets/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
# genes.stress <- genes.stress[1]
# colnames(genes.stress) <- "scDWS.Stress"
# results <- scScore(sc10x,score="Stress",geneset=as.list(genes.stress),cut.pt=0.75,anchor=c("EGR1","FOS","JUN"))
# sc10x.preStress <- results[[1]]
# sc10x <- results[[2]]
# rm(results)
# gc()
# sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nFeature_RNA","percent.mito","S.Score","G2M.Score"),do.par=TRUE,num.cores=45,verbose=FALSE)
# sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nFeature_RNA","percent.mito"),do.par=TRUE,num.cores=45,verbose=FALSE)
# gc()
# results <- scPC(sc10x,pc=50,hpc=0.9,file="post.stress",print="2",cca=TRUE)
# sc10x <- results[[1]]
# pc.use.poststress <- results[[2]]
# rm(results)
res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))
sc10x <- scCluster(sc10x,res=res,red="pca",dim=pc.use.prestress,print="2",folder="ALL")
save(sc10x,file=paste0("./analysis/sc10x.raw.rda"))
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