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

多時(shí)間點(diǎn)數(shù)據(jù)的聚類

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

過濾

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

因此，我們建議對(duì)k-means進(jìn)行以下調(diào)整：
1.樣條函數(shù)擬合：對(duì)于每個(gè)基因，計(jì)算(如moanin對(duì)象中包含的那樣)擬合樣條函數(shù)。

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

3.K-means：對(duì)樣條函數(shù)重新調(diào)整后的擬合值應(yīng)用k-means來估計(jì)聚類。

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

聚類步驟由moanin中的splines_kmeans函數(shù)執(zhí)行。我們先將cluster數(shù)設(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ù)返回一個(gè)命名列表，包含：

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

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

然后，我們使用plot_splines_data函數(shù)來可視化通過樣條k-means模型獲得的每個(gè)cluster的質(zhì)心(圖12)。

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

圖12，K-means centroids

由于我們對(duì)樣條擬合進(jìn)行了重新調(diào)整，這些質(zhì)心的范圍為 0-1，并且并不代表實(shí)際的基因表達(dá)水平。在圖13中，我們繪制了一些聚類到cluster2的實(shí)際基因：

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質(zhì)心的模式有一些相似性，但從匹配質(zhì)心的意義上來說，這些特定基因和cluster并沒有達(dá)到最佳擬合。由于基因與cluster的擬合程度存在差異，我們希望能夠?qū)蚺ccluster的擬合程度進(jìn)行評(píng)分。

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

因此，聚類分析的下一步部分是評(píng)分scoring和標(biāo)記label。每個(gè)基因都會(huì)獲得一個(gè)分?jǐn)?shù)，該分?jǐn)?shù)對(duì)應(yīng)于每個(gè)基因和每個(gè)cluster之間的擬合優(yōu)度。評(píng)分和標(biāo)簽是通過splines_kmeans_score_and_label函數(shù)實(shí)現(xiàn)。該函數(shù)計(jì)算基因與聚類質(zhì)心的擬合度，并在得分足夠高的情況下為基因提供聚類標(biāo)簽。

# 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列表包含三個(gè)元素：

• socres：維度為n_clusterxn_genes的矩陣，包含每個(gè)基因和每個(gè)cluster的擬合優(yōu)度，如上所述。

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

• Score_cutoff：socres上用于確定是否分配標(biāo)簽的閾值

分配cluster標(biāo)簽：我們可以根據(jù)哪個(gè)cluster給出最低分?jǐn)?shù)將基因分配到cluster。默認(rèn)情況下，splines_kmeans_score_and_label需要足夠低的擬合優(yōu)度分?jǐn)?shù)才會(huì)自動(dòng)分配?！白銐虻汀钡臉?biāo)準(zhǔn)基于所有簇上所有基因的分?jǐn)?shù)分布(即splines_kmeans_score_and_label返回的整個(gè)分?jǐn)?shù)矩陣)。然后，僅當(dāng)基因的最佳得分高于50%時(shí)，基因才會(huì)被分配到一個(gè)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])

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

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

# 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.

對(duì)于某些基因，它們?cè)谒衏luster中的得分為1。此類基因不適合所有cluster，這強(qiáng)調(diào)了過濾不適合聚類基因的重要性。

我們還可以根據(jù)分配給每個(gè)cluster的基因數(shù)量來研究樣條k-means模型提供的標(biāo)簽與我們的評(píng)分和標(biāo)記步驟之間的差異(圖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個(gè)基因中，有4946個(gè)基因根據(jù)其擬合優(yōu)度得分分配到新的聚類中。我們可以根據(jù)如下所示的混合對(duì)比矩陣來比較它們是否仍然分配到相同的cluster(圖16)。k-manes聚類分配顯示在行坐標(biāo)，擬合優(yōu)度分配顯示在列坐標(biāo)，每個(gè)單元格中的數(shù)字表示兩個(gè)聚類交叉點(diǎn)中的基因數(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

在我們之前繪制的四個(gè)基因中，其中兩個(gè)在我們?cè)u(píng)分后仍然分配到聚類2。我們可以再次查看聚類2中的一些基因，現(xiàn)在只選擇最佳評(píng)分基因(圖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

我們看到，正如預(yù)期的那樣，這些基因與cluster質(zhì)心更好地匹配。

詳細(xì)查看特定簇cluster

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

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

熱圖可用于研究特定基因的表達(dá)模式。在這里，我們將并排繪制標(biāo)準(zhǔn)化基因表達(dá)模式和重新縮放基因表達(dá)模式的熱圖。

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

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

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

這兩個(gè)熱圖表明，聚類方法成功地對(duì)不同尺度的基因進(jìn)行了聚類，但對(duì)不同處理具有相同的動(dòng)態(tài)響應(yīng)。

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

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

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

• Silhouette法。與Elbow法類似，Silhouette法是指驗(yàn)證cluster內(nèi)一致性的方法，并提供可視化來選擇cluster的數(shù)量。

• 穩(wěn)定性法。穩(wěn)定性法方法比任何其他方法的計(jì)算量更大，因?yàn)樗鼈円蕾囉谠u(píng)估每個(gè)k在隨機(jī)抽樣數(shù)據(jù)的聚類穩(wěn)定性。然后，讓用戶根據(jù)相似性度量來選擇聚類的數(shù)量。

首先，我們對(duì)所有感興趣的k值運(yùn)行聚類。對(duì)于每個(gè)聚類，我們將保留(1)聚類內(nèi)平方(WCSS)；(2)每個(gè)基因的聚類分配（或標(biāo)簽）。

下面我們對(duì)k等于2–20進(jìn)行聚類。splines_kmeans函數(shù)返回每個(gè)聚類的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內(nèi)平方和(WCSS)。在某個(gè)K值之后，WCSS將開始更緩慢地減小，從而在圖表中給出一個(gè)角度或“肘部”。cluster的數(shù)量是在這個(gè)“肘點(diǎn)”處選擇的。

我們?cè)诖死L制k=2–20時(shí)的WCSS(圖19)。我們看到，正如預(yù)期的那樣，WCSS繼續(xù)下降，但除了k值非常小(3-4 個(gè)cluster)之外，下降幅度沒有明顯下降。然而，考慮到復(fù)雜性，只有3-4個(gè)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ù)點(diǎn)與其自身聚類(凝聚力)與其他聚類(分離度)相比的相似程度的指標(biāo)，如圖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í)數(shù)據(jù)上，聚類的數(shù)量不僅是未知的，而且是模糊的：它將取決于用戶所需的聚類分辨率。然而，就生物數(shù)據(jù)而言，結(jié)果的穩(wěn)定性和再現(xiàn)性對(duì)于確保當(dāng)數(shù)據(jù)或模型受到合理擾動(dòng)時(shí)結(jié)果的生物學(xué)解釋成立是必要的。

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

