SEACells Tutorial（一）：基于SEACells的scATAC-seq经典下游分析

写在前面：
目前，单细胞基因组学数据的细胞分辨率与大多数分析的集群级分辨率之间存在着根本的脱节，简单来说就是在聚类分析的过程中，通常会导致集群级分辨率，这意味着一组细胞被归为同一类别，而忽略了它们内部的细微差异。
为了克服数据的稀疏性和噪声的问题，我们通常会把metacells作为一个有效的解决办法。通俗来说，metacells是单细胞测序数据的细胞分组，即在不同深度下对单细胞进行聚合，其细胞内的变异都是由技术因素引起的，而非生物学因素，这样就使得细胞的均匀性大大提高。
2023年3月27日，在nature biotechnology上发布了一边文章，作者给出了一种计算metacells的算法，可从单细胞基因组学数据中推断出转录和表观基因组的细胞状态。在这里我们对作者发布的SEACells进行复现，期望能够在实践中加以应用。

备注：本教程代码部分在Google Colab运行环境中完成，Google Colab是Google提供的一个Jupyter Notebook式的交互环境。

---

---

概述

本文章主要为利用SEAcells进行scATAC-seq的下游分析提供一个教程，包含了四部分的内容：gene-peak associations的计算、highly regulated genes的识别、ATAC gene 得分计算以及基因可及性得分（gene accessibility scores），这大概包含了ATAC分析下游的大部分分析流程。

---

一、准备工作

#载入包
import numpy as np
import pandas as pd
import scanpy as sc
import SEACells
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

#画图设置
%matplotlib inline

sns.set_style('ticks')
matplotlib.rcParams['figure.figsize'] = [4, 4]
matplotlib.rcParams['figure.dpi'] = 100

二、加载数据

我们使用 scanpy Anndata 对象作为加载和过滤数据的首选模式。

# !mkdir data/
# !wget https://dp-lab-data-public.s3.amazonaws.com/SEACells-multiome/cd34_multiome_rna.h5ad -O data/cd34_multiome_rna.h5ad # RNA data
# !wget https://dp-lab-data-public.s3.amazonaws.com/SEACells-multiome/cd34_multiome_atac.h5ad -O data/cd34_multiome_atac.h5ad # ATAC data

该数据集包含 RNA 和 ATAC 模式作为两个不同的 Anndata 对象。ATAC 数据集包含预先计算的 SEACell 元细胞
（如何计算seacell后续SEACells系列文章会提到）

数据集介绍
这是过滤后的由 CD34+ 分选出来的的骨髓细胞数据集，未进行标准化，用于分析人类造血功能，详细信息可根据上述下载链接进行查看。

# 用scanpy加载数据集
rna_ad = sc.read('data/cd34_multiome_rna.h5ad')
atac_ad = sc.read('data/cd34_multiome_atac.h5ad')

# ATAC的参考细胞类型
sc.pl.scatter(atac_ad, basis='umap', color='celltype', frameon=False)

# SEACells的metacells计算结果
SEACells.plot.plot_2D(atac_ad, key='X_umap', colour_metacells=False)

三、提取元细胞

ATAC 和 RNA 单细胞 Anndata 对象应包含同一组细胞。仅相同细胞将用于下游分析。

atac_meta_ad, rna_meta_ad = SEACells.genescores.prepare_multiome_anndata(atac_ad, rna_ad, SEACells_label='SEACell')

atac_meta_ad,rna_meta_ad

AnnData object with n_obs × n_vars = 99 × 246113
    obs: 'n_counts'
    var: 'GC_bin', 'counts_bin', 'n_cells'
    uns: 'log1p'
    obsm: 'X_svd'
    layers: 'raw'
AnnData object with n_obs × n_vars = 99 × 12464
    obs: 'n_counts'
    uns: 'log1p'
    layers: 'raw'

四、Gene-peak correlations计算

要计算基因表达和基因附近峰的可及性的相关性，需要：
1、带有基因注释的 GTF 文件及基因周围的基因组跨度信息以测试相关性
2、染色体名称应编号为 1：、2：形式

#基因组文件获取
! wget https://dp-lab-data-public.s3.amazonaws.com/SEACells-multiome/hg38.gtf -O data/hg38.gtf

# 在此示例中，我们计算自此以来的前十个基因的基因峰值相关性，方便进行演示
# 这个过程是计算密集型的
gene_set = rna_meta_ad.var_names[:10]
gene_peak_cors = SEACells.genescores.get_gene_peak_correlations(atac_meta_ad, rna_meta_ad, 
                                           path_to_gtf='data/hg38.gtf', 
                                           span=100000, 
                                           n_jobs=1,
                                           gene_set=gene_set)

#查看运行结果
gene_peak_cors['FAM41C'].head()

五、高度调控基因的计算

对于下游分析，完整的基因峰值相关性结果可作为 pickle 文件提供，网址为 https://dp-lab-data-public.s3.amazonaws.com/SEACells-multiome/cd34_multiome_gene_peak_cors.p

