接下来我们继续多时间点样本实战分析流程的第二部分：聚类和富集分析。第一部分的完整流程请参照：RNA-Seq 分析流程：多时间点样本分析实战（一）

多时间点数据的聚类

前面我们发现70% 的基因是差异表达的，几乎所有通路都受到处理的影响。因此，分析流程的下一步是根据基因表达对处理的动态反应进行聚类。

过滤

在对基因进行聚类之前，我们首先将感兴趣的基因集浓缩为(1)显著差异表达的基因；(2) 处理之间存在较大倍数的基因。减少执行聚类的基因集可以增强聚类的可信度。为此，我们首先使用Fisher方法将时序差异表达步骤中获得的所有p值聚合为单个p值(pvalues_fisher_method)。接下来，我们选择Fisher方法调整的p值低于0.05且至少在一种处理和一个时间点的对数倍数变化至少为2的所有基因。

# Then rank by fisher's p-value and take max the number of genes of interest
# Filter out q-values for the pvalues table
fishers_pval = pvalues_fisher_method(pvalues)
qvalues = apply(pvalues, 2, p.adjust)
fishers_qval = p.adjust(fishers_pval)

genes_to_keep = row.names(
    log_fold_change_max[
    (rowSums(log_fold_change_max > 2) > 0) &
    (fishers_qval < 0.05), ])
# Keep the data corresponding to the genes of interest in another variable.
# by subsetting the `moanin_model`, which contains the data.
de_moanin_model = moanin_model[genes_to_keep,]

基于样条拟合的聚类

过滤后，我们剩下5521个基因。然后我们就可以开始聚类的常规分析流程。通过观察差异表达的基因发现，许多基因具有相似的基因表达模式，但上下调的幅度不同。

因此，我们建议对k-means进行以下调整：
1.样条函数拟合：对于每个基因，计算(如moanin对象中包含的那样)拟合样条函数。

2.重新归一化样条函数：对于每个基因，重新归一化估计的样条函数，使得值限制在0和1之间，从而在基因之间具有可比性。

3.K-means：对样条函数重新调整后的拟合值应用k-means来估计聚类。

样条函数的重新归一化旨在使每个基因达到可比较的尺度，类似于通常在没有时间维度的基因表达研究中对基因表达进行的归一化和缩放。

聚类步骤由moanin中的splines_kmeans函数执行。我们先将cluster数设置为8，不过后面我们将讨论选择最佳cluster数的问题。

# First fit the kmeans clusters
kmeans_clusters = splines_kmeans(de_moanin_model, n_clusters=8,
    random_seed=42,
    n_init=20)

splines_kmeans函数返回一个命名列表，包含：

• centroids：包含cluster质心的矩阵。矩阵的维度为(n_centroid, n_samples)。

• cluster：大小为n_genes的向量，包含每个基因的kmeans聚类被分配到哪一个cluster的信息。

然后，我们使用plot_splines_data函数来可视化通过样条k-means模型获得的每个cluster的质心(图12)。

plot_splines_data(de_moanin_model,
    data=kmeans_clusters$centroids, 
    colors=ann_colors$Group,
    smooth=TRUE)

图12，K-means centroids

由于我们对样条拟合进行了重新调整，这些质心的范围为 0-1，并且并不代表实际的基因表达水平。在图13中，我们绘制了一些聚类到cluster2的实际基因：

cluster_of_interest = 2
cluster2Genes = names(
    kmeans_clusters$clusters[kmeans_clusters$clusters==cluster_of_interest])
plot_splines_data(de_moanin_model,  
    centroid=kmeans_clusters$centroids[cluster_of_interest,], 
    colors=ann_colors$Group, smooth=TRUE,simpleY =FALSE,
    subset_data=cluster2Genes[3:6],
    mar=c(1.5,2.5,2,0.1))

图13，Genes from cluster 2

将基因分配给cluster

正如我们所看到的，虽然这些基因与cluster质心的模式有一些相似性，但从匹配质心的意义上来说，这些特定基因和cluster并没有达到最佳拟合。由于基因与cluster的拟合程度存在差异，我们希望能够对基因与cluster的拟合程度进行评分。

此外，前面我们只计算了根据过滤得到的基因子集，我们希望有一种方法将所有基因分配到各个cluster里面。

因此，聚类分析的下一步部分是评分scoring和标记label。每个基因都会获得一个分数，该分数对应于每个基因和每个cluster之间的拟合优度。评分和标签是通过splines_kmeans_score_and_label函数实现。该函数计算基因与聚类质心的拟合度，并在得分足够高的情况下为基因提供聚类标签。

