-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathKalinski_reanalysis.R
131 lines (96 loc) · 4.84 KB
/
Kalinski_reanalysis.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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
suppressMessages(library(ggbeeswarm))
suppressMessages(library(dplyr))
suppressMessages(library(Seurat))
suppressMessages(library(Matrix))
suppressMessages(library(gplots))
suppressMessages(library(ggplot2))
suppressMessages(library(openxlsx))
suppressMessages(library(cowplot))
suppressMessages(library(patchwork))
library(scales)
library(viridis)
library(SeuratWrappers)
library(slingshot)
require(BiocStyle)
library(data.table)
library(SingleCellExperiment)
sessionInfo()
###############################################
# Examine Uninjured post injury samples from kalinski et al 2020.
# Analysis of the immune response to sciatic nerve injury identifies efferocytosis as a key mechanism of nerve debridement
###############################################
rep1_path='/Users/maurizio.aurora/Dropbox (HSR Global)/WORKSPACE/Bonanomi/Bonanomi_1287_scRNA_injury/7_bioinfo/public_data/Kalinski/GSE153762_RAW/rep1'
list.files(rep1_path)
rep1 <- Read10X(data.dir = rep1_path)
rep1 <- CreateSeuratObject(counts = rep1, project = "rep1", min.cells = 5, min.features = 200)
rep1[["percent.mt"]] <- PercentageFeatureSet(rep1, pattern = "^mt-")
rep1[["percent.Rpl"]] <- PercentageFeatureSet(rep1, pattern = "Rpl")
rep1$stim <- "3D_kalinski_rep1"
options(repr.plot.width=8, repr.plot.height=6)
VlnPlot(rep1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
min_nFeature_RNA = 200
max_nFeature_RNA = 7500
max_percent_MT = 20
rep2_path='/Users/maurizio.aurora/Dropbox (HSR Global)/WORKSPACE/Bonanomi/Bonanomi_1287_scRNA_injury/7_bioinfo/public_data/Kalinski/GSE153762_RAW/rep2'
list.files(rep2_path)
rep2 <- Read10X(data.dir = rep2_path)
rep2 <- CreateSeuratObject(counts = rep2, project = "rep2", min.cells = 5, min.features = 200)
rep2[["percent.mt"]] <- PercentageFeatureSet(rep2, pattern = "^mt-")
rep2[["percent.Rpl"]] <- PercentageFeatureSet(rep2, pattern = "Rpl")
rep2$stim <- "3D_kalinski_rep2"
options(repr.plot.width=8, repr.plot.height=6)
VlnPlot(rep2, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
min_nFeature_RNA = 200
max_nFeature_RNA = 7500
max_percent_MT = 20
rep3_path='/Users/maurizio.aurora/Dropbox (HSR Global)/WORKSPACE/Bonanomi/Bonanomi_1287_scRNA_injury/7_bioinfo/public_data/Kalinski/GSE153762_RAW/rep3'
list.files(rep3_path)
rep3 <- Read10X(data.dir = rep3_path)
rep3 <- CreateSeuratObject(counts = rep3, project = "rep3", min.cells = 5, min.features = 200)
rep3[["percent.mt"]] <- PercentageFeatureSet(rep3, pattern = "^mt-")
rep3[["percent.Rpl"]] <- PercentageFeatureSet(rep3, pattern = "Rpl")
rep3$stim <- "3D_kalinski_rep3"
options(repr.plot.width=8, repr.plot.height=6)
VlnPlot(rep3, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
min_nFeature_RNA = 200
max_nFeature_RNA = 7500
max_percent_MT = 20
# normalize and find varible features in the objects
object_clean_new.list <- lapply(X = c(rep1, rep2, rep3 ), FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
# find integration anchors
integrated <- FindIntegrationAnchors(object.list = object_clean_new.list, dims = 1:35)
features.to.integrate = integrated@anchor.features
integrated <- IntegrateData(anchorset = integrated, dims = 1:35, features.to.integrate = features.to.integrate)
DefaultAssay(integrated) <- "integrated"
table(integrated$stim)
res = 1
# Run the standard workflow for visualization and clustering
integrated <- ScaleData(integrated, verbose = FALSE)
integrated <- RunPCA(integrated, npcs = 30, verbose = FALSE)
# t-SNE and Clustering
integrated <- RunUMAP(integrated, reduction = "pca", dims = 1:30)
integrated <- RunTSNE(integrated, reduction = "pca", dims = 1:30)
integrated <- FindNeighbors(integrated, reduction = "pca", dims = 1:30)
integrated <- FindClusters(integrated, resolution = res)
DimPlot(integrated, reduction = "umap", group.by = "stim")
DimPlot(integrated, reduction = "umap", split.by = "stim")
DefaultAssay(integrated) = "RNA"
saveRDS(integrated, "integrated_all_CT_Kalinski.Rds")
FeaturePlot(integrated, "Rgs5")
integrated_sub = subset(integrated, idents = c("20","10","32","25"))
integrated_sub <- NormalizeData(integrated_sub)
integrated_sub<- FindVariableFeatures(integrated_sub, selection.method = "vst", nfeatures = 2000)
res = 1
# Run the standard workflow for visualization and clustering
integrated_sub <- ScaleData(integrated_sub, verbose = FALSE)
integrated_sub <- RunPCA(integrated_sub, npcs = 30, verbose = FALSE)
# t-SNE and Clustering
integrated_sub <- RunUMAP(integrated_sub, reduction = "pca", dims = 1:30)
integrated_sub <- FindNeighbors(integrated_sub, reduction = "pca", dims = 1:30)
integrated_sub <- FindClusters(integrated_sub, resolution = res)
#identify EC after removing pericytes
integrated_sub_new = subset(integrated_sub, idents = c("0","1","2","3","4","5","6","7","8"))
saveRDS(integrated_sub_new, "integrated_EC_Kalinski_no_pericytes.Rds")