deconstructSigs-mutation signature看一下你的數據是什么“氣質”的？

本文首發於“生信補給站” https://mp.weixin.qq.com/s/k7yzk9hPX3Bi-ohAo83ZYw

還有其他 R統計繪圖生信的干貨，也許有需要的呢？

Mutational Signatures 首次出現在2013年的nature文章Signatures of mutational processes in human cancer中（https://www.nature.com/articles/nature12477）。**將mutation位置加上前后一個鹼基，構成三鹼基模式，然后統計96（6 * 4 * 4）種突變組合的情況。

好奇為什么是96種的，可以查一下文獻。

本文介紹如何利用deconstructSigs-R包進行mutation signature分析。

一准備R包，數據

#install.packages("deconstructSigs") 
library(deconstructSigs)
#讀入數據
head(sample.mut.ref)
 Sample  chr      pos ref alt
1      1 chr1   905907   A   T
2      1 chr1  1192480   C   A
3      1 chr1  1854885   G   C
4      1 chr1  9713992   G   A
5      1 chr1 12908093   C   A
6      1 chr1 17257855   C   T

class(sample.mut.ref)
## [1] "data.frame"

只需要將自己的數據整理成以上五列（ID，chr，pos，ref，alt ）信息即可，如果是TCGA中的MAF文件也是很好提取的。

二 mut.to.sigs.input構建輸入文件

使用 mut.to.sigs.input 函數，構建計算signature的輸入文件，得到每個樣本的96種三鹼基類型。

# Convert to deconstructSigs input
sigs.input <- mut.to.sigs.input(mut.ref = sample.mut.ref, 
                                sample.id = "Sample", 
                                chr = "chr", 
                                pos = "pos", 
                                ref = "ref", 
                                alt = "alt")

注：這一步也許會提示沒有XX包，按照要求下載指定R包即可（也許是數據庫，耐心安裝）。

#查看結果信息
dim(sigs.input)
#[1] 2 96  
head(t(sigs.input)) #只有兩個sample：“1”和“2”
         1 2
A[C>A]A  9 1
A[C>A]C  7 1
A[C>A]G  5 0
A[C>A]T  7 0
C[C>A]A 10 3
C[C>A]C 18 2

以上就得到了sample.mut.ref文件中的每一個sample的96種三鹼基類型的結果了。

三推斷signature的組成

# Determine the signatures contributing to the two example samples
sample_1 = whichSignatures(tumor.ref = sigs.input, 
                           signatures.ref = signatures.cosmic, 
                           sample.id = 1, 
                           contexts.needed = TRUE,
                           tri.counts.method = 'default')

其中：

tumor.ref：每個sample的96種三鹼基突變序列 signatures.ref：已知的signatures參考文件，可選signatures.nature2013和signatures.cosmic sample.id：對應tumor.ref文件中的樣本名 contexts.needed ：是否需要突變上下文 tri.counts.method：三核酸序列標准化方式，默認“default” 不進行標准化；或者選擇exome,genome,exome2genome,genome2exome 來限定區域。

3.2 查看返回結果

#查看結果
class(sample_1)
#查看權重結果
sample_1$weights

#輸出tumor的三鹼基序列百分比
sample_1$tumor

#三鹼基序列百分比 * 權重
sample_1$product

whichSignatures會輸出5個元素的list文件：

weights -- data frame containing the weights assigned to each of the k signatures of the input signatures matrix
tumor -- matrix of the trinucleotide contexts for the tumor sample used as input
product -- matrix obtained when the tumor matrix is multiplied by the assigned weights
diff -- matrix representing the difference between the tumor matrix and product matrix
unknown -- numeric weight not assigned to any of the input signatures

3.3 指定signature權重

通過associated參數指定參與計算的signature

sample_1.associate = whichSignatures(tumor.ref = sigs.input, 
                           signatures.ref = signatures.cosmic, 
                           sample.id = 1, 
                           associated = c("Signature.1","Signature.22"),
                           contexts.needed = TRUE,
                           tri.counts.method = 'default')
sample_1.associate$weights

3.4 設定signature的閾值

通過signature.cutoff設定閾值，小於此值的為0

sample_1.cutoff = whichSignatures(tumor.ref = sigs.input, 
                           signatures.ref = signatures.cosmic, 
                           sample.id = 1, 
                           contexts.needed = TRUE,
                           signature.cutoff = 0.08 ,
                           tri.counts.method = 'default')

sample_1.cutoff$weights

四 plotSignatures 可視化

# Plot example
plot_example <- whichSignatures(tumor.ref = sigs.input,
                      signatures.ref = signatures.cosmic,
                      sample.id = 1)

# Plot output
plotSignatures(plot_example, sub = 'example')