Seurat的scATAC-seq分析流程

Seurat 3.X版本能夠整合scRNA-seq和scATAC-seq, 主要體現(xiàn)在：
  1. 基于scRNA-seq的聚類結(jié)果對scATAC-seq的細(xì)胞進(jìn)行聚類
  2. scRNA-seq和scATAC-seq共嵌入(co-embed)分析
整合步驟包括如下步驟:
  1. 從ATAC-seq中估計RNA-seq表達(dá)水平，即從ATAC-seq reads定量基因表達(dá)活躍度
  2. 使用LSI學(xué)習(xí)ATAC-seq數(shù)據(jù)的內(nèi)部結(jié)構(gòu)
  3. 鑒定ATAC-seq和RNA-seq數(shù)據(jù)集的錨點
  4. 數(shù)據(jù)集間進(jìn)行轉(zhuǎn)移，包括聚類的標(biāo)簽，在ATAC-seq數(shù)據(jù)中推測RNA水平用于共嵌入分析
    數(shù)據(jù)下載
    可以復(fù)制下載鏈接到瀏覽器下載，也可以直接在R語言用download.file中進(jìn)行下載。
# 下載peak
atac_peak <- "http://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5"
atac_peak_file <- "atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5"
download.file(atac_peak, atac_peak_file)

# 下載singlecell.csv
singlecell <- "http://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_singlecell.csv"
singlecell_file <- "atac_v1_pbmc_10k_singlecell.csv"
download.file(singlecell, singlecell_file)

# 下載rna-seq
rna_seq <- "http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_10k_v3/pbmc_10k_v3_filtered_feature_bc_matrix.h5"
rna_seq_file <- "pbmc_10k_v3_filtered_feature_bc_matrix.h5"
download.file(rna_seq, rna_seq_file)

# 下載GTF
gtf <- "ftp://ftp.ensembl.org/pub/grch37/release-82/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz"
gtf_file <- "Homo_sapiens.GRCh37.82.gtf.gz"
download.file(gtf, gtf_file)
基因活躍度定量

首先，先將peak矩陣轉(zhuǎn)成基因活躍度矩陣。Seurat做了一個簡單的假設(shè)，基因活躍度可以通過簡單的將落在基因區(qū)和其上游2kb的count相加得到基因活躍度，并且這個結(jié)果Cicero等工具返回gene-by-cell矩陣是類似的。

library(Seurat)
library(ggplot2)
# 讀取peak
peaks <- Read10X_h5(atac_peak_file)
activity.matrix <- CreateGeneActivityMatrix(peak.matrix = peaks, 
                                            annotation.file = gtf_file, 
                                            seq.levels = c(1:22, "X", "Y"), 
                                            upstream = 2000, 
                                            verbose = TRUE)

activity.matrix是一個dgCMatrix對象，其中行為基因，列為細(xì)胞。因此如果對于Cicero的輸出結(jié)果，只要提供相應(yīng)的矩陣數(shù)據(jù)結(jié)構(gòu)即可。

設(shè)置對象

下一步，我們要來設(shè)置Seurat對象，將原始的peak counts保存到assay中，命名為"ATAC"

pbmc.atac <- CreateSeuratObject(counts = peaks, assay = "ATAC", project = "10x_ATAC")

增加基因活躍度矩陣，命名為"ACTIVITY".

pbmc.atac[["ACTIVITY"]] <- CreateAssayObject(counts = activity.matrix)

增加細(xì)胞的元信息，該信息來自于scATAC-seq的CellRanger處理的singlecell.csv

meta <- read.table(singlecell_file, sep = ",", header = TRUE, row.names = 1, 
    stringsAsFactors = FALSE)
meta <- meta[colnames(pbmc.atac), ]
pbmc.atac <- AddMetaData(pbmc.atac, metadata = meta)

過濾掉scATAC-seq數(shù)據(jù)中總count數(shù)低于5000的細(xì)胞，這個閾值需要根據(jù)具體實驗設(shè)置

pbmc.atac <- subset(pbmc.atac, subset = nCount_ATAC > 5000)
pbmc.atac$tech <- "atac"
數(shù)據(jù)預(yù)處理

為了找到scATAC-seq數(shù)據(jù)集和scRNA-seq數(shù)據(jù)集之間的錨定點，我們需要對基因活躍度矩陣進(jìn)行預(yù)處理

設(shè)置pbmc.atac的默認(rèn)Assay為"ACTIVITY"， 然后尋找高表達(dá)的基因，對基因活躍度矩陣進(jìn)行標(biāo)準(zhǔn)化和Scale。

DefaultAssay(pbmc.atac) <- "ACTIVITY"
pbmc.atac <- FindVariableFeatures(pbmc.atac)
pbmc.atac <- NormalizeData(pbmc.atac)
pbmc.atac <- ScaleData(pbmc.atac)

