引言
本文將介紹如何利用Signac和Seurat這兩個(gè)工具私杜，對(duì)一個(gè)同時(shí)記錄了DNA可接觸性和基因表達(dá)的單細(xì)胞數(shù)據(jù)集進(jìn)行綜合分析救欧。我們將以一個(gè)公開的10x Genomics Multiome數(shù)據(jù)集為例笆怠，該數(shù)據(jù)集針對(duì)的是人體的外周血單核細(xì)胞誊爹。
數(shù)據(jù)準(zhǔn)備
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
# load the RNA and ATAC data
counts <- Read10X_h5("../vignette_data/multiomic/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5")
fragpath <- "../vignette_data/multiomic/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz"
# get gene annotations for hg38
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevels(annotation) <- paste0('chr', seqlevels(annotation))
# create a Seurat object containing the RNA adata
pbmc <- CreateSeuratObject(
counts = counts$`Gene Expression`,
assay = "RNA"
)
# create ATAC assay and add it to the object
pbmc[["ATAC"]] <- CreateChromatinAssay(
counts = counts$Peaks,
sep = c(":", "-"),
fragments = fragpath,
annotation = annotation
)
pbmc
質(zhì)控
我們可以通過DNA可及性數(shù)據(jù)來評(píng)估每個(gè)細(xì)胞的質(zhì)量控制指標(biāo)，并排除那些指標(biāo)異常的細(xì)胞办成。此外，對(duì)于那些在RNA或ATAC檢測中計(jì)數(shù)特別低或特別高的細(xì)胞某弦，我們也會(huì)進(jìn)行剔除而克。
DefaultAssay(pbmc) <- "ATAC"
pbmc <- NucleosomeSignal(pbmc)
pbmc <- TSSEnrichment(pbmc)
對(duì)象數(shù)據(jù)中變量之間的相互關(guān)系可以通過 DensityScatter() 函數(shù)來直觀展示员萍。此外，設(shè)置 quantiles=TRUE 選項(xiàng)螃壤，可以幫助我們迅速確定不同質(zhì)量控制指標(biāo)的適宜閾值筋帖。
DensityScatter(pbmc, x = 'nCount_ATAC', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)
VlnPlot(
object = pbmc,
features = c("nCount_RNA", "nCount_ATAC", "TSS.enrichment", "nucleosome_signal"),
ncol = 4,
pt.size = 0
)
# filter out low quality cells
pbmc <- subset(
x = pbmc,
subset = nCount_ATAC < 100000 &
nCount_RNA < 25000 &
nCount_ATAC > 1800 &
nCount_RNA > 1000 &
nucleosome_signal < 2 &
TSS.enrichment > 1
)
pbmc
基因表達(dá)數(shù)據(jù)處理
我們可以使用 SCTransform 對(duì)基因表達(dá)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化蚁滋，并使用 PCA 降低維度赘淮。
DefaultAssay(pbmc) <- "RNA"
pbmc <- SCTransform(pbmc)
pbmc <- RunPCA(pbmc)
DNA可及性數(shù)據(jù)處理
在這里，我們通過執(zhí)行潛在語義索引 ( LSI )走诞，以處理 scATAC-seq 數(shù)據(jù)集的相同方式處理 DNA 可及性檢測蛤高。
DefaultAssay(pbmc) <- "ATAC"
pbmc <- FindTopFeatures(pbmc, min.cutoff = 5)
pbmc <- RunTFIDF(pbmc)
pbmc <- RunSVD(pbmc)
注釋細(xì)胞類型
為了注釋數(shù)據(jù)集中的細(xì)胞類型，我們可以使用 Seurat 包中的工具塞绿，將細(xì)胞標(biāo)簽從現(xiàn)有的 PBMC 參考數(shù)據(jù)集中轉(zhuǎn)移過來恤批。
我們將使用 Hao 等人（2020年）的注釋 PBMC 參考數(shù)據(jù)集，可以從這里下載：https://atlas.fredhutch.org/data/nygc/multimodal/pbmc_multimodal.h5seurat
請(qǐng)注意诀浪，加載參考數(shù)據(jù)集需要安裝 SeuratDisk 包。
library(SeuratDisk)
# load PBMC reference
reference <- LoadH5Seurat("../vignette_data/multiomic/pbmc_multimodal.h5seurat", assays = list("SCT" = "counts"), reductions = 'spca')
reference <- UpdateSeuratObject(reference)
DefaultAssay(pbmc) <- "SCT"
# transfer cell type labels from reference to query
transfer_anchors <- FindTransferAnchors(
reference = reference,
query = pbmc,
normalization.method = "SCT",
reference.reduction = "spca",
recompute.residuals = FALSE,
dims = 1:50
)
predictions <- TransferData(
anchorset = transfer_anchors,
refdata = reference$celltype.l2,
weight.reduction = pbmc[['pca']],
dims = 1:50
)
pbmc <- AddMetaData(
object = pbmc,
metadata = predictions
)
# set the cell identities to the cell type predictions
Idents(pbmc) <- "predicted.id"
# remove low-quality predictions
pbmc <- pbmc[, pbmc$prediction.score.max > 0.5]
聯(lián)合 UMAP 可視化
使用 Seurat v4 中的加權(quán)最近鄰方法睛竣，我們可以計(jì)算代表基因表達(dá)和 DNA 可及性測量的UMAP圖射沟。
# build a joint neighbor graph using both assays
pbmc <- FindMultiModalNeighbors(
object = pbmc,
reduction.list = list("pca", "lsi"),
dims.list = list(1:50, 2:40),
modality.weight.name = "RNA.weight",
verbose = TRUE
)
# build a joint UMAP visualization
pbmc <- RunUMAP(
object = pbmc,
nn.name = "weighted.nn",
assay = "RNA",
verbose = TRUE
)
DimPlot(pbmc, label = TRUE, repel = TRUE, reduction = "umap") + NoLegend()
將峰與基因聯(lián)系起來
為了找到可能調(diào)控每個(gè)基因的峰值集合躏惋，我們可以計(jì)算基因表達(dá)與其附近峰值可及性之間的相關(guān)性嚷辅，并校正由于 GC 含量、整體可及性和峰值大小引起的偏差簸搞。
在整個(gè)基因組上執(zhí)行這一步驟可能非常耗時(shí)，因此我們在這里以部分基因?yàn)槔虺穑故痉?基因鏈接寺擂。通過省略 genes.use 參數(shù)，可以使用相同的函數(shù)來找到所有基因的鏈接：
DefaultAssay(pbmc) <- "ATAC"
# first compute the GC content for each peak
pbmc <- RegionStats(pbmc, genome = BSgenome.Hsapiens.UCSC.hg38)
# link peaks to genes
pbmc <- LinkPeaks(
object = pbmc,
peak.assay = "ATAC",
expression.assay = "SCT",
genes.use = c("LYZ", "MS4A1")
)
我們可以使用 CoveragePlot() 函數(shù)可視化這些鏈接垦细，或者我們可以在交互式分析中使用 CoverageBrowser() 函數(shù)：
idents.plot <- c("B naive", "B intermediate", "B memory",
"CD14 Mono", "CD16 Mono", "CD8 TEM", "CD8 Naive")
p1 <- CoveragePlot(
object = pbmc,
region = "MS4A1",
features = "MS4A1",
expression.assay = "SCT",
idents = idents.plot,
extend.upstream = 500,
extend.downstream = 10000
)
p2 <- CoveragePlot(
object = pbmc,
region = "LYZ",
features = "LYZ",
expression.assay = "SCT",
idents = idents.plot,
extend.upstream = 8000,
extend.downstream = 5000
)
patchwork::wrap_plots(p1, p2, ncol = 1)
本文由mdnice多平臺(tái)發(fā)布