# Then assign scores and labels to *all* the data, using a goodness-of-fit
# scoring function.
scores_and_labels = splines_kmeans_score_and_label(
    moanin_model, kmeans_clusters)

Scores_and_labels列表包含三个元素：

• socres：维度为n_clusterxn_genes的矩阵，包含每个基因和每个cluster的拟合优度，如上所述。

• labels：具有足够好的拟合优度得分的所有基因的标签。

• Score_cutoff：socres上用于确定是否分配标签的阈值

分配cluster标签：我们可以根据哪个cluster给出最低分数将基因分配到cluster。默认情况下，splines_kmeans_score_and_label需要足够低的拟合优度分数才会自动分配。“足够低”的标准基于所有簇上所有基因的分数分布(即splines_kmeans_score_and_label返回的整个分数矩阵)。然后，仅当基因的最佳得分高于50%时，基因才会被分配到一个cluster，其余基因则将NA作为其分配(此选择可以通过ratio_genes_to_label更改)。

labels = scores_and_labels$labels

# Let's keep only the list of genes that have a label.
labels = unlist(labels[!is.na(labels)])

# And also keep track of all the genes in cluster 2.
genes_in_cluster2 = names(labels[labels==cluster_of_interest])

对数据运行Scores_and_labels后，我们现在有19772个基因被分配到cluster。

如果我们没有设置splines_kmeans_score_and_label过滤阈值，我们可以通过考虑每个cluster的拟合优度分数的分布来可视化此过滤过程的影响，基于它们最低分数，将每个基因分配给一个cluster。我们可视化了分配给每个cluster的基因的分数(图14)，并与我们的过滤后该分数必须达到的阈值进行比较。我们可以通过重新运行splines_kmeans_score_and_label并设置过滤percentage_genes_to_label=1实现(并且我们可以通过参数previous_scores给出之前计算的分数来加速计算)。

# Get the best score and best label for all of the genes
# This time without filtering labels
# We can give the previous calculated scores to `previous_scores` to save time
unfiltered_scores = splines_kmeans_score_and_label(
    moanin_model, kmeans_clusters,
    proportion_genes_to_label=1,
    previous_scores=scores_and_labels$scores)
best_score = rowMin(scores_and_labels$scores)
best_label = unfiltered_scores$labels

par(mfrow=c(3, 3))
n_clusters = dim(kmeans_clusters$centroids)[1]
for(cluster_id in 1:n_clusters){
    hist(best_score[best_label==cluster_id],
     breaks=(1:50/50), xlim=c(0, 1),
     col="black", main=paste("C", cluster_id, sep=""),
     xlab="score", ylab="Num. genes")
    abline(v=scores_and_labels$score_cutoff, col="red", lwd=3, lty=2)
}

图14，Distribution of goodness-of-fit scores for each cluster. The dashed red line indicates the scoring threshold: all genes with a score above this threshold will not be labeled.

对于某些基因，它们在所有cluster中的得分为1。此类基因不适合所有cluster，这强调了过滤不适合聚类基因的重要性。

我们还可以根据分配给每个cluster的基因数量来研究样条k-means模型提供的标签与我们的评分和标记步骤之间的差异(图15)。

图15，Number of genes assigned to each cluster based on kmeans criteria of nearest centroid (top) versus our goodness-of-fit filtering strategy (bottom)

在聚类中使用的原始5521个基因中，有4946个基因根据其拟合优度得分分配到新的聚类中。我们可以根据如下所示的混合对比矩阵来比较它们是否仍然分配到相同的cluster(图16)。k-manes聚类分配显示在行坐标，拟合优度分配显示在列坐标，每个单元格中的数字表示两个聚类交叉点中的基因数量。

genes_scored = names(labels)

# of those used in clustering, keep only those 
# with high goodness of fit
kmeans_labels = kmeans_clusters$clusters
kmeans_labels = kmeans_labels[genes_scored]

labels = labels[names(kmeans_labels)]

both_labels = data.frame("kmeans"=kmeans_labels, "scored"=labels)
# Create another cluster.
#both_labels[is.na(both_labels)] = 21
confusion_matrix = table(both_labels)
text = as.matrix(as.character(confusion_matrix))
dim(text) = dim(confusion_matrix)
aheatmap(confusion_matrix, treeheight=0, main="Confusion matrix", 
     border_color="black", txt=text, legend=FALSE, Rowv=NA,
     Colv=NA, sub="Goodness-of-fit assignments")

