使用R包deconstructSigs根据已知的signature进行比例推断

但是更多的时候，研究者不会去自定义signature的，而是会直接使用sanger研究所科学家【1】提出来了肿瘤somatic突变的signature概念 ，把96突变频谱的非负矩阵分解后的30个特征，在cosmic数据库可可以查询到的30个特征。不同的特征有不同的生物学含义【2】，比如文章【3】 就是使用了 这些signature区分生存！主要是R包deconstructSigs可以把自己的96突变频谱对应到cosmic数据库的30个突变特征。

• 【1】https://software.broadinstitute.org/cancer/cga/msp
• 【2】https://en.wikipedia.org/wiki/Mutational_signatures
• 【3】https://www.nature.com/articles/s41586-019-1056-z

而且我在教程：比较不同的肿瘤somatic突变的signature 也分享了如何比较不同方法拿到的signature，这样它们的生物学意义就可以联系起来了。

首先R包deconstructSigs可以把自己的96突变频谱对应到cosmic数据库的30个突变特征

在文章主页下载；https://static-content.springer.com/esm/art%3A10.1038%2Fs41422-020-0333-6/MediaObjects/41422_2020_333_MOESM23_ESM.csv

这个是大于500M的CSV文件，下载后修改名字，然后**读入R，**筛选拿到SNV突变位点，去对应的参考基因组里面获取突变上下文，就是  mut.to.sigs.input 函数：

library(data.table)
a=fread('../maf.csv',data.table = F)
a[1:4,1:3]
colnames(a)
mut=a
table(mut$Variant_Type)
mut=mut[mut$Variant_Type=='SNP',] 
a=mut[,c(10,2,3,8,9)]
colnames(a)=c( 'Sample','chr', 'pos','ref',  'alt')
# 下面的 mut.to.sigs.input 函数仅仅是需要 SNV的5列信息即可
sigs.input <- mut.to.sigs.input(mut.ref = a, 
                                sample.id = 'Sample', 
                                chr = 'chr', 
                                pos = 'pos', 
                                ref = 'ref', 
                                alt = 'alt',
                                bsg = BSgenome.Hsapiens.UCSC.hg19)


然后 对每个样本，循环运行 whichSignatures 函数，判断每个样本的signature组成：

w=lapply(unique(a$Sample) , function(i){
  ## signatures.cosmic signatures.nature2013
  sample_1 = whichSignatures(tumor.ref = sigs.input[,], 
                             signatures.ref = signatures.cosmic, 
                             sample.id =  i, 
                             contexts.needed = TRUE,
                             tri.counts.method = 'default')
  print(i)
  return(sample_1$weights)
})
w=do.call(rbind,w)
library(pheatmap)
pheatmap(t(w),cluster_rows = F,cluster_cols = T)
pheatmap(w,cluster_rows = T,cluster_cols = F)
mut.wt=w
save(mut.wt,file = 'wgs-mut.wt.Rdata')

这个时候，可以选择signatures.cosmic 和 signatures.nature2013，这两个内置的signature，比如signatures.cosmic，其实在网络文件，signatures_probabilities.txt 可以查看。

R包SomaticSignatures进行denovo的signature推断后的11个signature

就是前面对每个样本，循环运行 whichSignatures 函数，判断每个样本的signature组成的时候，替换掉内置的signatures.cosmic 和 signatures.nature2013，代码如下：

signatures.cosmic
rowSums(signatures.cosmic)
colnames(signatures.cosmic)
load(file = 'escc_denovo_results.Rdata')
str(sigs_nmf) 
# sp signatures_probabilities
sp=sigs_nmf@signatures
head(sp)
colSums(sp)
sp=apply(sp,2,function(x){
  x/sum(x)
})
head(sp)
colSums(sp)
sp=t(sp)
chr=colnames(sp)
colnames(sp)=gsub(' ','',
                  paste(substring(chr,4,4),
                        '[',substring(chr,1,1),'>',substring(chr,2,2),']',
                        substring(chr,6,6)
                  ))
colnames(signatures.cosmic)
sp=sp[,colnames(signatures.cosmic)]
sc=signatures.cosmic


其中代码里面的escc_denovo_results.Rdata文件，来自于教程：使用R包SomaticSignatures进行denovo的signature推断。

把自己的11个signature制作成为R包内置的signatures.cosmic 和 signatures.nature2013样式，这个代码非常复杂，需要大家自行认真的理解。

接下来对每个样本，循环运行 whichSignatures 函数的代码如下：

b=a[a$Sample %in% head(s,20),]
w=lapply(unique(b$Sample) , function(i){
  ## signatures.cosmic signatures.nature2013
  sample_1 = whichSignatures(tumor.ref = sigs.input[,], 
                             signatures.ref = as.data.frame(sp), 
                             sample.id =  i, 
                             contexts.needed = TRUE,
                             tri.counts.method = 'default')
  print(i)
  return(sample_1$weights)
})
w=do.call(rbind,w)
library(pheatmap)
pheatmap(t(w),cluster_rows = F,cluster_cols = T)
pheatmap(w,cluster_rows = T,cluster_cols = F)
sp=sigs_nmf@samples
sp=sp[rownames(sp) %in% head(s,20),]
pheatmap(sp,cluster_rows = F,cluster_cols = F)
pheatmap(w,cluster_rows = F,cluster_cols = F)

可以看到，自己制作好的 as.data.frame(sp) 就替代了 R包内置的signatures.cosmic 和 signatures.nature2013。

这个时候，就会根据自己的11个signature进行分解，而不是原来的R包内置的signatures.cosmic 和 signatures.nature2013两种分解模式。

但是可以对比两次的11个signature分解的差异。

首先看看教程：使用R包deconstructSigs根据已知的signature进行比例推断，的比例情况：

然后看看教程：使用R包SomaticSignatures进行denovo的signature推断，的比例情况；