chromVAR

艾丽西亚Schep

chromVAR是一个用于稀疏染色质可达性分析的R包。chromVAR将来自ATAC-seq或DNAse-seq实验的对齐片段(过滤重复和低质量)以及基因组注释(如motif位置)作为输入。chromVAR为每个注释和每个单元格或样本计算一个基于所有单元格或样本平均值的预期可访问性的偏差校正“偏差”。

这个小插图涵盖了chromVAR标准输入的基本用法。有关涵盖各种输入和其他应用程序的不同选项的更详细文档,请查看文档的网站.有关该方法和应用的更多描述,请参阅出版pdf补充).

##加载包

使用库或require加载包和有用的附加包。

图书馆(chromVAR)图书馆(motifmatchr)图书馆(矩阵)图书馆(SummarizedExperiment)图书馆(BiocParallel)set.seed2017

设置多处理选项包使用BiocParallel进行多处理。查看BiocParallel的文档以查看可用的选项。可以通过寄存器功能进行设置。例如,使用8核multicorepam:

注册MulticoreParam8))

若要为多处理任务启用进度条,请使用

注册MulticoreParam8progressbar =真正的))

对于Windows,MulticoreParam不会工作,但你可以使用SnowParam:

注册SnowParam工人=1类型=“袜子”))

即使你不想使用多个核心,最好使用SerialParam显式注册该选择:

注册SerialParam())

请参阅相关文档BiocParallel有关的更多信息注册函数和多处理的各种选项。

##读取输入

chromVAR以落在开放染色质峰的片段计数表作为输入。有许多软件包可以从表观基因组数据中创建计数表;chromVAR还提供了一种针对单细胞ATAC-seq计数的方法。有关计数表的选项以及如何格式化通过其他软件计算的计数表的详细信息,请参阅文档网站的计数部分

# # #的山峰

对于使用chromVAR,建议使用固定宽度、不重叠的峰值。该方法对峰值宽度的精确选择具有相当的鲁棒性,但建议宽度为250 ~ 500 bp。看到补充对论文的部分内容进行了讨论,并补充了相关图峰宽度的选择。

如果分析单细胞数据,使用来自同一种群或类似种群(或可能来自公共资源,如路线图表观基因组学项目)的大量ATAC或DNAse-seq数据的峰值是有意义的。

如果组合来自不同种群的多个峰值文件,建议将这些峰值组合在一起。的filterPeaks函数(在这个小插图的后面一点演示)将把峰值减少到一个非重叠的集合,基于这个集合,重叠的峰值在所有数据中具有更强的信号。

这个函数getPeaks读取峰值作为一个GenomicRanges对象。我们将通过阅读一个小样本床文件来展示它的使用。我们将使用sort_peaks参数表示要对峰进行排序。

peakfile < -执行“extdata / test_bed.txt”包=“chromVAR”山峰< -getPeaks(peakfilesort_peaks =真正的
##山峰排序

对于窄峰格式的峰值读取,chromVAR包含一个函数,readNarrowpeaks,它将读入文件,将峰值大小调整为基于宽度参数,并删除与较强信号的峰值重叠的峰值(如果non_overlapping设置为TRUE -默认值)。

# # #计数

这个函数getCounts返回一个chromVARCounts对象,其中包含每个检测峰的每个样品/细胞的片段计数矩阵。可以访问此数据计数(fragment_counts).使用Matrix包,因此如果矩阵是稀疏的,矩阵将被存储为稀疏矩阵。我们将用包中包含的一个非常小的读取集来演示这个函数:

bamfile < -执行“extdata / test_RG.bam”包=“chromVAR”fragment_counts < -getCounts(bamfile山峰,成对的=真正的by_rg =真正的格式=“砰”colData =DataFramecelltype =“通用汽车”))
##读取文件:/tmp/RtmpY8E6dV/Rinst16684f68f03d9/chromVAR/extdata/test_RG.bam

这里我们只传递了一个bam文件,但是我们也可以传递一个bam文件名的向量。在这种情况下,对于列数据,我们将为每个文件指定适当的值,例如DataFrame(celltype = c("GM","K562"))如果我们传递一个GM细胞文件和一个K562细胞文件。

如果RG标记不用于在一个文件中组合多个样本,请使用by_rg = FALSE.有关计数阅读的更多信息,请参见文档网站的计数部分