大多數(shù)通過穩(wěn)定性度量來查找聚類數(shù)量的方法僅提供可視化輔助來幫助選擇。通?？梢暬囊粋€(gè)方法是共識(shí)矩陣consensus matrix：共識(shí)矩陣是一個(gè)n×n矩陣，存儲(chǔ)兩個(gè)項(xiàng)目聚集在一起的聚類比例。一個(gè)完美的共識(shí)矩陣，其排序使得屬于同一cluster的所有元素彼此相鄰，沿對(duì)角線顯示的塊接近1。

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

Bootstrap法：為了獲得resampled的數(shù)據(jù)，我們通過隨機(jī)替換獲得與原始聚類大小相同的樣本(即5521個(gè)基因)：

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，我們獲取基因的獨(dú)特子集，保留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, ]

我們對(duì)差異表達(dá)且對(duì)數(shù)倍數(shù)變化高于2(使用lfc_max方法計(jì)算)的所有基因運(yùn)行bootstrap方法，并對(duì)每個(gè)k=2–10執(zhí)行B = 20次重復(fù)。我們顯示下面的代碼，但由于計(jì)算所花費(fèi)的時(shí)間很長，我們單獨(dú)計(jì)算了這些值，并的得到每一個(gè)k的結(jié)果。

# 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的結(jié)果。

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

每列對(duì)應(yīng)一個(gè)bootstrap樣本，每行對(duì)應(yīng)一個(gè)基因，每個(gè)條目對(duì)應(yīng)為該特定聚類找到的cluster標(biāo)簽。因此，對(duì)于每個(gè)基因，我們?cè)?5次重采樣運(yùn)行中分配到一個(gè)cluster。

現(xiàn)在，我們使用moanin包中的函數(shù)consensus_matrix來計(jì)算每對(duì)樣本在25次重采樣運(yùn)行中聚集在一起的次數(shù)比例，并繪制k = 5和k = 20時(shí)這些矩陣的熱圖。

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)，在重采樣運(yùn)行中，K = 5的選擇似乎比K = 20的選擇更穩(wěn)定。

Cluster的下游分析。

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

首先，過濾需要使用的基因，只選擇我們要在富集分析中使用的基因?？梢钥紤](1)不進(jìn)行任何進(jìn)一步過濾的聚類結(jié)果，(2)僅考慮一組每個(gè)cluster中差異表達(dá)的基因，或(3)合適的cluster基因子集(基于某些標(biāo)準(zhǔn))。

使用biomaRt尋找GO terms富集通路

Gene Ontology(GO)數(shù)據(jù)庫將基因分類為有意義的biological ontologies，并可通過topGO包進(jìn)行富集分析。我們使用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結(jié)果是一個(gè)包含兩列的矩陣：基因名稱和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提供了一個(gè)簡(jiǎn)單的函數(shù)(create_go_term_mapping)來進(jìn)行此轉(zhuǎn)換：

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

一旦創(chuàng)建了基因ID到GO映射列表，moanin就會(huì)為topGO提供一個(gè)接口來確定富集的GO terms。它還會(huì)執(zhí)行p值校正，并且以data.frame形式返回顯著的GO terms富集。在這里展示了一個(gè)在“Biological process”上運(yùn)行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")

結(jié)語

此分析流程提供了使用包moanin分析多時(shí)間點(diǎn)基因表達(dá)數(shù)據(jù)的實(shí)戰(zhàn)教程。分析流程由三個(gè)主要步驟組成，需要在質(zhì)量控制和標(biāo)準(zhǔn)化之后執(zhí)行：(1)差異表達(dá)分析；(2)時(shí)序基因表達(dá)數(shù)據(jù)的聚類；(3)聚類的下游分析。