生信软件工具箱

chromVAR:从单细胞表观数据中推断转录因子相关开放区域

chromVAR是一个用于分析稀疏染色质开放的R包。chromVAR的输入文件包括，ATAC-seq处理后的fragments文件（过滤重复和低质量数据）, DNAse-seq实验结果，以及基因组注释（例如motif位置）

chromVAR先根据所有细胞或者样本的平均情况来计算期望开放性， 然后用它来计算每个注释，每个细胞或样本的偏差，最后对开放进行纠正。

安装加载

安装R包

1
2
3
4
5
BiocManager::install("GreenleafLab/chromVAR")
BiocManager::install("motifmatchr")
BiocManager::install(SummarizedExperiment)
BiocManager::install(BiocParallel)
BiocManager::install("JASPAR2018") # JASPAR 2018数据库

chromVAR的运行需要先加载下面这些R包

1
2
3
4
5
6
library(chromVAR)
library(motifmatchr)
library(Matrix)
library(SummarizedExperiment)
library(BiocParallel)
set.seed(2019)

chromVAR能够使用多核进行并行运算，调用方法如下

1
2
3
4
# 全部用户
register(MulticoreParam(8, progressbar = TRUE))
# Windows用户，MulticoreParam不好用就用SnowParam
register(SnowParam(workers = 1, type = "SOCK"))

注意: 不运行的话，后续代码可能会报错

读取输入

chromVAR接受的输入是落入开放区域的read数统计表。有许多软件可以做到，chromVAR也提供了相应的方法。

首先要提供一个peak文件，文档建议这个peak文件存放的peak为等宽非重叠，建议宽度在250-500 bp之间。peak文件可以利用已有的bulk ATAC-seq或DNAse-seq数据来获取。对于来自于多个样本的peak，需要先用filterPeaks函数保证peak之间不重合。

使用getPeaks读取peak

1
2
peakfile <- system.file("extdata/test_bed.txt", package = "chromVAR")
peaks <- getPeaks(peakfile, sort_peaks = TRUE)

MACS2分析结果里会提供narrowpeak格式文件，chromVAR提供了readNarrowpeaks函数进行读取。

随后用getCounts函数基于BAM文件和加载的peak获取count

1
2
3
4
5
6
7
bamfile <- system.file("extdata/test_RG.bam", package = "chromVAR")
fragment_counts <- getCounts(bamfile, 
                             peaks, 
                             paired =  TRUE, 
                             by_rg = TRUE, 
                             format = "bam", 
                             colData = DataFrame(celltype = "GM"))

bamfile可以有多个bam文件路径，bam文件的类型和colData对应。此外by_rg定义是否要根据BAM文件中的RG对输入分组。

实际演示的时候，我们官方提供的示例数据

1
data(example_counts, package = "chromVAR")

因为这是一个*SummarizedExperiment 对象，因此我们就能用assay,colData和rowData 了解这个数据集

主体是一个存放count的dgCMatrix

1
2
3
4
5
6
7
8
9
10
assay(example_counts)[1:5,1:2]

5 x 2 sparse Matrix of class "dgCMatrix"
     singles-GM12878-140905-1 singles-GM12878-140905-2
[1,]                        .                        1
[2,]                        .                        .
[3,]                        .                        .
[4,]                        .                        .
[5,]                        .                        .

行是特征(feature), 这里就是peak

1
2
3
4
5
6
7
head(rowData(example_counts), n=2)

DataFrame with 2 rows and 3 columns
      score      qval        name
  <integer> <numeric> <character>
1       259  25.99629   GM_peak_6
2        82   8.21494   H1_peak_

列表示样本，

1
2
3
4
5
6
7
head(colData(example_counts), n=2)

DataFrame with 2 rows and 2 columns
                           Cell_Type     depth
                         <character> <integer>
singles-GM12878-140905-1          GM      9220
singles-GM12878-140905-2          GM      9401

前期准备

分析peak中的GC含量

GC含量信息用于确定哪些peak可能是背景， addGCBias返回更新后的SummarizedExperiment 。其中genome支持BSgenome, FaFile和DNAStringSet对象。

1
2
3
4
library(BSgenome.Hsapiens.UCSC.hg19)
example_counts <- addGCBias(example_counts, 
                            genome = BSgenome.Hsapiens.UCSC.hg19)
