介绍

背景介绍

细胞的基因组不断受到内源性和环境性dna损伤的威胁，例如紫外线和自发反应。为了维护它们的基因组完整性，细胞利用各种机制修复受损的dna。无论是在复制之前不正确地修复还是未修复，这些都会导致突变被整合到基因组中。每一个突变过程都会留下一个不同的基因组标记。例如，紫外光优先诱导cc>tt。二核苷酸替换，而5-甲基胞嘧啶的自发脱氨导致cpg位点的c>t取代，因此，突变模式可以用来推断哪些变异。

在过去的几年里，对不同人类癌症类型的肿瘤基因组数据的大规模分析揭示了30种突变模式，这些所谓的“突变信号”的特征是碱基替代类型的特定贡献。在一定的序列背景下，每个突变特征被认为反映了单一的突变机制。然而，大多数突变特征的病因目前尚不清楚。为了在功能上将突变特征与生物过程联系起来，评估这些突变特征在暴露于特定诱变剂或细胞的细胞中的贡献，Mutational Pattern的R包提供了一套易于使用的工具集，用于在肿瘤样本或DNA修复缺陷细胞的碱基替换目录中描述和可视化突变模式。软件包涵盖广泛的模式，包括：突变特征、转录链偏倚、基因组分布和与基因组特征的关联，这对于研究突变过程的活动具有共同的意义。重新提取突变特征，并推断先前识别的突变特征的贡献。

方法介绍

该软件包包涵：

(1)新突变特征的提取
(2)对用户指定的突变特征的贡献进行了量化

虽然第一种方法可以用于识别新的突变特征，但这只是有意义的。对于具有大量突变谱样本的数据集，由于它依赖于非负矩阵分解的降维方法。第二种方法可以用来研究单个样本中的突变过程，并通过评估它们在不同系统或不同条件下的贡献来进一步表征先前识别的突变特征。用于探讨其他类型的模式，如转录链不对称、基因组分布以及与染色质组织等(可公开获得的)注释的关联。这些特征对于识别活跃的突变诱导过程和参与特定的DNA修复途径。例如，基因区域存在转录链偏差，这可能意味着活性。

任何一组基本替换调用都可以从VCF文件中导入，基因组构建一个突变矩阵，计数所有96个可能的三核苷酸变化。此外，还包括转录链等其他特征，形成192个特征计数矩阵(96个三核苷酸*2个链)。为此，可以从ucsc中检索到的基因定义用于确定基因中的碱基替换是位于转录的链上还是位于未转录的链上。

下载安装

下载地址： https://github.com/CuppenResearch/MutationalPatterns

数据

要执行突变模式分析，需要加载一个或多个vcf文件，其中包含单核苷酸变异调用和相应的参考基因组。

列出参考基因组

 library(BSgenome)
 head(available.genomes())
[1] "BSgenome.Alyrata.JGI.v1" "BSgenome.Amellifera.BeeBase.assembly4"
[3] "BSgenome.Amellifera.UCSC.apiMel2" "BSgenome.Amellifera.UCSC.apiMel2.masked"
[5] "BSgenome.Athaliana.TAIR.04232008" "BSgenome.Athaliana.TAIR.TAIR9"

#Download and load your reference genome of interest
ref_genome <- "BSgenome.Hsapiens.UCSC.hg19"
library(ref_genome, character.only = TRUE)

加载样本数据

library(MutationalPatterns)
vcf_files <- list.files(path="./data",pattern = ".samtools.snp.reformated.vcf", full.names = TRUE)
sample_names <- c( "YDY019_OA","YDY019_PC","YDY022_OA","YDY022_PC","YDY069_OA","YDY069_PC","YDY106_OA", "YDY106_PC","YDY124_OA","YDY124_PC","YDY125_OA","YDY125_PC")
vcfs <- read_vcfs_as_granges(vcf_files, sample_names, ref_genome)

#定义样本上的相关元数据
tissue <- c("YDY019_OA","YDY019_PC","YDY022_OA","YDY022_PC","YDY069_OA","YDY069_PC","YDY106_OA", "YDY106_PC","YDY124_OA","YDY124_PC","YDY125_OA","YDY125_PC")

画图

突变谱显示了碱基替换目录中每个突变类型的相对贡献。图的谱函数绘制了6个碱基替换类型中的每一个在所有样品上的平均相对贡献。误差条表示所有样品的标准偏差。指示突变的总数

type_occurrences <- mut_type_occurrences(vcfs, ref_genome)

p1 <- plot_spectrum(type_occurrences)
p2 <- plot_spectrum(type_occurrences, CT = TRUE)
p3 <- plot_spectrum(type_occurrences, CT = TRUE, legend = FALSE)
library("gridExtra")
grid.arrange(p1, p2, p3, ncol=3, widths=c(3,3,1.75))

划分每个样本组，例如分别绘制每个组织的光谱

p4 <- plot _ spectrum(type _ occurrences, by = tissue, CT = TRUE, legend = TRUE)
#自定义颜色 
palette <- c("pink", "orange", "blue", "lightblue", "green", "red", "purple")
p5 <- plot _ spectrum(type _ occurrences, CT=TRUE, legend=TRUE, colors=palette)
grid.arrange(p4, p5, ncol=2, widths=c(4,2.3))