同樣的，我們還需要對peak矩陣進(jìn)行處理。這里我們用的隱語義(Latent semantic indexing, LSI)方法對scATAC-seq數(shù)據(jù)進(jìn)行降維。該步驟學(xué)習(xí)scRNA-seq數(shù)據(jù)的內(nèi)部結(jié)構(gòu)，并且在轉(zhuǎn)換信息時對錨點恰當(dāng)權(quán)重的決定很重要。

根據(jù) Cusanovich et al, Science 2015提出的LSI方法，他們搞了一個RunLSI函數(shù)。LSI的計算方法為TF-IDF加SVD。

我們使用在所有細(xì)胞中至少有100個read的peak，然后降維到50。該參數(shù)的選擇受之前的scATAC-seq研究啟發(fā)，所以沒有更改，當(dāng)然你你可以把它改了。

DefaultAssay(pbmc.atac) <- "ATAC"
VariableFeatures(pbmc.atac) <- names(which(Matrix::rowSums(pbmc.atac) > 100))
pbmc.atac <- RunLSI(pbmc.atac, n = 50, scale.max = NULL)
#pbmc.atac <- RunUMAP(pbmc.atac, reduction = "lsi", dims = 1:50)
pbmc.atac <- RunTSNE(pbmc.atac, reduction = "lsi", dims = 1:50)

注: 要將pbmc.atac的默認(rèn)assay切換成"ATAC", 非線性降維可以選擇UMAP或者t-SNE。

我們之前使用過Seurat對scRNA-seq數(shù)據(jù)進(jìn)行預(yù)處理和聚類，下載地址為Dropbox

pbmc.rna <- readRDS("pbmc_10k_v3.rds")
pbmc.rna$tech <- "rna"

將scRNA-seq和scATAC-seq共同展示，對一些骨髓(myeloid)和淋巴(lymphoid)PBMC中比較常見的聚類，其實是能從圖中看出來。

p1 <- DimPlot(pbmc.atac, reduction = "tsne") + 
        NoLegend() + 
        ggtitle("scATAC-seq")
p2 <- DimPlot(pbmc.rna, group.by = "celltype", reduction = "tsne", label = TRUE, repel = TRUE) + 
        NoLegend() + 
        ggtitle("scRNA-seq")
CombinePlots(plots = list(p1, p2))
image.png

現(xiàn)在，我們用FindTransferAnchors鑒定scATAC-seq數(shù)據(jù)集和scRNA-seq數(shù)據(jù)集的錨點，然后根據(jù)這些錨點將 10K scRNA-seq數(shù)據(jù)中鑒定的細(xì)胞類型標(biāo)記轉(zhuǎn)移到scATAC-seq細(xì)胞中。

transfer.anchors <- FindTransferAnchors(reference = pbmc.rna, 
                                        query = pbmc.atac, 
                                        features = VariableFeatures(object = pbmc.rna), 
                                        reference.assay = "RNA", 
                                        query.assay = "ACTIVITY", 
                                        reduction = "cca")

Seurat比較的是scRNA-seq表達(dá)量矩陣和scATAC-seq中基因活躍度矩陣，利用CCA降維方法比較兩者在scRNA-seq中的高變異基因的關(guān)系

為了轉(zhuǎn)移細(xì)胞類群的編號，我們需要一組之前注釋過的細(xì)胞類型，作為TransferData的 refdata 參數(shù)輸入。TransferData本質(zhì)上采用的是KNN算法，利用距離未知類型細(xì)胞的最近N個已知細(xì)胞所屬的類型來定義該細(xì)胞。weight.reduction參數(shù)是用來選擇設(shè)置權(quán)重的降維。

celltype.predictions <- TransferData(anchorset = transfer.anchors, 
                                     refdata = pbmc.rna$celltype, 
                                     weight.reduction = pbmc.atac[["lsi"]])
pbmc.atac <- AddMetaData(pbmc.atac, metadata = celltype.predictions)

我們可以檢查每個細(xì)胞預(yù)測得分的分布情況，選擇性的過濾哪些得分比較低的細(xì)胞。我們發(fā)現(xiàn)超過95%的細(xì)胞得分大于或等于0.5.

hist(pbmc.atac$prediction.score.max)
abline(v = 0.5, col = "red")
image.png
table(pbmc.atac$prediction.score.max > 0.5)
# FALSE  TRUE 
#   383  7483

之后，我們就可以在UMAP上檢查預(yù)測的細(xì)胞類型的分布，檢查細(xì)胞類型在scRNA-seq和scATAC-seq中的相對位置

pbmc.atac.filtered <- subset(pbmc.atac, 
                             subset = prediction.score.max > 0.5)
