scATAC-seq| motif 分析

scATAC-seq| motif 分析 by 生信技能树

在单细胞转录组冒尖的时候，单细胞其它组学却一直不温不火，比如单细胞ATAC以及单细胞免疫组库。所以虽然很早之前我们就推荐了github代码：https://github.com/GreenleafLab/10x-scATAC-2019

01_Filter_Cells_v2.R
02_Get_Peak_Set_hg19_v2.R
03_Run_chromVAR_v2.R
04_Run_Cicero_v2.R
05_Cluster_Unique_Peaks_v2.R
06_Analyze_UMAP_Trajectory.R
07_ChromVAR_For_GWAS_w_CoAccessbility_v2.R
08_Run_scCNV_v2.R

但是一直没有深入解读它，在2021的时候《单细胞天地》的一个小伙伴写了个开头：

但是没有能坚持下来，其实文章给的配套github代码非常齐全了，就是需要花时间钻研和解读。

恰好2022我们又看到了一个带有全套github代码以及配套数据的文章，见：10x的单细胞ATAC上游流程之cellranger-atac，在等待的期间，交流群的学徒就先拿 atac_v1_pbmc_10k 这样的示例数据开始练习了。所以转发学徒的系列教程给大家。

现在继续：

上回的scRNA文件pbmc已经分好群了，接下来对motif进行分析

对motif 简单的介绍：单细胞还能拉上转录因子？？？

流程

载入数据
数据库 JASPAR2020 获取motif矩阵：getMatrixSet
添加motif矩阵：AddMotifs
计算motif活性分数:RunChromVAR
寻找显著差异的peaks：FindMarkers
peaks对应的motif :FindMotifs
可视化：MotifPlot and FeaturePlot

1. 载入数据

注意：DefaultAssay(pbmc)<- 'peaks'要改成peaks，已经不是分群时候用到的RNA了

library(Signac)
library(Seurat)
#BiocManager::install('TFBSTools')
library(JASPAR2020)
library(TFBSTools)
library(motifmatchr)
library(BSgenome.Hsapiens.UCSC.hg19)
library(chromVAR)
library(ggseqlogo)
set.seed(1234)

# 载入数据
pbmc <- readRDS("/home/data/yulan/scATAC/data/1.4scATAC_fina.rds")
pbmc
# An object of class Seurat 
# 107571 features across 7060 samples within 2 assays 
# Active assay: RNA (20010 features, 0 variable features)
# 1 other assay present: peaks
# 2 dimensional reductions calculated: lsi, umap

# 选择peaks
pbmc@assays
# $peaks
# ChromatinAssay data with 87561 features for 7060 cells
# Variable features: 43781 
# Genome: hg19 
# Annotation present: TRUE 
# Motifs present: FALSE 
# Fragment files: 1 
# 
# $RNA
# Assay data with 20010 features for 7060 cells
# First 10 features:
#   PLCXD1, GTPBP6, PPP2R3B, SHOX, CRLF2, CSF2RA, IL3RA, SLC25A6, ASMTL, P2RY8
DefaultAssay(pbmc) <- 'peaks'

2. 获取motif矩阵

从数据库 JASPAR2020 获取motif矩阵：getMatrixSet

# 从公共数据库中获取motif矩阵
# Get a list of motif position frequency matrices from the JASPAR database
# http://www.bioconductor.org/packages/release/data/annotation/vignettes/JASPAR2020/inst/doc/JASPAR2020.html
# https://jaspar.genereg.net/
pfm <- getMatrixSet(
  x = JASPAR2020,
  opts = list(collection = "CORE", 
              tax_group = 'vertebrates', # 脊椎动物
              all_versions = FALSE)
)
pfm
# PFMatrixList of length 746
# names(746): MA0004.1 MA0006.1 MA0019.1 MA0029.1 MA0030.1 ... MA0526.3 MA0748.2 MA0528.2 MA0609.2

JASPAR是一个免费公开的转录因子数据库，在该数据库中收录了转录因子(TF)的motif信息，可以用来预测转录因子与序列的结合区域。存储六个分类组中多个物种的TF的位置频率矩阵（PFM）。在JASPAR的第8版中，CORE系列已扩展为245个新的PFM（脊椎动物169个，植物42个，线虫17个，昆虫10个，真菌7个），更新了156个PFM（脊椎动物125个，植物28个，昆虫3个）

添加motif矩阵：AddMotifs

# add motif information
pbmc <- AddMotifs(
  object = pbmc,
  genome = BSgenome.Hsapiens.UCSC.hg19,
  pfm = pfm
)
# Building motif matrix
# Finding motif positions
# Creating Motif object

# 得到motifs矩阵
GetAssayData(object = pbmc, slot = "motifs")
# A Motif object containing 746 motifs in 87561 regions

计算motif活性分数:RunChromVAR

# 计算 motif deviation score
pbmc <- RunChromVAR(
  object = pbmc,
  genome = BSgenome.Hsapiens.UCSC.hg19
)
# Computing GC bias per region
# Selecting background regions
# Computing deviations from background
# Constructing chromVAR assay

