Extrahieren k-bedeutet Cluster-spezifische Funktionen

Question

Ich bin Clustering einige Genexpressionsdaten mit k-Mittel in der reduzierten, PCA-Raum und jetzt möchte ich verschiedene Features extrahieren, die am besten beschreiben jeden Cluster. Dies sind Merkmale, die in jedem Cluster stark ausgeprägt sind.

Ich habe unten ein reproduzierbares Beispiel veröffentlicht, um meine Logik zu zeigen und wo ich aufgehört habe.

# Create test matrix 
test = matrix(rnorm(200), 20, 10) 
test[1:10, seq(1, 10, 2)] = test[1:10, seq(1, 10, 2)] + 3 
test[11:20, seq(2, 10, 2)] = test[11:20, seq(2, 10, 2)] + 2 
test[15:20, seq(2, 10, 2)] = test[15:20, seq(2, 10, 2)] + 4 
colnames(test) = paste("Cell", 1:10, sep = "") 
rownames(test) = paste("Gene", 1:20, sep = "") 

# plot the inital heatmap 
library(pheatmap) 
pheatmap(t(test)) 

# preform PCA 
pca = prcomp(t(test), center=TRUE, scale=TRUE) 
rotation = data.frame(pca$x) 
plot(rotation[1:3], pch=16, cex=0.6, cex.main=0.9) 

# preform Kmeans in PCA space 
wss = (nrow(rotation)-1)*sum(apply(rotation,2,var)) 
for (i in 2:9) wss[i] <- sum(kmeans(rotation, centers=i)$withinss) 
plot(1:9, wss, type="b", xlab="Number of Clusters", ylab="Within groups sum of squares") 
km = kmeans(rotation, 2) 
cluster_assignment = as.factor(km$cluster) 

# plot k-means cluster assignment in PCA space 
library(ggplot2) 
ggplot(rotation, aes(rotation$PC1, rotation$PC2, color=cluster_assignment, label=rownames(rotation))) + geom_point() + geom_text() 

# create a cluster annotation data.frame 
n_clusters = length(unique(km$cluster)) 
temp_cluster_design_list = list() 
for(n in 1:n_clusters){ 
    temp_cluster_design = data.frame(row.names=row.names(t(test)[cluster_assignment %in% n,])) 
    temp_cluster_design$cluster = n 
    temp_cluster_design_list[[n+1]] <- temp_cluster_design 
} 
cluster_design = do.call(rbind, temp_cluster_design_list) 

# cluster_design looks something like this: 
#  cluster 
# Cell1  1 
# Cell3  1 
# Cell5  1 
# Cell7  1 
# Cell9  1 
# Cell2  2 
# Cell4  2 
# Cell6  2 
# Cell8  2 
# Cell10  2 

# heatmap with cell cluster annotation 
PCA_heatmap_data = t(test)[row.names(cluster_design),] 
cluster_design$cluster = as.factor(cluster_design$cluster) 
pheatmap(PCA_heatmap_data, annotation_row=cluster_design) 


# extract out expressed genes from each cluster 
for(n in 1:n_clusters) { 
    temp_cluster = t(test)[cluster_assignment %in% n,] 
#  ??? somehow extract out the expressed gene names specific to each cluster 
} 

Der obige Code schließlich lässt mich mit einem Heatmap, die so etwas wie this aussieht. Nun, was ich mit jedem Cluster in der Heatmap machen möchte, ist, Gen-Namen, die stark exprimiert sind, zu extrahieren. Ich würde schließlich gerne einen Tisch schreiben, die etwa wie folgt aussieht:

GENE  CLUSTER 
Gene20 cluster2 
Gene19 cluster2 
Gene15 cluster2 
Gene18 cluster2 
Gene16 cluster2 
Gene17 cluster2 
Gene9  cluster1 
Gene8  cluster1 
Gene4  cluster1 
Gene3  cluster1 
...  ... 

Ich bin nicht sicher, der beste und effizienteste Weg, um dies zu realisieren. Ich würde jede Hilfe schätzen, die Sie meinen Weg werfen konnten! Vielen Dank!

EDIT

Wie definieren wir stark exprimiert? Das bin ich mir nicht ganz sicher und hoffe auf einen Einblick. Vielleicht den Durchschnitt für jedes Gen über alle Zellen hinweg vergleichen und das mit dem Durchschnitt im Cluster vergleichen? Das mag funktionieren, aber ich denke, dass dies durch Ausreißer verzerrt wird. Ein anderer Gedanke wäre, jeden Cluster zu skalieren und Gene zu nehmen, die signifikant stark exprimiert werden?

Wie würde ich bei der Identifizierung von stark exprimierten Genen unter sehr ähnlichen Clustern vorgehen? Zum Beispiel zeigt die Heatmap this drei Cluster, aber die roten und grünen Cluster sind sehr ähnlich. Es scheint, dass Gene9 cluster2-spezifisch ist?

Answer 1

threshold <- 1 
nms  <- as.character() 
for(n in 1:n_clusters) { 

    for(i in 1:ncol(t(test)[row.names(t(test)) %in% names(cluster_assignment[cluster_assignment == n]),])){ 
    if(mean(t(test)[row.names(t(test)) %in% names(cluster_assignment[cluster_assignment == n]),i]) > threshold){ 
      nms <- c(nms, colnames(t(test)[row.names(t(test)) %in% names(cluster_assignment[cluster_assignment == n]),])[i]) 
    } 
    } 
    if(n == 1) result <- data.frame(GENE = nms, CLUSTER = rep(n,length(nms))); rm(nms); nms <- as.character() 
    if(n > 1) result <- rbind(result, data.frame(GENE = nms, CLUSTER = rep(n,length(nms)))) 
} 

result 

 GENE CLUSTER 
1 Gene7  1 
2 Gene11  1 
3 Gene12  1 
4 Gene13  1 
5 Gene14  1 
6 Gene15  1 
7 Gene16  1 
8 Gene17  1 
9 Gene18  1 
10 Gene19  1 
11 Gene20  1 
12 Gene1  2 
13 Gene2  2 
14 Gene3  2 
15 Gene4  2 
16 Gene5  2 
17 Gene6  2 
18 Gene7  2 
19 Gene8  2 
20 Gene9  2 
21 Gene10  2 

verließ ich die threshold als Parameter, so dass Sie es jedoch definieren, die Sie mögen. In diesem Beispiel habe ich 1 als Schwellenwert verwendet.