Signac｜成年小鼠大脑单细胞ATAC分析（2）

数据科学工厂

发布于 2024-06-18 14:17:24

730

发布于 2024-06-18 14:17:24

文章被收录于专栏：数据科学（冷冻工厂）

引言

在本教程中，我们将探讨由10x Genomics公司提供的成年小鼠大脑细胞的单细胞ATAC-seq数据集。本教程中使用的所有相关文件均可在10x Genomics官方网站上获取。

本教程复现了之前在人类外周血单核细胞（PBMC）的Signac入门教程中执行的命令。我们通过在不同的系统上进行相同的分析，来展示其性能以及对不同组织类型的适用性，并提供了一个来自不同物种的示例。

创建基因活动矩阵

# compute gene activities
gene.activities <- GeneActivity(brain)

# add the gene activity matrix to the Seurat object as a new assay
brain[['RNA']] <- CreateAssayObject(counts = gene.activities)
brain <- NormalizeData(
  object = brain,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(brain$nCount_RNA)
)

DefaultAssay(brain) <- 'RNA'
FeaturePlot(
  object = brain,
  features = c('Sst','Pvalb',"Gad2","Neurod6","Rorb","Syt6"),
  pt.size = 0.1,
  max.cutoff = 'q95',
  ncol = 3
)

与 scRNA-seq 数据整合

为了更好地解读单细胞ATAC-seq数据，我们可以根据来自相同生物体系（即成年小鼠大脑）的单细胞RNA测序（scRNA-seq）实验结果，

# Load the pre-processed scRNA-seq data
allen_rna <- readRDS("../vignette_data/allen_brain.rds")
allen_rna <- UpdateSeuratObject(allen_rna)
allen_rna <- FindVariableFeatures(
  object = allen_rna,
  nfeatures = 5000
)

transfer.anchors <- FindTransferAnchors(
  reference = allen_rna,
  query = brain,
  reduction = 'cca',
  dims = 1:30
)

predicted.labels <- TransferData(
  anchorset = transfer.anchors,
  refdata = allen_rna$subclass,
  weight.reduction = brain[['lsi']],
  dims = 2:30
)

brain <- AddMetaData(object = brain, metadata = predicted.labels)

plot1 <- DimPlot(allen_rna, group.by = 'subclass', label = TRUE, repel = TRUE) + NoLegend() + ggtitle('scRNA-seq')
plot2 <- DimPlot(brain, group.by = 'predicted.id', label = TRUE, repel = TRUE) + NoLegend() + ggtitle('scATAC-seq')
plot1 + plot2

您可以看到基于 RNA 的分类与 UMAP 可视化一致，仅根据 ATAC-seq 数据计算。

查找簇之间可差异访问的峰值

在这里，我们发现皮层不同层的兴奋性神经元之间的可访问区域存在差异。

#switch back to working with peaks instead of gene activities
DefaultAssay(brain) <- 'peaks'
Idents(brain) <- "predicted.id"

da_peaks <- FindMarkers(
  object = brain,
  ident.1 = c("L2/3 IT"), 
  ident.2 = c("L4", "L5 IT", "L6 IT"),
  test.use = 'LR',
  latent.vars = 'nCount_peaks'
)

head(da_peaks)
##                                  p_val avg_log2FC pct.1 pct.2    p_val_adj
## chr4-86523678-86525285    3.266647e-69   3.691294 0.426 0.037 5.135267e-64
## chr2-118700082-118704897  8.553383e-61   2.092487 0.648 0.182 1.344617e-55
## chr15-87605281-87607659   3.864918e-55   2.450827 0.499 0.097 6.075767e-50
## chr10-107751762-107753240 1.534485e-52   1.801355 0.632 0.192 2.412257e-47
## chr4-101303935-101305131  5.949521e-51   3.427059 0.356 0.031 9.352825e-46
## chr13-69329933-69331707   1.604991e-49  -2.254722 0.140 0.435 2.523094e-44

plot1 <- VlnPlot(
  object = brain,
  features = rownames(da_peaks)[1],
  pt.size = 0.1,
  idents = c("L4","L5 IT","L2/3 IT")
)
plot2 <- FeaturePlot(
  object = brain,
  features = rownames(da_peaks)[1],
  pt.size = 0.1,
  max.cutoff = 'q95'
)
plot1 | plot2