# 新增了一个assay：chromvar
pbmc
# An object of class Seurat 
# 108317 features across 7060 samples within 3 assays 
# Active assay: peaks (87561 features, 43781 variable features)
# 2 other assays present: RNA, chromvar
# 2 dimensional reductions calculated: lsi, umap

寻找显著差异的 peaks：FindMarkers

# 显著差异的peaks
DimPlot(pbmc) + DimPlot(pbmc,group.by = 'predicted.id')
table(pbmc$predicted.id)
table(Idents(pbmc))
Idents(pbmc) <- pbmc$predicted.id

# 对其中两个细胞亚群差异分析
da_peaks <- Seurat::FindMarkers(
  object = pbmc,
  ident.1 = 'CD4 Naive',
  ident.2 = 'CD8 Naive',
  only.pos = TRUE,
  test.use = 'LR',# 推荐用logistic regression
  min.pct = 0.05,# 该基因表达数目占该类细胞总数的比例
  latent.vars = 'nCount_peaks'
)
# get top differentially accessible peaks
top.da.peak <- rownames(da_peaks[da_peaks$p_val < 0.005, ])

> head(da_peaks)
                                p_val avg_log2FC pct.1 pct.2    p_val_adj
chr1-234544345-234545516 8.454563e-16  0.5295198 0.191 0.022 7.402900e-11
chr6-167362837-167366950 1.648728e-13  0.2943871 0.285 0.085 1.443643e-08
chr12-6895592-6896474    3.370162e-13  0.4385442 0.123 0.006 2.950948e-08
chr12-6898094-6899126    8.623449e-13  0.2885348 0.106 0.003 7.550778e-08
chr9-90340551-90341713   1.854913e-12  0.2595704 0.142 0.016 1.624180e-07
chr12-6899696-6902624    1.489885e-11  0.2842456 0.206 0.041 1.304558e-06

求peaks对应的 motifs: FindMotifs

# 求peaks对应的motifs
# test enrichment
enriched.motifs <- FindMotifs(
  object = pbmc,
  features = top.da.peak
)
head(enriched.motifs)

> head(enriched.motifs)
            motif observed background percent.observed percent.background fold.enrichment
MA0769.2 MA0769.2       19       5146         41.30435            12.8650        3.210598
MA0523.1 MA0523.1       20       6377         43.47826            15.9425        2.727192
MA0162.4 MA0162.4       19       5961         41.30435            14.9025        2.771639
MA0733.1 MA0733.1       14       3953         30.43478             9.8825        3.079664
MA1650.1 MA1650.1       13       3442         28.26087             8.6050        3.284238
MA0732.1 MA0732.1       13       4031         28.26087            10.0775        2.804353
               pvalue motif.name    p.adjust
MA0769.2 1.475866e-06       TCF7 0.001100996
MA0523.1 8.828251e-06     TCF7L2 0.003292938
MA0162.4 1.336985e-05       EGR1 0.003324637
MA0733.1 9.317482e-05       EGR4 0.013940194
MA1650.1 9.343293e-05     ZBTB14 0.013940194
MA0732.1 4.502527e-04       EGR3 0.048616281

motif 可视化：MotifPlot and FeaturePlot

head(rownames(enriched.motifs)
# [1] "MA0769.2" "MA0523.1" "MA0162.4" "MA0733.1" "MA1650.1" "MA0732.1"
MotifPlot(
  object = pbmc,
  motifs = head(rownames(enriched.motifs)),
  assay = 'peaks'
)
# 默认以转录因子为作图标题

FeaturePlot(
  object = pbmc,
  features = head(rownames(enriched.motifs)),
  min.cutoff = 'q10',
  max.cutoff = 'q90',
  pt.size = 0.1
)

同样的方法，寻找显著差异的 motifs

# 显著差异的 motifs
DefaultAssay(pbmc) <- 'chromvar'
da_motifs <- Seurat::FindMarkers(
  object = pbmc,
  ident.1 = 'CD4 Naive',
  ident.2 = 'CD8 Naive',
  only.pos = TRUE,
  test.use = 'LR', # 推荐用logistic regression
  min.pct = 0.05, # 该基因表达数目占该类细胞总数的比例
  latent.vars = 'nCount_peaks'
)
head(da_motifs)
#                 p_val avg_log2FC pct.1 pct.2  p_val_adj
# MA0050.2 2.459215e-05  0.2994300 0.186 0.110 0.01834574
# MA0038.2 2.275283e-04  0.4947875 0.428 0.312 0.16973610
# MA0836.2 3.334793e-04  2.0343057 0.062 0.025 0.24877553
# MA1418.1 3.555421e-04  0.3582367 0.367 0.274 0.26523439
# MA0826.1 3.669158e-04  0.5638397 0.484 0.394 0.27371917
# MA0466.2 4.664331e-04  1.6827524 0.054 0.032 0.34795909

参考:https://satijalab.org/signac/articles/motif_vignette.html

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

生物信息学马拉松授课（买一得五），你的生物信息学入门课

ixxmu / mp_duty

scATAC-seq| motif 分析 #3045

scATAC-seq| motif 分析 by 生信技能树

流程

1. 载入数据

2. 获取motif矩阵

文末友情宣传