接下來我們繼續(xù)多時間點樣本實戰(zhàn)分析流程的第二部分：聚類和富集分析摧阅。第一部分的完整流程請參照：RNA-Seq 分析流程：多時間點樣本分析實戰(zhàn)（一）

多時間點數(shù)據(jù)的聚類

前面我們發(fā)現(xiàn)70% 的基因是差異表達的，幾乎所有通路都受到處理的影響珠洗。因此诺核，分析流程的下一步是根據(jù)基因表達對處理的動態(tài)反應進行聚類抄肖。

過濾

在對基因進行聚類之前，我們首先將感興趣的基因集濃縮為(1)顯著差異表達的基因窖杀；(2) 處理之間存在較大倍數(shù)的基因漓摩。減少執(zhí)行聚類的基因集可以增強聚類的可信度。為此入客，我們首先使用Fisher方法將時序差異表達步驟中獲得的所有p值聚合為單個p值(pvalues_fisher_method)幌甘。接下來，我們選擇Fisher方法調整的p值低于0.05且至少在一種處理和一個時間點的對數(shù)倍數(shù)變化至少為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個基因。然后我們就可以開始聚類的常規(guī)分析流程鞍泉。通過觀察差異表達的基因發(fā)現(xiàn)皱埠，許多基因具有相似的基因表達模式，但上下調的幅度不同咖驮。

因此边器，我們建議對k-means進行以下調整：
1.樣條函數(shù)擬合：對于每個基因，計算(如moanin對象中包含的那樣)擬合樣條函數(shù)托修。

2.重新歸一化樣條函數(shù)：對于每個基因忘巧，重新歸一化估計的樣條函數(shù)，使得值限制在0和1之間睦刃，從而在基因之間具有可比性砚嘴。

3.K-means：對樣條函數(shù)重新調整后的擬合值應用k-means來估計聚類。

樣條函數(shù)的重新歸一化旨在使每個基因達到可比較的尺度，類似于通常在沒有時間維度的基因表達研究中對基因表達進行的歸一化和縮放际长。

聚類步驟由moanin中的splines_kmeans函數(shù)執(zhí)行耸采。我們先將cluster數(shù)設置為8，不過后面我們將討論選擇最佳cluster數(shù)的問題工育。

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

splines_kmeans函數(shù)返回一個命名列表虾宇，包含：

• centroids：包含cluster質心的矩陣。矩陣的維度為(n_centroid, n_samples)如绸。

• cluster：大小為n_genes的向量嘱朽，包含每個基因的kmeans聚類被分配到哪一個cluster的信息。

然后怔接，我們使用plot_splines_data函數(shù)來可視化通過樣條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的擬合程度進行評分。

此外趁蕊，前面我們只計算了根據(jù)過濾得到的基因子集坞生，我們希望有一種方法將所有基因分配到各個cluster里面。

因此掷伙，聚類分析的下一步部分是評分scoring和標記label是己。每個基因都會獲得一個分數(shù)，該分數(shù)對應于每個基因和每個cluster之間的擬合優(yōu)度任柜。評分和標簽是通過splines_kmeans_score_and_label函數(shù)實現(xiàn)卒废。該函數(shù)計算基因與聚類質心的擬合度，并在得分足夠高的情況下為基因提供聚類標簽宙地。

# 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的擬合優(yōu)度，如上所述宅粥。

• labels：具有足夠好的擬合優(yōu)度得分的所有基因的標簽参袱。

• Score_cutoff：socres上用于確定是否分配標簽的閾值

分配cluster標簽：我們可以根據(jù)哪個cluster給出最低分數(shù)將基因分配到cluster。默認情況下，splines_kmeans_score_and_label需要足夠低的擬合優(yōu)度分數(shù)才會自動分配蓖柔〕狡螅“足夠低”的標準基于所有簇上所有基因的分數(shù)分布(即splines_kmeans_score_and_label返回的整個分數(shù)矩陣)。然后况鸣，僅當基因的最佳得分高于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])

對數(shù)據(jù)運行Scores_and_labels后潜索，我們現(xiàn)在有19772個基因被分配到cluster。