open_l23 <- rownames(da_peaks[da_peaks$avg_log2FC > 3, ])
open_l456 <- rownames(da_peaks[da_peaks$avg_log2FC < 3, ])
closest_l23 <- ClosestFeature(brain, open_l23)
closest_l456 <- ClosestFeature(brain, open_l456)
head(closest_l23)

##                                 tx_id gene_name            gene_id
## ENSMUST00000151481 ENSMUST00000151481   Fam154a ENSMUSG00000028492
## ENSMUST00000131864 ENSMUST00000131864   Gm12796 ENSMUSG00000085721
## ENSMUST00000139527 ENSMUST00000139527     Yipf1 ENSMUSG00000057375
## ENSMUSE00001329193 ENSMUST00000185379   Gm29414 ENSMUSG00000099392
## ENSMUSE00000514286 ENSMUST00000077353      Hmbs ENSMUSG00000032126
## ENSMUST00000161356 ENSMUST00000161356      Reln ENSMUSG00000042453
##                      gene_biotype type           closest_region
## ENSMUST00000151481 protein_coding  gap   chr4-86487920-86538964
## ENSMUST00000131864        lincRNA  gap chr4-101292521-101318425
## ENSMUST00000139527 protein_coding  cds chr4-107345009-107345191
## ENSMUSE00001329193        lincRNA exon   chr1-25026581-25026779
## ENSMUSE00000514286 protein_coding exon   chr9-44344010-44344228
## ENSMUST00000161356 protein_coding  gap   chr5-21891568-21895988
##                                query_region distance
## ENSMUST00000151481   chr4-86523678-86525285        0
## ENSMUST00000131864 chr4-101303935-101305131        0
## ENSMUST00000139527 chr4-107344435-107345145        0
## ENSMUSE00001329193   chr1-25008426-25009334    17246
## ENSMUSE00000514286   chr9-44345250-44346015     1021
## ENSMUST00000161356   chr5-21894051-21894682        0

head(closest_l456)
##                                 tx_id gene_name            gene_id
## ENSMUST00000104937 ENSMUST00000104937   Ankrd63 ENSMUSG00000078137
## ENSMUSE00000647021 ENSMUST00000068088   Fam19a5 ENSMUSG00000054863
## ENSMUST00000165341 ENSMUST00000165341     Otogl ENSMUSG00000091455
## ENSMUST00000044081 ENSMUST00000044081     Papd7 ENSMUSG00000034575
## ENSMUST00000070198 ENSMUST00000070198    Ppp3ca ENSMUSG00000028161
## ENSMUST00000084628 ENSMUST00000084628    Hs3st2 ENSMUSG00000046321
##                      gene_biotype type            closest_region
## ENSMUST00000104937 protein_coding  cds  chr2-118702266-118703438
## ENSMUSE00000647021 protein_coding exon   chr15-87625230-87625486
## ENSMUST00000165341 protein_coding  utr chr10-107762223-107762309
## ENSMUST00000044081 protein_coding  utr   chr13-69497959-69499915
## ENSMUST00000070198 protein_coding  utr  chr3-136935226-136937727
## ENSMUST00000084628 protein_coding  cds  chr7-121392730-121393214
##                                 query_region distance
## ENSMUST00000104937  chr2-118700082-118704897        0
## ENSMUSE00000647021   chr15-87605281-87607659    17570
## ENSMUST00000165341 chr10-107751762-107753240     8982
## ENSMUST00000044081   chr13-69329933-69331707   166251
## ENSMUST00000070198  chr3-137056475-137058371   118747
## ENSMUST00000084628  chr7-121391215-121395519        0

绘制基因组区域

我们同样可以利用CoveragePlot()函数，根据不同的细胞聚类、细胞类型或对象中存储的其他任何元数据信息，为特定的基因组区域绘制出分组的覆盖度图。这些覆盖度图实际上是伪批量的可访问性轨迹图，通过将同一组内所有细胞的信号进行平均，从而在视觉上展示出特定区域内DNA的可访问性情况。

# show cell types with at least 50 cells
idents.plot <- names(which(table(Idents(brain)) > 50))

CoveragePlot(
  object = brain,
  region = c("Neurod6", "Gad2"),
  idents = idents.plot,
  extend.upstream = 1000,
  extend.downstream = 1000,
  ncol = 1
)