-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathQuickAnalysis100kcells.R
77 lines (62 loc) · 3.4 KB
/
QuickAnalysis100kcells.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
require("Seurat")
liver.integrated <- readRDS("basicintegration.rds")
set.seed(7742)
require("sctransform")
liver.integrated <- ScaleData(liver.integrated);
liver.integrated <- RunPCA(liver.integrated, features = VariableFeatures(object = liver.integrated))
ElbowPlot(liver.integrated)
set.seed(8291)
npcs <- 10
res <- 1.2
nkNN <- 20
liver.integrated <- FindNeighbors(liver.integrated, dims = 1:npcs, k.param=nkNN)
liver.integrated <- FindClusters(liver.integrated, resolution = res, k.param=nkNN)
liver.integrated <- RunTSNE(liver.integrated, dims = 1:npcs)
liver.integrated <- RunUMAP(liver.integrated, dims = 1:npcs, parallel=FALSE, n.neighbour=nkNN)
png("Highres_integrated_umap.png", width=6, height=6, units="in", res=50)
DimPlot(liver.integrated, reduction = "umap")
dev.off()
png("Highres_integrated_tsne.png", width=6, height=6, units="in", res=50)
DimPlot(liver.integrated, reduction = "tsne")
dev.off()
png("Donor_integrated_umap.png", width=6, height=6, units="in", res=50)
DimPlot(liver.integrated, reduction = "umap", group.by="orig.ident")
dev.off()
png("Donor_integrated_tsne.png", width=6, height=6, units="in", res=50)
DimPlot(liver.integrated, reduction = "tsne", group.by="orig.ident")
dev.off()
png("Integrated_clusters_by_donor.png", width=6, height=6, units="in", res=50)
barplot(table([email protected]$orig.ident, [email protected]$seurat_clusters), col=rainbow(20))
dev.off()
#saveRDS(liver.integrated, "basic_integration_analysis.rds");
# Automated annotation? - doesn't work as too few genes left, but does work on original datasets.
# automated annotation at individual dataset level then add it here.
source("Setup_autoannotation.R")
#liver.integrated <- readRDS("basic_integration_analysis.rds")
#test <- Use_markers_for_anno(liver.integrated@assays$RNA$data, [email protected]$seurat_clusters)
all_anno <- readRDS("All20_automatedannotation.rds");
#cell_barcodes <- unlist(lapply(strsplit(rownames([email protected]), "_"), function(x){return(x[[1]])}))
#[email protected]$cell_barcode <- cell_barcodes
[email protected]$scmap_anno <- rep("unknown", ncol(liver.integrated));
[email protected]$scmap_anno2 <- rep("unknown", ncol(liver.integrated));
for (donor in unique([email protected]$orig.ident)) {
cell_ids <- [email protected]$cell_barcode[[email protected]$donor == donor]
anno <- all_anno[[donor]];
anno <- anno[anno$cell_barcode %in% cell_ids,]
anno <- anno[match(cell_ids, anno$cell_barcode),]
[email protected]$scmap_anno[[email protected]$orig.ident == donor] <- as.character(anno$scmap_cluster_anno$lm1)
[email protected]$scmap_anno2[[email protected]$orig.ident == donor] <- as.character(anno$scmap_cell_anno)
}
png("AutoAnno_integrated_umap.png", width=8, height=6, units="in", res=50)
DimPlot(liver.integrated, reduction = "umap", group.by="scmap_anno")
dev.off()
png("AutoAnno_integrated_tsne.png", width=8, height=6, units="in", res=50)
DimPlot(liver.integrated, reduction = "tsne", group.by="scmap_anno")
dev.off()
png("AutoAnno2_integrated_umap.png", width=8, height=6, units="in", res=50)
DimPlot(liver.integrated, reduction = "umap", group.by="scmap_anno2")
dev.off()
png("AutoAnno2_integrated_tsne.png", width=8, height=6, units="in", res=50)
DimPlot(liver.integrated, reduction = "tsne", group.by="scmap_anno2")
dev.off()
saveRDS(liver.integrated, "basic_integration_analysis.rds");