pycisTopic | 新手友好的scATAC分析文档

scATAC数据的分析软件也有不少，Signac、snapATAC2、ArchR等都是比较常见的软件，其中ArchR的功能最为全面。那pycisTopic有什么优势呢？好处是分析生态对接，使用pycisTopic分析scATAC数据生成的结果可以直接作为SCENIC+的输入，后续可以分析增强子驱动的转录因子调控网络。

pycisTopic的软件文档也是相当的详细，分析过程的一些注意事项都写的很清楚，即使是新手跟着流程来一遍也可以完成分析。文档使用的数据是human cerebellum单细胞ATAC和配套的单细胞RNA多组学数据。

import pycisTopic
import warnings
warnings.filterwarnings('ignore')

下载数据

本教程使用的数据可自由获取，可从https://www.10xgenomics.com/resources/datasets/frozen-human-healthy-brain-tissue-3-k-1-standard-1-0-0下载，处理后的loom文件可在https://scope.aertslab.org/#/scenic-v2上找到，数据也可在UCSC基因组浏览器(https://genome-euro.ucsc.edu/s/cbravo/SCENIC%2B_cerebellum)上探索。

!wget -O data/human_brain_3k_atac/fragments.tsv.gz https://cf.10xgenomics.com/samples/cell-arc/1.0.0/human_brain_3k/human_brain_3k_atac_fragments.tsv.gz
!wget -O data/human_brain_3k_atac/fragments.tsv.gz.tbi https://cf.10xgenomics.com/samples/cell-arc/1.0.0/human_brain_3k/human_brain_3k_atac_fragments.tsv.gz.tbi
!wget -O data/human_brain_3k_atac/cell_data.tsv https://raw.githubusercontent.com/aertslab/pycisTopic/polars/data/cell_data_human_cerebellum.tsv

环境准备

创建一些目录用于存储pycisTopic的输出结果。

import os

out_dir = "data"
os.makedirs(out_dir, exist_ok = True)

定义一个字典，将样本ID ("10x_multiome_brain") 映射到片段文件 ("fragments.tsv.gz")。

如果有多个样本，可以向该字典中添加多个条目。

注意： pycisTopic 会自动在细胞条形码后附加样本ID，以避免样本间条形码冲突！

fragments_dict = {"10x_multiome_brain": "data/human_brain_3k_atac/fragments.tsv.gz"}

细胞注释

在本教程中，我们假设分析的是来自多组学数据集的scATAC-seq数据，这样可以方便地从scRNA-seq分析中获得细胞注释。对于独立的scATAC-seq数据，细胞注释也可以通过其他方法获得，例如无标注/初步聚类分析（使用预定义的区域，例如针对小鼠和人类的SCREEN）。在后一种情况下，您可以跳过本节，直接将整体区域作为QC步骤的输入。

首先，将条形码到细胞类型的注释读取为 pd.DataFrame。

数据框的索引应对应于 fragments.tsv.gz 文件中的条形码列（可以使用 split_pattern（本例中为 "-"）选择性地将样本ID附加到条形码），并且必须包含一个包含样本ID的列。这些样本ID应与上面定义的片段字典中的样本ID匹配。

或者，可以在元数据中添加一个名为 "barcode" 的列，其中包含细胞条形码（不带任何后缀！）。在这种情况下，将不会使用数据框的索引。

import pandas as pd
cell_data = pd.read_table("data/human_brain_3k_atac/cell_data.tsv", index_col = 0)

还需要每条染色体的大小，可以从UCSC数据库下载。UCSC的下载地址无法访问，即使翻墙也无济于事，自行准备吧。

chromsizes = pd.read_table("http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes", header=None, names=["Chromosome", "End"])
chromsizes.insert(1, "Start", 0)

现在可以生成每个细胞类型的伪批量ATAC-seq图谱。

该函数将为variable参数定义的每种细胞类型生成一个fragments.tsv.gz和一个 bigwig 文件。

对于每种细胞类型：

• fragments.tsv.gz包含该细胞类型所有条形码的片段。
• 每个fragments.tsv.gz文件会生成一个bigwig文件，您可以在任何基因组浏览器中查看。

from pycisTopic.pseudobulk_peak_calling import export_pseudobulk
os.makedirs(os.path.join(out_dir, "consensus_peak_calling"), exist_ok = True)
os.makedirs(os.path.join(out_dir, "consensus_peak_calling/pseudobulk_bed_files"), exist_ok = True)
os.makedirs(os.path.join(out_dir, "consensus_peak_calling/pseudobulk_bw_files"), exist_ok = True)