Mutational signatures

Mutational signatures突变特征被认为代表了突变过程，其特征是96种碱基替换类型对某一序列的特定贡献。突变特征可以从你的突变计数矩阵中提取出来，并使用非负矩阵因式分解(Nmf)。nmf中的一个关键参数是因式分解秩，即突变特征的数量。使用nmf包确定最优的因式分解等级。

mut_mat <- mut_mat + 0.0001
estimate <- nmf(mut_mat, rank=2:5, method="brunet", nrun=10, seed=123456)
plot(estimate)

使用extract _ signatures从具有ExtractSignals的突变计数矩阵中提取2个突变特征

#rank值指定特征数量
#对于较大的数据集，通过更改nrun参数以实现稳定性和避免局部极小值来执行更多的迭代是明智的
nmf_res <- extract_signatures(mut_mat, rank = 2, nrun = 10)
colnames(nmf_res$signatures) <- c("Signature A", "Signature B")
rownames(nmf_res$contribution) <- c("Signature A", "Signature B")
plot_96_profile(nmf_res$signatures, condensed = TRUE)
pc1 <- plot_contribution(nmf_res$contribution, nmf_res$signature,mode = "relative")
pc2 <- plot_contribution(nmf_res$contribution, nmf_res$signature,mode = "absolute")
grid.arrange(pc1, pc2)

#X和Y轴翻转
plot_contribution(nmf_res$contribution, nmf_res$signature,mode = "absolute", coord_flip = TRUE)

每个样本的每一个特征的相对贡献也可以被绘制为一个热图，它可能比堆叠的树刺图更容易解释和比较。这些样本可以根据它们的欧几里得dis-tance进行分层聚类。这些特征可以按照用户指定的顺序绘制。

#将特征贡献绘制为具有样本聚类树状图和指定特征顺序的热图
pch1 <-plot_contribution_heatmap(nmf_res$contribution,sig_order = c("Signature B", "Signature A"))
pch2 <- plot_contribution_heatmap(nmf_res$contribution, cluster_samples=FALSE)
grid.arrange(pch1, pch2, ncol = 2, widths = c(2,1.6))

将重构的突变剖面与原始突变剖面进行比较

plot_compare_profiles(mut_mat[,1],nmf_res$reconstructed[,1],profile_names = c("Original", "Reconstructed"),condensed = TRUE)

根据COSMIC特征与平均链接的相似性对COSMIC特征进行分级聚类

sp_url <- paste("http://cancer.sanger.ac.uk/cancergenome/assets/","signatures_probabilities.txt", sep = "")
cancer_signatures = read.table(sp_url, sep = "\t", header = TRUE)
#将突变类型的顺序与变异模式标准相匹配
new_order = match(row.names(mut_mat), cancer_signatures$Somatic.Mutation.Type)
# Reorder cancer signatures dataframe
cancer_signatures = cancer_signatures[as.vector(new_order),]
# Add trinucletiode changes names as row.names
row.names(cancer_signatures) = cancer_signatures$Somatic.Mutation.Type
# Keep only 96 contributions of the signatures in matrix
cancer_signatures = as.matrix(cancer_signatures[,4:33])
hclust_cosmic = cluster_signatures(cancer_signatures, method = "average")
# store signatures in new order
cosmic_order = colnames(cancer_signatures)[hclust_cosmic$order]
plot(hclust_cosmic)

计算突变剖面与 COSMIC特征之间的成对余弦相似性

cos_sim_samples_signatures = cos_sim_matrix(mut_mat, cancer_signatures)
# Plot heatmap with specified signature order
plot_cosine_heatmap(cos_sim_samples_signatures,col_order = cosmic_order,cluster_rows = TRUE)

除了重新提取特征外，还可以量化任何一组特征对样本突变轮廓的贡献。这种独特的特性特别适用于小群体或单个样本的突变特征分析，但也可以将自己的发现与已知的签名和已发表的发现联系起来。FIT_to_Signals函数可以找到突变签名的最佳线性组合，这是大多数突变签名的最佳线性组合。通过求解一个非负最小二乘约束问题来构造变异矩阵。

fit_res <-fit_to_signatures(mut_mat, cancer_signatures)

# Select signatures with some contribution
select <- which(rowSums(fit_res$contribution) > 10)
# Plot contribution barplot
plot_contribution(fit_res$contribution[select,],cancer_signatures[,select],coord_flip = FALSE,mode = "absolute")

用样本聚类绘制样本中癌症特征的相对贡献图

plot_contribution_heatmap(fit_res$contribution,cluster_samples = TRUE,method = "complete")

参考

http://bioconductor.org/packages/release/bioc/vignettes/MutationalPatterns/inst/doc/Introduction_to_MutationalPatterns.pdf
https://www.biorxiv.org/content/biorxiv/early/2016/08/30/071761.full.pdf

转载请注明出处
简书作者：ODDXIX

微信公众号：oddxix