示例数据

对于插图的其余部分,我们将使用一个非常小的(但比前面的示例略大)数据集,其中包含10个GM细胞和10个H1细胞,已将其作为summarize实验对象读入。

数据(example_counts包=“chromVAR”(example_counts)
##类:rangedsummarize实验## dim: 6 50 ##元数据(0):## assays(1):计数## rownames: NULL ## rowData names(3): score qval name ## colnames(50): single - gm12878 -140905-1 single - gm12878 -140905-2…## colData name (2): Cell_Type depth

##获取峰值的GC内容

GC含量将用于确定背景峰。这个函数addGCBias返回一个更新后的summarizeexperiment,包含一个名为“bias”的新rowData列。该函数需要输入基因组序列,该序列可以作为BSgenome、FaFile或DNAStringSet对象提供。

图书馆(BSgenome.Hsapiens.UCSC.hg19)example_counts < -addGCBias(example_counts基因组=BSgenome.Hsapiens.UCSC.hg19)rowData(example_counts))
##数据帧与6行4列##分数qval名称偏置## <整数> <数字> <字符> <数字> ## 1 259 25.99629 GM_peak_6 0.652 ## 2 82 8.21494 H1_peak_7 0.680 ## 3 156 15.65456 GM_peak_17 0.788 ## 4 82 8.21494 H1_peak_13 0.674 ## 5 140 14.03065 GM_peak_27 0.600 ## 6 189 18.99529 GM_peak_29a 0.742

看看available.genomes从BSgenome包中获取哪些基因组可作为BSgenome包。要创建自己的BSgenome对象,请查看BSgenomeForge

# #过滤输入

如果处理单细胞数据,建议过滤掉读取量不足或峰值读取量比例低的样本,因为这些可能代表空孔或死细胞。有两个参数用于过滤——min_in_peaks和min_depth。如果没有提供(如上所述),这些截止值是根据数据的中位数估计的。Min_in_peaks设置为峰值中位数比例的0.5倍。Min_depth设置为最大值500或库中值大小的10%。

除非plot = FALSE作为函数的参数filterSamples,将生成一个图。

#查找要保留的样本索引counts_filtered < -filterSamples(example_countsmin_depth =1500min_in_peaks =0.15闪亮的=

如果shiny参数被设置为TRUE(默认值),一个闪亮的小工具将弹出,允许您使用过滤参数,并查看哪些单元格通过过滤器或没有。

要获得过滤后的图像,请使用filterSamplesPlot.默认情况下,情节是交互式的-将其设置为非交互式使用use_plotly = FALSE

#查找要保留的样本索引filtering_plot < -filterSamplesPlot(example_countsmin_depth =1500min_in_peaks =0.15use_plotly =filtering_plot

要返回要保留的样本的索引,而不是一个新的summarizeexperiment对象,请使用ix_return = TRUE。