bw_paths, bed_paths = export_pseudobulk(
    input_data = cell_data,
    variable = "VSN_cell_type",
    sample_id_col = "VSN_sample_id",
    chromsizes = chromsizes,
    bed_path = os.path.join(out_dir, "consensus_peak_calling/pseudobulk_bed_files"),
    bigwig_path = os.path.join(out_dir, "consensus_peak_calling/pseudobulk_bw_files"),
    path_to_fragments = fragments_dict,
    n_cpu = 10,
    normalize_bigwig = True,
    temp_dir = "/tmp/ray_spill",
    split_pattern = "-")

稍后需要这些bed文件的路径，所以现在将它们保存到磁盘。

with open(os.path.join(out_dir, "consensus_peak_calling/bw_paths.tsv"), "wt") as f:
    for v in bw_paths:
        _ = f.write(f"{v}\t{bw_paths[v]}\n")
        
with open(os.path.join(out_dir, "consensus_peak_calling/bed_paths.tsv"), "wt") as f:
    for v in bed_paths:
        _ = f.write(f"{v}\t{bed_paths[v]}\n")

推断共有峰

接下来，我们将使用MACS对每个伪批量fragments.tsv.gz文件进行峰calling。

from pycisTopic.pseudobulk_peak_calling import peak_calling

os.makedirs(os.path.join(out_dir, "consensus_peak_calling/MACS"), exist_ok=True)

narrow_peak_dict = peak_calling(
    macs_path = 'macs2',
    bed_paths = bed_paths,
    outdir = os.path.join(os.path.join(out_dir, "consensus_peak_calling/MACS")),
    genome_size = 'hs',
    n_cpu = 10,
    input_format = 'BEDPE',
    shift = 73,
    ext_size = 146,
    keep_dup = 'all',
    q_value = 0.05,
    _temp_dir = '/tmp/ray_spill'
)

最后，是时候推导共有峰了。为此，我们使用 TGCA 迭代峰过滤方法。
首先，每个峰中心 (summit) 向两侧各扩展peak_half_width，然后我们迭代地过滤掉那些与更显著的峰重叠的较不显著峰。在此过程中，峰会被合并，根据合并中包含的峰数量，将发生不同的处理：

• 1个峰：保留原始峰区域
• 2个峰：保留得分较高的原始峰区域
• 3个或更多峰：取得分最显著的原始峰区域，并移除与该显著峰区域重叠的所有原始峰区域。然后对下一个最显著的峰 (如果未被移除) 重复该过程，直到所有峰都被处理。

该过程将进行两次，首先对每个伪批量峰分别进行；然后在峰得分归一化后，对所有峰一起进行处理。

from pycisTopic.iterative_peak_calling import get_consensus_peaks

consensus_peaks = get_consensus_peaks(
    narrow_peaks_dict = narrow_peak_dict,
    peak_half_width = 250,
    chromsizes = chromsizes,
    path_to_blacklist = 'hg38.blacklist.v2.encode.bed')

consensus_peaks.to_bed(
    path = os.path.join(out_dir, "consensus_peak_calling/consensus_regions.bed"),
    keep =True,
    compression = 'infer',
    chain = False)

质量控制

下一步是对scATAC-seq样本 (本例中只有一个样本) 进行质量控制。此步骤包含多项指标和可视化：

• 条形码秩图（Barcode rank plot）
• 重复率（Duplication rate）
• 插入片段大小（Insertion size）
• TSS富集（TSS enrichment）
• 峰内读段比例（FRIP）

要计算TSS富集，我们需要提供TSS注释。您可以通过 pycistopic tss get_tss 命令轻松下载。

如果不确定Ensembl数据库中用于指定基因名称的列名，请运行pycistopic tss gene_annotation_list并搜索您的物种。

!pycistopic tss gene_annotation_list | grep Human

确定好了，然后下载：

!mkdir -p data/qc
!pycistopic tss get_tss \
    --output data/qc/tss.bed \
    --name "hsapiens_gene_ensembl" \
    --to-chrom-source ucsc \
    --ucsc hg38

接下来，使用pycistopic qc命令计算QC指标。

!pycistopic qc \
    --fragments data/human_brain_3k_atac/fragments.tsv.gz \
    --regions data/consensus_peak_calling/consensus_regions.bed \
    --tss data/qc/tss.bed \
    --output data/qc/10x_multiome_brain

如果有多个样本，可以按如下方式并行运行QC步骤。

regions_bed_filename = os.path.join(out_dir, "consensus_peak_calling/consensus_regions.bed")
tss_bed_filename = os.path.join(out_dir, "qc", "tss.bed")

pycistopic_qc_commands_filename = "pycistopic_qc_commands.txt"

