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

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

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

安装加载

安装R包

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

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

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

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

  1. # 全部用户
  2. register(MulticoreParam(8, progressbar = TRUE))
  3. # Windows用户,MulticoreParam不好用就用SnowParam
  4. 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. peakfile <- system.file("extdata/test_bed.txt", package = "chromVAR")
  2. peaks <- getPeaks(peakfile, sort_peaks = TRUE)

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

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

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

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

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

  1. data(example_counts, package = "chromVAR")

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

主体是一个存放count的dgCMatrix

  1. assay(example_counts)[1:5,1:2]
  2.  
  3. 5 x 2 sparse Matrix of class "dgCMatrix"
  4. singles-GM12878-140905-1 singles-GM12878-140905-2
  5. [1,] . 1
  6. [2,] . .
  7. [3,] . .
  8. [4,] . .
  9. [5,] . .

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

  1. head(rowData(example_counts), n=2)
  2.  
  3. DataFrame with 2 rows and 3 columns
  4. score qval name
  5. 1 259 25.99629 GM_peak_6
  6. 2 82 8.21494 H1_peak_

列表示样本,

  1. head(colData(example_counts), n=2)
  2.  
  3. DataFrame with 2 rows and 2 columns
  4. Cell_Type depth
  5. singles-GM12878-140905-1 GM 9220
  6. singles-GM12878-140905-2 GM 9401

前期准备

分析peak中的GC含量

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

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

过滤输入(单细胞)

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

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

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

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

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

确认标准后,就可以过滤

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

获取Motifs和分析包含motif的peak

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

  1. # BiocManager::install("JASPAR2018")
  2. library(JASPAR2018)
  3. species <- "Homo sapiens"
  4. collection <- "CORE"
  5. opts <- list()
  6. opts["species"] <- species
  7. opts["collection"] <- collection
  8. out <- TFBSTools::getMatrixSet(JASPAR2018, opts)
  9.  
  10. if (!isTRUE(all.equal(TFBSTools::name(out), names(out))))
  11. names(out) <- paste(names(out), TFBSTools::name(out),
  12. sep = "_")
  13. 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. library(motifmatchr)
  2. motif_ix <- matchMotifs(motifs, counts_filtered,
  3. 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. bg <- getBackgroundPeaks(object = counts_filtered)
  2. dev <- computeDeviations(object = counts_filtered, annotations = motif_ix,
  3. background_peaks = bg)

变异度

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

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

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

偏离度可视化

可以利用tSNE可视化偏离度

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

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

参考资料

发表评论

匿名网友

拖动滑块以完成验证
加载失败