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
computeVariabilityreturns 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
