R语言实现突变信号（Mutational Signatures）分析
R语言实现突变信号（Mutational Signatures）分析
突变信号（Mutational Signatures）首次2013年在《nature》进行报道。并做了相关的定义：细胞在成长过程中，基因组不断受到内源性和外源性DNA损伤的威胁，正是由于这些威胁，使得细胞基因组不断发生变化，并最终发生一些突变的积累。每一个突变过程都会留下一个不同的基因组标记，也就称为突变信号。
今天为大家介绍一个R包MutationalPatterns，可以用于在肿瘤样本或DNA修复缺陷细胞的碱基替换目录中描述和可视化突变模式。包括：突变特征、转录链偏倚、基因组分布和与基因组特征的关联。当然在他之前还有另外一个R包deconstructSigs也可以解决类似的问题，这两个包的联合使用将会把突变信号从个体到群体的状态全部描述清楚。
接下来我们看下他们的安装：
- deconstructSigs安装很简单，它已存在于CRAN上，直接install.packages()就行。
- MutationalPatterns安装相对复杂点，一种是bioconductor安装，不再赘述；另一种需要先安装devtools包，基于这个包利用GitHub源进行安装。首先打开GitHub中包的界面：
安装需要红框中的信息（用户/包名）：
install_github('UMCUGenetics/MutationalPatterns')。至此两个包的安装就结束了，然后我们看下其数据的导入，支持VCF，txt等多格式导入：
1. VCF文件的读取。需要利用MutationalPatterns包中的read_vcfs_as_granges函数进行读取转化：
其中主要的参数：
Vcf_files vcf文件的列表。主要是每一个vcf 的路径信息（eg：colon1-sample.vcf）。
Sample_names 主要是每一vcf文件对应的样本名称（eg: colon1）。
Genome 参考基因组名称（eg：BSgenome.Hsapiens.UCSC.hg19）。
实例：
sample_names <- c ( "colon1","colon2", "colon3","intestine1","intestine2", "intestine3","liver1","liver2", "liver3" )
# We assemble a list of files we want toload. These files match the
# sample names defined above.
vcf_files <-list.files(system.file("extdata", package="MutationalPatterns"),pattern= ".vcf", full.names = TRUE)
# Get a reference genome BSgenome object.
ref_genome <-"BSgenome.Hsapiens.UCSC.hg19"
library("BSgenome")
library(ref_genome, character.only = TRUE)
# This function loads the files as GRangesobjects
vcfs <- read_vcfs_as_granges(vcf_files,sample_names, ref_genome)
2. TXT文件的读取。需要利用deconstructSigs包中的mut.to.sigs.input可以按照一定的结构进行读取数据，同时计算出96种突变情况百分比。
其中主要的参数：
Mut.ref 主要是文件的路径，其中主要包括chr，pos，ref，alt信息的文件。
Bsg 如果参考基因不是BSgenome.Hsapiens.UCSC.hg19，需要自己设定。
实例：
sigs.input = mut.to.sigs.input(mut.ref =sample.mut.ref, sample.id = "Sample",
chr = "chr", pos ="pos", ref = "ref", alt = "alt", bsg = BSgenome.Hsapiens.UCSC.hg19)
3. RDS格式数据读取源自base包readRDS函数，其实就是R语言自己的数据格式，此格式在可以输出到文件。
实例：
vcfs <-readRDS(system.file("states/read_vcfs_as_granges_output.rds",package="MutationalPatterns")) 接下来是突变信号的计算，两个包都有自己的计算函数，但是核心都一样，突变信号数量是一致的。deconstructSigs包中的mut.to.sigs.input可以直接计算输入的数据的突变信号分布。MutationalPatterns 包中mutations_from_vcf也可以对每一个样本进行突变信号的统计。
实例：
vcfs <-readRDS(system.file("states/read_vcfs_as_granges_output.rds",package="MutationalPatterns"))
muts = mutations_from_vcf(vcfs[[1]])
当然MutationalPatterns 包中还提供了更高级的统计方式，可以直接把每一个样本的突变信号一次性都展示出来，需要用到函数mut_type_occurrences（vcf_list,ref_genome）。
实例：
type_occurrences =mut_type_occurrences(vcfs, ref_genome)
这样两个包形成了互补，对输入的不同数据都可以得到所有样本的突变信号。那么接下来三核苷酸的突变信息需要mut_matrix {MutationalPatterns}函数计算。
实例：
mut_mat <- mut_matrix(vcf_list = vcfs,ref_genome = ref_genome)
以上是我们得到的96个突变信号，那么我们如果将染色体的正反链都引入，我们的突变信号数量将会加倍。为此此包作者还设计了其192个信号的查看函数mut_matrix_stranded。
实例：
## Replication strand analysis:
## Read example bed file with replicationdirection annotation
repli_file =system.file("extdata/ReplicationDirectionRegions.bed", package ="MutationalPatterns")
repli_strand = read.table(repli_file,header = TRUE)
repli_strand_granges = GRanges(seqnames =repli_strand$Chr,
ranges =IRanges(start = repli_strand$Start + 1,
end =repli_strand$Stop),
strand_info =repli_strand$Class)
## UCSC seqlevelsstyle
seqlevelsStyle(repli_strand_granges) ="UCSC"
# The levels determine the order in whichthe features
# will be countend and plotted in thedownstream analyses
# You can specify your preferred order ofthe levels:
repli_strand_granges$strand_info =factor(repli_strand_granges$strand_info, levels = c("left","right"))
mut_mat_s_rep = mut_matrix_stranded(vcfs,ref_genome, repli_strand_granges,
mode ="replication")
接下来就是利用NMF模型对上面的分布情况进行进一步的数据标准化，或者说重构。首先我们看下extract_signatures{MutationalPatterns}主要是直接通过上面的96/192突变信号进行数据重构。其中主要的参数rank，指的需要拆分后矩阵的大小。
实例：
nmf_res_s <- extract_signatures(mut_mat_s_rep,rank = 2)
nmf_res <- extract_signatures(mut_mat,rank = 2)
计算后获得NMF 的权重值（contribution），特征矩阵（signatures）以及重构后的矩阵（reconstructed）。
然后，我们重构了我们的突变信号，那么接下来就需要把我们的信号与已经存在数据进行比对，看看每个样本的在别人提出的典型信号之间的关系，就需要用到函数whichSignatures {deconstructSigs}。或者也可以把原始的信号分布数据直接放进去。
主要参数：
Tumor.ref 三核苷酸的分布矩阵行为样本名称，列为三核苷酸信号。
Sample.id 样本的编号名称，也就是行名。
signatures.ref 主要是你要参看的数据集，在deconstructSigs中提供的数据集是signatures.nature2013。
contexts.needed 如果数据是已经标准化的矩阵，那么false；如果是突变信号的分布数那么为TRUE。
tri.counts.method 突变信号分布的数据标准化的方法。默认不标准化。
实例：
test = whichSignatures(tumor.ref =randomly.generated.tumors,sample.id = "2", contexts.needed = FALSE)
接下来就是这些突变信号的可视化，首先我们看下plotSignatures{deconstructSigs}函数，就可以将上面的结果绘制相应的bar图。
实例：
plotSignatures(test, sub ="example")
然后makePie {deconstructSigs}函数可以对获得权重值进行饼图的绘制：
实例：
接下来我们看下在MutationalPatterns包中更加丰富的可视化函数：
1. plot_96_profile
实例：
plot_96_profile(nmf_res$signatures)
2. plot_192_profile
实例：
plot_192_profile(nmf_res_s$signatures)
3. plot_compare_profiles 多组数据之间的对比情况可视化
实例：
plot_compare_profiles(mut_mat[,1],
nmf_res$reconstructed[,1],
profile_names =c("Original", "Reconstructed"))
4. plot_contribution 突变信号贡献度的可视化。
具体不再一个一个参数细讲，我们直接看实例：
plot_contribution(nmf_res$contribution,nmf_res$signature,mode= "relative")
plot_contribution(nmf_res$contribution,nmf_res$signature,mode= "absolute")
plot_contribution(nmf_res$contribution,nmf_res$signature,mode= "absolute",index = c(1,2))
plot_contribution(nmf_res$contribution,nmf_res$signature,mode= "absolute",coord_flip = TRUE)
4. plot_contribution_heatmap 突变信号贡献的热图绘制。
实例：
rownames(nmf_res$contribution) =c("Signature A", "Signature B")
## Define signature order for plotting
sig_order = c("Signature B","Signature A")
## Contribution heatmap with signatures indefined order
plot_contribution_heatmap(nmf_res$contribution,
sig_order =c("Signature B", "Signature A"))
## Contribution heatmap without sampleclustering
plot_contribution_heatmap(nmf_res$contribution,
sig_order =c("Signature B", "Signature A"),
cluster_samples =FALSE, method = "complete")
5. plot_spectrum突变分布的可视化。
实例：
type_occurrences =mut_type_occurrences(vcfs, ref_genome)
## Plot the point mutation spectrum overall samples
plot_spectrum(type_occurrences)
## Or with distinction of C>T at CpGsites
plot_spectrum(type_occurrences, CT = TRUE)
此包还包含其他一些可视化的图形，但是感觉都没啥用徒增感官度而已。就不一一介绍了。
欢迎大家学习交流！