关于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开始的。

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

sigminer/R/generate_matrices.R

Lines 21 to 33 in d401688

	ss <- BSgenome::getSeq(
	x = ref_genome,
	names = extract.tbl$Chromosome,
	start = extract.tbl$Start - 2,
	end = extract.tbl$End + 2,
	as.character = TRUE
	)
	send_info("Extracting +/- 20bp around mutated bases for background C>T estimation.")
	updwn <- BSgenome::getSeq(
	x = ref_genome, names = extract.tbl$Chromosome, start = extract.tbl$upstream,
	end = extract.tbl$downstream, as.character = FALSE
	)

你能否具体看下是哪一行报错，内部用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中的问题。