# recreate, so in same order, etc as before
labels = scores_and_labels$labels
labels = unlist(labels[!is.na(labels)])

图16，Confusion matrix between the two sets of labels

在我们之前绘制的四个基因中，其中两个在我们评分后仍然分配到聚类2。我们可以再次查看聚类2中的一些基因，现在只选择最佳评分基因(图17)。

# order them by their score in cluster
ord = order(scores_and_labels$scores[genes_in_cluster2,
                     cluster_of_interest])
cluster2Score = genes_in_cluster2[ord]
plot_splines_data(
    moanin_model, subset_data=cluster2Score[1:4], 
    centroid=kmeans_clusters$centroids[cluster_of_interest,], 
    colors=ann_colors$Group, smooth=TRUE, simpleY=FALSE,
    mar=c(1.5, 2.5, 2, 0.1))

图17，Top scoring genes in Cluster 2

我们看到，正如预期的那样，这些基因与cluster质心更好地匹配。

详细查看特定簇cluster

现在，让我们更详细地了解一些特定的cluster。Cluster6似乎特别有趣：它捕获了不同流感病毒处理和对照之间存在显著差异的基因，而对照组则表达相对稳定。

我们已经展示了如何绘制一些示例基因，但是，考虑到噪声量，很难理解单个基因，也很难根据一些基因得出结论。

热图可用于研究特定基因的表达模式。在这里，我们将并排绘制标准化基因表达模式和重新缩放基因表达模式的热图。

首先，我们根据拟合优度分配选择感兴趣的基因，即cluster6中的基因。

cluster_to_plot = 6
genes_to_plot =  names(labels[labels == cluster_to_plot])

然后，创建这些基因的热图(3746个基因，图18)。

layout(matrix(c(1,2),nrow=1), widths=c(1.5, 2))

# order the observations by Group, then time, then replicate
orderedMeta<-data.frame(colData(moanin_model))
ord = order(
  orderedMeta$Group,
  orderedMeta$Timepoint,
  orderedMeta$Replicate)
orderedMeta$Timepoint = as.factor(orderedMeta$Timepoint)

orderedMeta = orderedMeta[ord, ]

res = aheatmap(
    as.matrix(assay(moanin_model[genes_to_plot,ord])),
    Colv=NA,
    color="YlGnBu",
    annCol=orderedMeta[,c("Group", "Timepoint")],
    annLegend=FALSE,
    annColors=ann_colors,
    main=paste("Cluster", cluster_to_plot, "(raw)"),
    treeheight=0, legend=FALSE)
# Now use the results of the previous call to aheatmap to reorder the genes.
aheatmap(
    rescale_values(moanin_model[genes_to_plot,ord])[res$rowInd,],
    Colv=NA,
    Rowv=NA,
    annCol=orderedMeta[, c("Group", "Timepoint")],
    annLegend=TRUE,
    annColors=ann_colors,
    main=paste("Cluster", cluster_to_plot, "(rescaled)"),
    treeheight=0)

图18，Heatmap of genes in Cluster 6

这两个热图表明，聚类方法成功地对不同尺度的基因进行了聚类，但对不同处理具有相同的动态响应。

如何选择cluster的数量。

这一部分的内容其实非常多，我做了很多省略，感兴趣的可以私聊我，问得人多就单独出一期。 执行聚类时出现的一个常见问题是如何选择聚类数量。cluster数K的选择取决于分析目的。cluster的数量不应超过我们想要解释的基因集的数量。首先我们假设这个数字是20。一旦设置了最大cluster数，就有几种策略可以识别cluster的数量：

• Elbow法。 Elbow法于1953年由Thorndike首次提出，该方法利用cluster内总平方和(WCSS)决定cluster数。当增加cluster没有将WCSS降低足够的量时，可以考虑增加。因此，该方法为用户选择cluster的数量提供了可视化的图形，但通常在实际数据上很难看到“肘部”，因此cluster数没有明确的数值。

• Silhouette法。与Elbow法类似，Silhouette法是指验证cluster内一致性的方法，并提供可视化来选择cluster的数量。

• 稳定性法。稳定性法方法比任何其他方法的计算量更大，因为它们依赖于评估每个k在随机抽样数据的聚类稳定性。然后，让用户根据相似性度量来选择聚类的数量。

