关于sig fit 和 sig denovo的问题

[KirsieMin]

大神好！最近进行sigminer包的相关操作时，又发现了一个问题：我用某PCAWG样本进行sig_fit_bootstrap（n=100）后得到的optimal exposure与 官网该样本的sig denovo exposure 完全不一致==，望能解答~ fit数据来源：https://xenabrowser.net/datapages/?dataset=October_2016_all_patients_2778.snv_mnv_indel.maf.coding.xena&host=https%3A%2F%2Fpcawg.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443 fit代码： 选取SP10084：

[KirsieMin]

补充：denovo exposure来源：https://dcc.icgc.org/releases/PCAWG/mutational_signatures/Signatures_in_Samples/SP_Signatures_in_Samples

[ShixiangWang]

应该是你这里的数据不一致导致的，Xena 的数据只包含 coding 区域，要少很多的。

download maf file, then cut -f 1,2,3,4,6,8,10,13,29,38,39 October_2016_all_patients_2778.snv_mnv_indel.maf, mapped to specimen ids, averaged VAF (MNVs have multiple VAF, one for each base call), selected only coding mutations, about 0.8% of all mutations calls.

如果你想要进一步验证的话，可以下载原始的文件 https://dcc.icgc.org/api/v1/download?fn=/PCAWG/consensus_snv_indel/final_consensus_passonly.snv_mnv_indel.icgc.public.maf.gz

然后选择对应的样本进行查看。

---

另外可以使用模拟数据确定这个函数的有效性：

library(sigminer)

data("simulated_catalogs")
mat = t(simulated_catalogs$set1) # 由10个COSMIC signature随机组合模拟生成的矩阵

sig = sig_extract(mat, n_sig = 10, nrun = 30, cores = 4)
sig$Exposure

fit_expo = sig_fit(t(mat), sig = sig$Signature)

fit_expo

for (i in 1:ncol(fit_expo)) {
  print(cor(sig$Exposure[, i], fit_expo[, i]))
}

注：sig_fit 是 sig_fit_bootstrap_*提取optimal exposure的核心。

输出：

> for (i in 1:ncol(fit_expo)) {
+   print(cor(sig$Exposure[, i], fit_expo[, i]))
+ }
[1] 0.9999712
[1] 0.9999382
[1] 0.9999875
[1] 0.9999996
[1] 0.9999992
[1] 0.9990981
[1] 0.9999549
[1] 0.9997031
[1] 0.9999474
[1] 0.9999878
[1] 0.9999834
[1] 0.9999935
[1] 0.9999623
[1] 0.9999662
[1] 0.9999982
[1] 0.9999883
[1] 0.9999464
[1] 0.9999765
[1] 0.9999803
[1] 0.9999906
[1] 0.999859
[1] 0.9999982
[1] 0.9999897
[1] 0.9999991
[1] 0.9999441
[1] 0.9999846
[1] 0.999997
[1] 0.9999212
[1] 0.9998311
[1] 0.999982

[ShixiangWang]

@KirsieMin 有重复和解决吗？

[ShixiangWang]

我关闭了。如果有问题请重新开启。

[KirsieMin]

大神新年好！之前在假期所以咩有及时回复感到抱歉！目前模拟数据已经成功尝试，但是根据大神提供的原始数据进行sig_tally的时候遇到之前未知的报错： 代码： 报错：

这边进行debug的时候发现具体出错语句是：

所以我进行了以下检查，发现并没有不妥之处：

故向您请教，希望解决该问题～

[ShixiangWang]

@KirsieMin 你能否尝试一部分少量数据是否会报错？

这个如果输入数据没有问题的话，原因可能是基因组版本不完全匹配。例如，有一个突变是chr1:100，而基因组chr1是从200开始的。

我猜测可能是下面这段代码哪里出问题了

https://github.com/ShixiangWang/sigminer/blob/d40168821143b5c7174702324d5ef4a61f5dbe29/R/generate_matrices.R#L21-L33

你能否具体看下是哪一行报错，内部用for循环找出一个具体报错的变异？这样我这边再看看能否调试和解决下。

[KirsieMin]

目前尝试了随机取了十个样本的突变数据以及之前报错样本的突变数据进行了重新tally，发现成功运行。==说明不是报错数据的问题，这还需要进行基因组版本的测试嘛？或许可以尝试将基因组版本换成非UCSC的版本进行尝试？

[ShixiangWang]

@KirsieMin 最细致地处理是使用icgc/文献使用的参考基因组文件。

[KirsieMin]

