肿瘤突变特征推断工具SomaticSignatures
之前我们介绍了如何根据体细胞突变结果推断已知的突变特征，这里的已知突变特征主要是来自COSMIC网站，从分析结果中我们发现有一部分特征无法跟COSMIC数据匹配，这一部分我们一般归纳为未知的突变特征（Unknow Signature），研究者完全可以使用自己整理的数据进行突变特征推断和研究，今天我们就一起来学习使用R包SomaticSignatures进行样本未知突变特征推断。
首先是软件的下载安装，这里推荐用conda安装，简单方便。支持最新版本R，小编就是在R 4.3.3下安装的。
conda install r-devtools
conda install r-biocmanager
conda install bioconductor-somaticsignatures
BiocManager::install("BSgenome.Hsapiens.1000genomes.hs37d5")
PS：如果体细胞突变检测用的UCSC参考基因组，需要安装对应的参考基因组文件。
conda install bioconductor-somaticcanceralterations
PS：该测试数据来自8个不同癌症的外显子测序的项目。
library(SomaticCancerAlterations)
sca_data = unlist(scaLoadDatasets()) #unlist()函数将嵌套的数据结构转化为普通向量
sca_data$study = factor(gsub("(.*)_(.*)", "\\1", toupper(names(sca_data))))
sca_data = unname(subset(sca_data, Variant_Type %in% "SNP"))
sca_data = keepSeqlevels(sca_data, hsAutosomes(), pruning.mode = "coarse")
head(sca_data)
说明：这里我们使用测试数据进行演示分析，也可以使用自有数据，只需要提供c( "Sample","chr", "pos","ref", "alt") 这几列数据，就可以构建输入文件。
library(SomaticSignatures)
library(BSgenome.Hsapiens.1000genomes.hs37d5)
sca_vr = VRanges(
seqnames = a$chr,
ranges = IRanges(start = a$pos,end = a$pos+1),
ref = a$ref,
alt = a$alt,
sampleNames = as.character(a$Sample),
study = as.character(a$study))
sca_vr
sca_motifs = mutationContext(sca_vr, BSgenome.Hsapiens.1000genomes.hs37d5)
head(sca_motifs)
说明：这个对象跟上面sca_vr区别在于末尾增加两列alteration和context，最后一列context为上下游碱基。
这一步完成，我们可以绘制每个样本突变频谱图。
plotMutationSpectrum(sca_motifs, "study")
说明：这里仅仅是利用统计的突变频数，下一步我们会对频数进行标准化，方便进行sginature推断。
escc_sca_mm = motifMatrix(sca_motifs, group = "study", normalize = TRUE)
head(escc_sca_mm)
说明：这里按照癌种统计的96中突变模式，并对统计结果进行标准化；行是突变模式（96种），列是各个癌种数据。
colSums(escc_sca_mm)
接下来就是Signature推断了，但是别急，在这之前需要确定Signature合适的分解数目，这里使用
assessNumberSignatures()函数来评估找到最合适的signature分解数目。
n_sigs = 2:8
gof_nmf = assessNumberSignatures(escc_sca_mm , n_sigs, nReplicates = 5)
head(gof_nmf)
结果说明：
NumberSignatures：指定的突变特征数；
RSS：残差平方和，用于评估模型的拟合程度，越小表示模型拟合程度越好；
ExplainedVariance：可解释方差得分（EVS），用以衡量从波动性角度解释模型对数据的契合程度，取值范围通常为[0,1]，值越接近于1表示模型性能越好；
Replicate：模型重复次数。
为了更直观的观察结果，绘制在不同Signature数目下RSS和ExplainedVariance分布图。
plotNumberSignatures(gof_nmf)
说明：看图可知，当signature数到7时，对应的RSS和ExplainedVariance接近水平，故我们选择signature合适分解数为7，进行下一步signature推断。
sigs_nmf = identifySignatures(escc_sca_mm, 7, nmfDecomposition)
#sigs_nmf = identifySignatures(escc_sca_mm, 7, pcaDecomposition)
说明：这个包提供了两种方法，分别是NMF和PCA，通常我们使用NMF方法。
head(sigs_nmf@signatures)
head(sigs_nmf@samples)
save(escc_sca_mm, sigs_nmf, file = 'escc_denovo_results.Rata')
load(file = 'escc_denovo_results.Rata')
str(sigs_nmf)
library(ggplot2)
plotSignatureMap(sigs_nmf) +
ggtitle("Somatic Signatures: NMF - Heatmap")
说明：横坐标是确定的7种denovo signature，纵坐标是96种突变模式，热图用来展示不同signature在哪些突变模式上具有共同的特征。
plotSignatures(sigs_nmf, normalize =T) +
ggtitle("Somatic Signatures: NMF - Barchart") +
facet_grid(signature ~ alteration,scales = "free_y")
说明：横坐标是96种突变模式，且将其归纳为6种特征模式，纵坐标是确定的7种denovo signature，频谱图用来展示每个signature的6中特征模式分布情况，从图上可以看出signature之间的相似性。
到这里我们就已经完成了未知突变特征（denovo signature）推断，到这里就结束了吗？NO！获得未知突变特征只是第一步，接下来需要对突变特征进一步分析，特别是生物学意义研究。
还记得开头提到的COSMIC突变特征数据吗？我们可以将获得的denovo signature跟cosmic signature进行相关性分析，这样就能给一些的denovo signature快速注释生物学意义。
load(file = 'escc_denovo_results.Rata')
str(sigs_nmf)
sp <- sigs_nmf@signatures
sp <- apply(sp,2,function(x){x/sum(x)})
denovo <- sp
#https://cancer.sanger.ac.uk/cancergenome/assets/signatures_probabilities.txt
cosmic <- read.table('signatures_probabilities.txt',header = T,sep = '\t')[,1:33]
head(cosmic[,1:3])
tmp=cosmic[,2];substr(tmp,2,2) <- '.'
rownames(cosmic)=paste(gsub('>','',cosmic[,1]),tmp)
PS：进行简单的转换，保证两个signature的矩阵行名是一样的。
cosmic <- cosmic[,4:33]
comp <- cbind(denovo[rownames(cosmic),], cosmic)
colSums(comp)
corr <- cor(comp)
library(pheatmap)
pheatmap::pheatmap(corr)
说明：通过热图可以发现哪些denovo signature跟COSMIC中已知的signature具有相似性，从而将某些signature归纳，对没有进行归纳的signature继续研究其生物学意义。
