Spatial Transcriptomics Analysis Pipeline for 10x Visium Platform Using Snakemake

public public 1yr ago 0 bookmarks
View Workflow
Report
spatial-workflow
spatial-workflow
View Workflow

Introduction

This is a spatial transcriptomics analysis pipeline for Visium 10x platform. The pipeline is built in Snakemake and can be run on different platforms and high performance computing (HPC) systems.

Quick Start Example

The workflow excepts Visum 10x samples under data folder in this format:

"data/{sample}/outs/filtered_feature_bc_matrix.h5"

This will register the directory name as sample name for later processing.

You can start the pipeline by calling, 

snakemake -j 20

which will create a 20 threads job.

Code Snippets

20
21
shell:
    "workflow/scripts/spatial-rds.R {wildcards.sample}"
SnakeMake From line 20 of rules/spatial.smk 
29
30
shell:
    "workflow/scripts/clusteringTree.R {wildcards.sample}"
SnakeMake From line 29 of rules/spatial.smk 
38
39
40
41
42
shell:
    """
    #convert {input} -colorspace HCL -channel R -evaluate set 67% +channel -colorspace sRGB {output}
    workflow/scripts/saturation 0.4 {input} {output}
    """
SnakeMake From line 38 of rules/spatial.smk 
49
50
51
52
53
54
shell:
    """
    #convert {input} -colorspace HCL -channel R -evaluate set 67% +channel -colorspace sRGB {output}
    workflow/scripts/saturation 0.4 {input} {output}

    """
SnakeMake From line 49 of rules/spatial.smk 
64
65
66
67
shell:
    """
    workflow/scripts/spatial-qc.R {wildcards.sample}
    """
SnakeMake From line 64 of rules/spatial.smk 
75
76
77
78
shell:
    """
    workflow/scripts/spatial-umap.R {wildcards.sample} {wildcards.res}
    """
SnakeMake From line 75 of rules/spatial.smk 
86
87
88
89
shell:
    """
    workflow/scripts/spatial-metrics.R {wildcards.sample} {wildcards.res}
    """
SnakeMake From line 86 of rules/spatial.smk 
 97
 98
 99
100
shell:
    """
    workflow/scripts/spatial-cluster-markers.R {wildcards.sample} {wildcards.res}
    """
SnakeMake From line 97 of rules/spatial.smk 
109
110
111
112
shell:
    """
    workflow/scripts/spatial-cluster-markerplots.R {wildcards.sample} {wildcards.res}
    """
SnakeMake From line 109 of rules/spatial.smk 
121
122
123
124
125
shell:
    """
    mkdir -p {output[1]}
    workflow/scripts/spatial-spatial-markers.R {wildcards.sample}
    """
SnakeMake From line 121 of rules/spatial.smk 
133
134
135
136
137
shell:
    """
    mkdir -p {output}
    workflow/scripts/spatial-selected-markers.R {wildcards.sample}
    """
SnakeMake From line 133 of rules/spatial.smk 
145
146
147
148
shell:
    """
    workflow/scripts/spatial-sc-cluster.R {wildcards.datafile}
    """
SnakeMake From line 145 of rules/spatial.smk 
160
161
162
163
shell:
    """
    workflow/scripts/spatial-spotlight.R {wildcards.sample} {wildcards.datafile}
    """
SnakeMake From line 160 of rules/spatial.smk 
176
177
178
179
shell:
    """
    workflow/scripts/spatial-rctd.R {wildcards.sample} {wildcards.datafile}
    """
SnakeMake From line 176 of rules/spatial.smk 
188
189
190
191
shell:
    """
    workflow/scripts/spatial-rctd-pdf.R {wildcards.sample} {wildcards.datafile}
    """
SnakeMake From line 188 of rules/spatial.smk 
201
202
203
204
shell:
    """
    workflow/scripts/spatial-spotlight-pdf.R {wildcards.sample} {wildcards.datafile}
    """
SnakeMake From line 201 of rules/spatial.smk 
213
214
215
216
shell:
    """
    workflow/scripts/spatial-seurat-decon.R {wildcards.sample} {wildcards.datafile}
    """
SnakeMake From line 213 of rules/spatial.smk 
224
225
226
227
shell:
    """
    workflow/scripts/spatial-giotto.R {wildcards.sample} {wildcards.datafile}
    """
SnakeMake From line 224 of rules/spatial.smk 
235
236
237
238
shell:
    """
    workflow/scripts/spatial-convert.R {wildcards.datafile}
    """