#基因组文件获取
! wget https://dp-lab-data-public.s3.amazonaws.com/SEACells-multiome/cd34_multiome_gene_peak_cors.p -O data/cd34_multiome_gene_peak_cors.p
gene_peak_cors = pd.read_pickle('data/cd34_multiome_gene_peak_cors.p') 

# 高度调控的基因，即与多个峰相关的基因，可以使用 get_gene_peak_assocations 函数来识别。get_gene_peak_assocations 返回与每个基因相关的显着峰的数量
peak_counts = SEACells.genescores.get_gene_peak_assocations(gene_peak_cors, 
                                                           pval_cutoff=1e-1,
                                                           cor_cutoff=0.1)

参数说明：
pval_cutoff=1e-1: 这个参数设置了一个p值截止值。相关性的p值大于这个截止值（在这种情况下是0.1）的基因峰值相关性可能被过滤掉。这意味着只有p值小于或等于0.1的相关性被认为是显著的。
cor_cutoff=0.1: 这个参数设置了一个相关系数截止值。相关性系数小于这个截止值（在这种情况下是0.1）的基因峰值相关性可能被过滤掉。这意味着只有相关系数大于或等于0.1的相关性被认为是显著的。

#绘制分布图以识别具有较高调控程度的基因
plt.scatter(np.arange(len(peak_counts)), 
           np.sort(peak_counts), s=20)
sns.despine()
plt.xlabel('Gene rank')
plt.ylabel('No. of correlated peaks')

六、基因得分

基因分数计算为相关峰的可及性的加权和，并且可以使用 get_gene_scores 计算。

gene_scores = SEACells.genescores.get_gene_scores(atac_meta_ad, 
                                                  gene_peak_cors,
                                                  pval_cutoff=1e-1,
                                                  cor_cutoff=0.1)

gene_scores 是一个 pandas.Data Frame，其中元细胞为行，基因为列。
这可用于任何下游分析，例如聚类、可视化等。

七、Gene-accessibility计算

元细胞中的开放峰识别

基因分数计算为相关峰的可及性的加权和，并且可以使用 get_gene_scores 计算。

该函数还需要低维嵌入，例如 X_svd、X_umap等等。
此处我们应用scATAC 的 SVD 来进行下面的分析

# 用原始计数矩阵创建一个元细胞的anndata数据
atac_meta_ad = SEACells.core.summarize_by_SEACell(atac_ad, SEACells_label='SEACell')
atac_meta_ad.obs['n_counts'] = np.ravel(atac_meta_ad.X.sum(axis=1))

# 加入svd低维嵌入
seacell_label = 'SEACell'
sc_svd = pd.DataFrame(atac_ad.obsm['X_svd'], index=atac_ad.obs_names)
atac_meta_ad.obsm['X_svd'] = sc_svd.groupby(atac_ad.obs[seacell_label]).mean().loc[atac_meta_ad.obs_names, :]

# 确定每个元细胞中的开放峰
SEACells.accessibility.determine_metacell_open_peaks(atac_meta_ad, peak_set=None, low_dim_embedding='X_svd', pval_cutoff=1e-2,
                                  read_len=147, n_neighbors=3, n_jobs=1)
# 该函数会将“Open Peaks”添加到 Anndata 层，并且是一个二进制矩阵，指示元细胞中峰是开放的还是封闭的

基因可及性指标

开放峰用于计算基因可及性度量，该度量表示相关开放峰的分数。
注意：只有当每个基因都有足够数量的开放峰时，此指标才是可靠的。
建议仅用于高度调控的基因，即我们第五部分计算得到的结果。

# 根据肘点选择基因，此处的参考图为高度调控基因的分布图
high_reg_genes = peak_counts.index[peak_counts > 9]

# 计算基因可及性
SEACells.accessibility.get_gene_accessibility(atac_meta_ad, gene_peak_cors, 
                                              gene_set=high_reg_genes, pval_cutoff=1e-1, cor_cutoff=0.1)
# 此函数会将“基因可访问性”添加到 Anndata `.obsm` 字段

可视化

# 首先加入umap以可视化基因可及性
rna_umap = pd.DataFrame(rna_ad.obsm['X_umap'], index=rna_ad.obs_names)
rna_meta_ad.obsm['X_umap'] = rna_umap.groupby(atac_ad.obs[seacell_label]).mean().loc[rna_meta_ad.obs_names, :].values

genes = ['KLF1', 'GATA1', 'SPI1']
# Copy accessibility to RNA meta anndata
temp = rna_meta_ad[:, genes]
temp.layers['GeneAccessibility'] = atac_meta_ad[rna_meta_ad.obs_names].obsm['GeneAccessibility'][genes].values

# Plot expression
sc.pl.scatter(rna_meta_ad, basis='umap', color=genes)

# Plot accessibility
sc.pl.scatter(temp, basis='umap', color=genes, layers='GeneAccessibility')

Reference

[1]https://github.com/dpeerlab/SEACells
[2]Persad, S., Choo, ZN., Dien, C. et al. SEACells infers transcriptional and epigenomic cellular states from single-cell genomics data. Nat Biotechnol (2023).