pbmc.atac.filtered$predicted.id <- factor(pbmc.atac.filtered$predicted.id, 
                                          levels = levels(pbmc.rna))  # to make the colors match
p1 <- DimPlot(pbmc.atac.filtered, 
              group.by = "predicted.id",
              label = TRUE, 
              repel = TRUE) + 
        ggtitle("scATAC-seq cells") + 
        NoLegend() + 
        scale_colour_hue(drop = FALSE)
p2 <- DimPlot(pbmc.rna, group.by = "celltype", label = TRUE, repel = TRUE) +
        ggtitle("scRNA-seq cells") + 
        NoLegend()
CombinePlots(plots = list(p1, p2))
image.png

在轉(zhuǎn)移細(xì)胞類型標(biāo)記之后，我們就可以進(jìn)行細(xì)胞特意水平上的下有分析。舉個例子，我們可以去找一些某些細(xì)胞類型特異的增強(qiáng)子，尋找富集的motif。目前這些分析Seurat還不直接支持，還在調(diào)試中。

共嵌入(co-embedding)

最后，如果你想將所有的細(xì)胞一同展示，可以將scRNA-seq和scATAC-seq數(shù)據(jù)嵌入到相同的低維空間。

我們使用之前的錨點從scATAC-seq細(xì)胞中推斷RNA-seq的值，后續(xù)分析就相當(dāng)于兩個單細(xì)胞數(shù)據(jù)的分析流程。

注意: 這一步只是為了可視化，其實不做也行。

選擇用于估計的基因，可以是高變動基因，也可以是所有基因。

# 只選擇高變動的基因作為參考
genes.use <- VariableFeatures(pbmc.rna)
data <- GetAssayData(pbmc.rna, assay = "RNA", slot = "data")[genes.use, ]

之后利用TransferData推斷scATAC-seq在這些基因中的可能值，輸出結(jié)果就是ATAC細(xì)胞對應(yīng)的scRNA-seq矩陣

imputation <- TransferData(anchorset = transfer.anchors, 
                           refdata = refdata, 
                           weight.reduction = pbmc.atac[["lsi"]])
# this line adds the imputed data matrix to the pbmc.atac object
pbmc.atac[["RNA"]] <- imputation

合并兩個的結(jié)果，然后就是scRNA-seq的常規(guī)分析。

coembed <- merge(x = pbmc.rna, y = pbmc.atac)
# 標(biāo)準(zhǔn)化分析
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)
coembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
#coembed <- RunUMAP(coembed, dims = 1:30)
coembed <- RunTSNE(coembed, dims = 1:30)
coembed$celltype <- ifelse(!is.na(coembed$celltype), coembed$celltype, coembed$predicted.id)

在t-SNE上繪制結(jié)果

p1 <- DimPlot(coembed, reduction="tsne", group.by = "tech")
p2 <- DimPlot(coembed, reduction="tsne", group.by = "celltype", label = TRUE, repel = TRUE)
CombinePlots(list(p1, p2))
image.png

從上面的結(jié)果中，你可能會發(fā)現(xiàn)某些細(xì)胞只有在一類技術(shù)中存在。舉個例子，從巨噬細(xì)胞(megakaryocytes)到成熟的血小板細(xì)胞(pletelet)由于沒有細(xì)胞核，無法被scATAC-seq技術(shù)檢測到。我們可以單獨(dú)展示每個技術(shù)，進(jìn)行檢查

DimPlot(coembed, reduction="tsne", split.by = "tech", group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend()
image.png

此外，你還可以發(fā)現(xiàn)B細(xì)胞前體類型附近有一類細(xì)胞只由scATAC-seq的細(xì)胞構(gòu)成。通過展示這些細(xì)胞在CellRanger分析結(jié)果提供的黑名單位置所占read數(shù)，可以推測出這類細(xì)胞可能是死細(xì)胞，或者是其他scRNA-seq無法檢測的人為因素。

coembed$blacklist_region_fragments[is.na(coembed$blacklist_region_fragments)] <- 0
FeaturePlot(coembed, features = "blacklist_region_fragments", max.cutoff = 500)
image.png

Seurat這種基于基因活躍度得分進(jìn)行細(xì)胞類型預(yù)測，是否靠譜，開發(fā)者提供了如下幾個證據(jù)

總體預(yù)測得分(pbmc.atac$prediction.score.max)高，意味著用scRNA-seq定義細(xì)胞類型比較可靠
我們可以在scATC-seq降維結(jié)果中
利用相同錨點的貢嵌入分析，發(fā)現(xiàn)兩類形態(tài)能很好的混合
將ATAC-seq數(shù)據(jù)根據(jù)聚類結(jié)果構(gòu)建pseduo bulk, 發(fā)現(xiàn)和真實的bulk數(shù)據(jù)近似

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容