前往小程序，Get更优阅读体验！
立即前往
文档建议反馈控制台
首页
学习
活动
专区
工具
TVP
最新优惠活动
文章/答案/技术大牛
发布
首页
学习
活动
专区
工具
TVP最新优惠活动
返回腾讯云官网
社区首页 >专栏 >使用 sigminer 进行突变模式分析

使用 sigminer 进行突变模式分析

作者头像
生信菜鸟团
发布于 2020-06-10 11:24:58
2K0
发布于 2020-06-10 11:24:58
举报
文章被收录于专栏：生信菜鸟团

突变模式分析（Mutual Signature Analysis）已经逐步成为变异检测后一个通用分析，本文简单介绍如何使用sigminer进行突变模式分析，以解决2大分析任务：

  • 从头发现签名
  • 已知一些参考Signatures，我们想要定量Signature的Exposure（或称为Activity）

支持SBS，DBS，INDEL以及副本号（研究中）。

安装

代码语言：javascript
复制
BiocManager::install("sigminer", dependencies = TRUE)

加载包

代码语言：javascript
复制
library(sigminer)
library(NMF)

数据输入

sigminer的开发与maftools很有渊源，使用MAF对象作为数据的存储结构。如果你会使用maftools读入突变数据，那么就会使用sigminer读入突变数据，支持 data.frame 和MAF文件。

这里我们使用maftools内置数据集，maftools其实本身也可以做signature分析（但不是它的核心）。

代码语言：javascript
复制
laml.maf <- system.file("extdata", "tcga_laml.maf.gz", package = "maftools", mustWork = TRUE)
laml <- read_maf(maf = laml.maf)
#> -Reading
#> -Validating
#> -Silent variants: 475
#> -Summarizing
#> -Processing clinical data
#> --Missing clinical data
#> -Finished in 0.290s elapsed (0.269s cpu)
# 与 maftools::read.maf(maf = laml.maf) 结果是一样的
# read_maf 是对 read.maf 的封装

生成突变分类矩阵

使用 sig_tally() 对突变进行归类整理，针对MAF对象，支持设置 mode 为'SBS'，'DBS'，'ID'以及'ALL'。

mode 为'ALL`时会试图生成所有矩阵：

代码语言：javascript
复制
mats <- mt_tally <- sig_tally(
  laml,
  ref_genome = "BSgenome.Hsapiens.UCSC.hg19",
  useSyn = TRUE,
  mode = "ALL"
)

str(mats, max.level = 1)
#> List of 6
#>  $ SBS_6        : int [1:193, 1:6] 1 0 0 2 1 0 1 1 2 2 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ SBS_96       : int [1:193, 1:96] 0 0 0 0 0 0 1 0 1 0 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ SBS_1536     : int [1:193, 1:1536] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ ID_28        : int [1:193, 1:28] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ ID_83        : int [1:193, 1:83] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ APOBEC_scores:Classes 'data.table' and 'data.frame':  182 obs. of  44 variables:
#>   ..- attr(*, ".internal.selfref")=<externalptr>
#>   ..- attr(*, "index")= int(0)
#>   .. ..- attr(*, "__APOBEC_Enriched")= int [1:182] 106 147 5 6 8 9 10 11 12 13 ...

这个数据集没有双碱基替换，所以没有相应的矩阵。

代码语言：javascript
复制
mt_tally$SBS_96[1:5, 1:5]
#>              A[T>C]A C[T>C]A G[T>C]A T[T>C]A A[C>T]A
#> TCGA-AB-2802       0       0       0       0       0
#> TCGA-AB-2803       0       0       0       0       2
#> TCGA-AB-2804       0       0       0       0       2
#> TCGA-AB-2805       0       0       0       0       0
#> TCGA-AB-2806       0       0       0       0       0

针对整个数据集的分类就可以画图，签名其实就是它的分解。

代码语言：javascript
复制
show_catalogue(mt_tally$SBS_96 %>% t(), mode = "SBS", style = "cosmic")
代码语言：javascript
复制
show_catalogue(mt_tally$SBS_6 %>% t(), mode = "SBS", style = "cosmic")

