scATAC分析实战（step3）

Step3-Run chromVAR on se

Running chromVAR on single-cell summarized experiment

准备工作

chromVARmotifs包需要在Github上下载：https://github.com/GreenleafLab/chromVARmotifs

PS:这是我见过的为数不多的还算比较稳定的包，一般都能下下来

The TFBSTools package from Bioconductor is the only direct dependency; the PWMatrixList object from that package is used to store the motifs.

这个包里存储了很多东西，包括motif Reference PWM，主要是人和鼠的motif数据：

data("human_pwms_v1") #curated collection of human motifs from cisBP database
data("mouse_pwms_v1") #curated collection of mouse motifs from cisBP database
data("human_pwms_v2") #filtered collection of human motifs from cisBP database
data("mouse_pwms_v2") #filtered collection of mouse motifs from cisBP database
data("homer_pwms") #motifs from HOMER
data("encode_pwms") #motifs from ENCODE

函数说明

matchMotifs

可以实现peaks和motif数据库的比对，除此之外该函数还可以返回其他信息：每个peak有多少motif，每个peak的最大motif score（对应参数out = scores），每个匹配motif的具体位置（对应参数 out = positions），然后返回值out = matches or out = scores 都可以传给computeDeviations函数进行计算

computeDeviations

返回两个assay，第一个矩阵（可以通过deviations(dev) or assays(dev)$deviations访问）是每个样本（或细胞）的每个peaks可及性的bias corrected "deviation" 分数（原理是反映accessibility的分布，计算每个样本相对bulk样本平均accessibility的“偏离”分数，类似多重t检验的差异基因分析），考虑了GC偏倚和background peaks level。第二个矩阵（可通过(deviationScores(dev) or assays(deviations)$z)访问）给得是Z-score，这个分数可以一定程度反映“偏离”分数的置信度。

dev <- computeDeviations(object = counts_filtered, annotations = motif_ix)

computeVariability

对上一步的结果进一步计算整体peaks的变异程度，应用了bootstrap方法

The function computeVariability returns a data.frame that contains the variability (standard deviation of the z scores computed above across all cell/samples for a set of peaks), bootstrap confidence intervals for that variability (by resampling cells/samples), and a p-value for the variability being greater than the null hypothesis of 1.

variability <- computeVariability(dev)

plotVariability(variability, use_plotly = FALSE) 

完整代码

library(chromVAR)
library(SummarizedExperiment)
library(chromVARmotifs)
library(motifmatchr)
library(BiocParallel)
library(BSgenome.Hsapiens.UCSC.hg19)
register(SerialParam())
set.seed(1)

#-----------------
# Read Inputs
#-----------------
genome <- BSgenome.Hsapiens.UCSC.hg19
se <- readRDS("results/scATAC-Summarized-Experiment.rds") # single-cell summarized experiment rowRanges as peaks  上一步的se，union peaks的counts
se <- addGCBias(se, genome = genome)#Computes GC content for peaks
data("human_pwms_v1")#这个数据在chromVARmotifs这个包里
matches <- matchMotifs(human_pwms_v1, rowRanges(se), genome = "BSgenome.Hsapiens.UCSC.hg19")#这个计算需要一段时间，比对union peaks和motif

#compute deviations
dev <- computeDeviations(object = se, annotations = matches)#Computes deviations in chromatin accessibility across sets of annotations  这个也很耗时

#compute variability
metadata(dev)$Variability <- computeVariability(dev)

#add se
metadata(dev)$SummarizedExperiment <- se

#add matches
metadata(dev)$motifMatches <- matches

saveRDS(dev, "results/chromVAR-Summarized-Experiment.rds")

引用：http://www.bioconductor.org/packages/release/bioc/vignettes/chromVAR/inst/doc/Introduction.html

看一下dev的结构：每一行应该就是motif TFgene

dev
class: chromVARDeviations 
dim: 1764 3770 
metadata(3): Variability SummarizedExperiment motifMatches
assays(2): deviations z
rownames(1764): ENSG00000008196_LINE2_TFAP2B_D_N1
  ENSG00000008196_LINE3_TFAP2B_D_N1 ...
  ENSG00000122145_LINE20002_TBX22_I_N1
  ENSG00000122145_LINE20003_TBX22_I_N1
rowData names(3): name fractionMatches fractionBackgroundOverlap
colnames(3770): PBMC#GAGGTCCGTCTCTGCT-1 PBMC#CCCTGATGTCTTAGCA-1 ...
  PBMC#GGAACTTCATTTGGCA-1 PBMC#AAACTGCCATTCCCGT-1
colData names(5): FRIP uniqueFrags Clusters TSNE1 TSNE2

引用：[https://doi.org/10.1038/s41587-019-0206-z]