Week4— chromVAR:预测染色质可及性相关的转录因子

第四周 2018— 06.11-06.17
原文： chromVAR: Inferring transcription factor-associated accessibility from single-cell epigenomic data
期刊：Nat Methods ，14(10): 975–978.
发表日期：2017 October ;
doi: 10.1038/nmeth.4401.
关键词：

概述

ChromVAR是一个R包，可用于鉴定不同细胞类型的转录因子基序。它通过估计染色质可及性的激活和关闭来分析稀疏染色质可及性数据，能够对scATAC-seq图谱进行精确地聚类，并鉴定与染色质可及性变化相关的已知的和新的序列基序。不仅可用于单细胞也可用于bulk细胞的表观数据，从而鉴定细胞类型和细胞特异性的反式调控因子。另外chromVAR对测序深度低和多种细胞类型或样本的数据具有较高的灵敏性。

输入数据

ChromVAR 包的输入数据包括3个：
1）比对后的reads;
2）chromatin accessibility peaks；
3）一组代表motif位置的权重矩阵（PWM）或基因组注释的染色质特征数据集

他们从cisBP 数据库提取了一系列人和鼠的PWMs，cisBP数据库包含多样而全面的已知TF基序集合。用户也可以自己提供TF 基序或其他类型的基因组注释，如增强子模块，ChIP-seq peaks, GWAS 疾病注释等。

工作流程

• 首先计算Raw deviation in accessibility
• 计算Raw deviation for background peaks
• Bias corrected deviations and z-scores
• 下游应用

具体使用方法

安装

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("chromVAR")
# BiocManager::install("GO.db")

• 并行
  建议使用BiocParallel的register函数指定首选的并行化方法

# unix
library(BiocParallel)
register(MulticoreParam(8)) # Use 8 cores
# windows
register(SnowParam(SnowParam(workers=1, type = "SOCK")))
# 如果不想用多核，建议用SerialParam注册该选项
register(SerialParam())

应用

chromVAR的应用包括以下几个方面：

• 首先计算偏差

library(chromVAR)
library(motifmatchr)
library(SummarizedExperiment)
library(Matrix)
library(ggplot2)
library(BiocParallel)
library(BSgenome.Hsapiens.UCSC.hg19)
register(SerialParam())
data(example_counts, package = "chromVAR")
set.seed(2017)
example_counts <- addGCBias(example_counts, genome = BSgenome.Hsapiens.UCSC.hg19)
counts_filtered <- filterSamples(example_counts, min_depth = 1500, 
                                 min_in_peaks = 0.15, shiny = FALSE)
counts_filtered <- filterPeaks(counts_filtered)
motifs <- getJasparMotifs()
motif_ix <- matchMotifs(motifs, counts_filtered, genome = BSgenome.Hsapiens.UCSC.hg19)
dev <- computeDeviations(object = counts_filtered, 
                                 annotations = motif_ix)

• 1. 计算感兴趣的cell或样本之间每个motif或注释的变异性
  用函数plotVariability

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

Variability

• 2.聚类
  也可以用偏差校正偏差（bias corrected deviations）对样本聚类，函数getSampleCorrelation首先剔除高度相关注释和低可变性注释，然后细胞间的相关性。

sample_cor <- getSampleCorrelation(dev)

library(pheatmap)
pheatmap(as.dist(sample_cor), 
         annotation_row = colData(dev), 
         clustering_distance_rows = as.dist(1-sample_cor), 
         clustering_distance_cols = as.dist(1-sample_cor))

Clustering

• 3. 细胞或样本间的相似性
  也可以用tSNE看细胞的相似性，函数deviationsTsne可以实现此功能。参数shiny=TRUE可以首先交互界面，perplexity根据细胞的实际数目调整，通常情况下，100个细胞时，perplexity在30-50都是可以的。

tsne_results <- deviationsTsne(dev, threshold = 1.5, perplexity = 10, 
                               shiny = FALSE)
# plotDeviationsTsne 绘图
tsne_plots <- plotDeviationsTsne(dev, tsne_results, annotation = "TEAD3", 
                                   sample_column = "Cell_Type", shiny = FALSE)