注意上面对矩阵进行了转置。

估计签名数

这一步实际上是多次运行NMF，查看一些指标的变化，用于后续确定提取多少个签名。

代码语言：javascript
复制
est <- sig_estimate(mt_tally$SBS_96, range = 2:5, nrun = 2, pConstant = 1e-9, verbose = TRUE)
#> Compute NMF rank= 2  ... + measures ... OK
#> Compute NMF rank= 3  ... + measures ... OK
#> Compute NMF rank= 4  ... + measures ... OK
#> Compute NMF rank= 5  ... + measures ... OK
#> Estimation of rank based on observed data.
#>   method   seed rng metric rank sparseness.basis sparseness.coef  rss  evar
#> 2 brunet random   1     KL    2            0.588           0.653 1924 0.411
#> 3 brunet random   1     KL    3            0.648           0.623 1875 0.426
#> 4 brunet random   1     KL    4            0.685           0.633 1802 0.448
#> 5 brunet random   1     KL    5            0.698           0.687 1723 0.472
#>   silhouette.coef silhouette.basis residuals niter   cpu cpu.all nrun
#> 2           1.000            1.000      2880  1040 0.496    4.86    2
#> 3           0.769            0.830      2718   800 0.409    4.72    2
#> 4           0.656            0.767      2543  1790 0.946    1.82    2
#> 5           0.609            0.712      2415  1750 1.044    5.23    2
#>   cophenetic dispersion silhouette.consensus
#> 2      0.888      0.576                0.712
#> 3      0.802      0.577                0.494
#> 4      0.784      0.655                0.473
#> 5      0.748      0.700                0.398

这里加入了一个 pConstant = 1e-9，是因为NMF包调用矩阵分解存在bug，一个分类的和不能为0。报错就加一个伪小数，不报错就可以不管。

代码语言：javascript
复制
show_sig_number_survey2(est$survey)

cophenetic是一个主要指标，我们看到一直在下降。不过我们观察到残差，稀疏在往好的方向变化，这里可以选择4个尝试（上面运行最好30-50次可以得到稳定结果）。

提取签名和可视化

代码语言：javascript
复制
sigs <- sig_extract(mt_tally$SBS_96, n_sig = 4, nrun = 2, pConstant = 1e-9)

生成的是一个带 Signature 类信息的列表：

代码语言：javascript
复制
str(sigs, max.level = 1)
#> List of 6
#>  $ Signature     : num [1:96, 1:4] 2.67e-16 7.00 6.00 1.00 1.06e-07 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ Signature.norm: num [1:96, 1:4] 4.86e-19 1.27e-02 1.09e-02 1.82e-03 1.93e-10 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ Exposure      : num [1:4, 1:193] 3.687 4.047 0.243 2.023 2.402 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ Exposure.norm : num [1:4, 1:193] 0.3687 0.4047 0.0243 0.2023 0.1716 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ K             : int 4
#>  $ Raw           :List of 3
#>  - attr(*, "class")= chr "Signature"
#>  - attr(*, "nrun")= num 2
#>  - attr(*, "method")= chr "brunet"
#>  - attr(*, "pConstant")= num 1e-09
#>  - attr(*, "seed")= num 123456
#>  - attr(*, "call_method")= chr "NMF"

sigminer也有很多函数专门用户获取一些信息，可视化等等。

我们先看一个最常见的突变模式图谱：

代码语言：javascript
复制
p <- show_sig_profile(sigs, mode = "SBS", style = "cosmic")
p

接下来可以计算下它与COSMIC签名的相似性，评估病因，对于SBS有2个COSMIC数据库版本legacy（30个，目前最常用的）与SBS v3。