如果我們沒有設置splines_kmeans_score_and_label過濾閾值懂酱，我們可以通過考慮每個cluster的擬合優(yōu)度分數(shù)的分布來可視化此過濾過程的影響竹习，基于它們最低分數(shù)，將每個基因分配給一個cluster列牺。我們可視化了分配給每個cluster的基因的分數(shù)(圖14)整陌，并與我們的過濾后該分數(shù)必須達到的閾值進行比較。我們可以通過重新運行splines_kmeans_score_and_label并設置過濾percentage_genes_to_label=1實現(xiàn)(并且我們可以通過參數(shù)previous_scores給出之前計算的分數(shù)來加速計算)瞎领。

# 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震放，這強調了過濾不適合聚類基因的重要性。

我們還可以根據(jù)分配給每個cluster的基因數(shù)量來研究樣條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個基因根據(jù)其擬合優(yōu)度得分分配到新的聚類中乙各。我們可以根據(jù)如下所示的混合對比矩陣來比較它們是否仍然分配到相同的cluster(圖16)墨礁。k-manes聚類分配顯示在行坐標，擬合優(yōu)度分配顯示在列坐標耳峦，每個單元格中的數(shù)字表示兩個聚類交叉點中的基因數(shù)量恩静。

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中的一些基因，現(xiàn)在只選擇最佳評分基因(圖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

現(xiàn)在罕扎，讓我們更詳細地了解一些特定的cluster。Cluster6似乎特別有趣：它捕獲了不同流感病毒處理和對照之間存在顯著差異的基因丐重，而對照組則表達相對穩(wěn)定腔召。

我們已經(jīng)展示了如何繪制一些示例基因，但是扮惦，考慮到噪聲量臀蛛，很難理解單個基因，也很難根據(jù)一些基因得出結論崖蜜。

熱圖可用于研究特定基因的表達模式浊仆。在這里，我們將并排繪制標準化基因表達模式和重新縮放基因表達模式的熱圖豫领。

首先抡柿，我們根據(jù)擬合優(yōu)度分配選擇感興趣的基因，即cluster6中的基因等恐。

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

然后乒疏，創(chuàng)建這些基因的熱圖(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

這兩個熱圖表明斧账，聚類方法成功地對不同尺度的基因進行了聚類跋炕，但對不同處理具有相同的動態(tài)響應剃浇。

如何選擇cluster的數(shù)量耘柱。

這一部分的內容其實非常多虚循，我做了很多省略同欠，感興趣的可以私聊我，問得人多就單獨出一期横缔。 執(zhí)行聚類時出現(xiàn)的一個常見問題是如何選擇聚類數(shù)量铺遂。cluster數(shù)K的選擇取決于分析目的。cluster的數(shù)量不應超過我們想要解釋的基因集的數(shù)量茎刚。首先我們假設這個數(shù)字是20襟锐。一旦設置了最大cluster數(shù)，就有幾種策略可以識別cluster的數(shù)量：

• Elbow法膛锭。 Elbow法于1953年由Thorndike首次提出粮坞，該方法利用cluster內總平方和(WCSS)決定cluster數(shù)。當增加cluster沒有將WCSS降低足夠的量時初狰，可以考慮增加莫杈。因此，該方法為用戶選擇cluster的數(shù)量提供了可視化的圖形奢入，但通常在實際數(shù)據(jù)上很難看到“肘部”筝闹，因此cluster數(shù)沒有明確的數(shù)值。

• Silhouette法。與Elbow法類似关顷，Silhouette法是指驗證cluster內一致性的方法糊秆，并提供可視化來選擇cluster的數(shù)量。

• 穩(wěn)定性法议双。穩(wěn)定性法方法比任何其他方法的計算量更大痘番，因為它們依賴于評估每個k在隨機抽樣數(shù)據(jù)的聚類穩(wěn)定性。然后平痰，讓用戶根據(jù)相似性度量來選擇聚類的數(shù)量夫偶。

首先，我們對所有感興趣的k值運行聚類觉增。對于每個聚類兵拢，我們將保留(1)聚類內平方(WCSS)；(2)每個基因的聚類分配（或標簽）逾礁。

下面我們對k等于2–20進行聚類说铃。splines_kmeans函數(shù)返回每個聚類的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法

選擇聚類數(shù)量的Elbow法依賴于可視化來選擇聚類數(shù)量腻扇。該方法依賴于繪制cluster內平方和(WCSS)。在某個K值之后砾嫉，WCSS將開始更緩慢地減小幼苛，從而在圖表中給出一個角度或“肘部”。cluster的數(shù)量是在這個“肘點”處選擇的焕刮。

我們在此繪制k=2–20時的WCSS(圖19)舶沿。我們看到，正如預期的那樣配并，WCSS繼續(xù)下降括荡，但除了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輪廓值是衡量數(shù)據(jù)點與其自身聚類(凝聚力)與其他聚類(分離度)相比的相似程度的指標，如圖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

考察聚類的穩(wěn)定性

在真實數(shù)據(jù)上，聚類的數(shù)量不僅是未知的婶溯，而且是模糊的：它將取決于用戶所需的聚類分辨率鲸阔。然而偷霉，就生物數(shù)據(jù)而言，結果的穩(wěn)定性和再現(xiàn)性對于確保當數(shù)據(jù)或模型受到合理擾動時結果的生物學解釋成立是必要的褐筛。

依賴聚類結果的穩(wěn)定性來選擇k的方法可確保聚類的生物學解釋在數(shù)據(jù)受到擾動后仍然成立类少。此外，使用明確定義的k生成數(shù)據(jù)的模擬表明渔扎，對于正確數(shù)量的cluster硫狞，聚類更加穩(wěn)定。

大多數(shù)通過穩(wěn)定性度量來查找聚類數(shù)量的方法僅提供可視化輔助來幫助選擇晃痴。通巢蟹裕可視化的一個方法是共識矩陣consensus matrix：共識矩陣是一個n×n矩陣，存儲兩個項目聚集在一起的聚類比例倘核。一個完美的共識矩陣泣侮，其排序使得屬于同一cluster的所有元素彼此相鄰，沿對角線顯示的塊接近1紧唱。

要執(zhí)行此類分析活尊，第一步是要獲得resampled的數(shù)據(jù)集，有兩種方法：bootstrap和subsampling漏益。

Bootstrap法：為了獲得resampled的數(shù)據(jù)蛹锰，我們通過隨機替換獲得與原始聚類大小相同的樣本(即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, ]

我們對差異表達且對數(shù)倍數(shù)變化高于2(使用lfc_max方法計算)的所有基因運行bootstrap方法铜犬，并對每個k=2–10執(zhí)行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")

現(xiàn)在煎谍，我們引入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唆垃。

現(xiàn)在五芝，我們使用moanin包中的函數(shù)consensus_matrix來計算每對樣本在25次重采樣運行中聚集在一起的次數(shù)比例，并繪制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的選擇更穩(wěn)定醉途。

Cluster的下游分析矾瑰。

一旦獲得了好的聚類，下一步就是利用聚類來解釋定義隘擎，可以對每個聚類定義的基因集進行經(jīng)典的富集分析殴穴，例如GO term富集分析和基序motif富集分析(KEGG本來也想做，但是網(wǎng)站太不穩(wěn)定了货葬，舍棄）采幌。

首先，過濾需要使用的基因震桶，只選擇我們要在富集分析中使用的基因休傍。可以考慮(1)不進行任何進一步過濾的聚類結果蹲姐，(2)僅考慮一組每個cluster中差異表達的基因尊残，或(3)合適的cluster基因子集(基于某些標準)。

使用biomaRt尋找GO terms富集通路

Gene Ontology(GO)數(shù)據(jù)庫將基因分類為有意義的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提供了一個簡單的函數(shù)(create_go_term_mapping)來進行此轉換：

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

一旦創(chuàng)建了基因ID到GO映射列表慰毅，moanin就會為topGO提供一個接口來確定富集的GO terms。它還會執(zhí)行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分析多時間點基因表達數(shù)據(jù)的實戰(zhàn)教程东臀。分析流程由三個主要步驟組成着饥，需要在質量控制和標準化之后執(zhí)行：(1)差異表達分析；(2)時序基因表達數(shù)據(jù)的聚類惰赋；(3)聚類的下游分析宰掉。