head(rowData(example_counts))

过滤输入(单细胞)

如果是处理单细胞数据，建议过滤测序深度不够，或者是信噪比低不够(也就是peak中所占read比例低)的数据，主要靠两个参数,min_in_peaks和min_depth.

1
2
counts_filtered <- filterSamples(example_counts, min_depth = 1500, 
                                 min_in_peaks = 0.15, shiny = FALSE)

然后我们用filterSamplesPlot对过滤情况进行可视化

1
2
3
filtering_plot <- filterSamplesPlot(example_counts, min_depth = 1500, 
                                    min_in_peaks = 0.15, use_plotly = FALSE)
filtering_plot

确认标准后，就可以过滤

1
counts_filtered <- filterPeaks(counts_filtered, non_overlapping = TRUE)

获取Motifs和分析包含motif的peak

可以用getJasparMotifs在JASPAR数据库中提取motif信息，但是其中getJasparMotifs只是TFBSTools::getMatrixSet一个简单的封装而已, 默认用的是JASPAR2016, 建议通过下面的代码获取最新版本的数据。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# BiocManager::install("JASPAR2018")
library(JASPAR2018)
species <- "Homo sapiens"
collection <- "CORE"
opts <- list()
opts["species"] <- species
opts["collection"] <- collection
out  <- TFBSTools::getMatrixSet(JASPAR2018, opts)

if (!isTRUE(all.equal(TFBSTools::name(out), names(out)))) 
  names(out) <- paste(names(out), TFBSTools::name(out), 
                      sep = "_")
motif <- out

collection参数中接受: “CORE”, “CNE”, “PHYLOFACTS”, “SPLICE”, “POLII”, “FAM”, “PBM”, “PBM_HOMEO”, “PBM_HLH” 选项。

获取Motif的另外一种方式是用chromVARmotifs, 这个R包的主要功能就是整合了一些motif, 通过devtools::install_github("GreenleafLab/chromVARmotifs")进行安装。

之后用motifmatchr中的macthMotifs去分析peak中motif包含情况。默认会返回一个SummarizedExperiment 对象，就是一个稀疏矩阵，来提示是不是匹配了motif

1
2
3
library(motifmatchr)
motif_ix <- matchMotifs(motifs, counts_filtered, 
                        genome = BSgenome.Hsapiens.UCSC.hg19)

关键参数是p.cutoff用于设置严格度，默认是0.00005其实能够返回比较合理的结果。如果需要一些额外信息，可以通过参数out来调整，例如out=positions表示返回实际的匹配位置

偏离度和变异度分析(Deviations and Variability )

上面都是准备阶段，计算deviations和Variability才是这个软件的重点。

计算偏离度

函数computeDeviations返回的SummarizedExperiment 对象, 包含两个”assays”,

  • 偏差纠正后的开放偏离度:assays(dev)$deviations
  • 偏离度的Z-score: assays(deviations)$z: 考虑到了当从基因组上随机抽取一些GC含量类似的片段分析时，有多大概率会得到该结果。
1
dev <- computeDeviations(object = counts_filtered, annotations = motif_ix)

compuateDeviations支持背景peak的输入，用于的偏移度得分进行标准化，默认会自动计算不会返回结果。你可以选择用getBackgroundPeaks分析背景，然后传递给computeDeviations

1
2
3
bg <- getBackgroundPeaks(object = counts_filtered)
dev <- computeDeviations(object = counts_filtered, annotations = motif_ix,
                         background_peaks = bg)

变异度

computeVariability会返回一个数据框data.frame，里面有变异度, 该变异度的bootstrap置信区间和衡量拒绝原假设的概率(即 > 1)

1
2
variability <- computeVariability(dev)
plotVariability(variability, use_plotly = FALSE)

偏离度可视化

可以利用tSNE可视化偏离度

1
2
3
4
5
6
7
tsne_results <- deviationsTsne(dev, threshold = 1.5, perplexity = 10)
tsne_plots <- plotDeviationsTsne(dev, tsne_results, 
                                 annotation_name = "TEAD3", 
                                 sample_column = "Cell_Type", 
                                 shiny = FALSE)
cowplot::plot_grid(tsne_plots[[1]],tsne_plots[[1]] )

参考资料