SnakeMake From line 235 of rules/spatial.smk 
251
252
253
254
shell:
    """
    workflow/scripts/spatial-tangram.py data/{wildcards.sample}/outs {input[0]} {output}
    """
SnakeMake From line 251 of rules/spatial.smk 
262
263
264
265
shell:
    """
    workflow/scripts/spatial-tangram-pdf.R {wildcards.sample} {wildcards.datafile}
    """
SnakeMake From line 262 of rules/spatial.smk 
277
278
279
280
shell:
    """
    workflow/scripts/spatial-tangram-gene.py data/{wildcards.sample}/outs {input[0]} {output}
    """
SnakeMake From line 277 of rules/spatial.smk 
288
289
290
291
shell:
    """
    workflow/scripts/spatial-giotto-pdf.R {wildcards.sample} {wildcards.datafile}
    """
SnakeMake From line 288 of rules/spatial.smk 
301
302
303
304
shell:
    """
    workflow/scripts/spatial-gbmtest.R {wildcards.sample} {wildcards.modelfile}
    """
SnakeMake From line 301 of rules/spatial.smk 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
library(Seurat)
require(tidyverse)
require(clustree)
arguments=commandArgs(TRUE)

sampleID=arguments[1]


Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))

Spatial_Data <- FindClusters(Spatial_Data, verbose = FALSE,resolution = seq(0.1,2.5,0.1))

clustree(Spatial_Data) -> p1

ggsave(paste0("results/",sampleID,"/clusteringTree/clusteringTree-",sampleID,".pdf"),p1,width=8,height=15)
R tidyverse Seurat clustree From line 3 of scripts/clusteringTree.R 
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
library(Seurat)
require(tidyverse)
require(patchwork)
source("workflow/scripts/spatial-functions.R")

arguments=commandArgs(TRUE)

sampleID=arguments[1]
res=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))

function_image_fixer(Spatial_Data) -> Spatial_Data

SCT_=paste0("SCT_snn_res.",res)


Idents(object = Spatial_Data) <- Spatial_Data@meta.data[[SCT_]]


Positive_Features=openxlsx::read.xlsx(paste0(sampleID,"/resolution-",res,"/",sampleID,".positive-markers-forAllClusters.xlsx")) %>% select(cluster,gene) 


for(d in (Positive_Features %>% distinct(cluster) %>% pull())) {

    dir.create(paste0(sampleID,"/resolution-",res,"/markers/","cluster",d,"/"),recursive=TRUE)
}

options(warn=-1)
suppressMessages(for (i in 1:nrow(Positive_Features)) {

    gene=Positive_Features[i,]$gene
    cluster=Positive_Features[i,]$cluster

p1 <- FeaturePlot(Spatial_Data,assay = "SCT", reduction = "umap", features=gene) + scale_colour_gradientn(colours = rev(RColorBrewer::brewer.pal(n = 11, name = "RdYlGn")))
p2 <- SpatialFeaturePlot(Spatial_Data,assay = "SCT", features=gene,images=paste0("image"),alpha=c(0.1,1),pt.size.factor=1.2) + scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(n = 11, name = "RdYlGn")))
p3 <- DotPlot(Spatial_Data,assay = "SCT", features=gene)
p4 <- VlnPlot(Spatial_Data,assay = "SCT",features=gene)

suppressWarnings(wrap_plots(p1,p2,p3,p4,ncol=2) -> wp)

ggsave(paste0("results/",sampleID,"/resolution-",res,"/markers/","cluster",cluster,"/",gene,".pdf"),wp,height=10,width=10)

})
R tidyverse Seurat patcHwork From line 2 of scripts/spatial-cluster-markerplots.R 
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
library(Seurat)
require(tidyverse)


arguments=commandArgs(TRUE)

sampleID=arguments[1]
res=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))

SCT_=paste0("SCT_snn_res.",res)


Idents(object = Spatial_Data) <- Spatial_Data@meta.data[[SCT_]]


all_markers=FindAllMarkers(Spatial_Data)

openxlsx::write.xlsx(all_markers,file=paste0("results/",sampleID,"/resolution-",res,"/",sampleID,".all-markers-forAllClusters",".xlsx"))

openxlsx::write.xlsx(all_markers %>% filter(avg_log2FC > 0),file=paste0("results/",sampleID,"/resolution-",res,"/",sampleID,".positive-markers-forAllClusters",".xlsx"))
R tidyverse Seurat From line 2 of scripts/spatial-cluster-markers.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
require(Seurat)
require(SeuratDisk)
require(dplyr)

source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

scrnaID=arguments[1]

scrna_data=readRDS(paste0("scrna/",scrnaID,".rds"))

