老規(guī)矩，先奉上學(xué)習(xí)資料鏈接：
這一次學(xué)習(xí)Mutational signatures的構(gòu)建。
Mutational signatures被認(rèn)為代表了突變過(guò)程，并以具有特定序列上下文的突變類(lèi)型的特定貢獻(xiàn)為特征。
NMF：使用非負(fù)矩陣分解(NMF)，可以從突變計(jì)數(shù)矩陣中重新提取Mutational signatures。
signature refitting：還可以確定突變計(jì)數(shù)矩陣對(duì)先前定義的Mutational signatures的暴露（exposure）。
NMF對(duì)大樣本最有用，而signature refitting可用于單個(gè)樣本。
接下來(lái)首先討論NMF，然后再討論signature refitting，最后，我們將討論如何直接分析mutational profile和signatures之間的相似性。
NMF提取De novo mutational signature
1）NMF
NMF中的一個(gè)關(guān)鍵參數(shù)是分解秩(the factorization rank)，它是提取的mutational signature的數(shù)量?？梢允褂肗MF包確定最優(yōu)分解秩。
先跟上一個(gè)教程一樣，處理vcf文件得到mutaional矩陣：
ref_genome <- "BSgenome.Hsapiens.UCSC.hg19"
library(ref_genome, character.only = TRUE)
vcf_files <- list.files(system.file("extdata", package = "MutationalPatterns"), pattern = "sample.vcf", full.names = TRUE)
sample_names <- c(
"colon1", "colon2", "colon3",
"intestine1", "intestine2", "intestine3",
"liver1", "liver2", "liver3")
grl <- read_vcfs_as_granges(vcf_files, sample_names, ref_genome)
mut_mat <- mut_matrix(vcf_list = grl, ref_genome = ref_genome)
head(mut_mat)
一般來(lái)說(shuō)，數(shù)據(jù)集越大，rank會(huì)更高。這里我們展示SNV的NMF。其他突變類(lèi)型執(zhí)行NMF的工作方式與此相同。
使用NMF包生成估計(jì)秩圖，需要很長(zhǎng)時(shí)間，最好是最后寫(xiě)成傳參腳本后臺(tái)運(yùn)行。
# First add a small pseudocount to your mutation count matrix:
mut_mat <- mut_mat + 0.0001
library("NMF")
estimate <- nmf(mut_mat, rank = 2:5, method = "brunet", nrun = 10, seed = 123456, .opt = "v-p")
# And plot it:
p <- plot(estimate)
rank plot圖：
使用extract_signatures從突變計(jì)數(shù)矩陣中提取mutational signatures。在本例中，rank 為2，提取2個(gè)mutational signatures(對(duì)于較大的數(shù)據(jù)集，明智的做法是通過(guò)改變nrun參數(shù)來(lái)執(zhí)行更多的迭代，以獲得穩(wěn)定性并避免不良的局部最小值)。
上圖繪制了rank2-5的各種評(píng)價(jià)指標(biāo)，所以到底是因?yàn)槭裁粗颠x取的rank 2？？？
想起來(lái)前面還有段話(huà)：
The most common approach is to choose the smallest rank for which cophenetic correlation coefficient starts decreasing.
也就是上圖左上角的那個(gè)cophenetic，rank2-5的時(shí)候一直往下降，最小的rank就是2。
確定好rank之后，提取mutational signatures：
nmf_res <- extract_signatures(mut_mat, rank = 2, nrun = 10, single_core = TRUE)
這樣得到了兩個(gè)mutational signatures：
2）Bayesian NMF
也可以使用變分貝葉斯NMF（variational bayes NMF）。這樣可以更容易地確定正確的rank。需要安裝ccfindR包，然后可以確定最優(yōu)signatures數(shù)，這同樣需要很長(zhǎng)時(shí)間。
# BiocManager::install("ccfindR")
library(ccfindR)
sc <- scNMFSet(count = mut_mat)
set.seed(4)
estimate_bayes <- vb_factorize(sc, ranks = 1:3, nrun = 1, progress.bar = FALSE, verbose = 0)
png(filename = paste0(opt$od, "/NMF_estimate_bayes.png"),width = 1000, height = 700, res=200)
plot(estimate_bayes)
dev.off()
結(jié)果圖如下：
然后提取signatures：
nmf_res_bayes <- extract_signatures(mut_mat, rank = 2, nrun = 10, nmf_type = "variational_bayes")
head(nmf_res_bayes)
3）Changing the names of the extracted signatures
我們可以對(duì)提取出來(lái)的signatures的名字更改為習(xí)慣用名：
## 3）Changing the names of the extracted signatures
colnames(nmf_res$signatures) <- c("Signature A", "Signature B")
rownames(nmf_res$contribution) <- c("Signature A", "Signature B")
head(nmf_res[["signatures"]])
Signature A Signature B
A[C>A]A 82.30033 101.35159
A[C>A]C 41.87138 36.46964
A[C>A]G 40.09448 21.64503
A[C>A]T 59.09189 42.90475
C[C>A]A 76.61626 55.09600
C[C>A]C 77.84636 24.82175
此外，NMF 提取的一些signatures可能與已知的signatures非常相似，因此，將NMF signatures的名稱(chēng)更改為這些已知signatures可能會(huì)很有用。這通常使你更容易解釋你的結(jié)果。
要做到這一點(diǎn)，首先需要讀取一些已經(jīng)存在的signatures。在這里，我們將使用COSMIC (v3.2) (Alexandrov et al. 2020)的signatures。(稍后我們將討論如何使用其他signatures矩陣。)
現(xiàn)在可以更改NMF提取的signatures的名稱(chēng)。在本例中，如果signatures的名稱(chēng)與現(xiàn)有COSMIC signatures的余弦相似度大于0.85，則更改該signatures的名稱(chēng)。
signatures = get_known_signatures()
nmf_res <- rename_nmf_signatures(nmf_res, signatures, cutoff = 0.85)
colnames(nmf_res$signatures)
[1] "SBS5-like" "SBS1-like"
我們現(xiàn)在看到，我們提取的signatures與COSMIC signatures SBS1和SBS5非常相似。這有助于解釋?zhuān)驗(yàn)镾BS1的病因是已知的。這也告訴我們，我們沒(méi)有發(fā)現(xiàn)任何全新的過(guò)程。與任何先前定義的signatures不相似的提取signatures也不能證明是“新”signatures。提取的signatures可能是已知signatures的組合，不能被NMF分割。例如，當(dāng)用于NMF的樣本太少時(shí)，就會(huì)發(fā)生這種情況。
4）Visualizing NMF results
我們可以繪制96-profile的signatures(在查看SNV時(shí))
## 4）Visualizing NMF results
p <- plot_96_profile(nmf_res$signatures, condensed = TRUE)
ggsave(filename = paste0(opt$od, "/NMF_96_profile.png"), width = 9, height = 6, plot = p)
也可以使用barplot可視化signatures的貢獻(xiàn)：
p <- plot_contribution(nmf_res$contribution, nmf_res$signature, mode = "relative")
ggsave(filename = paste0(opt$od, "/NMF_contribution.png"), width = 7, height = 3.5, plot = p)
每個(gè)樣本中的signature的相對(duì)貢獻(xiàn)對(duì)也可以使用熱圖進(jìn)行展示：
p <- plot_contribution_heatmap(nmf_res$contribution, cluster_samples = TRUE, cluster_sigs = TRUE)
ggsave(filename = paste0(opt$od, "/NMF_contribution_heatmap.png"), width = 10, height = 7, plot = p)
我們也可以事先對(duì)樣本和signature進(jìn)行聚類(lèi)然后繪制熱圖：
# cluster the signatures
hclust_signatures <- cluster_signatures(nmf_res$signatures, method = "average")
signatures_order <- colnames(nmf_res$signatures)[hclust_signatures$order]
signatures_order
# cluster the samples
hclust_samples <- cluster_signatures(mut_mat, method = "average")
samples_order <- colnames(mut_mat)[hclust_samples$order]
samples_order
# plot
p <- plot_contribution_heatmap(nmf_res$contribution,
sig_order = signatures_order, sample_order = samples_order,
cluster_sigs = FALSE, cluster_samples = FALSE)
ggsave(filename = paste0(opt$od, "/NMF_contribution_heatmap_1.png"), width = 6, height = 4, plot = p)
基于the signatures and their contribution，NMF對(duì)每個(gè)樣本重新構(gòu)建了一個(gè)突變譜，NMF工作得越好，重建的profile就越接近原始profile，將單個(gè)樣本的重建突變譜與原始突變譜進(jìn)行比較：
p <- plot_compare_profiles(mut_mat[, 1],
nmf_res$reconstructed[, 1],
profile_names = c("Original", "Reconstructed"),
condensed = TRUE)
ggsave(filename = paste0(opt$od, "/NMF_compare_profiles.png"), width = 7, height = 4.5, plot = p)
還可以繪制所有樣本的原始矩陣和重構(gòu)矩陣之間的余弦相似度。當(dāng)重建profile與原始profile的余弦相似度大于0.95時(shí)，認(rèn)為重建profile非常好。
p <- plot_original_vs_reconstructed(mut_mat, nmf_res$reconstructed, y_intercept = 0.95)
ggsave(filename = paste0(opt$od, "/NMF_original_vs_reconstructed.png"), width = 7, height = 4.5, plot = p)
下次學(xué)習(xí)Signature refitting~