首先，我们对所有感兴趣的k值运行聚类。对于每个聚类，我们将保留(1)聚类内平方(WCSS)；(2)每个基因的聚类分配（或标签）。

下面我们对k等于2–20进行聚类。splines_kmeans函数返回每个聚类的WCSS，我们将其求和以获得总WCSS。

all_possible_n_clusters = c(2:20)
all_clustering = list()
wss_values = list()

i = 1
for(n_cluster in all_possible_n_clusters){
    clustering_results = splines_kmeans(de_moanin_model,
    n_clusters=n_cluster, random_seed=42,
    n_init=10)
    wss_values[i] = sum(clustering_results$WCSS_per_cluster)
    all_clustering[[i]] = clustering_results$clusters
    i = i + 1
}

Elbow法

选择聚类数量的Elbow法依赖于可视化来选择聚类数量。该方法依赖于绘制cluster内平方和(WCSS)。在某个K值之后，WCSS将开始更缓慢地减小，从而在图表中给出一个角度或“肘部”。cluster的数量是在这个“肘点”处选择的。

我们在此绘制k=2–20时的WCSS(图19)。我们看到，正如预期的那样，WCSS继续下降，但除了k值非常小(3-4 个cluster)之外，下降幅度没有明显下降。然而，考虑到复杂性，只有3-4个cluster聚类又似乎太少。

plot(all_possible_n_clusters, wss_values,
     type="b", pch=19, frame=FALSE, 
     xlab="Number of clusters K",
     ylab="Total within-clusters sum of squares")

图19，Plot of WCSS as a function of k

平均Silhouette法

Silhouette value轮廓值是衡量数据点与其自身聚类(凝聚力)与其他聚类(分离度)相比的相似程度的指标，如图20所示。

# function to compute average silhouette for k clusters
average_silhouette = function(labels, y) {
    silhouette_results = silhouette(unlist(labels[1]),dist(y))
    return(mean(silhouette_results[, 3]))
}

# extract the average silhouette
average_silhouette_values = list()
i = 1
for(i in 1:length(all_clustering)){
    clustering_results = all_clustering[i]
    average_silhouette_values[i] = average_silhouette(clustering_results,
        assay(de_moanin_model))
    i = i + 1
}

plot(all_possible_n_clusters, average_silhouette_values,
     type="b", pch=19, frame=FALSE,
     xlab="Number of clusters K",
     ylab="Average Silhouettes")

图20，Plot of silhouette values as a function of k

考察聚类的稳定性

在真实数据上，聚类的数量不仅是未知的，而且是模糊的：它将取决于用户所需的聚类分辨率。然而，就生物数据而言，结果的稳定性和再现性对于确保当数据或模型受到合理扰动时结果的生物学解释成立是必要的。

依赖聚类结果的稳定性来选择k的方法可确保聚类的生物学解释在数据受到扰动后仍然成立。此外，使用明确定义的k生成数据的模拟表明，对于正确数量的cluster，聚类更加稳定。

大多数通过稳定性度量来查找聚类数量的方法仅提供可视化辅助来帮助选择。通常可视化的一个方法是共识矩阵consensus matrix：共识矩阵是一个n×n矩阵，存储两个项目聚集在一起的聚类比例。一个完美的共识矩阵，其排序使得属于同一cluster的所有元素彼此相邻，沿对角线显示的块接近1。

要执行此类分析，第一步是要获得resampled的数据集，有两种方法：bootstrap和subsampling。

Bootstrap法：为了获得resampled的数据，我们通过随机替换获得与原始聚类大小相同的样本(即5521个基因)：

n_genes = dim(de_moanin_model)[1]
indices = sample(1:dim(de_moanin_model)[1], n_genes, replace=TRUE)

bootstrapped_y = de_moanin_model[indices, ]

第二种方法是使用二次采样策略subsampling，我们获取基因的独特子集，保留80%的基因：

subsample_proportion = 0.8
indices = sample(1:dim(de_moanin_model)[1],
         n_genes * subsample_proportion,
         replace=FALSE)

subsampled_y = de_moanin_model[indices, ]

我们对差异表达且对数倍数变化高于2(使用lfc_max方法计算)的所有基因运行bootstrap方法，并对每个k=2–10执行B = 20次重复。我们显示下面的代码，但由于计算所花费的时间很长，我们单独计算了这些值，并的得到每一个k的结果。

# You may want to set the random seed of the experiment so that the results
# don't vary if you rerun the experiment.
# set.seed(random_seed)
n_genes = dim(de_moanin_model)[1] * subsample_proportion
indices = sample(1:dim(de_moanin_model)[1], n_genes, replace=TRUE)