# 创建包含所有 pycistopic qc 命令行的文本文件。
with open(pycistopic_qc_commands_filename, "w") as fh:
    for sample, fragment_filename in fragments_dict.items():
        print("pycistopic qc",
            f"--fragments {fragment_filename}",
            f"--regions {regions_bed_filename}",
            f"--tss {tss_bed_filename}",
            f"--output {os.path.join(out_dir, 'qc')}/{sample}",
            sep=" ", file=fh)

然后在命令行环境中运行以下命令：

cat pycistopic_qc_commands.txt | parallel -j 4 {}

最后，可以可视化样本级别的统计信息。包括：

• 条形码秩图：显示非重复读段的分布以及哪些条形码被推断为与细胞相关。陡峭的下降（“knee”）表明细胞相关条形码与空分区相关条形码之间有良好的分离。
• 插入片段大小：ATAC-seq需要在DNA末端有Tn5转座酶的正确配对切割事件。在无核小体的开放染色质区域，许多Tn5分子可以作用并将DNA切割成小片段；在核小体占据的区域，Tn5只能进入连接区域。因此，在良好的ATAC-seq文库中，您应该期望在<100 bp区域（开放染色质）看到尖锐的峰，在约200 bp区域 (单核小体) 看到峰，以及其他更大的峰（多核小体）。清晰的核小体模式表明实验质量良好。
• 样本TSS富集：TSS富集计算是一种信噪比计算。收集参考TSS集合周围的读段，形成以TSS为中心、向两侧各延伸1000 bp (总共2000 bp) 的汇总分布。然后通过取分布两端各100 bp (总共200 bp) 的平均读段深度进行归一化，并计算每个位置相对于该平均读段深度的倍数变化。这意味着两侧翼应从1开始，如果在转录起始位点有高读段信号 (基因组高度开放区域)，则信号应增加直至中间达到峰值。

from pycisTopic.plotting.qc_plot import plot_sample_stats, plot_barcode_stats
import matplotlib.pyplot as plt

for sample_id in fragments_dict:
    fig = plot_sample_stats(sample_id=sample_id, pycistopic_qc_output_dir="data/qc")

可以可视化条形码水平的统计信息。这些统计信息可用于过滤细胞条形码，仅保留高质量细胞。

注意：pycistopic qc命令将自动确定最小唯一片段数和最小TSS富集的阈值。如果想更改这些阈值或希望基于FRIP进行阈值筛选，可以使用以下参数手动定义阈值：

• unique_fragments_threshold
• tss_enrichment_threshold
• frip_threshold

在本例中，将使用自动确定的阈值，请务必手动检查质量指标以确保这些阈值有效！

条形码水平统计信息包括：

• (唯一)片段总数
• TSS富集：每个条形码在TSS富集得分中位置0 (即TSS) 的得分。噪声细胞将具有较低的TSS富集。
• FRIP：每个条形码的峰内读段比例。噪声细胞的FRIP值较低。然而，该过滤应谨慎使用，因为它取决于原始峰的质量。例如，如果样本中存在罕见群体，其特异性峰可能被峰 calling 算法遗漏，导致其FRIP值降低。

from pycisTopic.qc import get_barcodes_passing_qc_for_sample

sample_id_to_barcodes_passing_filters = {}
sample_id_to_thresholds = {}
for sample_id in fragments_dict:
    (sample_id_to_barcodes_passing_filters[sample_id],
     sample_id_to_thresholds[sample_id]) = get_barcodes_passing_qc_for_sample(
            sample_id = sample_id,
            pycistopic_qc_output_dir = "data/qc",
            unique_fragments_threshold = None, 
            tss_enrichment_threshold = None, 
            frip_threshold = 0,
            use_automatic_thresholds = True)

for sample_id in fragments_dict:
    fig = plot_barcode_stats(sample_id = sample_id,
        pycistopic_qc_output_dir = "data/qc",
        bc_passing_filters = sample_id_to_barcodes_passing_filters[sample_id],
        detailed_title = False,
        **sample_id_to_thresholds[sample_id])

创建对象

在此步骤中，将创建一个cisTopic对象。这包括生成一个计数矩阵，其中包含在共有峰 (见上文) 上的片段计数，针对通过上述QC指标的每个细胞条形码。