UpdateSeuratObject(scrna_data) -> scrna_data
DefaultAssay(scrna_data) <- "RNA"
scrna_data[["SCT"]] <- NULL

DietSeurat(scrna_data) -> scrna_data

scrna_data@meta.data %>% dplyr::mutate(dplyr::across(where(is.factor), as.character)) -> scrna_data@meta.data

#AddMetaData(scrna_data,make.names(as.character(scrna_data$seurat_clusters)),"seurat_clusters_tangram") -> scrna_data


SaveH5Seurat(scrna_data,paste0("scrna/",scrnaID,".h5Seurat"))
SeuratDisk::Convert(paste0("scrna/",scrnaID,".h5Seurat"), dest = "h5ad")
R dplyr Seurat SeuratDisk From line 3 of scripts/spatial-convert.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
require(Seurat)
require(tidyverse)
#require(caret)


source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

sampleID=arguments[1]
modelID=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))
model=readRDS(paste0("models/",modelID,".rds"))





function_add_missing=function(df,vector) {

  for (i in vector) {
    df <- df %>% dplyr::mutate("{i}":=0)

  }
  return(df)
}

function_image_fixer(Spatial_Data) -> Spatial_Data


function_gbm_test=function(Seurat_Object) {
  GetAssayData(Seurat_Object,slot = "scale.data") %>% as.matrix() %>% t() %>% as.data.frame() %>% rownames_to_column("cell") %>% janitor::clean_names()-> testdf2


  testdf2 <-function_add_missing(testdf2,setdiff(c(model$coefnames,"cell"),colnames(testdf2)))

  predict(model,testdf2,type="prob") %>% janitor::clean_names() %>% bind_cols(testdf2 %>% select(barcodes=cell)) -> gbm_decon_df


cell_types_all=gbm_decon_df %>% select(-barcodes) %>% colnames()

Seurat_Object@meta.data <- Seurat_Object@meta.data %>%
  tibble::rownames_to_column("barcodes") %>%
  dplyr::left_join(gbm_decon_df, by = "barcodes") %>%
  tibble::column_to_rownames("barcodes")

wp=Seurat::SpatialFeaturePlot(
  object = Seurat_Object,
  features = cell_types_all,
  alpha = c(0.1, 1),pt.size.factor = 1,ncol=2,images=paste0("image")) & scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(n = 11, name = "RdYlGn")),limits=c(0,1)) & theme(legend.title = element_text(size=4),legend.key.size = unit(0.2,"cm"),legend.text = element_text(size=3),legend.margin=margin(t = 0,b = 0, unit='cm'),plot.margin = margin(0.1, 0.1, 0.1, 0.1, "cm"))


ggsave(paste0("results/",sampleID,"/deconvolution/gbm/",sampleID,"-",modelID,"-gbm.pdf"),wp,height=20,width=7)

return(Seurat_Object)
}

function_gbm_test(Spatial_Data) -> Spatial_Data
R tidyverse Seurat From line 3 of scripts/spatial-gbmtest.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
require(Seurat)
require(tidyverse)
require(viridis)



source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

sampleID=arguments[1]
scrnaID=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))
dwls_data=readRDS(paste0("DWLS_assay/",scrnaID,"/",sampleID,".rds"))


function_image_fixer(Spatial_Data) -> Spatial_Data



Spatial_Data[["DWLS"]] <- dwls_data


DefaultAssay(Spatial_Data) <- "DWLS"

cell_types_all=rownames(Spatial_Data)

wp=seurat_plotting()

ggsave(paste0("results/",sampleID,"/deconvolution/dwls/",sampleID,"-",scrnaID,"-dwls.pdf"),wp,height=18,width=8)
R tidyverse Seurat viridis From line 3 of scripts/spatial-giotto-pdf.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
require(Seurat)
require(tidyverse)


require(Giotto)


arguments=commandArgs(TRUE)

sampleID=arguments[1]
scrnaID=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))
scrna_data=readRDS(paste0("scrna/",scrnaID,".rds"))



python_path=system("which python",intern=T)

instrs = createGiottoInstructions(python_path = python_path)

