-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathMap2.2_Select_Clustering.R
124 lines (98 loc) · 5.21 KB
/
Map2.2_Select_Clustering.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
proj_name="harmony_integrated2_20livers"
require(Seurat)
liver.integrated <- readRDS(paste(proj_name, "harmony_plus_clusters.rds", sep="_"))
# Compare all these clusterings
require(igraph)
require(gplots)
clust_table <- [email protected][, grepl("^knn_", names([email protected]))]
clust_table <- as.matrix(apply(clust_table,2,as.numeric))
require("proxy")
clust_dists <- proxy::dist(clust_table, method=function(x,y){igraph::compare(x,y,method="vi")}, by_rows=FALSE)
clust_similr1 <- proxy::simil(clust_table, method=function(x,y){igraph::compare(x,y,method="nmi")}, by_rows=FALSE)
clust_similr2 <- proxy::simil(clust_table, method=function(x,y){igraph::compare(x,y,method="adjusted.rand")}, by_rows=FALSE)
# Find robust exemplar clustering(s)
require("apcluster")
require("gplots")
set.seed(18371)
res1 <- apcluster(-1*as.matrix(clust_dists), p=-2.5)
coarse_lvl <- names(res1@exemplars)[2]
fine_lvl <- names(res1@exemplars)[3]
#manually select which exemplar to use
[email protected]$Coarse_clusters <- [email protected][[coarse_lvl]]
[email protected]$Fine_clusters <- [email protected][[fine_lvl]]
#apcluster::heatmap(res1, -1*as.matrix(clust_dists))
png(paste(proj_name,"compare_clusterings_heatmap.png",sep="_"), width=6, height=6, units="in", res=300)
lab <- matrix("", ncol=ncol(clust_table), nrow=ncol(clust_table))
lab[colnames(clust_table)==fine_lvl, colnames(clust_table)==fine_lvl] <- "F"
lab[colnames(clust_table)==coarse_lvl, colnames(clust_table)==coarse_lvl] <- "C"
heatmap.2(as.matrix(clust_dists), trace="none", distfun=function(x){return(as.dist(clust_dists))}, cellnote=lab)
dev.off()
# Visualize the Chosen clusterings
print("plotting1")
png(paste(proj_name,"coarse_umap.png", sep="_"), width=6, height=6, units="in", res=150)
DimPlot(liver.integrated, reduction = "umap", group.by="Coarse_clusters")
dev.off()
png(paste(proj_name,"coarse_tsne.png",sep="_"), width=6, height=6, units="in", res=150)
DimPlot(liver.integrated, reduction = "tsne", group.by="Coarse_clusters")
dev.off()
png(paste(proj_name,"fine_umap.png", sep="_"), width=6, height=6, units="in", res=150)
DimPlot(liver.integrated, reduction = "umap", group.by="Fine_clusters")
dev.off()
png(paste(proj_name,"fine_tsne.png", sep="_"), width=6, height=6, units="in", res=150)
DimPlot(liver.integrated, reduction = "tsne", group.by="Fine_clusters")
dev.off()
saveRDS(liver.integrated, paste(proj_name, "harmony_plus_analysis.rds", sep="_"))
print("cluster anno")
#Use autoannotation to label clusters
cluster_assign <- function(scAnno) {
scAnno <- factor(scAnno)
freqs <- table(scAnno)/length(scAnno);
# specific labels - majority rule
if (max(freqs) > 0.5) {
return(levels(scAnno)[freqs == max(freqs)])
}
# general labels - majority rule
scAnno <- as.character(scAnno);
scAnno[grepl("Tcell", scAnno)] <- "Tcell"
scAnno[grepl("Bcell", scAnno)] <- "Bcell"
scAnno[grepl("LSEC", scAnno)] <- "LSEC"
scAnno[grepl("Macrophage", scAnno)] <- "Macrophage"
scAnno[grepl("Hep", scAnno)] <- "Hepatocyte"
scAnno <- factor(scAnno)
freqs <- table(scAnno)/length(scAnno);
if (max(freqs) > 0.5) {
return(levels(scAnno)[freqs == max(freqs)])
} else {
return("ambiguous")
}
}
consistent_cluster_labs <- sapply(split([email protected]$consistent_labs, [email protected]$Fine_clusters), cluster_assign)
general_cluster_labs <- sapply(split([email protected]$general_labs, [email protected]$Fine_clusters), cluster_assign)
[email protected]$cluster_quickanno <- consistent_cluster_labs[[email protected]$Fine_clusters];
# Cluster ID labelled figures
agg_coord_by_cluster <- function(coords, clusters) {
x <- split(seq(nrow(coords)), clusters)
result <- sapply(x, function(a) apply(coords[a,],2,median))
return(result)
}
tsne_lab_pos <- agg_coord_by_cluster(liver.integrated@[email protected], [email protected]$Fine_clusters)
umap_lab_pos <- agg_coord_by_cluster(liver.integrated@[email protected], [email protected]$Fine_clusters)
lab_id <- colnames(tsne_lab_pos)
require("ggplot2")
source("~/scripts/LiverMap2.0/Colour_Scheme.R")
new_colour_scheme <- Cell_type_colours[order(Cell_type_colours[,1]),]
[email protected]$short_cluster_anno <- factor(map_cell_types([email protected]$cluster_quickanno), levels=new_colour_scheme[,1]);
new_colour_scheme <- new_colour_scheme[new_colour_scheme[,1] %in% [email protected]$shoart_cluster_anno,]
print("plotting")
png("Autoanno_label_harmony_umap.png", width=7.5, height=6, units="in", res=300)
DimPlot(liver.integrated, reduction="umap", group.by="short_cluster_anno", pt.size=.1)+scale_color_manual(values=new_colour_scheme[,2])+annotate("text", x=umap_lab_pos[1,], y=umap_lab_pos[2,], label=lab_id, colour="grey35")
dev.off()
png("Autoanno_label_harmony_tsne.png", width=7.5, height=6, units="in", res=300)
DimPlot(liver.integrated, reduction="tsne", group.by="short_cluster_anno", pt.size=.1)+scale_color_manual(values=new_colour_scheme[,2])+annotate("text", x=tsne_lab_pos[1,], y=tsne_lab_pos[2,], label=lab_id, colour="grey35")
dev.off()
# T/NK = 10, 5, + 20
# Endo = 15, 12, 8, 19, 21
# Mac = 9, 11, 13, 25
# Hep Traj = 4, 1, 2, 24, 17, 3, 6, 7, 26
# Doublets? 20
# Hep / Erythroid = 0, 14, 27