第九< -filterSamples(example_countsix_return =真正的闪亮的=
min_in_peaks设置为0.105
min_depth设置为1120.65

对于大容量和单单元数据,峰值应该基于至少一定数量的片段进行过滤。至少,每个峰值应该在所有样本中至少有一个片段(由于使用由其他数据定义的峰值集,有可能具有零读取的峰值)。否则,下游功能将无法工作。这个函数filterPeaks如果non_overlap参数设置为TRUE(这是默认值),也将减少峰值设置为不重叠的峰值(对于重叠的峰值保持计数较高的峰值)。

counts_filtered < -filterPeaks(counts_filterednon_overlapping =真正的

获得母题和包含母题的峰

这个函数getJasparMotifs从JASPAR数据库中获取图案。

主题< -getJasparMotifs()

getJasparMotifs()函数默认从JASPAR核心数据库获取人体图案。对于其他物种母题,改变物种参数。

yeast_motifs < -getJasparMotifs物种=“酿酒酵母”

若要使用内核以外的集合,请使用集合论点。选项包括:" CORE ", " CNE ", " PHYLOFACTS ", " SPLICE ", " POLII ", " FAM ", " PBM ", " PBM_HOMEO ", " PBM_HLH "。

getJasparMotifs函数只是一个简单的包装getMatrixSetfrom TFBSTools -如果你喜欢,你也可以使用该函数从JASPAR中获取图案,和/或查看该函数的文档以获得更多信息。

这个函数matchMotifs从motifmatchr包中找到哪些峰值包含哪些motif。默认情况下,它返回一个summarizeexperiment对象,该对象包含一个稀疏矩阵,表示motif匹配与否。这个函数requires an input of a genome sequence, which can be provided as a BSgenome, FaFile, or DNAStringSet object.

图书馆(motifmatchr)motif_ix < -matchMotifs(主题、counts_filtered基因组=BSgenome.Hsapiens.UCSC.hg19)

一种选择是用p.c otoff来决定motif调用的严格程度。默认值是0.00005,它倾向于为人类图案提供合理的图案匹配数量。

除了返回motif匹配,该函数还可以返回每个峰值的motif匹配数量和每个峰值的最大motif得分的附加矩阵(存储为分析)。有关此附加信息,请使用Out =分数.若要返回主题匹配的实际位置,请使用Out =位置.要么是输出Out =匹配Out =分数可以传递给computedeviation函数。

如果不是使用已知的图案,您想使用一定长度的所有kmer,即matchKmers函数可以使用。有关使用kmer作为输入的更多信息,请参见文档网站的注释部分

kmer_ix < -matchKmers6counts_filtered,基因组=BSgenome.Hsapiens.UCSC.hg19)

计算偏差

这个函数computeDeviations返回一个包含两个“assays”的摘要实验。第一个矩阵(通过偏差(dev)化验(dev)美元偏差)将给出每个单元格或样本(列)的每一组峰值(行)的可访问性的偏差校正“偏差”。这一指标代表了峰集相对于基于跨细胞/样本的相等染色质可访问性轮廓的期望的可访问性,由一组与GC和平均可访问性匹配的背景峰集标准化。第二个矩阵(deviationScores (dev)化验(偏差)美元z)给出了偏差z分数,该分数考虑了如果随机采样具有相似GC含量和平均可及性的喙集,出现这种分数的可能性。

dev < -computeDeviations对象=counts_filtered,注释=motif_ix)

背景的山峰

函数computedeviation将使用一组背景峰值来规范化偏差分数。默认情况下,此计算在内部完成,不返回—为了更好地控制此步骤,用户可以运行getBackgroundPeaks函数本身,并将结果传递给background_peaks参数下的computedeviation。

背景峰值是指与GC内容和平均可访问性中的峰值相似的峰值。

bg < -getBackgroundPeaks对象=counts_filtered)

结果是getBackgroundPeaks是一个索引矩阵,其中每一列表示作为背景峰值的峰值的索引。

要使用计算出的背景峰值,只需将它们添加到computedeviation调用中:

dev < -computeDeviations对象=counts_filtered,注释=motif_ix,background_peaks =bg)

可变性

这个函数computeVariability返回一个data.frame,其中包含可变性(对一组峰值在所有单元格/样本上计算的z分数的标准偏差),该可变性的自举置信区间(通过重新采样单元格/样本),以及可变性大于零假设1的p值。