st_data <- createGiottoObject(
    raw_exprs = Spatial_Data@assays$Spatial@counts,
    instructions = instrs
)
# st_data <- filterGiotto(gobject = st_data,
#                              expression_threshold = 1,
#                              gene_det_in_min_cells = 50,
#                              min_det_genes_per_cell = 1000,
#                              expression_values = c('raw'),
#                              verbose = T)
st_data <- normalizeGiotto(gobject = st_data)
st_data <- calculateHVG(gobject = st_data)
gene_metadata = fDataDT(st_data)
featgenes = gene_metadata[hvg == 'yes']$gene_ID
st_data <- runPCA(gobject = st_data, genes_to_use = featgenes, scale_unit = F)
signPCA(st_data, genes_to_use = featgenes, scale_unit = F)
st_data <- runUMAP(st_data, dimensions_to_use = 1:10)
st_data <- createNearestNetwork(gobject = st_data, dimensions_to_use = 1:10, k = 15)
st_data <- doLeidenCluster(gobject = st_data, resolution = 0.4, n_iterations = 1000)
sc_data <- createGiottoObject(
    raw_exprs = scrna_data@assays$RNA@counts,
    instructions = instrs
)
sc_data <- normalizeGiotto(gobject = sc_data)
sc_data <- calculateHVG(gobject = sc_data)
gene_metadata = fDataDT(sc_data)
featgenes = gene_metadata[hvg == 'yes']$gene_ID
sc_data <- runPCA(gobject = sc_data, genes_to_use = featgenes, scale_unit = F)
signPCA(sc_data, genes_to_use = featgenes, scale_unit = F)


sc_data@cell_metadata$leiden_clus <- as.character(scrna_data@meta.data[,"seurat_clusters"])
scran_markers_subclusters = findMarkers_one_vs_all(gobject = sc_data,
                                                   method = 'scran',
                                                   expression_values = 'normalized',
                                                   cluster_column = 'leiden_clus')
Sig_scran <- unique(scran_markers_subclusters$genes[which(scran_markers_subclusters$ranking <= 100)])
norm_exp<-2^(sc_data@norm_expr)-1
id<-sc_data@cell_metadata$leiden_clus
ExprSubset<-norm_exp[Sig_scran,]
Sig_exp<-NULL
for (i in unique(id)){
  Sig_exp<-cbind(Sig_exp,(apply(ExprSubset,1,function(y) mean(y[which(id==i)]))))
}
colnames(Sig_exp)<-unique(id)
st_data <- runDWLSDeconv(st_data,sign_matrix = Sig_exp, n_cell = 20)



DWLS <- CreateAssayObject(st_data@spatial_enrichment$DWLS %>% column_to_rownames("cell_ID") %>% t())

saveRDS(DWLS,file = paste0("DWLS_assay/",scrnaID,"/",sampleID,".rds"))
R tidyverse Seurat From line 3 of scripts/spatial-giotto.R 
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
library(Seurat)
require(tidyverse)


arguments=commandArgs(TRUE)

sampleID=arguments[1]
res=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))





SCT_=paste0("SCT_snn_res.",res)

metrics=table(Spatial_Data@meta.data[[SCT_]], Spatial_Data@meta.data$orig.ident)

openxlsx::write.xlsx(metrics %>% as.data.frame() %>% select(sampleID=1),file=paste0("results/",sampleID,"/resolution-",res,"/",sampleID,".number-of-cells-per-cluster",".xlsx"))
R tidyverse Seurat From line 2 of scripts/spatial-metrics.R 
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
library(Seurat)
require(tidyverse)
require(patchwork)
source("workflow/scripts/spatial-functions.R")

arguments=commandArgs(TRUE)

sampleID=arguments[1]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))

Spatial_Data[["percent.mt"]] <- PercentageFeatureSet(Spatial_Data, pattern = "^MT-")
Spatial_Data[["percent.rp"]] <- PercentageFeatureSet(Spatial_Data, pattern = "^RP[SL]")



function_image_fixer(Spatial_Data) -> Spatial_Data


plot1 <- VlnPlot(Spatial_Data, features = "nCount_Spatial", pt.size = 0.1,assay="Spatial",group.by="orig.ident") + NoLegend()
plot2 <- SpatialFeaturePlot(Spatial_Data, features = "nCount_Spatial",images="image") + theme(legend.position = "right")

wrap_plots(plot1, plot2,ncol=2)

ggsave(filename=paste0("results/",sampleID,"/technicals/",sampleID,".n_counts.pdf"),width=13,height=7)


Spatial_Data <- NormalizeData(Spatial_Data, verbose = FALSE, assay = "Spatial")

Spatial_Data <- GroupCorrelation(Spatial_Data, group.assay = "Spatial", assay = "Spatial", slot = "data", do.plot = FALSE)
Spatial_Data <- GroupCorrelation(Spatial_Data, group.assay = "Spatial", assay = "SCT", slot = "scale.data", do.plot = FALSE)

