我順便把LouvainKRange 也寫一下，方便后來人復(fù)現(xiàn)。感謝樓主
LouvainKRange <- function(data.obj, kmin = 5, kmax = 50, kstep = 5) {
require(igraph, quietly = TRUE)
# check if seurat object
if (class(data.obj)[1] == "Seurat") {
dist.mat <- data.obj@assays[[data.obj@active.assay]]@misc$dist.mat
} else {
dist.mat <- data.obj
}
# create lists
cluster.list <- list()
sil.list <- list()
# setup iteration
k <- kmin
while (k <= kmax) {
print(paste("Clustering with k = ", k, "...", sep = ''))
# generate clustering and silhouette score
clust.vec <- LouvainClust(dist.mat, k)
sil.score <- cluster::silhouette(clust.vec, dist.mat)
# add to list
k.ind <- paste('k', k, sep = '.')
cluster.list[[k.ind]] <- clust.vec
sil.list[[k.ind]] <- mean(sil.score[,3])
# iterate
k <- k + kstep
}
# identify optimal cluster
opt.clust <- cluster.list[[which.max(sil.list)]]
# add to data.object
if (class(data.obj)[1] == "Seurat") {
data.obj@assays[[data.obj@active.assay]]@misc[['pisces.cluster']] <- opt.clust
data.obj@assays[[data.obj@active.assay]]@misc[['clustering.obj']] <- list('clusterings' = cluster.list, 'sils' = sil.list)
return(data.obj)
} else {
return(list('pisces.cluster' = opt.clust, 'clustering.obj' = list('clusterings' = cluster.list, 'sils' = sil.list)))
}
}
基于單細(xì)胞測序的蛋白活性推斷（PISCES）分析流程筆記，凎！寫在前面 感謝califano-lab[https://github.com/califano-lab]團隊將代碼無私的分享出來=。=！ 本流程需要一點點R語言與linux基...