目前找到icgc参考基因组信息网页如下： https://dcc.icgc.org/releases/PCAWG/reference_data/pcawg-bwa-mem sig_tally目前貌似只可以直接使用BSgenome中的参考基因组进行tally，我可以将sig_tally中ref_genome输入为网页的参考基因组直接tally吗？

[ShixiangWang]

@KirsieMin 目前是不行，只能尝试自建有一个BSGenome对象（包）看看。

• https://www.jianshu.com/p/1fa3c7eb850e
• https://bioconductor.org/packages/devel/bioc/vignettes/BSgenome/inst/doc/BSgenomeForge.pdf

后面我有空看看能不能直接支持fasta文件作为输入。

[KirsieMin]

大神好！目前进行了两方面的尝试： 1.自建pcawg的BSGenome对象。在构建对象的时候需要建立一个seed file，之前获得的信息pcawg用的ref_genome是个组装genome，所以需要填写参数： 在这里我突然意识到pcawg的ref genome可能已经在BSGenome中了，由此我进行了探索发现，icgc参考基因组信息网页中：https://dcc.icgc.org/releases/PCAWG/reference_data/pcawg-bwa-mem/中的参考基因组是hs37d5： 网页描述与搜索结果相符，故选择BSgenome.Hsapiens.1000genomes.hs37d5进行tally，但又出现了之前一样的报错： 所以我理解为并非ref genome的问题导致的报错==，这一方面暂无其他思路，希望得到帮助~

2.我将之前的报错样本突变单独tally（并无报错），与其他样本tally数据合并，继续进行fit。对比denovo sigs和使用三种不同方法（SA还在处理）的fit sigs ，使用了包中cosine进行了重构矩阵相似性比较： 不同方法之间只有SBS和CNS表现出明显差异，DBS和ID则没有；另外可以看到ID的fit很差，SBS、DBS的fit较denovo也不太一致。这里判断SBS、DBS、ID的fit较denovo不一致的原因可能还是上一方面tally中的问题。

[ShixiangWang]

关于第1点，我目前稍微有点糊涂。如果你报错的样本sig_tally成功的话，但全部数据又报错，那很可能报错的样本找错了。 目前的错误是输入数据和参考数据的区间对不上，存在2种可能的原因：

1. 参考基因组版本还是不正确，有2种可能：一个是版本是对的，但你的输入数据里包含了参考基因组版本中没有的contigs，这种contigs的区间可能对不上。另一个可能是本身数据的基因组版本信息存在问题，你简单试下hg38是不是会成功？
2. tally中取区间存在问题，不过这个code是从maftools拉过来的，本身没有问题，但有可能取区间数据的时候因为加了上下游+/-20bp，所以出界了。你看得到报错前正在进行提取context的区间数据（maftools很多issue提到类似的报错）。

关于第2点我看到的是本身跟这个问题没有关系。这里你说明的是不同的方法结果不太一样，我不清楚你的操作，但这种对比方法发现结果不一样是合理的。而这里的问题最核心的关注点在于你最开始说的：

用某PCAWG样本进行sig_fit_bootstrap（n=100）后得到的optimal exposure与 官网该样本的sig denovo exposure 完全不一致

我还不清楚你是否确定该问题存在，还是说可以重复，已经解决了？

再回到第一个问题，你能否抽取能重复该问题的最小数据（比如一个样本的突变数据，你用它运行sig_tally，发现能够重复第一个问题的报错），然后放到这里或者邮件给我w_shixiang@163.com，我这边有空可以仔细检查和调试下源代码。

[KirsieMin]

1.我尝试抽取了该问题的最小数据（1，10，20，100，500，1000，1500，1950）发现1950个样本（全部数据）的时候才报错==，目前找的报错样本依旧是"a08ec059-7592-4698-bb45-25a9c3680c23"。样本数据为https://www.jianguoyun.com/p/DfZ5xP8Qs6rSCBiM6K8E(有效期为五天)。另外我同样尝试了hg38，发现其他报错： 2.目前我想通过cosine的计算来看denovo sigs和fit sigs的exposure是否一致，从上面的结果来看，我分别将【denovo sigs的exposure和fit exposure】 与 【相同的cosmic sigs 和 sigminer CNS】相乘得到的 两个重构矩阵进行cosine计算获得similarity，通过对similarity的值来判断两者的一致性，发现了SBS、ID、DBS的similarity都不是百分之百，说明还是不可重复。所以问题还是没有得到解决。

[ShixiangWang]

第1点我有空深入研究下再回你哈，最近有点忙，这周不一定能弄好。

