前言 Seurat 3.X版本能够整合scRNA-seq和scATAC-seq, 主要体现在:
整合步骤包括如下步骤:
从ATAC-seq中估计RNA-seq表达水平,即从ATAC-seq reads定量基因表达活跃度
使用LSI学习ATAC-seq数据的内部结构
鉴定ATAC-seq和RNA-seq数据集的锚点
数据集间进行转移,包括聚类的标签,在ATAC-seq数据中推测RNA水平用于共嵌入分析
数据下载 测试数据下载地址:
可以复制下载链接到浏览器下载,也可以直接在R语言用download.file
中进行下载。
1 # 下载peak 2 atac_peak <- "http://cf./samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5" 3 atac_peak_file <- "atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5" 4 download.file(atac_peak, atac_peak_file) 5 6 # 下载singlecell.csv 7 singlecell <- "http://cf./samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_singlecell.csv" 8 singlecell_file <- "atac_v1_pbmc_10k_singlecell.csv" 9 download.file(singlecell, singlecell_file)10 11 # 下载rna-seq 12 rna_seq <- "http://cf./samples/cell-exp/3.0.0/pbmc_10k_v3/pbmc_10k_v3_filtered_feature_bc_matrix.h5" 13 rna_seq_file <- "pbmc_10k_v3_filtered_feature_bc_matrix.h5" 14 download.file(rna_seq, rna_seq_file)15 16 # 下载GTF 17 gtf <- "ftp://ftp.ensembl.org/pub/grch37/release-82/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz" 18 gtf_file <- "Homo_sapiens.GRCh37.82.gtf.gz" 19 download.file(gtf, gtf_file)
基因活跃度定量 首先,先将peak矩阵转成基因活跃度矩阵。Seurat做了一个简单的假设,基因活跃度可以通过简单的将落在基因区和其上游2kb的count相加得到基因活跃度,并且这个结果Cicero等工具返回gene-by-cell矩阵是类似的。
1 library (Seurat)2 library (ggplot2)3 # 读取peak 4 peaks <- Read10X_h5(atac_peak_file)5 activity.matrix <- CreateGeneActivityMatrix(peak.matrix = peaks, 6 annotation.file = gtf_file, 7 seq.levels = c(1 :22 , "X" , "Y" ), 8 upstream = 2000 , 9 verbose = TRUE )
activity.matrix是一个dgCMatrix
对象,其中行为基因,列为细胞。因此如果对于Cicero的输出结果,只要提供相应的矩阵数据结构即可。
设置对象 下一步,我们要来设置Seurat
对象,将原始的peak counts保存到assay中,命名为"ATAC"
1 pbmc.atac <- CreateSeuratObject(counts = peaks, assay = "ATAC" , project = "10x_ATAC" )
增加基因活跃度矩阵,命名为"ACTIVITY".
1 pbmc.atac[["ACTIVITY" ]] <- CreateAssayObject(counts = activity.matrix)
增加细胞的元信息,该信息来自于scATAC-seq的CellRanger处理的singlecell.csv
1 meta <- read.table(singlecell_file, sep = "," , header = TRUE , row.names = 1 , 2 stringsAsFactors = FALSE )3 meta <- meta[colnames(pbmc.atac), ]4 pbmc.atac <- AddMetaData(pbmc.atac, metadata = meta)
过滤掉scATAC-seq数据中总count数低于5000的细胞,这个阈值需要根据具体实验设置
1 pbmc.atac <- subset(pbmc.atac, subset = nCount_ATAC > 5000 )2 pbmc.atac$tech <- "atac"
数据预处理 为了找到scATAC-seq数据集和scRNA-seq数据集之间的锚定点,我们需要对基因活跃度矩阵进行预处理
设置pbmc.atac的默认Assay为"ACTIVITY", 然后寻找高表达的基因,对基因活跃度矩阵进行标准化和Scale。
1 DefaultAssay(pbmc.atac) <- "ACTIVITY" 2 pbmc.atac <- FindVariableFeatures(pbmc.atac)3 pbmc.atac <- NormalizeData(pbmc.atac)4 pbmc.atac <- ScaleData(pbmc.atac)
同样的,我们还需要对peak矩阵进行处理。这里我们用的隐语义(Latent semantic indexing, LSI)方法对scATAC-seq数据进行降维。该步骤学习scRNA-seq数据的内部结构,并且在转换信息时对锚点恰当权重的决定很重要。
根据 Cusanovich et al, Science 2015提出的LSI方法,他们搞了一个RunLSI
函数。LSI的计算方法为TF-IDF加SVD。
我们使用在所有细胞中至少有100个read的peak,然后降维到50。该参数的选择受之前的scATAC-seq研究启发,所以没有更改,当然你你可以把它改了。
1 DefaultAssay(pbmc.atac) <- "ATAC" 2 VariableFeatures(pbmc.atac) <- names(which(Matrix::rowSums(pbmc.atac) > 100 ))3 pbmc.atac <- RunLSI(pbmc.atac, n = 50 , scale.max = NULL )4 #pbmc.atac <- RunUMAP(pbmc.atac, reduction = "lsi", dims = 1:50) 5 pbmc.atac <- RunTSNE(pbmc.atac, reduction = "lsi" , dims = 1 :50 )
注 : 要将pbmc.atac的默认assay切换成"ATAC", 非线性降维可以选择UMAP或者t-SNE。
我们之前使用过Seurat对scRNA-seq数据进行预处理和聚类,下载地址为Dropbox。(https://www./s/3f3p5nxrn5b3y4y/pbmc_10k_v3.rds?dl=1 )
1 pbmc.rna <- readRDS("pbmc_10k_v3.rds" )2 pbmc.rna$tech <- "rna"
将scRNA-seq和scATAC-seq共同展示,对一些骨髓(myeloid)和淋巴(lymphoid)PBMC中比较常见的聚类,其实是能从图中看出来。
1 p1 <- DimPlot(pbmc.atac, reduction = "tsne" ) + 2 NoLegend() + 3 ggtitle("scATAC-seq" )4 p2 <- DimPlot(pbmc.rna, group.by = "celltype" , reduction = "tsne" , label = TRUE , repel = TRUE ) + 5 NoLegend() + 6 ggtitle("scRNA-seq" )7 CombinePlots(plots = list(p1, p2))
标签转移前 现在,我们用FindTransferAnchors
鉴定scATAC-seq数据集和scRNA-seq数据集的锚点,然后根据这些锚点将 10K scRNA-seq数据中鉴定的细胞类型标记转移到scATAC-seq细胞中。
1 transfer.anchors <- FindTransferAnchors(reference = pbmc.rna, 2 query = pbmc.atac, 3 features = VariableFeatures(object = pbmc.rna), 4 reference.assay = "RNA" , 5 query.assay = "ACTIVITY" , 6 reduction = "cca" )
Seurat比较的是scRNA-seq表达量矩阵和scATAC-seq中基因活跃度矩阵,利用CCA降维方法比较两者在scRNA-seq中的高变异基因的关系
为了转移细胞类群的编号,我们需要一组之前注释过的细胞类型,作为TransferData
的 refdata 参数输入。TransferData
本质上采用的是KNN算法,利用距离未知类型细胞的最近N个已知细胞所属的类型来定义该细胞。weight.reduction
参数是用来选择设置权重的降维。
1 celltype.predictions <- TransferData(anchorset = transfer.anchors, 2 refdata = pbmc.rna$celltype, 3 weight.reduction = pbmc.atac[["lsi" ]])4 pbmc.atac <- AddMetaData(pbmc.atac, metadata = celltype.predictions)
我们可以检查每个细胞预测得分的分布情况,选择性的过滤哪些得分比较低的细胞。我们发现超过95%的细胞得分大于或等于0.5.
1 hist(pbmc.atac$prediction.score.max)2 abline(v = 0.5 , col = "red" )
得分分布 1 table(pbmc.atac$prediction.score.max > 0.5 )2 # FALSE TRUE 3 # 383 7483
之后,我们就可以在UMAP上检查预测的细胞类型的分布,检查细胞类型在scRNA-seq和scATAC-seq中的相对位置
1 pbmc.atac.filtered <- subset(pbmc.atac, 2 subset = prediction.score.max > 0.5 ) 3 pbmc.atac.filtered$predicted.id <- factor(pbmc.atac.filtered$predicted.id, 4 levels = levels(pbmc.rna)) # to make the colors match 5 p1 <- DimPlot(pbmc.atac.filtered, 6 group.by = "predicted.id" , 7 label = TRUE , 8 repel = TRUE ) + 9 ggtitle("scATAC-seq cells" ) + 10 NoLegend() + 11 scale_colour_hue(drop = FALSE )12 p2 <- DimPlot(pbmc.rna, group.by = "celltype" , label = TRUE , repel = TRUE ) +13 ggtitle("scRNA-seq cells" ) + 14 NoLegend()15 CombinePlots(plots = list(p1, p2))
预测后结果 在转移细胞类型标记之后,我们就可以进行细胞特意水平上的下游分析。举个例子,我们可以去找一些某些细胞类型特异的增强子,寻找富集的motif。目前这些分析Seurat还不直接支持,还在调试中。
共嵌入(co-embedding) 最后,如果你想将所有的细胞一同展示,可以将scRNA-seq和scATAC-seq数据嵌入到相同的低维空间。
我们使用之前的锚点从scATAC-seq细胞中推断RNA-seq的值,后续分析就相当于两个单细胞数据的分析流程。
注意: 这一步只是为了可视化,其实不做也行。
选择用于估计的基因,可以是高变动基因,也可以是所有基因。
1 # 只选择高变动的基因作为参考 2 genes.use <- VariableFeatures(pbmc.rna)3 refdata <- GetAssayData(pbmc.rna, assay = "RNA" , slot = "data" )[genes.use, ]
之后利用TransferData
推断scATAC-seq在这些基因中的可能值,输出结果就是ATAC细胞对应的scRNA-seq矩阵
1 imputation <- TransferData(anchorset = transfer.anchors, 2 refdata = refdata, 3 weight.reduction = pbmc.atac[["lsi" ]])4 # this line adds the imputed data matrix to the pbmc.atac object 5 pbmc.atac[["RNA" ]] <- imputation
合并两个的结果,然后就是scRNA-seq的常规分析。
1 coembed <- merge(x = pbmc.rna, y = pbmc.atac)2 # 标准化分析 3 coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE )4 coembed <- RunPCA(coembed, features = genes.use, verbose = FALSE )5 #coembed <- RunUMAP(coembed, dims = 1:30) 6 coembed <- RunTSNE(coembed, dims = 1 :30 )7 coembed$celltype <- ifelse(!is.na(coembed$celltype), coembed$celltype, coembed$predicted.id)
在t-SNE上绘制结果
1 p1 <- DimPlot(coembed, reduction="tsne" , group.by = "tech" )2 p2 <- DimPlot(coembed, reduction="tsne" , group.by = "celltype" , label = TRUE , repel = TRUE )3 CombinePlots(list(p1, p2))
tSNE plot 从上面的结果中,你可能会发现某些细胞只有在一类技术中存在。举个例子,从巨噬细胞(megakaryocytes)到成熟的血小板细胞(pletelet)由于没有细胞核,无法被scATAC-seq技术检测到。我们可以单独展示每个技术,进行检查
1 DimPlot(coembed, reduction="tsne" , split.by = "tech" , group.by = "celltype" , label = TRUE , repel = TRUE ) + NoLegend()
分别展示 此外,你还可以发现B细胞前体类型附近有一类细胞只由scATAC-seq的细胞构成。通过展示这些细胞在CellRanger分析结果提供的黑名单位置所占read数,可以推测出这类细胞可能是死细胞,或者是其他scRNA-seq无法检测的人为因素。
1 coembed$blacklist_region_fragments[is.na(coembed$blacklist_region_fragments)] <- 0 2 FeaturePlot(coembed, features = "blacklist_region_fragments" , max.cutoff = 500 )
黑名单区 Seurat这种基于基因活跃度得分进行细胞类型预测,是否靠谱,开发者提供了如下几个证据
总体预测得分(pbmc.atac$prediction.score.max
)高,意味着用scRNA-seq定义细胞类型比较可靠
我们可以在scATC-seq降维结果中
利用相同锚点的贡嵌入分析,发现两类形态能很好的混合
将ATAC-seq数据根据聚类结果构建pseduo bulk, 发现和真实的bulk数据近似
参考资料: