-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathReading_in_pre_processing_Seurat.R
57 lines (38 loc) · 1.88 KB
/
Reading_in_pre_processing_Seurat.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
setwd("your_path") #must contain dirs (one for each sample) with barcodes.tsv.gz, matrix.mtx.gz, features.tsv.gz
data.10x = list()
dirs <- list.dirs(".", recursive = FALSE)
for (i in 1:length(dirs)) {
data.10x[[i]] <- Read10X(data.dir = dirs[[i]])
}
names(data.10x) <- sub('./', '', dirs)
scrna.list = list()
samples<-sub('./', '', dirs)
for (i in 1:length(data.10x)) {
scrna.list[[i]] = CreateSeuratObject(counts = data.10x[[i]], min.cells=3, min.features=200, project=samples[i]);
scrna.list[[i]] =NormalizeData(object = scrna.list[[i]]);
scrna.list[[i]] =ScaleData(object = scrna.list[[i]]);
scrna.list[[i]] =FindVariableFeatures(object = scrna.list[[i]]);
scrna.list[[i]] =RunPCA(object = scrna.list[[i]], verbose = FALSE)
scrna.list[[i]][["percent.mt"]] = PercentageFeatureSet(object=scrna.list[[i]], pattern = "^MT-")
}
names(scrna.list) <- sub('./', '', dirs)
anchors <- FindIntegrationAnchors(object.list = scrna.list, reduction = "rpca", dims = 1:50)
aggr <- IntegrateData(anchorset = anchors, dims = 1:50)
VlnPlot(aggr, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
aggr <- subset(aggr, subset = nFeature_RNA > 500 & percent.mt < 10)
aggr <- FindVariableFeatures(aggr)
aggr <- ScaleData(aggr, verbose = FALSE)
aggr <- RunPCA(aggr, verbose = FALSE)
aggr <- RunUMAP(aggr, dims = 1:50)
aggr <- FindNeighbors(aggr, dims = 1:20)
aggr <- FindClusters(aggr, resolution = c(0.01, 0.05, 0.1, 0.2, 0.3), graph.name = 'integrated_snn')
Idents(aggr)<-'integrated_snn_res.0.01'
res0.01markers<-FindAllMarkers(aggr, only.pos = T)
Idents(aggr)<-'integrated_snn_res.0.05'
res0.05markers<-FindAllMarkers(aggr, only.pos = T)
Idents(aggr)<-'integrated_snn_res.0.1'
res0.1markers<-FindAllMarkers(aggr, only.pos = T)
Idents(aggr)<-'integrated_snn_res.0.2'
res0.2markers<-FindAllMarkers(aggr, only.pos = T)
dev.off()
save.image("you_path/analysis.RData")