老规矩，先奉上学习资料链接：
这一次学习Genomic distribution章节部分的内容。
突变不是随机分布在整个基因组中。通过MutationalPatterns，您可以可视化突变是如何在整个基因组中分布的。您还可以查看特定的基因组区域，如启动子、CTCF结合位点和转录因子结合位点。在这些区域内，你可以寻找突变的富集/耗尽，你可以寻找它们之间突变谱的差异。
一、 Rainfall plot
降雨图显示突变类型和突变间隔。降雨图可以用来可视化突变沿基因组或染色体子集的分布。
y轴对应于突变与前一个突变的距离，并进行log10变换。图中的下拉框表示突变的集群或“热点”。可以为snv、indeds、DBSs和mbs制作降雨图。
在这个例子中，我们对一个样本的常染色体做了一个降雨图：
library(BSgenome)
ref_genome <- "BSgenome.Hsapiens.UCSC.hg19"
library(ref_genome, character.only = TRUE)
vcf_files <- list.files(system.file("extdata", package = "MutationalPatterns"),
pattern = "sample.vcf", full.names = TRUE)
sample_names <- c("colon1", "colon2", "colon3","intestine1", "intestine2", "intestine3","liver1", "liver2", "liver3")
tissue <- c(rep("colon", 3), rep("intestine", 3), rep("liver", 3))
grl <- read_vcfs_as_granges(vcf_files, sample_names, ref_genome)
## Rainfall plot
# Define autosomal chromosomes
chromosomes <- seqnames(get(ref_genome))[1:22]
# Make a rainfall plot
p <- plot_rainfall(grl[[1]], title = names(grl[1]), chromosomes = chromosomes, cex = 1.5, ylim = 1e+09 )
ggsave(filename = paste0(opt$od,"/rainfall.png"), width = 8, height = 5, plot = p)
二、Define genomic regions
要查看特定类型的基因组区域，首先需要在一个名为GRangesList的列表中定义它们。您可以使用自己的基因组区域定义(例如基于ChipSeq实验)，也可以使用公开可用的基因组注释数据，如下例所示。
下面的示例展示了如何使用Biocpkg(“biomaRt”)从Ensembl下载基因组构建hg19的启动子、CTCF结合位点和转录因子结合位点区域。有关其他数据集，请参阅biomaRt文档。
下载基因组区域：
注意:这里我们通过加载示例数据的结果采取了一些捷径。下载这些数据的相应代码可以在我们运行的命令上面找到：
library(biomaRt)
# regulatory <- useEnsembl(biomart="regulation",
# dataset="hsapiens_regulatory_feature",
# GRCh = 37)
## Download the regulatory CTCF binding sites and convert them to
## a GRanges object.
# CTCF <- getBM(attributes = c('chromosome_name',
# 'chromosome_start',
# 'chromosome_end',
# 'feature_type_name'),
# filters = "regulatory_feature_type_name",
# values = "CTCF Binding Site",
# mart = regulatory)
#
# CTCF_g <- reduce(GRanges(CTCF$chromosome_name,
# IRanges(CTCF$chromosome_start,
# CTCF$chromosome_end)))
CTCF_g <- readRDS(system.file("states/CTCF_g_data.rds", package = "MutationalPatterns"))
CTCF_g
## Download the promoter regions and convert them to a GRanges object.
# promoter = getBM(attributes = c('chromosome_name', 'chromosome_start',
# 'chromosome_end', 'feature_type_name'),
# filters = "regulatory_feature_type_name",
# values = "Promoter",
# mart = regulatory)
# promoter_g = reduce(GRanges(promoter$chromosome_name,
# IRanges(promoter$chromosome_start,
# promoter$chromosome_end)))
promoter_g <- readRDS(system.file("states/promoter_g_data.rds",package = "MutationalPatterns"))
promoter_g
## Download the promoter flanking regions and convert them to a GRanges object.
# flanking = getBM(attributes = c('chromosome_name',
# 'chromosome_start',
# 'chromosome_end',
# 'feature_type_name'),
# filters = "regulatory_feature_type_name",
# values = "Promoter Flanking Region",
# mart = regulatory)
# flanking_g = reduce(GRanges(
# flanking$chromosome_name,
# IRanges(flanking$chromosome_start,
# flanking$chromosome_end)))
flanking_g <- readRDS(system.file("states/promoter_flanking_g_data.rds", package = "MutationalPatterns"))
flanking_g
将所有基因组区域(GRanges对象)组合在一个命名的GrangesList中
regions <- GRangesList(promoter_g, flanking_g, CTCF_g)
names(regions) <- c("Promoter", "Promoter flanking", "CTCF")
regions
确保这些区域使用与突变数据相同的染色体命名约定:
seqlevelsStyle(regions) <- "UCSC"
三、 Enrichment or depletion of mutations in genomic regions
有必要在每个样本的分析中包括一个区域范围的列表，即基因组中你有足够高质量的reads来支持识别的突变。这可以使用例如GATK的CallableLoci来确定。
如果你不把调查区域包括在你的分析中，你可能会看到一个特定基因组区域的突变减少，这仅仅是该区域低覆盖率的结果，因此并不代表实际的突变减少。
我们在包中提供了一个示例调查区域数据文件。为简单起见，这里我们对每个示例使用相同的调查文件。为了进行正确的分析，确定每个样本的调查区域并在分析中使用这些区域。
加载示例调查区域数据：
## Get the filename with surveyed/callable regions
surveyed_file <- system.file("extdata/callableloci-sample.bed", package = "MutationalPatterns")
## Import the file using rtracklayer and use the UCSC naming standard
library(rtracklayer)
surveyed <- import(surveyed_file)
seqlevelsStyle(surveyed) <- "UCSC"
## For this example we use the same surveyed file for each sample.
surveyed_list <- rep(list(surveyed), 9)
首先，您需要计算每个样本的每个基因组区域中观察到的突变数量和预期突变数量。
distr <- genomic_distribution(grl, surveyed_list, regions)
接下来，可以使用双侧二项检验来测试定义的基因组区域中突变的富集或减少。在这个测试中，观察到突变的概率被计算为突变的总数除以被调查碱基的总数，并进行多次测试校正。fdr和p值的显著性截止值可以用与strand_bias_test相同的方式更改。在本例中，我们按组织类型执行富集/减少试验。
distr_test <- enrichment_depletion_test(distr, by = tissue)
head(distr_test)
最后，您可以绘制结果。星号表示显著的突变富集/减少。这里我们使用p值来绘制星号。
p <- plot_enrichment_depletion(distr_test, sig_type = "p")
ggsave(filename = paste0(opt$od,"/enrichment_depletion.png"), width = 8, height = 5, plot = p)
四、Mutational patterns of genomic regions
1）根据基因组区域分裂突变
你也可以看看基因组区域的突变模式。然而，请记住，突变很少的区域将导致不太可靠的结果。
首先，可以根据定义的基因组区域拆分包含突变的GRangesList。
grl_region <- split_muts_region(grl, regions)
names(grl_region)
现在，您可以将这些样本/区域组合视为完全独立的样本。例如，你可以对这些进行NMF，试图识别特定于某些基因组区域的特征。
mut_mat_region <- mut_matrix(grl_region, ref_genome)
nmf_res_region <- extract_signatures(mut_mat_region, rank = 2, nrun = 10, single_core = TRUE)
signatures = get_known_signatures()
nmf_res_region <- rename_nmf_signatures(nmf_res_region, signatures, cutoff = 0.85)
p <- plot_contribution_heatmap(nmf_res_region$contribution, cluster_samples = TRUE, cluster_sigs = TRUE)
ggsave(filename = paste0(opt$od,"/enrichment_depletion_nmf.png"), width = 8, height = 5, plot = p)
2）Mutation Spectrum
也可以使用plot_spectrum_region函数绘制每个基因组区域的谱，而不是将样品/区域组合作为单独的样本处理。plot_spectrum的参数也可用于此函数。默认情况下，y轴表示变异数除以该样本和基因组区域中的变异总数。这样，突变很少的区域的谱可以更容易地与突变多的区域进行比较。
type_occurrences_region <- mut_type_occurrences(grl_region, ref_genome)
p <- plot_spectrum_region(type_occurrences_region)
ggsave(filename = paste0(opt$od,"/enrichment_depletion_spectrum.png"), width = 8, height = 5, plot = p)
还可以在y轴上绘制变异数除以该样本中变异总数的图。在这种情况下，你不会对每个基因组区域的变异数量进行标准化。如下图所示，本例中绝大多数突变发生在“其他”区域。
p <- plot_spectrum_region(type_occurrences_region, mode = "relative_sample")
ggsave(filename = paste0(opt$od,"/enrichment_depletion_spectrum.relative_sample.png"), width = 8, height = 5, plot = p)
3）Mutational profiles
除了绘制spectra谱外，还可以绘制突变profile。要做到这一点，首先需要制作一个“长”突变矩阵。在这个矩阵中，不同的基因组区域被认为是不同的突变类型，而不是像以前那样作为不同的样本。
mut_mat_region <- mut_matrix(grl_region, ref_genome)
mut_mat_long <- lengthen_mut_matrix(mut_mat_region)
mut_mat_long[1:5, 1:5]
现在可以使用plot_profile_region绘制它。plot_96_profile的参数也可用于此函数。y轴的选项与plot_spectrum_region相同。然而，默认情况下，没有对每个基因组区域的变体数量进行标准化，因为每种突变类型的突变数量通常是有限的。
p <- plot_profile_region(mut_mat_long[, c(1, 4, 7)])
ggsave(filename = paste0(opt$od,"/enrichment_depletion_profile.png"), width = 8, height = 5, plot = p)
4）Mutation density
在上面的例子中，我们使用了已知的功能，如区域的启动子。也可以根据突变密度来定义区域。你可以把基因组分成3个不同突变密度的箱子，像这样：
regions_density <- bin_mutation_density(grl, ref_genome, nrbins = 3)
names(regions_density) <- c("Low", "Medium", "High")
这些区域可以像前面的区域一样使用。例如，将kataegis区域的spectrum与基因组其余部分的spectrum进行比较，这可能是有用的。
grl_region_dens <- split_muts_region(grl, regions_density, include_other = FALSE)
五、无监督的局部突变模式
区域突变模式也可以使用带有determine_regional_similarity函数的无监督方法进行研究。该函数使用滑动窗口方法来计算全局突变谱和较小基因组窗口的突变谱之间的余弦相似度，从而允许对具有突变谱的区域进行无偏识别，这些区域与基因组的其余部分不同。由于这个函数的无偏方法，它在包含至少100,000个substitutions的大型数据集上工作得最好。
首先我们把所有的样品组合在一起。通常，您只会对来自相同癌症类型/组织的样本进行此操作，但由于示例数据中的substitutions数量有限，因此在这里我们将所有内容结合起来。
gr = unlist(grl)
接下来，与基因组其他部分不同的突变模式区域被识别出来。这里我们使用一个较小的窗口大小，因为示例数据的大小较小。在实践中，窗口大小为100或更多效果更好。
regional_sims <- determine_regional_similarity(gr,
ref_genome,
chromosomes = c("chr1", "chr2", "chr3", "chr4", "chr5", "chr6"),
window_size = 40,
stepsize = 10,
max_window_size_gen = 40000000 )
determine_regional_similarity的结果可以可视化。每个点表示单个窗口的突变谱与基因组其余部分之间的余弦相似度。具有不同突变谱的区域将具有较低的余弦相似度。这些点是根据窗户的大小来着色的。这个大小是窗口中第一个和最后一个突变之间的距离。
p <- plot_regional_similarity(regional_sims)
ggsave(filename = paste0(opt$od,"/regional_similarity.png"), width = 8, height = 5, plot = p)
这个包的内容很多，这次学习只能粗略过一遍了，后面还有一些，下次见~