关于第2点，我觉得你的处理可能有点问题。Similarity是用来看不同signature（组成）的相似情况，如果你想要计算exposure的一致性，在固定好signature，计算的exposure可以算残差或均方误差来对比结果的一致性。

[KirsieMin]

好的谢谢大神！你先忙你的事情，我这边也不急！第一个问题你可以先下载好数据放着备用，第二个问题我继续处理一下，处理之后给你反馈下结果，感谢！~

[KirsieMin]

大神好！我这边已经处理好了pcawg 相同sigs的denovo exposure和fit exposure的mse，结果确实是SBS，DBS，ID的fit结果很差，但是cns的挺好的，目前还是无法判断是基因组问题还是包内的问题，这边我下一步将使用sigprofile进行fit尝试看看效果：

[ShixiangWang]

再开方算RMSE好一点，可以直接得到平均多少个突变偏移了。

第一个问题我这几天有时间处理了，搞完再回复你。

[ShixiangWang]

通过下面的代码发现 b42d183c-bc9c-4652-9e56-10c54c5ee96e 这个样本出错了。调试发现为了取上游位点做减法，最小值小于0了。我不知道

library(sigminer)

obj = readRDS("~/../Downloads/pcawg_sbs_maf.rds")

# IDs = maftools::getSampleSummary(obj)$Tumor_Sample_Barcode
# for (i in IDs) {
#   test = maftools::subsetMaf(obj, i)
#   cli::cli_alert_success("process sample {i}")
#   obj_tally = suppressMessages(
#     sig_tally(test, mode = "SBS",
#               ref_genome = "BSgenome.Hsapiens.1000genomes.hs37d5",
#               genome_build = "hg19")
#   )
# }

# b42d183c-bc9c-4652-9e56-10c54c5ee96e
test = maftools::subsetMaf(obj, "b42d183c-bc9c-4652-9e56-10c54c5ee96e")
obj_tally = sig_tally(test, mode = "SBS",
            ref_genome = "BSgenome.Hsapiens.1000genomes.hs37d5",
            genome_build = "hg19")

#debug(sig_tally)

obj_tally = sig_tally(test, mode = "ALL",
                      ref_genome = "BSgenome.Hsapiens.1000genomes.hs37d5",
                      genome_build = "hg19")

重新安装github版本后就不会出现这个问题了。我目前的机器比较小，运行全部数据内存会爆掉，你有空再试试全部的数据tally还是出现问题不：

obj_tally = sig_tally(obj, mode = "ALL",
                      ref_genome = "BSgenome.Hsapiens.1000genomes.hs37d5",
                      genome_build = "hg19")

[ShixiangWang]

@KirsieMin

[ShixiangWang]

@KirsieMin 谢谢赞赏😉 不知道上面的问题是否已经解决呢？

[KirsieMin]

客气了~ 目前安装最新的sigminer以后 全部数据tally已经成功运行！之后我将根据最近tally数据进行fit操作看RMSE的情况，完成后再给大神反馈！

[ShixiangWang]

可以的。

关于 PCAWG 官网的 signature exposure结果对比，你可以随机抽十个样本的数据，然后用选一个fitting方法简单看下SBS exposure，对比下官网提供的结果，是不是非常接近。这样判定就可以了。如果你不是在做一个系统的对比研究，不需要处理全部的样本进行分析。

[ShixiangWang]

@KirsieMin 谢谢你这个 issue 的反馈，最新的修改版本接下来会推送到CRAN

[KirsieMin]

大神好~ 我最近处理了根据原始基因组hs37d5 tally的数据，进行了系统的fitting分析，发现结果比之前稍微好了一点点，但sbs和id的结果仍然较差，平均突变偏移都很高，尤其是sbs==。这里考虑是fitting方法本身的原因导致的，想请教一下这里是否还有可以优化fitting结果的方法？还请答疑解惑~

[ShixiangWang]

这个RMSE明显跟检出的变异记录数目成正比。

如果已知肿瘤类型，参考COSMIC中不同肿瘤类型的signature分布图谱，在fitting的时候指定用哪些signature，可能可以减少误差值。

你这使用的是WGS的数据？平均样本的变异有几万个吗？

[KirsieMin]

关于变异记录数目比较的话，我这里使用了pcawg （WGS）原始的数据预处理好的maf文件，分别对SNP，DNP，Indel（DEL+INS），CNV，以及总的突变（SNP+DNP+Indel+CNV） 这五类 进行了相应样本的平均突变数目计算，得到了下图：结果与RMSE的结果貌似并不成正比==，平均样本的突变 在这五类中均在200以内，未有上万==，所以目前对RMSE的结果抱有疑惑🤔。 另外对于已知肿瘤类型的fitting方法，略读过该文章。也就是说，在未知肿瘤类型的情况下 目前未有合适的方法对fitting结果进行优化对吗？