p1 <- GroupCorrelationPlot(Spatial_Data, assay = "Spatial", cor = "nCount_Spatial_cor") + ggtitle("Log Normalization") + theme(plot.title = element_text(hjust = 0.5))
p2 <- GroupCorrelationPlot(Spatial_Data, assay = "SCT", cor = "nCount_Spatial_cor") + ggtitle("SCTransform Normalization") + theme(plot.title = element_text(hjust = 0.5))


wrap_plots(p1, p2,ncol=2)

ggsave(filename=paste0("results/",sampleID,"/technicals/",sampleID,".normalization.pdf"),width=13,height=7)
R tidyverse Seurat patcHwork From line 2 of scripts/spatial-qc.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
library(Seurat)
require(tidyverse)
require(viridis)



source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

sampleID=arguments[1]
scrnaID=arguments[2]

Spatial_Data=readRDS(paste0("rds_rctd/",scrnaID,"/",sampleID,".rds"))
scrna_data=readRDS(paste0("scrna/",scrnaID,".rds"))


cell_types_all=Idents(scrna_data) %>% unique() %>% as.character() %>% make.names()

function_image_fixer(Spatial_Data) -> Spatial_Data


wp=seurat_plotting()

ggsave(paste0("results/",sampleID,"/deconvolution/rctd/",sampleID,"-",scrnaID,"-rctd.pdf"),wp,height=18,width=6)
R tidyverse Seurat viridis From line 3 of scripts/spatial-rctd-pdf.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
require(Seurat)
require(tidyverse)
require(spacexr)


source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

sampleID=arguments[1]
scrnaID=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))
scrna_data=readRDS(paste0("scrna/",scrnaID,".rds"))


cell_types_all=Idents(scrna_data) %>% unique() %>% as.character() %>% make.names()

function_image_fixer(Spatial_Data) -> Spatial_Data


counts <- scrna_data@assays$RNA@counts
colnames(counts) <- colnames(scrna_data)
meta_data <- data.frame(scrna_data@meta.data)
cell_types <- meta_data$seurat_clusters %>% make.names()
names(cell_types) <- rownames(scrna_data@meta.data)
cell_types <- as.factor(cell_types)
nUMI_df <- data.frame(colSums(scrna_data@assays$RNA@counts))
nUMI <- nUMI_df$colSums.scrna_data.assays.RNA
names(nUMI) <- rownames(nUMI_df)



reference <- Reference(counts, cell_types, nUMI)




coords <- data.frame(colnames(Spatial_Data))
colnames(coords) <- 'barcodes'
coords$xcoord <- seq_along(colnames(Spatial_Data))
coords$ycoord <- seq_along(colnames(Spatial_Data))
counts <- data.frame(Spatial_Data@assays$Spatial@counts) # load in counts matrix
colnames(counts) <- colnames(Spatial_Data)
coords <- data.frame(colnames(Spatial_Data))
colnames(coords) <- 'barcodes'
coords$xcoord <- seq_along(colnames(Spatial_Data))
coords$ycoord <- seq_along(colnames(Spatial_Data))
rownames(coords) <- coords$barcodes; coords$barcodes <- NULL # Move barcodes to rownames
nUMI <- colSums(counts) # In this case, total counts per pixel is nUMI

puck<- SpatialRNA(coords = coords,counts = counts,nUMI = nUMI)

myRCTD <- create.RCTD(puck,reference,max_cores = 5,UMI_min = 20)
myRCTD <- run.RCTD(myRCTD,doublet_mode = "full")


Spatial_Data@meta.data <- Spatial_Data@meta.data %>%
  tibble::rownames_to_column("barcodes") %>%
  dplyr::left_join(myRCTD@results$weights %>% as.data.frame() %>% rownames_to_column("barcodes"), by = "barcodes") %>%
  tibble::column_to_rownames("barcodes")


saveRDS(Spatial_Data,file = paste0("rds_rctd/",scrnaID,"/",sampleID,".rds"))
R tidyverse Seurat From line 3 of scripts/spatial-rctd.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
library(Seurat)
require(tidyverse)

source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

sampleID=arguments[1]
slice=make.names(arguments[1])
datadir=c(paste0("data/",sampleID,"/outs/"))

Spatial_Data=function_analyse_spatial(datadir = datadir,slice = slice)
saveRDS(Spatial_Data,file = paste0("rds/",sampleID,".rds"))
R tidyverse Seurat From line 3 of scripts/spatial-rds.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
library(Seurat)
require(tidyverse)
library(SPOTlight)


source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

scrnaID=arguments[1]

scrna_data=readRDS(paste0("scrna/",scrnaID,".rds"))





