cisTopic——从scATAC-seq数据中同时识别细胞状态和顺式调控主题的R包

简介

主页： https://github.com/aertslab/cisTopic
cisTopic——从scATAC-seq数据中同时识别细胞状态和顺式调控主题的R包。

cisTopic是一个用于从单细胞表观基因组学数据同时识别顺式调控主题和细胞状态的R Bioconductor package。它依赖于Latent Dirichlet Allocation (LDA)算法，这种算法是一种有效的用于文本挖掘的贝叶斯方法，将相似主题和相关单词的文档分组到主题中。

cisTopic包括四个主要步骤: (1)生成一个二进制可及性矩阵作为LDA的输入;(2) LDA和模型选择;(3)利用LDA中的主题细胞分布进行细胞状态识别，(4)探索区域主题分布。

推荐理由：联合优化了细胞聚类和增强子的分类，可以识别细胞亚型和可及性的增强子。

软件安装

建议使用R3.6以上的版本

devtools::install_github("aertslab/RcisTarget")
devtools::install_github("aertslab/AUCell")
devtools::install_github("aertslab/cisTopic")

安装过程或许没那么顺利，我的测试环境是linux，平时倾向于使用conda install 来安装R包，因为它可以帮你解决一些依赖。

数据

目前该软件仅支持人、小鼠和果蝇。这里使用cellranger atac 的outs结果作为测试数据。下面链接可下载练习数据（物种：人）

测试

读取数据

#读取上述测试数据
 pathTo10X <- '/public/analysis/xxx/Research/cisTopic/sample_hsa/outs/'
 data_folder <- paste0(pathTo10X, 'filtered_peak_bc_matrix')
 metrics <- paste0(pathTo10X,'singlecell.csv')
 cisTopicObject <- createcisTopicObjectFrom10Xmatrix(data_folder, metrics,  project.name='test')

读取成功后可以查看一下数据@cell.data and @region.data

head(cisTopicObject@cell.data) 

添加一些其他信息

使用addRegionMetadata可添加一些其他信息，比如graph-based clusters。

pathTographBasedClusters_CRA <- paste0(pathTo10X, 'analysis/clustering/graphclust/clusters.csv') 
graphBasedClusters_CRA <- read.table(pathTographBasedClusters_CRA, sep=',', header=TRUE, row.names = 1)
colnames(graphBasedClusters_CRA) <- 'graphBasedClusters_CRA'
graphBasedClusters_CRA[,1] <- as.factor(graphBasedClusters_CRA[,1])
cisTopicObject <- addCellMetadata(cisTopicObject, graphBasedClusters_CRA)

加上之后，cisTopicObject@cell.data 中多了一列graphBasedClusters_CRA

构建模型

cisTopicObject2 <- runWarpLDAModels(cisTopicObject, topic=30, seed=987, nCores=30, iterations = 500, addModels=FALSE)
#INFO  [17:44:45.757] early stopping at 50 iteration 迭代500次但实际只进行了50次

#查看一下所用的参数
head(cisTopicObject@calc.params)

构建模型这一步，运行约30min。在测试时若设置topic=c(20,30,40)这样子同时跑多个模型再去筛选最佳模型，运行失败，提示是Matrix的问题，暂时还没解决。

关于参数的理解可以使用？runWarpLDAModels 进行查看，其他函数同理。

选择最佳topic数

pdf(file="/public/analysis/xxx/Research/cisTopic/results/selectModel300707.pdf")
cisTopicObject <- selectModel(cisTopicObject, type='maximum') #有3中筛选方法，默认maximum
dev.off()

跑完selectModel之后，cisTopicObject2@log.lik 存储了model的信息。

对模型进行分析

上一步的LDA返回两个分布：
(1) the topic contributions per cell （cells为列，topics为行，贡献(contributions)为值。）
(2) the region contribution to a topic（topics为列，regions为行，贡献(contributions)为值。）

第一部分：对cisTopic-cell分布进行分析

针对每个分布，都可以使用不同的方法进行聚类和可视化。cisTopic包含包装器函数，可轻松运行Umap、tSNE、 diffussion maps 和PCA（结果保存在slot@dr$cell中），对于cell 和region两个部分，只需要通过指定 target=‘cell’ 就可以区分。下面是针对cell 部分进行降维与绘图：

cisTopicObject <- runtSNE(cisTopicObject, target='cell', seed=123, pca=FALSE, method='Probability')
cisTopicObject <- runUmap(cisTopicObject, target='cell', seed=123, method='Probability')

定义clusters

我们可以使用 cell-topic matrix来定义clusters 。

cellassign <- modelMatSelection(cisTopicObject, 'cell', 'Probability')
set.seed(123)
library(Rtsne)
DR <- Rtsne(t(cellassign), pca=F)
DRdist <- dist(DR$Y)