Blacklist区域 (https://www.nature.com/articles/s41598-019-45839-z) 将从该计数矩阵中移除。

path_to_regions = os.path.join(out_dir, "consensus_peak_calling/consensus_regions.bed")
path_to_blacklist = "hg38.blacklist.v2.encode.bed"
pycistopic_qc_output_dir = "data/qc"

from pycisTopic.cistopic_class import create_cistopic_object_from_fragments
import polars as pl

cistopic_obj_list = []
for sample_id in fragments_dict:
    sample_metrics = pl.read_parquet(os.path.join(pycistopic_qc_output_dir, f'{sample_id}.fragments_stats_per_cb.parquet')
    ).to_pandas().set_index("CB").loc[ sample_id_to_barcodes_passing_filters[sample_id] ]
    cistopic_obj = create_cistopic_object_from_fragments(
        path_to_fragments = fragments_dict[sample_id],
        path_to_regions = path_to_regions,
        path_to_blacklist = path_to_blacklist,
        metrics = sample_metrics,
        valid_bc = sample_id_to_barcodes_passing_filters[sample_id],
        n_cpu = 1,
        project = sample_id,
        split_pattern = '-')
    cistopic_obj_list.append(cistopic_obj)

本例中只有一个样本，因此只生成了一个cisTopic对象。如果您有多个样本，则需要使用 merge() 函数合并您的cisTopic对象列表。

import pickle

cistopic_obj = cistopic_obj_list[0]
pickle.dump(cistopic_obj, open(os.path.join(out_dir, "cistopic_obj.pkl"), "wb"))

添加元数据

可以向cisTopic对象添加额外的元数据 (针对区域或细胞)。例如，我们来添加scRNA-seq 数据的注释。缺失值将填充为Nan。

import pandas as pd
cell_data = pd.read_table("data/cell_data.tsv", index_col = 0)
cistopic_obj.add_cell_data(cell_data, split_pattern='-')
pickle.dump(cistopic_obj, open(os.path.join(out_dir, "cistopic_obj.pkl"), "wb"))

双胞检测

还可以选择在片段计数矩阵上运行scrublet，以从scATAC-seq数据中推断双细胞 (doublets)。

import scrublet as scr
scrub = scr.Scrublet(cistopic_obj.fragment_matrix.T, expected_doublet_rate=0.1)
doublet_scores, predicted_doublets = scrub.scrub_doublets()
scrub.plot_histogram()
scrub.call_doublets(threshold=0.22)
scrub.plot_histogram()
scrublet = pd.DataFrame([scrub.doublet_scores_obs_, scrub.predicted_doublets_], columns=cistopic_obj.cell_names, index=['Doublet_scores_fragments', 'Predicted_doublets_fragments']).T

cistopic_obj.add_cell_data(scrublet, split_pattern = '-')
sum(cistopic_obj.cell_data.Predicted_doublets_fragments == True)

pickle.dump(cistopic_obj, open(os.path.join(out_dir, "cistopic_obj.pkl"), "wb"))

singlets = cistopic_obj.cell_data[cistopic_obj.cell_data.Predicted_doublets_fragments == False].index.tolist()
# 子集 cisTopic 对象
cistopic_obj_raw = cistopic_obj
cistopic_obj = cistopic_obj.subset(singlets, copy=True, split_pattern='-')
pickle.dump(cistopic_obj, open(os.path.join(out_dir, "cistopic_obj_nodbl.pkl"), "wb"))

topic建模

接下来，我们将使用折叠吉布斯采样器 (Collapsed Gibbs Sampler) 执行实际的LDA topic建模。有两个函数可用于执行topic建模，两者产生相似的结果。

• 串行LDA：并行化是在topic之间进行，而不是在每个模型内部。适用于中小型数据集，其中需要测试多个不同topic数的模型。您可以使用runCGSModels()运行这些模型。
• 使用MALLET的并行LDA：并行化在每个模型内部进行。适用于大型数据集，其中只需测试少数几个不同topic数的模型。如果在集群中工作，我们建议为每个模型提交一个作业，以便它们可以同时运行。您可以使用runCGSModelsMallet()运行它。

在这里，使用的是Mallet。

注意： 为了能够运行Mallet，需要Mallet二进制文件，可以从https://github.com/mimno/Mallet/releases下载。也可以从源代码编译二进制文件，更多信息请访问Mallet Github：https://github.com/mimno/Mallet。

!wget https://github.com/mimno/Mallet/releases/download/v202108/Mallet-202108-bin.tar.gz
!tar -xf Mallet-202108-bin.tar.gz

!mkdir -p mkdir -p /tmp/ray_spill/mallet/tutorial

因为不知道对于我们的数据集最佳topic数是多少，可以运行多个topic模型，每个具有不同的topic数。

os.environ['MALLET_MEMORY'] = '100G'
from pycisTopic.lda_models import run_cgs_models_mallet

models = run_cgs_models_mallet(
    cistopic_obj,
    n_topics=[2, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50],
    n_cpu=12,
    n_iter=500,
    alpha=50,
    alpha_by_topic=True,
    eta=0.1,
    eta_by_topic=False,
    tmp_path="/tmp/ray_spill/mallet/tutorial",
    save_path="/tmp/ray_spill/mallet/tutorial",
    mallet_path='Mallet-202108/bin/mallet')
    
pickle.dump(models, open(os.path.join(out_dir, "models.pkl"), "wb"))

模型选择

选择具有最佳topic数的模型，没有最优的方法来进行此选择，然而获得精确的最优topic数也并非关键。为了简化选择，实现了几项指标：

• Minmo_2011：使用Mimno等人 (2011) 计算的平均模型一致性 (coherence)。为了减少topic数的影响，我们基于前n个选定平均值计算平均一致性。模型越好，一致性越高。
• Log-likelihood：使用Griffiths和Steyvers (2004) 最后迭代的对数似然。模型越好，对数似然越高。
• Arun_2010：使用Arun等人 (2010) 基于topic-region分布、cell-topic分布和细胞覆盖度的密度指标。模型越好，指标越低。
• Cao_Juan_2009：使用Cao Juan等人 (2009) 基于主题-区域分布的散度指标。模型越好，指标越低。

注意： 请注意，对于 Arun 和 Cao 指标，较低的分数对应较好的模型，为便于可视化，我们反转了这些分数。
在下图中，这些指标的较高分数因此对应更好的模型。

另请注意，并非所有指标都一致 (例如参见Arun指标)。

对于scATAC-seq数据模型，最有帮助的方法是Minmo和log-likelihood。

from pycisTopic.lda_models import evaluate_models

model = evaluate_models(models, select_model=40, return_model=True)

cistopic_obj.add_LDA_model(model)
pickle.dump(cistopic_obj, open(os.path.join(out_dir, "cistopic_obj.pkl"), "wb"))

聚类与可视化

可以使用Leiden算法对细胞 (或区域) 进行聚类，并使用UMAP和TSNE进行降维。在这些示例中，将重点关注细胞。在这些步骤中，将使用模型的细胞-主题贡献。

from pycisTopic.clust_vis import find_clusters, run_umap, run_tsne, plot_metadata, plot_topic, cell_topic_heatmap

find_clusters(cistopic_obj, target='cell', k=10, res=[0.6, 1.2, 3], prefix='pycisTopic_', scale=True, split_pattern='-')
run_umap(cistopic_obj, target='cell', scale=True)
run_tsne(cistopic_obj, target='cell', scale=True)

plot_metadata(
    cistopic_obj,
    reduction_name='UMAP',
    variables=['Seurat_cell_type', 'pycisTopic_leiden_10_0.6', 'pycisTopic_leiden_10_1.2', 'pycisTopic_leiden_10_3'], 
    target='cell', num_columns=4,
    text_size=10,
    dot_size=5)

根据与scRNA-seq注释的重叠来注释每个聚类。

for resolution in [0.6, 1.2, 3]:
    col = f'pycisTopic_leiden_10_{resolution}'
    # 计算各簇的众数细胞类型
    mode_series = cistopic_obj.cell_data.groupby(col)["Seurat_cell_type"].agg(lambda x: x.value_counts().idxmax())
    # 生成注释列
    new_col = f'{col}_anno'
    cistopic_obj.cell_data[new_col] = cistopic_obj.cell_data[col].map(lambda cluster: f"{mode_series[cluster]}({cluster})")

plot_metadata(
    cistopic_obj,
    reduction_name='UMAP',
    variables=['Seurat_cell_type', 'pycisTopic_leiden_10_0.6_anno', 'pycisTopic_leiden_10_1.2_anno', 'pycisTopic_leiden_10_3_anno'], 
    target='cell', num_columns=4,
    text_size=10,
    dot_size=5)

也可以绘制连续值：

plot_metadata(
    cistopic_obj,
    reduction_name='UMAP',
    variables=['log10_unique_fragments_count', 'tss_enrichment', 'Doublet_scores_fragments', 'fraction_of_fragments_in_peaks'],
    target='cell', num_columns=4,
    text_size=10,
    dot_size=5)

可视化cell-topic贡献：

plot_topic(cistopic_obj, reduction_name='UMAP', target='cell', num_columns=5)

或者也可以绘制topic贡献的热图：

cell_topic_heatmap(cistopic_obj, variables=['Seurat_cell_type'], scale=False, legend_loc_x=1, legend_loc_y=-1.3, figsize=(8,8))

topic二值化

接下来，可以对topic-region和cell-topic分布进行二值化。前者有助于使用其他处理区域集的工具 (例如GREAT、cisTarget) 探索主题；而后者有助于自动注释主题。

我们将首先对topic-region分布进行二值化。有几种方法可用于此：‘otsu’(Otsu, 1979)、‘yen’ (Yen et al., 1995)、‘li’ (Li & Lee, 1993)、‘aucell’ (Van de Sande et al.，2020) 或‘ntop’ (取每个主题的前n个区域)。Otsu 和 Yen 的方法对topic-region分布效果良好；但对于某些下游分析 (例如深度学习)，使用ntop以获得平衡的区域集可能更方便。

from pycisTopic.topic_binarization import binarize_topics

region_bin_topics_top_3k = binarize_topics(cistopic_obj, method='ntop', ntop=3000, plot=True, num_columns=5)

类似地，现在可以对cell-topic分布进行二值化。

binarized_cell_topic = binarize_topics(cistopic_obj, target='cell', method='li', plot=True, num_columns=5, nbins=100)

接下来，可以计算topic质量控制指标。包括：

• 分配数 (Number of assignments)
• topic一致性 (Mimno et al., 2011)：衡量主题中高评分区域在原始数据中实际共可及性的程度。如果较低，则表明topic较为随机。越高，主题越好。
• 边缘topic分布 (Marginal topic distribution)：指示每个topic对模型的贡献程度。越高，主题越好。
• 基尼指数 (Gini index)：介于0和1之间的值，指示主题的特异性 (0：一般，1：特异)
• 如果topic已被二值化，将添加每个topic的区域/细胞数量。

from pycisTopic.topic_qc import compute_topic_metrics, plot_topic_qc, topic_annotation
import matplotlib.pyplot as plt
from pycisTopic.utils import fig2img

topic_qc_metrics = compute_topic_metrics(cistopic_obj)

fig_dict={}
fig_dict['CoherenceVSAssignments'] = plot_topic_qc(topic_qc_metrics, var_x='Coherence', var_y='Log10_Assignments', var_color='Gini_index', plot=False, return_fig=True)
fig_dict['AssignmentsVSCells_in_bin'] = plot_topic_qc(topic_qc_metrics, var_x='Log10_Assignments', var_y='Cells_in_binarized_topic', var_color='Gini_index', plot=False, return_fig=True)
fig_dict['CoherenceVSCells_in_bin'] = plot_topic_qc(topic_qc_metrics, var_x='Coherence', var_y='Cells_in_binarized_topic', var_color='Gini_index', plot=False, return_fig=True)
fig_dict['CoherenceVSRegions_in_bin'] = plot_topic_qc(topic_qc_metrics, var_x='Coherence', var_y='Regions_in_binarized_topic', var_color='Gini_index', plot=False, return_fig=True)
fig_dict['CoherenceVSMarginal_dist'] = plot_topic_qc(topic_qc_metrics, var_x='Coherence', var_y='Marginal_topic_dist', var_color='Gini_index', plot=False, return_fig=True)
fig_dict['CoherenceVSGini_index'] = plot_topic_qc(topic_qc_metrics, var_x='Coherence', var_y='Gini_index', var_color='Gini_index', plot=False, return_fig=True)

fig=plt.figure(figsize=(40, 43))
i = 1
for fig_ in fig_dict.keys():
    plt.subplot(2, 3, i)
    img = fig2img(fig_dict[fig_])
    plt.imshow(img)
    plt.axis('off')
    i += 1
plt.subplots_adjust(wspace=0, hspace=-0.70)
plt.show()

接下来，可以自动注释topic，本例中按细胞类型注释。这里计算每个组中分配到二值化topic的细胞比例，与整个数据集中的比例进行比较。如果整个数据集中二值化topic的细胞比例与分配组中总细胞比例的差值大于0.2，这表明该topic是普遍的，如果topic在前景 (该组) 和背景 (整个数据集) 中都富集，比例检验可能失效：导致比例差异较大。

topic_annot = topic_annotation(cistopic_obj, annot_var='Seurat_cell_type', binarized_cell_topic=binarized_cell_topic, general_topic_thr=0.2)

差异可及区域

除了使用调控主题外，还可以识别细胞类型之间的差异可及区域 (DARs)。首先，将利用细胞-主题和主题-区域概率来推算区域可及性。为了将非常低的概率值收缩为0，使用一个缩放因子 (默认：10^6)。

from pycisTopic.diff_features import impute_accessibility, normalize_scores, find_highly_variable_features, find_diff_features
import numpy as np

imputed_acc_obj = impute_accessibility(cistopic_obj, scale_factor=10**6)
normalized_imputed_acc_obj = normalize_scores(imputed_acc_obj, scale_factor=10**4)

可选步骤，可以识别高度可变区域。这不是必须的，但会加速识别DARs的假设检验步骤。

variable_regions = find_highly_variable_features(normalized_imputed_acc_obj, plot=True)

现在可以识别组间的差异可及区域。默认情况下，该函数将使用指定的变量对每个组与其余组进行Wilcoxon秩和检验。或者，可以以列表形式提供指定的对比，包含前景组和背景组 (例如，组1与组2和组3，以及组2与组1和组3：[[[‘Group_1’], [‘Group_2, ‘Group_3’]], [[‘Group_2’], [‘Group_1, ‘Group_3’]]])。

markers_dict = find_diff_features(cistopic_obj, imputed_acc_obj, variable='Seurat_cell_type', var_features=variable_regions, log2fc_thr=np.log2(1.5), n_cpu=5,
    _temp_dir='/tmp/ray_spill', split_pattern='-')

还可以将区域可及性绘制到细胞-主题UMAP上。例如，查看一些细胞类型的最佳DARs。

from pycisTopic.clust_vis import plot_imputed_features

plot_imputed_features(
    cistopic_obj,
    reduction_name='UMAP',
    imputed_data=imputed_acc_obj,
    features=[markers_dict[x].index.tolist()[0] for x in ['BG', 'GC', 'INH_SST', 'COP']],
    scale=False,
    num_columns=4)

保存区域集

os.makedirs(os.path.join(out_dir, "region_sets"), exist_ok = True)
os.makedirs(os.path.join(out_dir, "region_sets", "topics_top_3k"), exist_ok = True)
os.makedirs(os.path.join(out_dir, "region_sets", "dars_cell_type"), exist_ok = True)

from pycisTopic.utils import region_names_to_coordinates

for topic in region_bin_topics_top_3k:
    region_names_to_coordinates(region_bin_topics_top_3k[topic].index).sort_values(["Chromosome", "Start", "End"]
    ).to_csv(os.path.join(out_dir, "region_sets", "topics_top_3k", f"{topic}.bed"), sep="\t", header=False, index=False)

for cell_type in markers_dict:
    region_names_to_coordinates(markers_dict[cell_type].index).sort_values(["Chromosome", "Start", "End"]
    ).to_csv(os.path.join(out_dir, "region_sets", "dars_cell_type", f"{cell_type}.bed"), sep="\t", header=False, index=False)

基因活性

推断基因活性，此函数中有几个选项可供评估：

• 搜索空间：用户可以选择搜索空间是否应包含其他基因 (use_gene_boundaries)，以及最小和最大距离 (上游和下游)
• 距离权重：与距离权重相关的参数以指数函数形式衡量距离对推断区域到基因权重的影响。用户可以控制是否使用此权重 (distance_weight) 以及距离的影响 (decay_rate)。
• 基因大小权重：大基因可能偶然有更多峰。用户可以选择基于每个基因的大小应用权重 (gene_size_weight)，默认是将每个基因的大小除以基因组中基因大小的中位数。或者，用户也可以使用average_scores，这将计算基因活性作为该基因所有关联区域加权区域可及性的平均值。
• 基尼权重：此权重将给予更特异的区域更多重要性 (gini_weight)。

import pyranges as pr
from pycisTopic.gene_activity import get_gene_activity

chromsizes = pr.PyRanges(chromsizes[["Chromosome", "Start", "End"]])
pr_annotation = pd.read_table(os.path.join(out_dir, "qc", "tss.bed")).rename({"Name": "Gene", "# Chromosome": "Chromosome"}, axis=1)
pr_annotation["Transcription_Start_Site"] = pr_annotation["Start"]
pr_annotation = pr.PyRanges(pr_annotation)

gene_act, weigths = get_gene_activity(
    imputed_acc_obj,
    pr_annotation,
    chromsizes,
    use_gene_boundaries=True,
    upstream=[1000, 100000],
    downstream=[1000,100000],
    distance_weight=True,
    decay_rate=1,
    extend_gene_body_upstream=10000,
    extend_gene_body_downstream=500,
    gene_size_weight=False,
    gene_size_scale_factor='median',
    remove_promoters=False,
    average_scores=True,
    scale_factor=1,
    extend_tss=[10,10],
    gini_weight = True,
    return_weights= True,
    project='Gene_activity')

正如之前对推算区域可及性所做的那样，现在也可以推断差异可及基因 (DAGs)。

dag_markers_dict= find_diff_features(cistopic_obj, gene_act, variable='Seurat_cell_type', log2fc_thr=np.log2(1.5), n_cpu=5,
    _temp_dir='/tmp/ray_spill', split_pattern='-')

plot_imputed_features(cistopic_obj, reduction_name='UMAP', imputed_data=gene_act,
    features=['PDGFRA', 'OLIG2', 'MBP', 'SOX10', # Olig differentiation
              'CTNNA3', 'ENPP6', 'OLIG1', # Olig differentiation
              'GAD2', 'VIP', 'SST', 'CTXN3', # Int
              'NFIB', 'SOX9',  # Ast
              'LEF1', # Endo
              'SPI1'], # Glia
    scale=True, num_columns=4)

标签转移

利用基因活性得分，可以从参考数据集 (如scRNA-seq) 转移标签。作为示例，我们将从该数据集的scRNA-seq层转移标签。可用于标签转移的方法有：

• ingest (scanpy)
• harmony (Korsunsky et al, 2019)
• bbknn (Polański et al, 2020)
• scanorama(Hie et al, 2019)
• cca

除ingest外，这些方法返回共同的联合嵌入，并使用查询细胞与参考细胞之间的距离作为权重来推断标签。

from pycisTopic.label_transfer import label_transfer
import scanpy as sc
rna_anndata = sc.read_h5ad("data/adata.h5ad").raw.to_adata()
atac_anndata = sc.AnnData(gene_act.mtx.T, obs=pd.DataFrame(index=gene_act.cell_names), var=pd.DataFrame(index=gene_act.feature_names))

atac_anndata.obs["sample_id"] = "10x_multiome_brain"
rna_anndata.obs["sample_id"] = "10x_multiome_brain"

label_dict = label_transfer(rna_anndata, atac_anndata,
    labels_to_transfer=['Seurat_cell_type'], variable_genes=True,
    methods = ['ingest', 'harmony', 'bbknn', 'scanorama', 'cca'],
    return_label_weights=False, _temp_dir= '/tmp/ray_spill')

现在可以将注释添加到cisTopic对象中，并在细胞-主题UMAP上可视化它们。Scanorama和harmony效果最好。

label_dict_x = [label_dict[key] for key in label_dict.keys()]
label_pd = pd.concat(label_dict_x, axis=1, sort=False)
label_pd.index = cistopic_obj.cell_names
label_pd.columns = ['pycisTopic_' + x for x in label_pd.columns]
cistopic_obj.add_cell_data(label_pd, split_pattern = '-')

plot_metadata(cistopic_obj, reduction_name='UMAP', variables=['Seurat_cell_type'] + label_pd.columns.to_list(), num_columns=3)

绘制混淆矩阵：

import seaborn as sns
fig, axs = plt.subplots(ncols = 3, nrows = 2, figsize = (3 * 5, 2 * 5))
for method, ax in zip(label_pd.columns.to_list(), axs.ravel()):
    conf_mat = pd.crosstab(cistopic_obj.cell_data["Seurat_cell_type"], cistopic_obj.cell_data[method])
    conf_mat = conf_mat / conf_mat.sum()
    sns.heatmap(conf_mat.loc[conf_mat.columns], ax = ax)
fig.tight_layout()
fig.show()

保持结果

可以保存为loom文件方便后续进一步探索结果。有两种类型的loom文件：

• 区域可及性：这些loom文件包括作为矩阵的推算可及性、作为调控子(regulons)的主题以及作为AUC矩阵的细胞-主题分布。推算值、用于推算的cistopic对象以及细胞-主题和主题-区域的二值化分布都是必需的。或者，我们也可以提供不同的聚类以及每个聚类每组的标记区域 (DARs)。
• 基因活性：这些loom文件包含作为矩阵的基因活性、基于scRNA-seq衍生的调控子及其基于基因活性的AUC值。基因活性值、cistopic对象和scRNA-seq衍生的调控子是必需的。或者，我们也可以提供不同的聚类以及每个聚类每组的标记基因 (DAGs)。

from pycisTopic.loom import export_region_accessibility_to_loom, export_gene_activity_to_loom

cluster_markers = {'Seurat_cell_type': markers_dict}
os.makedirs(os.path.join(out_dir, "loom"), exist_ok=True)

export_region_accessibility_to_loom(
    accessibility_matrix = imputed_acc_obj,
    cistopic_obj = cistopic_obj,
    binarized_topic_region = region_bin_topics_top_3k,
    binarized_cell_topic = binarized_cell_topic,
    selected_cells = cistopic_obj.projections['cell']['UMAP'].index.tolist(),
    out_fname = os.path.join(out_dir, "loom", "10x_multiome_brain_pycistopic_region_accessibility.loom"),
    cluster_annotation = ['Seurat_cell_type'],
    cluster_markers = cluster_markers,
    tree_structure = ('10x_multiome_brain', 'pycistopic', 'nodbl_all'),
    title = 'Region accessibility all',
    nomenclature = "hg38",
    split_pattern = '-')

export_gene_activity_to_loom(
    gene_activity_matrix = gene_act,
    cistopic_obj = cistopic_obj,
    out_fname = os.path.join(out_dir, "loom", "10x_multiome_brain_pycistopic_gene_activity.loom"),
    cluster_annotation = ['Seurat_cell_type'],
    cluster_markers = cluster_markers,
    tree_structure = ('10x_multiome_brain', 'pycistopic', 'atac'),
    title = 'Gene activity',
    nomenclature = "hg38",
    split_pattern = '-')