cluster_markers_all <- Seurat::FindAllMarkers(object = scrna_data, 
                                              assay = "SCT",
                                              slot = "data",
                                              verbose = TRUE, 
                                              only.pos = TRUE)


openxlsx::write.xlsx(cluster_markers_all,file=paste0("scrna/",scrnaID,".cluster_markers.xlsx"))
R tidyverse Seurat Spotlight From line 3 of scripts/spatial-sc-cluster.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
require(Seurat)
require(tidyverse)
require(viridis)



params=list(k.anchor=23,k.score=5,k.filter=100,n.trees=100)


source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

sampleID=arguments[1]
scrnaID=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))
scrna_data=readRDS(paste0("scrna/",scrnaID,".rds"))


function_image_fixer(Spatial_Data) -> Spatial_Data


function_decon_seurat = function(reference,query,anc){

anchors <- FindTransferAnchors(reference = reference, query = query, normalization.method = "SCT",
k.anchor = anc,
k.score = params$k.score,
k.filter=params$k.filter,
n.trees = params$n.trees)

predictions.assay <- TransferData(anchorset = anchors, refdata = reference$seurat_clusters, prediction.assay = TRUE)
query[["predictions"]] <- predictions.assay


return(query)
}



for (i in c(params$k.anchor,20,15,10,5)) {

try({
  function_decon_seurat(reference=scrna_data,query=Spatial_Data,anc=i) -> Spatial_Data
  break
  }
  )


}

DefaultAssay(Spatial_Data) <- "predictions"

cell_types_all=Idents(scrna_data) %>% unique() %>% as.character()

wp=seurat_plotting()

ggsave(paste0("results/",sampleID,"/deconvolution/seurat/",sampleID,"-",scrnaID,"-seurat.pdf"),wp,height=18,width=8)
R tidyverse Seurat viridis From line 3 of scripts/spatial-seurat-decon.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
library(Seurat)
require(tidyverse)
library(SPOTlight)
require(viridis)


source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

sampleID=arguments[1]
scrnaID=arguments[2]

Spatial_Data=readRDS(paste0("rds_decon/",scrnaID,"/",sampleID,".rds"))
scrna_data=readRDS(paste0("scrna/",scrnaID,".rds"))


cell_types_all=Idents(scrna_data) %>% unique() %>% as.character() %>% make.names()

function_image_fixer(Spatial_Data) -> Spatial_Data


wp=seurat_plotting()

ggsave(paste0("results/",sampleID,"/deconvolution/spotlight/",sampleID,"-",scrnaID,"-spotlight.pdf"),wp,height=18,width=6)
R tidyverse Seurat viridis Spotlight From line 3 of scripts/spatial-spotlight-pdf.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
library(Seurat)
require(tidyverse)
library(SPOTlight)


source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

sampleID=arguments[1]
scrnaID=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))
scrna_data=readRDS(paste0("scrna/",scrnaID,".rds"))



openxlsx::read.xlsx(paste0("scrna/",scrnaID,".cluster_markers.xlsx")) -> cluster_markers_all

function_spotlight=function(Seurat_Object,scrna_rds) {
set.seed(123)

spotlight_ls <- spotlight_deconvolution(
  se_sc = scrna_rds,
  counts_spatial = Seurat_Object@assays$Spatial@counts,
  clust_vr = "seurat_clusters", # Variable in sc_seu containing the cell-type annotation
  cluster_markers = cluster_markers_all, # Dataframe with the marker genes
  cl_n = 100, # number of cells per cell type to use
  hvg = 3000, # Number of HVG to use
  ntop = NULL, # How many of the marker genes to use (by default all)
  transf = "uv", # Perform unit-variance scaling per cell and spot prior to factorzation and NLS
  method = "nsNMF", # Factorization method
  min_cont = 0 # Remove those cells contributing to a spot below a certain threshold 
  )
nmf_mod <- spotlight_ls[[1]]
decon_mtrx <- spotlight_ls[[2]]

decon_mtrx_sub <- decon_mtrx[, colnames(decon_mtrx) != "res_ss"]
decon_mtrx_sub[decon_mtrx_sub < 0.08] <- 0
decon_mtrx <- cbind(decon_mtrx_sub, "res_ss" = decon_mtrx[, "res_ss"])
rownames(decon_mtrx) <- colnames(Seurat_Object)

decon_df <- decon_mtrx %>%
  data.frame() %>%
  tibble::rownames_to_column("barcodes")

Seurat_Object@meta.data <- Seurat_Object@meta.data %>%
  tibble::rownames_to_column("barcodes") %>%
  dplyr::left_join(decon_df, by = "barcodes") %>%
  tibble::column_to_rownames("barcodes")

return(Seurat_Object)
}