代码语言：javascript
复制
get_sig_similarity(sigs, sig_db = "legacy")
#> -Comparing against COSMIC signatures
#> ------------------------------------
#> --Found Sig1 most similar to COSMIC_6
#>    Aetiology: defective DNA mismatch repair [similarity: 0.797]
#> --Found Sig2 most similar to COSMIC_1
#>    Aetiology: spontaneous deamination of 5-methylcytosine [similarity: 0.8]
#> --Found Sig3 most similar to COSMIC_15
#>    Aetiology: defective DNA mismatch repair [similarity: 0.535]
#> --Found Sig4 most similar to COSMIC_1
#>    Aetiology: spontaneous deamination of 5-methylcytosine [similarity: 0.89]
#> ------------------------------------
#> Return result invisiblely.
代码语言：javascript
复制
sim <- get_sig_similarity(sigs, sig_db = "SBS")
#> -Comparing against COSMIC signatures
#> ------------------------------------
#> --Found Sig1 most similar to SBS6
#>    Aetiology: defective DNA mismatch repair [similarity: 0.803]
#> --Found Sig2 most similar to SBS1
#>    Aetiology: spontaneous or enzymatic deamination of 5-methylcytosine [similarity: 0.737]
#> --Found Sig3 most similar to SBS15
#>    Aetiology: Defective DNA mismatch repair [similarity: 0.567]
#> --Found Sig4 most similar to SBS1
#>    Aetiology: spontaneous or enzymatic deamination of 5-methylcytosine [similarity: 0.878]
#> ------------------------------------
#> Return result invisiblely.
代码语言：javascript
复制
add_labels(p, x = 0.72, y = 0.25, y_end = 0.9, labels = sim, n_label = 4)

这里的坐标位置需要自己细调一下。

自动提取签名

SignatureAnalyzer提供的变异NMF提供了自动提取Sigantures的功能，无需要自己判断提取的sign数量，这个可以通过 sig_auto_extract() 实现。该算法会生成更大的稀疏（相互之间相互）的签名，因此偏向于生成更多的从我多年研究签名的经验来看，它对于单点突变还是非常友好的。

代码语言：javascript
复制
sigs2 <- sig_auto_extract(mt_tally$SBS_96, nrun = 2)
#> Warning: [ONE-TIME WARNING] Forked processing ('multicore') is disabled
#> in future (>= 1.13.0) when running R from RStudio, because it is
#> considered unstable. Because of this, plan("multicore") will fall
#> back to plan("sequential"), and plan("multiprocess") will fall back to
#> plan("multisession") - not plan("multicore") as in the past. For more details,
#> how to control forked processing or not, and how to silence this warning in
#> future R sessions, see ?future::supportsMulticore
#>
#> 载入程辑包：'purrr'
#> The following objects are masked from 'package:foreach':
#>
#>     accumulate, when
#> The following object is masked from 'package:XVector':
#>
#>     compact
#> The following object is masked from 'package:GenomicRanges':
#>
#>     reduce
#> The following object is masked from 'package:IRanges':
#>
#>     reduce
#> Select Run 2, which K = 2 as best solution.

画图方式是完全一样的。

代码语言：javascript
复制
p <- show_sig_profile(sigs2, mode = "SBS", style = "cosmic")
p

虽然上面都是粗略的分析，但这种方法感觉更好。

实际研究时选择某些方法都需要根据数据还自己的需要决定，也可以比较上面的结果。

签名活动图谱

sigminer提供绝对和相对两种签名活动度值。

代码语言：javascript
复制
get_sig_exposure(sigs2) %>% head()
#>          sample  Sig1  Sig2
#> 1: TCGA-AB-2802 0.000  9.50
#> 2: TCGA-AB-2803 0.000 13.06
#> 3: TCGA-AB-2804 0.000  6.75
#> 4: TCGA-AB-2805 0.000 13.06
#> 5: TCGA-AB-2806 0.911 12.18
#> 6: TCGA-AB-2807 0.000 23.86
get_sig_exposure(sigs2, type = "relative") %>% head()
#> Filtering the samples with no signature exposure:
#> TCGA-AB-2823 TCGA-AB-2834 TCGA-AB-2840 TCGA-AB-2848 TCGA-AB-2892 TCGA-AB-2909 TCGA-AB-2942
#>          sample   Sig1 Sig2
#> 1: TCGA-AB-2802 0.0000 1.00
#> 2: TCGA-AB-2803 0.0000 1.00
#> 3: TCGA-AB-2804 0.0000 1.00
#> 4: TCGA-AB-2805 0.0000 1.00
#> 5: TCGA-AB-2806 0.0696 0.93
#> 6: TCGA-AB-2807 0.0000 1.00

