.Rhistory

library(Seurat)
InstallData("stxBrain")
??InstallData
library(seurat-data)
library(SeuratData)
devtools::install_github('satijalab/seurat-data')
# devtools::install_github('satijalab/seurat-data')
library(SeuratData)
InstallData("stxBrain")
brain <- LoadData("stxBrain", type = "anterior1")
brain
library(qs)
BiocManager::install("qdata")
BiocManager::install("qs2")
library(qs2)
qd_save(brain, file="visium/visium.qs")
qsave(brain, file="visium/visium.qs")
test <- qread("visium/visium.qs")
# This set up the working directory to this file so all files can be found
library(rstudioapi)
setwd(fs::path_dir(getSourceEditorContext()$path))
stopifnot(R.version$major>= 4) # requires R4
if (compareVersion(R.version$minor,"3.1")<0) warning("We recommend >= R4.3.1")
stopifnot(compareVersion(as.character(BiocManager::version()), "3.16")>=0)
stopifnot(compareVersion(as.character(packageVersion("Seurat")), "5.1")>=0)
library(knitr)
library(glue)
library(qs2)
library(dplyr)
library(purrr)
library(ggplot2)
library(ggprism)
library(grafify)
library(ggpubr)
library(gridExtra)
library(scales)
library(Seurat)
import::from(magrittr,set_colnames,set_rownames,"%<>%")
library(import)
BiocManager::install("import")
library(knitr)
library(import)
library(glue)
library(qs2)
library(dplyr)
library(purrr)
library(ggplot2)
library(ggprism)
library(grafify)
library(ggpubr)
library(gridExtra)
library(scales)
library(Seurat)
import::from(magrittr, set_colnames, set_rownames, "%<>%")
invisible(list2env(params,environment()))
source(project_file)
ggplot2::theme_set(theme_prism(base_size = 12))
# https://grafify-vignettes.netlify.app/colour_palettes.html
# NOTE change colors here if you wish
scale_colour_discrete <- function(...)
scale_colour_manual(...,
values = as.vector(grafify:::graf_palettes[["kelly"]]))
scale_fill_discrete <- function(...)
scale_fill_manual(...,
values = as.vector(grafify:::graf_palettes[["kelly"]]))
opts_chunk[["set"]](
cache = F,
cache.lazy = FALSE,
dev = c("png", "pdf"),
error = TRUE,
highlight = TRUE,
message = FALSE,
prompt = FALSE,
tidy = FALSE,
warning = FALSE,
echo = T,
fig.height = 4)
# set seed for reproducibility
set.seed(1234567890L)
# Metrics like `nCount` and `nfeature` are named with the suffix of default assay name, to make the variable usage more generalizable, we removed the suffix by pulling out the default assay of the visium `seurat` object.
visium <- qread(visiumHD_obj)
visium <- PercentageFeatureSet(visium, "^mt-", col.name = "percent_mito")
visium <- PercentageFeatureSet(visium, "^Hb.*-", col.name = "percent_hb")
metaD <- visium@meta.data
metaD$log10GenesPerUMI <- log10(metaD$nFeature)/log10(metaD$nCount)
colnames(metaD)%<>%gsub(pattern=glue("_{DefaultAssay(visium)}"),replacement="")
summary_metaD <- apply(metaD[,-1],2,mean)
metacol_label <- list("nFeature"="Genes","nCount"="UMI")
refs <- list("nFeature"=100,"nCount"=100)
dists_before <- imap(metacol_label,\(label,col)
ggdensity(metaD,
x = col,xscale="log10",add = "mean", rug = TRUE,
alpha = 0.2,fill = "lightgray",
xlab=glue("Number of {label} per bin(in log10 scale)"),
ylab="Cell density",
title=glue('Pre-QC {label}/Bin'))+
geom_vline(xintercept = refs[[col]],color="darkred",cex=rel(1.3),linetype="dashed")+
annotate("text",x=summary_metaD[col],y = Inf,
label = glue("Mean \n = {round(summary_metaD[col],0)}"),
vjust = 1,hjust=2)
)
dists_before[[1]] | dists_before[[2]]
col <- "log10GenesPerUMI"
ggdensity(metaD,x = col,add = "mean", rug = TRUE,
alpha = 0.2,fill = "lightgray",
xlab="complexity",ylab="Cell density",title=glue('Novelty score'))+
geom_vline(xintercept = 0.8,color="darkred",cex=rel(1.3),linetype="dashed")+
annotate("text",x=summary_metaD[col],y = Inf,
label = glue("Mean = {round(summary_metaD[col],0)}"),
vjust = 1,hjust=2)+
theme(plot.title = element_text(hjust=0.5, face="bold"))
knit_with_parameters("~/Science/nfcore/templates/spatial-reports/visium/01_quality_assessment/qc.Rmd")
ggplot(metaD %>%
select(orig.ident,starts_with("percent_")) %>%
tidyr::gather(class,percent_unexpected,-orig.ident),
aes_string(x = "orig.ident", y = "percent_unexpected")) +
geom_violin(position=position_dodge(1),alpha=1, na.rm=TRUE,trim=FALSE)+
ggbeeswarm::geom_quasirandom(na.rm=TRUE,dodge.width=0.5,
method='quasirandom',alpha=0.01)+
geom_boxplot(width=0.1,outliers = F)+
geom_hline(yintercept=20)+
facet_grid(~class,scales = "free")+
theme(
axis.text.x = element_text(size=rel(1),face="bold"),
plot.title = element_text(hjust = 0.5),
strip.text.x = element_text(size = rel(1.5), colour = "black"),
legend.position = "none"
)+
scale_y_log10(breaks=c(1,5,10,20,100))+
# ylim(c(0,100))+
labs(x="",y="% of contamination genes")
features2check <- c(glue('nCount_{DefaultAssay(visium)}'),
glue('nFeature_{DefaultAssay(visium)}'),
"percent_mito","percent_hb")
for(f in features2check){
cat("### ", f, "\n\n")
p1 <- SpatialFeaturePlot(visium,
feature = f,
pt.size.factor = 8)
print(p1)
cat("\n\n")
}
for(f in features2check){
cat("### ", f, "\n\n")
p1 <- SpatialFeaturePlot(visium,
feature = f,
pt.size.factor = 4)
print(p1)
cat("\n\n")
}
GeneVar <- glue('nFeature_{DefaultAssay(visium)}')
UMIVar <- glue('nCount_{DefaultAssay(visium)}')
cutoffs <- list("nFeature"=100,"nCount"=100,"hb"=20,"mito"=20)
Qced <-  visium@meta.data[,GeneVar] > cutoffs$nFeature &
visium@meta.data[,UMIVar] > cutoffs$nCount &
visium$percent_hb < cutoffs$hb &
visium$percent_mito < cutoffs$mito
visium <- visium[,Qced]
# Filter Mitocondrial
visium <- visium[!grepl("^mt-", rownames(visium)), ]
# Filter Hemoglobin gene (optional if that is a problem on your data)
visium <- visium[!grepl("^Hb.*-", rownames(visium)), ]
C <- GetAssayData(visium, slot = "counts")
C@x <- C@x / rep.int(colSums(C), diff(C@p))
most_expressed <- order(Matrix::rowSums(C), decreasing = T)[20:1]
exprD <- as.data.frame(t(C[most_expressed, ])) %>%
tibble::rownames_to_column("bin") %>%
tidyr::gather(gene,expr,-bin)
ggplot(exprD,aes(x=gene,y=expr,color=gene,fill=gene))+
geom_violin(position=position_dodge(1),alpha=0.5,
na.rm=TRUE,trim=FALSE)+
geom_boxplot(width=0.1,outliers = F,color="black")+
theme_minimal()+
theme(
axis.text.x = element_text(size=rel(1),face="bold"),
plot.title = element_text(hjust = 0.5),
legend.position = "none"
)+
scale_y_log10(breaks=c(0.001,.01,.1,1),labels=c(0.1,1,10,100))+
labs(x="",y="% of total UMIs/bin \n (log10 scaled)")+
coord_flip()
if(!dir.exists(results_dir)){
system(glue("mkdir -p {results_dir}"))
}
saveRDS(visium, file.path(results_dir, "01_qc.RDS"))
outputPath = file.path(results_dir, "01_qc.RDS")
if(!dir.exists(results_dir)){
system(glue("mkdir -p {results_dir}"))
}
qsave(visium, file.path(results_dir, "01_qc.qs"))
outputPath = file.path(results_dir, "01_qc.qs")
# This set up the working directory to this file so all files can be found
library(rstudioapi)
setwd(fs::path_dir(getSourceEditorContext()$path))
stopifnot(R.version$major>= 4) # requires R4
if (compareVersion(R.version$minor,"3.1")<0) warning("We recommend >= R4.3.1")
stopifnot(compareVersion(as.character(BiocManager::version()), "3.16")>=0)
stopifnot(compareVersion(as.character(packageVersion("Seurat")), "5.1")>=0)
# This set up the working directory to this file so all files can be found
library(rstudioapi)
library(Seurat)
# devtools::install_github('satijalab/seurat-data')
library(SeuratData)
InstallData("stxBrain")
brain <- LoadData("stxBrain", type = "anterior1")
# BiocManager::install("qs2")
library(qs2)
qs_save(brain, file="visium/visium.qs")
test <- qs_read("visium/visium.qs")
BiocManager::install("SeuratWrappers")
BiocManager::install("Banksy")
BiocManager::install("spacexr")
BiocManager::install("dmcable/spacexr")
BiocManager::install("prabhakarlab/Banksy")
??RunBanksy
library(SeuratWrappers)
BiocManager::install("satijalab/seurat-wrappers")
ref <- qs::qsave("visium/allen_scRNAseq_ref_subset.qs")
ref <- qs::qread("visium/allen_scRNAseq_ref_subset.qs")
qs::qs_save("visium/allen_ref_subset.qs")
qs2::qs_save("visium/allen_ref_subset.qs")
qs2::qs_save(ref, file="visium/allen_ref_subset.qs")
getwd()
download.file("https://zenodo.org/records/15784846/files/allen_ref_subset.qs?download=1")
download.file("https://zenodo.org/records/15784846/files/allen_ref_subset.qs?download=1", "allen_ref_subset.qs")
library(curl)
# Optional: increase timeout (default is 60 sec, here set to 300 sec = 5 min)
handle <- new_handle(timeout = 300)
url <- "https://zenodo.org/records/15784846/files/allen_ref_subset.qs?download=1"
destfile <- "allen_ref_subset.qs"
curl_download(url, destfile, handle = handle)
renv::dependencies(path = ".")[["Package"]]
renv::snapshot()
?renv::snapshot
renv::dependencies(path = ".")[["Package"]]
# Only scan "scripts/" and "analysis.R"
renv::snapshot(
prompt = FALSE,
packages = NULL,  # auto-detect
packages = renv::dependencies(path = ".")[["Package"]]
)
# Only scan "scripts/" and "analysis.R"
renv::snapshot(
prompt = FALSE,
packages = renv::dependencies(path = ".")[["Package"]]
)