function_spotlight(Spatial_Data,scrna_data) -> Spatial_Data

saveRDS(Spatial_Data,file = paste0("rds_decon/",scrnaID,"/",sampleID,".rds"))
R tidyverse Seurat Spotlight From line 3 of scripts/spatial-spotlight.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
import scanpy as sc
import tangram as tg
import pandas as pd
import sys
import matplotlib.pyplot as plt
import torch



print(sys.argv[1])
adata_st = sc.read_visium(path=sys.argv[1])
adata_sc= sc.read(filename=sys.argv[2])

adata_sc.X=adata_sc.raw.X.copy()


tg.pp_adatas(adata_sc,adata_st,genes=None)


if torch.cuda.is_available():
       ad_map = tg.map_cells_to_space(
                   adata_sc, 
                   adata_st,device="cuda:0")
else:
       ad_map = tg.map_cells_to_space(
                   adata_sc, 
                   adata_st)


ad_ge = tg.project_genes(ad_map, adata_sc)

genes=["LYVE1","F13A1","FOLR2","SELENOP","APOE","SLC40A1","C1QB","DAB2","PDK4","SPP1","ACP5","CD9","FCER1A","CD1C","CLEC10A","HSPA6","DNAJB1",
       "HSPA1B","S100A8","S100A9","S100A12","EREG","G0S2","FCN1","CCL20","IL1B","IL23A","CXCL10","CXCL9","GBP1","CDC1C","CCL3L1",
       "CCL3","CCL4L2","MT1X","MT1E","CTSL","RGS1","FOS","DNASE1L3","MMP9","LYZ","AREG","VCAN","HSPA1A","MARCO","COLEC12"]

genes=list(set(list(ad_ge.var["features"].to_numpy())) & set(genes))
genes=list(map(lambda x: x.lower(),genes))


expid=list(adata_st.uns['spatial'])[0]

scale_fac=adata_st.uns['spatial'][expid]['scalefactors']['tissue_hires_scalef']

pl=tg.plot_genes_sc(genes, adata_measured=adata_st, adata_predicted=ad_ge, perc=0.02,scale_factor=scale_fac,spot_size=40,return_figure=True)

plt.rcParams['figure.figsize'] = [15, 15]
plt.rcParams['figure.dpi'] = 150


pl.savefig(sys.argv[3], format="pdf", bbox_inches="tight")
Python Pandas matplotlib torch Tangram From line 3 of scripts/spatial-tangram-gene.py 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
require(Seurat)
require(tidyverse)
require(viridis)


source("workflow/scripts/spatial-functions.R")
arguments=commandArgs(TRUE)

sampleID=arguments[1]
scrnaID=arguments[2]


Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))
tangram_csv=read.csv(paste0("tangram/",scrnaID,"/",sampleID,".csv"))



function_image_fixer(Spatial_Data) -> Spatial_Data

Spatial_Data[["predictions"]] <- CreateAssayObject(tangram_csv  %>% column_to_rownames("X") %>% t())


DefaultAssay(Spatial_Data) <- "predictions"


cell_types_all=tangram_csv  %>% column_to_rownames("X") %>% colnames()


wp=seurat_plotting()


ggsave(paste0("results/",sampleID,"/deconvolution/tangram/",sampleID,"-",scrnaID,"-tangram.pdf"),wp,height=18,width=8)
R tidyverse Seurat viridis From line 3 of scripts/spatial-tangram-pdf.R 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
import scanpy as sc
import tangram as tg
import pandas as pd
import sys
import torch

print(sys.argv[1])
adata_st = sc.read_visium(path=sys.argv[1])
adata_sc= sc.read(filename=sys.argv[2])

adata_sc.X=adata_sc.raw.X.copy()


tg.pp_adatas(adata_sc,adata_st,genes=None)


if torch.cuda.is_available():
    ad_map = tg.map_cells_to_space(
                   adata_sc, 
                   adata_st,         
                   mode='clusters',
                   cluster_label='seurat_clusters',device="cuda:0")
else:
    ad_map = tg.map_cells_to_space(
                   adata_sc, 
                   adata_st,         
                   mode='clusters',
                   cluster_label='seurat_clusters')




tg.project_cell_annotations(ad_map, adata_st, annotation="seurat_clusters")


adata_st.obsm['tangram_ct_pred'].to_csv(sys.argv[3])
Python Pandas torch Tangram From line 3 of scripts/spatial-tangram.py 
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
library(Seurat)
require(tidyverse)
require(patchwork)

