老规矩，先奉上学习资料链接：
这一次学习Mutational signatures的构建。
Mutational signatures被认为代表了突变过程，并以具有特定序列上下文的突变类型的特定贡献为特征。
NMF：使用非负矩阵分解(NMF)，可以从突变计数矩阵中重新提取Mutational signatures。
signature refitting：还可以确定突变计数矩阵对先前定义的Mutational signatures的暴露（exposure）。
NMF对大样本最有用，而signature refitting可用于单个样本。
接下来首先讨论NMF，然后再讨论signature refitting，最后，我们将讨论如何直接分析mutational profile和signatures之间的相似性。
NMF提取De novo mutational signature
1）NMF
NMF中的一个关键参数是分解秩(the factorization rank)，它是提取的mutational signature的数量。可以使用NMF包确定最优分解秩。
先跟上一个教程一样，处理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)
一般来说，数据集越大，rank会更高。这里我们展示SNV的NMF。其他突变类型执行NMF的工作方式与此相同。
使用NMF包生成估计秩图，需要很长时间，最好是最后写成传参脚本后台运行。
# 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从突变计数矩阵中提取mutational signatures。在本例中，rank 为2，提取2个mutational signatures(对于较大的数据集，明智的做法是通过改变nrun参数来执行更多的迭代，以获得稳定性并避免不良的局部最小值)。
上图绘制了rank2-5的各种评价指标，所以到底是因为什么值选取的rank 2？？？
想起来前面还有段话：
The most common approach is to choose the smallest rank for which cophenetic correlation coefficient starts decreasing.
也就是上图左上角的那个cophenetic，rank2-5的时候一直往下降，最小的rank就是2。
确定好rank之后，提取mutational signatures：
nmf_res <- extract_signatures(mut_mat, rank = 2, nrun = 10, single_core = TRUE)
这样得到了两个mutational signatures：
2）Bayesian NMF
也可以使用变分贝叶斯NMF（variational bayes NMF）。这样可以更容易地确定正确的rank。需要安装ccfindR包，然后可以确定最优signatures数，这同样需要很长时间。
# 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()
结果图如下：
然后提取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
我们可以对提取出来的signatures的名字更改为习惯用名：
## 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的名称更改为这些已知signatures可能会很有用。这通常使你更容易解释你的结果。
要做到这一点，首先需要读取一些已经存在的signatures。在这里，我们将使用COSMIC (v3.2) (Alexandrov et al. 2020)的signatures。(稍后我们将讨论如何使用其他signatures矩阵。)
现在可以更改NMF提取的signatures的名称。在本例中，如果signatures的名称与现有COSMIC signatures的余弦相似度大于0.85，则更改该signatures的名称。
signatures = get_known_signatures()
nmf_res <- rename_nmf_signatures(nmf_res, signatures, cutoff = 0.85)
colnames(nmf_res$signatures)
[1] "SBS5-like" "SBS1-like"
我们现在看到，我们提取的signatures与COSMIC signatures SBS1和SBS5非常相似。这有助于解释，因为SBS1的病因是已知的。这也告诉我们，我们没有发现任何全新的过程。与任何先前定义的signatures不相似的提取signatures也不能证明是“新”signatures。提取的signatures可能是已知signatures的组合，不能被NMF分割。例如，当用于NMF的样本太少时，就会发生这种情况。
4）Visualizing NMF results
我们可以绘制96-profile的signatures(在查看SNV时)
## 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的贡献：
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)
每个样本中的signature的相对贡献对也可以使用热图进行展示：
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)
我们也可以事先对样本和signature进行聚类然后绘制热图：
# 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对每个样本重新构建了一个突变谱，NMF工作得越好，重建的profile就越接近原始profile，将单个样本的重建突变谱与原始突变谱进行比较：
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)
还可以绘制所有样本的原始矩阵和重构矩阵之间的余弦相似度。当重建profile与原始profile的余弦相似度大于0.95时，认为重建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)
下次学习Signature refitting~