可变性< -computeVariability(dev)plotVariability(可变性,use_plotly =

plotVariability取的输出computeVariability并返回按等级排序的注释集及其可变性的图形。默认情况下,情节将是交互式的,除非您设置use_plotly = FALSE

可视化的偏差

对于可视化单元格,使用TSNE将偏差值投影到二维中是很有用的。中提供了这样做的方便函数deviationsTsne.如果在交互式会话中运行,shiny可以设置为TRUE来加载一个闪亮的小工具来探索参数。

tsne_results < -deviationsTsne(dev,阈值=1.5困惑=10

为了绘制结果,plotDeviationsTsne可以使用。如果运行在交互式会话或交互式Rmarkdown文档中,shiny可以设置为TRUE以生成一个闪亮的小部件。这里我们将展示静态结果。

tsne_plots < -plotDeviationsTsne(dev tsne_resultsannotation_name =“TEAD3”sample_column =“Cell_Type”闪亮的=tsne_plots [[1]]

tsne_plots [[2]]

会话信息

Sys。日期()
##[1]“2023-01-19”
sessionInfo()
## R版本4.2.2(2022-10-31)##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 20.04.5 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack。所以## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_GB LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4 stats graphics grDevices utils datasets methods ##[8]基础## ##其他附加包:# # # # [1] BSgenome.Hsapiens.UCSC.hg19_1.4.3 BSgenome_1.66.2 [3] rtracklayer_1.58.0 Biostrings_2.66.0 # # [5] XVector_0.38.0 BiocParallel_1.32.5 # # [7] SummarizedExperiment_1.28.0 Biobase_2.58.0 # # [9] GenomicRanges_1.50.2 GenomeInfoDb_1.34.7 # # [11] IRanges_2.32.0 S4Vectors_0.36.1 # # [13] BiocGenerics_0.44.0 MatrixGenerics_1.10.0 # # [15] matrixStats_0.63.0 Matrix_1.5-3 # # [17] motifmatchr_1.20.0 chromVAR_1.20.2 # # # #通过加载一个名称空间(而不是附加):# # # # [1] Rtsne_0.16 colorspace_2.0-3 [3] rjson_0.2.21 ellipsis_0.3.2 # # [5] farver_2.1.1 DT_0.27 # # [7] bit64_4.0.5 AnnotationDbi_1.60.0 # # [9] fansi_1.0.3 codetools_0.2-18 # # [11] R.methodsS3_1.8.2 cachem_1.0.6 # # [13] knitr_1.41 jsonlite_1.8.4 # # [15] Rsamtools_2.14.0 seqLogo_1.64.0 # # [17] annotate_1.76.0 GO.db_3.16.0 # # [19] png_0.1-8 R.oo_1.25.0 # # [21] shiny_1.7.4 readr_2.1.3 # # [23] compiler_4.2.2 httr_1.4.4 # # [25] assertthat_0.2.1 fastmap_1.1.0 # # [27] lazyeval_0.2.2 cli_3.6.0 # # [29]later_1.3.0 htmltools_0.5.4 # # [31] tools_4.2.2 gtable_0.3.1 # # [33] glue_1.6.2 TFMPvalue_0.0.9 # # [35] GenomeInfoDbData_1.2.9 reshape2_1.4.4 # # [37] dplyr_1.0.10 Rcpp_1.0.9 # # [39] jquerylib_0.1.4 vctrs_0.5.1 # # [41] xfun_0.36 CNEr_1.34.0 # # [43] stringr_1.5.0 mime_0.12 # # [45] miniUI_0.1.1.1 lifecycle_1.0.3 # # [47] restfulr_0.0.15 poweRlaw_0.70.6 # # [49] gtools_3.9.4 xml_3.99 - 0.13 # # [51] JASPAR2016_1.26.0 zlibbioc_1.44.0 # # [53] scales_1.2.1 vroom_1.6.0 # # [55] hms_1.1.2 promises_1.2.0.1 # #[57] parallel_4.2.2 RColorBrewer_1.1-3 ## [59] yaml_2.3.6 memoise_2.0.1 ## [61] ggplot2_3.4.0 sass_0.4.4 ## [63] stringi_1.7.12 RSQLite_2.2.20 ## [69] pkgconfig_2.0.3 bitops_1.0-7 ## [71] archive_1.1.5 pracma_2.4.2 ## [73] evaluate_0.20 lattice_0.20-45 ## [75] purrr_1.0.1 labeling_0.4.2 ## [77] GenomicAlignments_1.34.0 htmlwidgets_1.6.1 ## [79] bit_4.0.5 tidyselect_1.2.0 ## [81] plyr_1.8.8 magrittr_2.0.3 ## [83] R6_2.5.1generics_0.1.3 ## [85] DelayedArray_0.24.0 DBI_1.1.3 ## [87] withr_2.5.0 pillar_1.8.1 ## [89] KEGGREST_1.38.0 RCurl_1.98-1.9 ## [91] tibble_3.1.8 crayar_1 .5.2 ## [95] utf8_1.2.2 plotly_4.10.1 ## [95] nabor_0.5.0 tzdb_0.3.0 ## [97] rmarkdown_2.19 TFBSTools_1.36.0 ## [99] grid_4.2.2 data.table_1.14.6 ## [101] blob_1.2.3 digest_0.6.31 ## [103] xtable_1. 1.8-4 tidyr_1.2.1 ## [105] httpuv_1.6.8 r.t uls_2.12.2 ## [107] munsell_0.5.0 dirichlet多国al_1.40.0 ## [109] viridisLite_0.4.1 bslib_0.4.2