source("workflow/scripts/spatial-functions.R")

arguments=commandArgs(TRUE)

sampleID=arguments[1]
res=arguments[2]

Spatial_Data=readRDS(paste0("rds/",sampleID,".rds"))

function_image_fixer(Spatial_Data) -> Spatial_Data

p1 <- DimPlot(Spatial_Data, reduction = "umap", label = TRUE,label.size = 10,group.by = paste0("SCT_snn_res.",res)) 
p2 <- SpatialDimPlot(Spatial_Data, label = TRUE, label.size = 6,group.by = paste0("SCT_snn_res.",res),images=paste0("image")) 

suppressWarnings(wrap_plots(p1,p2,ncol=2) -> wp)


ggsave(plot =wp,filename=paste0("results/",sampleID,"/resolution-",res,"/",sampleID,".umap.spatial",".pdf"),width=13,height=7)
R tidyverse Seurat patcHwork From line 2 of scripts/spatial-umap.R
ShowHide 38 more snippets with no or duplicated tags.

Login to post a comment if you would like to share your experience with this workflow.

Do you know this workflow well? If so, you can request seller status , and start supporting this workflow.

Free

Created: 1yr ago
Updated: 1yr ago
Maitainers: public
URL: https://github.com/sinanugur/spatial-workflow
Name: spatial-workflow
Version: 1
Badge:
workflow icon

Insert copied code into your website to add a link to this workflow.

Downloaded: 0
Copyright: Public Domain
License: None
  • Future updates

Related Workflows

psychip_snakemake
Show Details View Workflow
ENCODE pipeline for histone marks developed for the psychENCODE project
public

public

psychip pipeline is an improved version of the ENCODE pipeline for histone marks developed for the psychENCODE project. The o...
raw sequence reads Alignment Sequence alignment report macs2 ucsc-bedclip bedGraphToBigWig BEDTools BWA Picard SAMtools Snakemake
Free
74
ncov_2
Show Details View Workflow
Near-real time tracking of SARS-CoV-2 in Connecticut
public

public

Repository containing scripts to perform near-real time tracking of SARS-CoV-2 in Connecticut using genomic data. This pipeli...
JSON nextclade Augur Biopython FOCUS Pandas Snakemake bs4 epiweeks geopy matplotlib numpy pycountry pycountry-convert uszipcode
Free
28
cellranger-snakemake-gke
Show Details View Workflow
snakemake workflow to run cellranger on a given bucket using gke.
public

public

A Snakemake workflow for running cellranger on a given bucket using Google Kubernetes Engine. The usage of this workflow ...
macs2 ucsc-bedclip bedGraphToBigWig BEDTools BWA Picard SAMtools Snakemake
Free
74
atlas
Show Details View Workflow
ATLAS - Three commands to start analyzing your metagenome data
public

public

Metagenome-atlas is a easy-to-use metagenomic pipeline based on snakemake. It handles all steps from QC, Assembly, Binning, t...
raw sequence reads Genome assembly Annotation track checkm2 gunc prodigal snakemake-wrapper-utils MEGAHIT Atlas BBMap Biopython BioRuby Bwa-mem2 cd-hit CheckM DAS Diamond eggNOG-mapper v2 MetaBAT 2 Minimap2 MMseqs MultiQC Pandas Picard pyfastx SAMtools SemiBin Snakemake SPAdes SqueezeMeta TADpole VAMB CONCOCT ete3 gtdbtk h5py networkx numpy plotly psutil utils metagenomics
Free
175
rna-seq-star-deseq2
Show Details View Workflow
RNA-seq workflow using STAR and DESeq2
public

public

This workflow performs a differential gene expression analysis with STAR and Deseq2. The usage of this workflow is described ...
raw sequence reads Gene expression profiling ID list snakemake-wrapper-utils biomaRt PCAtools BiocParallel BWA Cutadapt DESeq2 GFFutils MultiQC Pandas SAMtools Snakemake STAR cli stringr tidyverse RNA-Seq
Free
41
dna-seq-gatk-variant-calling
Show Details View Workflow
This Snakemake pipeline implements the GATK best-practices workflow
public

public

This Snakemake pipeline implements the GATK best-practices workflow for calling small germline variants. The usage of thi...
VCF raw sequence reads Variant calling genetic variants gatk rust-bio-tools snakemake-wrapper-utils tabix BCFtools BWA FastQC MultiQC Pandas Picard SAMtools Snakemake Trimmomatic Variant Effect Predictor (VEP) common matplotlib numpy seaborn DNA
Free
57