search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
118 stars 30 forks source link

使用 sigminer 进行突变模式分析 #241

Closed ixxmu closed 4 years ago

ixxmu commented 4 years ago

https://mp.weixin.qq.com/s/D_v_hYwxs_OWFUsnFyqDGw

github-actions[bot] commented 4 years ago

使用 sigminer 进行突变模式分析 by 生信菜鸟团


上周一我们生信菜鸟团的周一 “肿瘤基因测序数据分析专栏” 在：嫌弃cosmic数据库的30个肿瘤突变signatures? 预告了如果对COSMIC数据库的30个signature不满意，有全新的R包工具可以帮助大家自定义signature，然后我在生信技能树介绍了一个：使用R包SomaticSignatures进行denovo的signature推断，马上收到了我们生信技能树早期元老，王诗翔老铁对他自己开发的工具《sigminer》的强烈推荐，所以立马在这里打call发布《sigminer》的教程哈！如下：

官方文档：https : //shixiangwang.github.io/sigminer-doc/

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

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

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

安装

BiocManager::install("sigminer", dependencies = TRUE)

加载包

library(sigminer)
library(NMF)

数据输入

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

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

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`时会试图生成所有矩阵：

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 ...

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

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

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

show_catalogue(mt_tally$SBS_96 %>% t(), mode = "SBS", style = "cosmic")
show_catalogue(mt_tally$SBS_6 %>% t(), mode = "SBS", style = "cosmic")

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

估计签名数

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

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。报错就加一个伪小数，不报错就可以不管。

show_sig_number_survey2(est$survey)

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

提取签名和可视化

sigs <- sig_extract(mt_tally$SBS_96, n_sig = 4, nrun = 2, pConstant = 1e-9)

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

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也有很多函数专门用户获取一些信息，可视化等等。

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

p <- show_sig_profile(sigs, mode = "SBS", style = "cosmic")
p

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

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.
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.
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() 实现。该算法会生成更大的稀疏（相互之间相互）的签名，因此偏向于生成更多的从我多年研究签名的经验来看，它对于单点突变还是非常友好的。

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.

画图方式是完全一样的。

p <- show_sig_profile(sigs2, mode = "SBS", style = "cosmic")
p

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

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

签名活动图谱

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

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

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

show_sig_exposure(sigs, rm_space = TRUE, style = "cosmic")

根据已知的签名

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

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

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

画图也很简单：

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.

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

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/进一步  学习，后续我会根据研究情况进一步开发。当然，读者完全可以基于上面的分析的结果值进行各种个性化分析。

要是觉得不错，欢迎去GitHub https://github.com/ShixiangWang/sigminer 点个Star。阅读原文可以直达哦

号外：如果有读者对肿瘤基因测序数据分析感兴趣，欢迎加入：0元，10小时教学视频直播《跟着百度李彦宏学习肿瘤基因组测序数据分析》