[ShixiangWang]

如果你一个样本平均的SBS才不到100，那么计算出来的RMSE怎么会达到800呢。你看下RMSE的公式，如果SBS平均才60的样子，它不会超过60这个范围的。

感觉你应该RMSE算的有很大问题喔

[ShixiangWang]

下面的示例代码供你参考：

library(sigminer)

?get_sig_exposure

# Load mutational signature
load(system.file("extdata", "toy_mutational_signature.RData",
                 package = "sigminer", mustWork = TRUE
))
# Get signature exposure
expo1 <- get_sig_exposure(sig2)

expo_mat = as.matrix(expo1[, -1])
# assume the actual data is 
expo_mat2 = round(expo_mat)

# The row are samples/cases, columns are signatures

# If you want to summarize by samples
sqrt(rowMeans((expo_mat - expo_mat2)^2))

# If you want to summarize by signatures
sqrt(colMeans((expo_mat - expo_mat2)^2))

[KirsieMin]

我检查了之前的代码和数据，我将所有fitting exposure（NNLS、QP、SA）和denovo exposure 根据sample名和sig名排序放在了同一个dataframe里： 之后根据RMSE的算法进行的计算： 代码如下： sqrt(mean((sbsExposure$denovo_exposure-sbsExposure$nnls_exposure)^2)) -> sbsnnls_rmse sqrt(mean((sbsExposure$denovo_exposure-sbsExposure$qp_exposure)^2)) -> sbsqp_rmse sqrt(mean((sbsExposure$denovo_exposure-sbsExposure$sa_exposure)^2))-> sbssa_rmse 经过重新检查运算，sbs的rmse值还是在700+......所有这里我考虑的问题是fitting结果的过于不一致性所导致的差距。

[ShixiangWang]

@KirsieMin 你有空把数据表格发我看一看，我看下数据的分布情况怎么样。

[KirsieMin]

好的大神！ 数据如下：https://www.jianguoyun.com/p/DZmTTfMQs6rSCBj1g7QE (访问密码：pHvZMP)

[ShixiangWang]

sbsExposure
summary(sbsExposure$nnls_exposure)
summary(sbsExposure$denovo_exposure)

# Error rate
rmse_by_sample = sbsExposure[, list(rmse = sqrt(mean((nnls_exposure - denovo_exposure)^2)),
                                    N = pmax(sum(nnls_exposure), sum(denovo_exposure))), by = sample][order(rmse)]
rmse_by_sample[, percent := rmse / N]
hist(rmse_by_sample$rmse, breaks = 100)
hist(rmse_by_sample$percent, breaks = 100)
summary(rmse_by_sample$percent)

需要从比例去分析这个误差。如果单纯从突变数目看，有不少样本存在极大的差异，因此整体的平均误差被拉高了。

如果我们考虑样本本身的总突变数目，发生分配错误的突变数据实际占样本总突变数的10%之内。

因此denovo和fitting的结果存在平均7%（上下5%）的差异。

[ShixiangWang]

大神好~ 我最近处理了根据原始基因组hs37d5 tally的数据，进行了系统的fitting分析，发现结果比之前稍微好了一点点，但sbs和id的结果仍然较差，平均突变偏移都很高，尤其是sbs==。这里考虑是fitting方法本身的原因导致的，想请教一下这里是否还有可以优化fitting结果的方法？还请答疑解惑~

建议这个图也从平均的差异率方面计算更新一下

[KirsieMin]

Thanks♪(･ω･)ﾉ根据平均的差异率计算更新如图： 从这个角度看的话，fitting 与denovo 的差异确实在相对小的范围内。但是我疑惑的一点是为什么sbs中有些样本的错配率极高==，下一步我想分析下错配率高的这些样本具体的信息是什么样的，有没有什么特殊的地方。

[ShixiangWang]

@KirsieMin 确实可以作为一个点去研究下哪些因素对这种错配贡献最大

我能想到的包括

• 测序和变异检测本身的错误
• signature图谱有点存在相似性，导致不正确的signature分配，进而exposure计算存在问题
• signature存在多重效应，本身突变呈现的是 导致突变 和 修复突变 两种作用的结合，因而有一些深入的理论因素没有考虑到
• 可能一些signature存在特殊的性质，不能用NMF类似的算法/或线性分解理论去合理解析

[KirsieMin]

收到~感谢大神的指点！我准备先从临床信息的角度观察一下信息的区别，有问题我会及时反馈！