kmeans_clusters = splines_kmeans(de_moanin_model[indices,],
    n_clusters=n_clusters,
    random_seed=42,
    n_init=20)

# Perform prediction on the whole set of data.
kmeans_clusters = splines_kmeans_predict(de_moanin_model, kmeans_clusters, 
    method="distance")

现在，我们引入k=5和k=20的结果。

stability_5 = read.table("results/stability_5.tsv", sep="\t")
stability_20 = read.table("results/stability_20.tsv", sep="\t")

每列对应一个bootstrap样本，每行对应一个基因，每个条目对应为该特定聚类找到的cluster标签。因此，对于每个基因，我们在25次重采样运行中分配到一个cluster。

现在，我们使用moanin包中的函数consensus_matrix来计算每对样本在25次重采样运行中聚集在一起的次数比例，并绘制k = 5和k = 20时这些矩阵的热图。

par(mfrow=c(1, 2))
consensus_matrix_stability_5 = consensus_matrix(stability_5,
                        scale=FALSE)
aheatmap(consensus_matrix_stability_5[1:1000, 1:1000], Rowv=FALSE,
     Colv=FALSE,
     treeheight=0, main="K=5")

consensus_matrix_stability_20 = consensus_matrix(stability_20,
                         scale=FALSE)
aheatmap(consensus_matrix_stability_20[1:1000, 1:1000], Rowv=FALSE,
     Colv=FALSE,
     treeheight=0, main="K=20")

图21，Consensus matrix of proportion of times samples were in the same cluster over bootstrap samples for k=5 and 20.

我们可以看到(图21)，在重采样运行中，K = 5的选择似乎比K = 20的选择更稳定。

Cluster的下游分析。

一旦获得了好的聚类，下一步就是利用聚类来解释定义，可以对每个聚类定义的基因集进行经典的富集分析，例如GO term富集分析和基序motif富集分析(KEGG本来也想做，但是网站太不稳定了，舍弃）。

首先，过滤需要使用的基因，只选择我们要在富集分析中使用的基因。可以考虑(1)不进行任何进一步过滤的聚类结果，(2)仅考虑一组每个cluster中差异表达的基因，或(3)合适的cluster基因子集(基于某些标准)。

使用biomaRt寻找GO terms富集通路

Gene Ontology(GO)数据库将基因分类为有意义的biological ontologies，并可通过topGO包进行富集分析。我们使用biomaRt来查找和基因集匹配的GO terms。

## Set up right ensembl database
ensembl = useMart("ensembl")
ensembl = useDataset("mmusculus_gene_ensembl", mart=ensembl)

## Get gene names of genes in Cluster 8
cluster = 8
gene_names = names(labels)
genes = gene_names[labels == cluster]

genes = getBM(attributes=c("go_id", "refseq_mrna"),
          values=gene_names,
          filters="refseq_mrna",
          mart=ensembl)

biomaRtquery结果是一个包含两列的矩阵：基因名称和GO terms ID。例如：

$NM_199153
[1] "GO:0016020" "GO:0016021" "GO:0007186" "GO:0004930" "GO:0007165"
[6] "GO:0050896" "GO:0050909"

$NM_201361
 [1] "GO:0016020" "GO:0016021" "GO:0003674" "GO:0008150" "GO:0005794"
 [6] "GO:0005829" "GO:0005737" "GO:0005856" "GO:0005874" "GO:0005739"
[11] "GO:0005819" "GO:0000922" "GO:0072686"

moanin提供了一个简单的函数(create_go_term_mapping)来进行此转换：

# Create gene to GO id mapping
gene_id_go_mapping = create_go_term_mapping(genes)

一旦创建了基因ID到GO映射列表，moanin就会为topGO提供一个接口来确定富集的GO terms。它还会执行p值校正，并且以data.frame形式返回显著的GO terms富集。在这里展示了一个在“Biological process”上运行GO terms富集的例子。

# Create logical vector of whether in cluster 6
assignments = labels == cluster

go_terms_enriched = find_enriched_go_terms(
    assignments,
    gene_id_go_mapping, ontology="BP")

结语

此分析流程提供了使用包moanin分析多时间点基因表达数据的实战教程。分析流程由三个主要步骤组成，需要在质量控制和标准化之后执行：(1)差异表达分析；(2)时序基因表达数据的聚类；(3)聚类的下游分析。