-
Notifications
You must be signed in to change notification settings - Fork 0
/
03_seurat_per_treatment.R
91 lines (56 loc) · 3.03 KB
/
03_seurat_per_treatment.R
1
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
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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
library(Seurat)
library(tidyverse)
library(Matrix)
library(SeuratDisk)
seurat <- readRDS("./data_output/seurat_raw_annotated_filtered.rds")
strrep = sub(pattern = "\\_2","", [email protected]$replicate_name)
[email protected]$treatment = strrep
samples = unique([email protected]$treatment)
for (i in 1:length(samples)){
condition = samples[i]
subseurat = subset(seurat, subset = treatment == condition)
# determine which samples are from which batch
mutate(batch =
ifelse(replicate_name == paste0(condition), 1, 2))
# scaling per condition
seurat_sct <- SCTransform(subseurat, vars.to.regress = "batch")
print(paste0("Scaling done for condition ",condition," which is ",i,"/16"))
# QC plots
Idents(seurat_sct) = [email protected]$treatment
p1 = VlnPlot(subseurat, features = c("nFeature_RNA", "nCount_RNA", "mitoPercent"),
ncol = 3, pt.size = 0.0001)
ggsave(plot = p1, paste0("./figures/per_treatment/",condition,"_violin_features.png"))
p2 = FeatureScatter(subseurat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") +
theme(legend.position="none")
ggsave(plot =p2, paste0("./figures/per_treatment/",condition,"_scatter_features.png"))
# pca
Idents(seurat_sct) = [email protected]$treatment
seurat_sct <- RunPCA(seurat_sct, features = VariableFeatures(object = seurat_sct))
p3 = DimPlot(seurat_sct, reduction = "pca")+
labs(title = "All data filtered and scaled")
ggsave(plot =p3, paste0("./figures/per_treatment/",condition,"_PCA.png"))
# umap
seurat_sct <- FindNeighbors(seurat_sct, dims = 1:40)
seurat_sct <- FindClusters(seurat_sct, resolution = 0.5)
seurat_sct <- RunUMAP(seurat_sct, dims = 1:40, seed.use = 054057, n.neighbors = 5)
p4 = DimPlot(seurat_sct, reduction = "umap")+
labs(title = "n neighbors = 5")
ggsave(plot =p4, paste0("./figures/per_treatment/",condition,"_UMAP_nn5.png"))
seurat_sct <- RunUMAP(seurat_sct, dims = 1:40, seed.use = 054057, n.neighbors = 10)
p5 = DimPlot(seurat_sct, reduction = "umap")+
labs(title = "n neighbors = 10")
ggsave(plot =p5, paste0("./figures/per_treatment/",condition,"_UMAP_nn10.png"))
print(paste0("UMAP done for condition ",condition," which is ",i,"/16"))
write_rds(seurat_sct, paste0("./data_output/per_treatment/",condition,"_seurat_scaled.rds"))
print(paste0("Saving final rds object of ",condition," which is ",i,"/16"))
seurat_sct[["RNA"]] <- as(object = seurat_sct[["RNA"]], Class = "Assay")
seurat_sct[["SCT"]] <- as(object = seurat_sct[["SCT"]], Class = "Assay")
SaveH5Seurat(seurat_sct, filename = paste0("./data_output/per_treatment/",
condition,"_seurat_sct.h5Seurat"),
overwrite = T)
Convert(paste0("./data_output/per_treatment/",
condition,"_seurat_sct.h5Seurat"), dest = "h5ad", overwrite = TRUE)
rm(seurat_sct)
}