Skip to content

Commit 3caa8bc

Browse files
lpantanolpantano
committed
change name
1 parent eb21da7 commit 3caa8bc

File tree

2 files changed

+247
-1
lines changed

2 files changed

+247
-1
lines changed

‎.Rhistory‎

Lines changed: 246 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,246 @@
1+
library(Seurat)
2+
InstallData("stxBrain")
3+
??InstallData
4+
library(seurat-data)
5+
library(SeuratData)
6+
devtools::install_github('satijalab/seurat-data')
7+
# devtools::install_github('satijalab/seurat-data')
8+
library(SeuratData)
9+
InstallData("stxBrain")
10+
brain <- LoadData("stxBrain", type = "anterior1")
11+
brain
12+
library(qs)
13+
BiocManager::install("qdata")
14+
BiocManager::install("qs2")
15+
library(qs2)
16+
qd_save(brain, file="visium/visium.qs")
17+
qsave(brain, file="visium/visium.qs")
18+
test <- qread("visium/visium.qs")
19+
# This set up the working directory to this file so all files can be found
20+
library(rstudioapi)
21+
setwd(fs::path_dir(getSourceEditorContext()$path))
22+
stopifnot(R.version$major>= 4) # requires R4
23+
if (compareVersion(R.version$minor,"3.1")<0) warning("We recommend >= R4.3.1")
24+
stopifnot(compareVersion(as.character(BiocManager::version()), "3.16")>=0)
25+
stopifnot(compareVersion(as.character(packageVersion("Seurat")), "5.1")>=0)
26+
library(knitr)
27+
library(glue)
28+
library(qs2)
29+
library(dplyr)
30+
library(purrr)
31+
library(ggplot2)
32+
library(ggprism)
33+
library(grafify)
34+
library(ggpubr)
35+
library(gridExtra)
36+
library(scales)
37+
library(Seurat)
38+
import::from(magrittr,set_colnames,set_rownames,"%<>%")
39+
library(import)
40+
BiocManager::install("import")
41+
library(knitr)
42+
library(import)
43+
library(glue)
44+
library(qs2)
45+
library(dplyr)
46+
library(purrr)
47+
library(ggplot2)
48+
library(ggprism)
49+
library(grafify)
50+
library(ggpubr)
51+
library(gridExtra)
52+
library(scales)
53+
library(Seurat)
54+
import::from(magrittr, set_colnames, set_rownames, "%<>%")
55+
invisible(list2env(params,environment()))
56+
source(project_file)
57+
ggplot2::theme_set(theme_prism(base_size = 12))
58+
# https://grafify-vignettes.netlify.app/colour_palettes.html
59+
# NOTE change colors here if you wish
60+
scale_colour_discrete <- function(...)
61+
scale_colour_manual(...,
62+
values = as.vector(grafify:::graf_palettes[["kelly"]]))
63+
scale_fill_discrete <- function(...)
64+
scale_fill_manual(...,
65+
values = as.vector(grafify:::graf_palettes[["kelly"]]))
66+
opts_chunk[["set"]](
67+
cache = F,
68+
cache.lazy = FALSE,
69+
dev = c("png", "pdf"),
70+
error = TRUE,
71+
highlight = TRUE,
72+
message = FALSE,
73+
prompt = FALSE,
74+
tidy = FALSE,
75+
warning = FALSE,
76+
echo = T,
77+
fig.height = 4)
78+
# set seed for reproducibility
79+
set.seed(1234567890L)
80+
# 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.
81+
visium <- qread(visiumHD_obj)
82+
visium <- PercentageFeatureSet(visium, "^mt-", col.name = "percent_mito")
83+
visium <- PercentageFeatureSet(visium, "^Hb.*-", col.name = "percent_hb")
84+
metaD <- [email protected]
85+
metaD$log10GenesPerUMI <- log10(metaD$nFeature)/log10(metaD$nCount)
86+
colnames(metaD)%<>%gsub(pattern=glue("_{DefaultAssay(visium)}"),replacement="")
87+
summary_metaD <- apply(metaD[,-1],2,mean)
88+
metacol_label <- list("nFeature"="Genes","nCount"="UMI")
89+
refs <- list("nFeature"=100,"nCount"=100)
90+
dists_before <- imap(metacol_label,\(label,col)
91+
ggdensity(metaD,
92+
x = col,xscale="log10",add = "mean", rug = TRUE,
93+
alpha = 0.2,fill = "lightgray",
94+
xlab=glue("Number of {label} per bin(in log10 scale)"),
95+
ylab="Cell density",
96+
title=glue('Pre-QC {label}/Bin'))+
97+
geom_vline(xintercept = refs[[col]],color="darkred",cex=rel(1.3),linetype="dashed")+
98+
annotate("text",x=summary_metaD[col],y = Inf,
99+
label = glue("Mean \n = {round(summary_metaD[col],0)}"),
100+
vjust = 1,hjust=2)
101+
)
102+
dists_before[[1]] | dists_before[[2]]
103+
col <- "log10GenesPerUMI"
104+
ggdensity(metaD,x = col,add = "mean", rug = TRUE,
105+
alpha = 0.2,fill = "lightgray",
106+
xlab="complexity",ylab="Cell density",title=glue('Novelty score'))+
107+
geom_vline(xintercept = 0.8,color="darkred",cex=rel(1.3),linetype="dashed")+
108+
annotate("text",x=summary_metaD[col],y = Inf,
109+
label = glue("Mean = {round(summary_metaD[col],0)}"),
110+
vjust = 1,hjust=2)+
111+
theme(plot.title = element_text(hjust=0.5, face="bold"))
112+
knit_with_parameters("~/Science/nfcore/templates/spatial-reports/visium/01_quality_assessment/qc.Rmd")
113+
ggplot(metaD %>%
114+
select(orig.ident,starts_with("percent_")) %>%
115+
tidyr::gather(class,percent_unexpected,-orig.ident),
116+
aes_string(x = "orig.ident", y = "percent_unexpected")) +
117+
geom_violin(position=position_dodge(1),alpha=1, na.rm=TRUE,trim=FALSE)+
118+
ggbeeswarm::geom_quasirandom(na.rm=TRUE,dodge.width=0.5,
119+
method='quasirandom',alpha=0.01)+
120+
geom_boxplot(width=0.1,outliers = F)+
121+
geom_hline(yintercept=20)+
122+
facet_grid(~class,scales = "free")+
123+
theme(
124+
axis.text.x = element_text(size=rel(1),face="bold"),
125+
plot.title = element_text(hjust = 0.5),
126+
strip.text.x = element_text(size = rel(1.5), colour = "black"),
127+
legend.position = "none"
128+
)+
129+
scale_y_log10(breaks=c(1,5,10,20,100))+
130+
# ylim(c(0,100))+
131+
labs(x="",y="% of contamination genes")
132+
features2check <- c(glue('nCount_{DefaultAssay(visium)}'),
133+
glue('nFeature_{DefaultAssay(visium)}'),
134+
"percent_mito","percent_hb")
135+
for(f in features2check){
136+
cat("### ", f, "\n\n")
137+
p1 <- SpatialFeaturePlot(visium,
138+
feature = f,
139+
pt.size.factor = 8)
140+
print(p1)
141+
cat("\n\n")
142+
}
143+
for(f in features2check){
144+
cat("### ", f, "\n\n")
145+
p1 <- SpatialFeaturePlot(visium,
146+
feature = f,
147+
pt.size.factor = 4)
148+
print(p1)
149+
cat("\n\n")
150+
}
151+
GeneVar <- glue('nFeature_{DefaultAssay(visium)}')
152+
UMIVar <- glue('nCount_{DefaultAssay(visium)}')
153+
cutoffs <- list("nFeature"=100,"nCount"=100,"hb"=20,"mito"=20)
154+
Qced <- [email protected][,GeneVar] > cutoffs$nFeature &
155+
[email protected][,UMIVar] > cutoffs$nCount &
156+
visium$percent_hb < cutoffs$hb &
157+
visium$percent_mito < cutoffs$mito
158+
visium <- visium[,Qced]
159+
# Filter Mitocondrial
160+
visium <- visium[!grepl("^mt-", rownames(visium)), ]
161+
# Filter Hemoglobin gene (optional if that is a problem on your data)
162+
visium <- visium[!grepl("^Hb.*-", rownames(visium)), ]
163+
C <- GetAssayData(visium, slot = "counts")
164+
C@x <- C@x / rep.int(colSums(C), diff(C@p))
165+
most_expressed <- order(Matrix::rowSums(C), decreasing = T)[20:1]
166+
exprD <- as.data.frame(t(C[most_expressed, ])) %>%
167+
tibble::rownames_to_column("bin") %>%
168+
tidyr::gather(gene,expr,-bin)
169+
ggplot(exprD,aes(x=gene,y=expr,color=gene,fill=gene))+
170+
geom_violin(position=position_dodge(1),alpha=0.5,
171+
na.rm=TRUE,trim=FALSE)+
172+
geom_boxplot(width=0.1,outliers = F,color="black")+
173+
theme_minimal()+
174+
theme(
175+
axis.text.x = element_text(size=rel(1),face="bold"),
176+
plot.title = element_text(hjust = 0.5),
177+
legend.position = "none"
178+
)+
179+
scale_y_log10(breaks=c(0.001,.01,.1,1),labels=c(0.1,1,10,100))+
180+
labs(x="",y="% of total UMIs/bin \n (log10 scaled)")+
181+
coord_flip()
182+
if(!dir.exists(results_dir)){
183+
system(glue("mkdir -p {results_dir}"))
184+
}
185+
saveRDS(visium, file.path(results_dir, "01_qc.RDS"))
186+
outputPath = file.path(results_dir, "01_qc.RDS")
187+
if(!dir.exists(results_dir)){
188+
system(glue("mkdir -p {results_dir}"))
189+
}
190+
qsave(visium, file.path(results_dir, "01_qc.qs"))
191+
outputPath = file.path(results_dir, "01_qc.qs")
192+
# This set up the working directory to this file so all files can be found
193+
library(rstudioapi)
194+
setwd(fs::path_dir(getSourceEditorContext()$path))
195+
stopifnot(R.version$major>= 4) # requires R4
196+
if (compareVersion(R.version$minor,"3.1")<0) warning("We recommend >= R4.3.1")
197+
stopifnot(compareVersion(as.character(BiocManager::version()), "3.16")>=0)
198+
stopifnot(compareVersion(as.character(packageVersion("Seurat")), "5.1")>=0)
199+
# This set up the working directory to this file so all files can be found
200+
library(rstudioapi)
201+
library(Seurat)
202+
# devtools::install_github('satijalab/seurat-data')
203+
library(SeuratData)
204+
InstallData("stxBrain")
205+
brain <- LoadData("stxBrain", type = "anterior1")
206+
# BiocManager::install("qs2")
207+
library(qs2)
208+
qs_save(brain, file="visium/visium.qs")
209+
test <- qs_read("visium/visium.qs")
210+
BiocManager::install("SeuratWrappers")
211+
BiocManager::install("Banksy")
212+
BiocManager::install("spacexr")
213+
BiocManager::install("dmcable/spacexr")
214+
BiocManager::install("prabhakarlab/Banksy")
215+
??RunBanksy
216+
library(SeuratWrappers)
217+
BiocManager::install("satijalab/seurat-wrappers")
218+
ref <- qs::qsave("visium/allen_scRNAseq_ref_subset.qs")
219+
ref <- qs::qread("visium/allen_scRNAseq_ref_subset.qs")
220+
qs::qs_save("visium/allen_ref_subset.qs")
221+
qs2::qs_save("visium/allen_ref_subset.qs")
222+
qs2::qs_save(ref, file="visium/allen_ref_subset.qs")
223+
getwd()
224+
download.file("https://zenodo.org/records/15784846/files/allen_ref_subset.qs?download=1")
225+
download.file("https://zenodo.org/records/15784846/files/allen_ref_subset.qs?download=1", "allen_ref_subset.qs")
226+
library(curl)
227+
# Optional: increase timeout (default is 60 sec, here set to 300 sec = 5 min)
228+
handle <- new_handle(timeout = 300)
229+
url <- "https://zenodo.org/records/15784846/files/allen_ref_subset.qs?download=1"
230+
destfile <- "allen_ref_subset.qs"
231+
curl_download(url, destfile, handle = handle)
232+
renv::dependencies(path = ".")[["Package"]]
233+
renv::snapshot()
234+
?renv::snapshot
235+
renv::dependencies(path = ".")[["Package"]]
236+
# Only scan "scripts/" and "analysis.R"
237+
renv::snapshot(
238+
prompt = FALSE,
239+
packages = NULL, # auto-detect
240+
packages = renv::dependencies(path = ".")[["Package"]]
241+
)
242+
# Only scan "scripts/" and "analysis.R"
243+
renv::snapshot(
244+
prompt = FALSE,
245+
packages = renv::dependencies(path = ".")[["Package"]]
246+
)

‎.github/workflows/run_visium_qc.yaml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ on:
77
branches: [main]
88

99
jobs:
10-
r-restore:
10+
r-visum-qc:
1111
runs-on: ubuntu-noble
1212

1313
env:

0 commit comments

Comments
 (0)