画图直接把对象扔进去就可以了。

代码语言：javascript
复制
show_sig_exposure(sigs, rm_space = TRUE, style = "cosmic")

根据已知的签名

这个比较有名的软件就是deconstructSigs了。sigminer也支持了这个功能，而且能够使用目前cosmic的所有图谱，也可以使用自己从头发现的签名 sig_fit()。

我们试着处理5个样本（支持单样本），使用COSMIC v2 30个签名作为参考。

代码语言：javascript
复制
examp_fit <- sig_fit(mt_tally$SBS_96[1:5, ] %>% t(), sig_index = "ALL", sig_db = "legacy")
head(examp_fit)
#>          TCGA-AB-2802 TCGA-AB-2803 TCGA-AB-2804 TCGA-AB-2805 TCGA-AB-2806
#> COSMIC_1        3.129        10.00            7         13.6         10.8
#> COSMIC_2        0.895         0.00            0          0.0          0.0
#> COSMIC_3        0.000         0.00            0          0.0          0.0
#> COSMIC_4        0.000         1.72            0          0.0          0.0
#> COSMIC_5        0.000         0.00            0          0.0          0.0
#> COSMIC_6        1.088         0.00            0          0.0          0.0

画图也很简单：

代码语言：javascript
复制
show_sig_fit(examp_fit, palette = NULL) + ggpubr::rotate_x_text()
#> ℹ [2020-06-07 16:06:33]: Started.
#> ℹ [2020-06-07 16:06:33]: Checking input format.
#> ✓ [2020-06-07 16:06:33]: Checked.
#> ℹ [2020-06-07 16:06:33]: Checking filters.
#> ℹ [2020-06-07 16:06:33]: Checked.
#> ℹ [2020-06-07 16:06:33]: Plotting.
#> ℹ [2020-06-07 16:06:33]: 0.048 secs elapsed.

如果设置散点图，可以把单样本结果画出来。

代码语言：javascript
复制
show_sig_fit(examp_fit,
  palette = NULL, plot_fun = "scatter",
  signatures = paste0("COSMIC_", c(1, 2, 4, 6, 19))
) + ggpubr::rotate_x_text()
#> ℹ [2020-06-07 16:06:34]: Started.
#> ℹ [2020-06-07 16:06:34]: Checking input format.
#> ✓ [2020-06-07 16:06:34]: Checked.
#> ℹ [2020-06-07 16:06:34]: Checking filters.
#> ℹ [2020-06-07 16:06:34]: Checked.
#> ℹ [2020-06-07 16:06:34]: Plotting.
#> ! [2020-06-07 16:06:34]: When plot_fun='scatter', setting legend='top' is recommended.
#> ℹ [2020-06-07 16:06:34]: 0.047 secs elapsed.

上面展示了最核心的分析和可视化功能，sigminer还支持很多功能，我就不再出现了。用户可以阅读官方文档 https://shixiangwang.github.io/sigminer-doc/进一步 学习，后续我会根据研究情况进一步开发。当然，读者完全可以基于上面的分析的结果值进行各种个性化分析。

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2020-06-08，如有侵权请联系 cloudcommunity@tencent.com 删除
analysis
auto
bug
data
signature

本文分享自 生信菜鸟团 微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ，欢迎热爱写作的你一起参与！

analysis
auto
bug
data
signature
评论
登录后参与评论
0 条评论
热度
最新
登录 后参与评论
推荐阅读
LV.
文章
0
获赞
0
目录
  • 安装
  • 加载包
  • 数据输入
  • 生成突变分类矩阵
  • 估计签名数
  • 提取签名和可视化
  • 自动提取签名
  • 签名活动图谱
  • 根据已知的签名
领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2024 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 深公网安备号 44030502008569

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号 | 京公网安备号11010802020287

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档

Copyright © 2013 - 2024 Tencent Cloud.

All Rights Reserved. 腾讯云 版权所有

登录 后参与评论