tsne_plots[[1]]
tsne_plots[[2]]

Cell / sample similarity

• 4.可及性和变异性的差异
  函数differentialDeviations和differentialVariability可以看各组间对某一注释的偏差校正偏差是否存在显著性差异和显著性变异

diff_acc <- differentialDeviations(dev, "Cell_Type")
head(diff_acc)
##                           p_value p_value_adjusted
## MA0025.1_NFIL3       6.667785e-02     1.175006e-01
## MA0030.1_FOXF2       1.802723e-02     4.166773e-02
## MA0031.1_FOXD1       3.825026e-03     1.077708e-02

diff_var <- differentialVariability(dev, "Cell_Type")
head(diff_var)
##                        p_value p_value_adjusted
## MA0025.1_NFIL3       0.4330854        0.9135025
## MA0030.1_FOXF2       0.7011659        0.9652404

• 5. motif/kmer的相似性
  设置参数what= "annotations"就可以用tsne看motif间的相似性

inv_tsne_results <- deviationsTsne(dev, threshold = 1.5, perplexity = 8, 
                                    what = "annotations", shiny = FALSE)
ggplot(inv_tsne_results, aes(x = Dim1, y = Dim2)) + geom_point() + 
  chromVAR_theme()

Motif / kmer similarity

• 6. Kmers and sequence specificity of variation
  利用kmers作为注释，我们可以使用kmers来识别染色质可及性变异所需的精确核苷酸。

kmer_ix <- matchKmers(6, counts_filtered, genome = BSgenome.Hsapiens.UCSC.hg19)
kmer_dev <- computeDeviations(counts_filtered, kmer_ix)
kmer_cov <- deviationsCovariability(kmer_dev)
plotKmerMismatch("CATTCC",kmer_cov)

Kmers and sequence specificity of variation

• 7. kmer从头组装
  根据kmer偏差结果，可以使用assembleKmers函数构建de novo motif。

de_novos <- assembleKmers(kmer_dev, progress = FALSE) #no progress bar
de_novos

pwmDistance函数可以看de novo motif与已知motif的匹配关系，最终会返回3个列表，dist是每对motifs间的距离，strandmotif匹配的链，offset 是motifs间的偏移。

dist_to_known <- pwmDistance(de_novos, motifs)
closest_match1 <- which.min(dist_to_known$dist[1,])
dist_to_known$strand[1,closest_match1]

library(ggmotif) # Package on github at AliciaSchep/ggmotif. Can use seqLogo alternatively
library(TFBSTools)
# De novo motif
ggmotif_plot(de_novos[[1]])
# Closest matching known
ggmotif_plot(toPWM(reverseComplement(motifs[[closest_match1]]),type = "prob"))

De novo kmer assembly

• 8. motif之间的协同作用
  chromVAR还包括计算注释/motif对之间“协同作用”的函数。

getAnnotationSynergy(counts_filtered, motif_ix[,c(83,24,20)])

• 9. 相关性
  一个极端的协同值表明转录因子间可能存在合作或竞争的作用。函数getAnnotationCorrelation可以计算motifs间的相关性。

getAnnotationCorrelation(counts_filtered, motif_ix[,c(83,24,20)])

##                      MA0139.1_CTCF MA0107.1_RELA MA0140.2_GATA1::TAL1
## MA0139.1_CTCF            1.0000000   -0.41950962          -0.27827092
## MA0107.1_RELA           -0.4195096    1.00000000          -0.08943439
## MA0140.2_GATA1::TAL1    -0.2782709   -0.08943439           1.00000000

---

相关资料

• chromVAR github：https://github.com/GreenleafLab/chromVAR
• https://github.com/GreenleafLab/chromVAR/blob/master/vignettes/Introduction.Rmd
• https://github.com/GreenleafLab/chromVAR/tree/master/vignettes/Articles
• https://greenleaflab.github.io/chromVAR/articles/Articles/Applications.html
• motifmatchr : https://github.com/GreenleafLab/motifmatchr
• chromVARmotifs:https://github.com/GreenleafLab/chromVARmotifs