[1]:

import scanpy as sc
import pandas as pd
from scipy.sparse import csr_matrix
from scipy import sparse
import scipy.io
import os
import anndata as ad # Anndata version must > 0.8
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme(style='white')

from scATAnno.SnapATAC2_spectral import *
from scATAnno.SnapATAC2_tools import *
from scATAnno.SnapATAC2_utils import *
from scATAnno import scATAnno_preprocess
from scATAnno import scATAnno_assignment
from scATAnno import scATAnno_integration
from scATAnno import scATAnno_plotting

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.simplefilter(action='ignore', category=DeprecationWarning)

default_28 = scATAnno_plotting.get_palettes("default28")
default_102 = scATAnno_plotting.get_palettes("default102")

TNBC sample

Download TIL reference data

[2]:

# !wget -O BCC_TIL_reference_atlas_final.h5ad "https://www.dropbox.com/s/ky4jezsj3pf2qwi/BCC_TIL_reference_atlas_final.h5ad?dl=1"

Download TNBC query data, whose features are aligned with BCC TIL reference peak

[3]:

# !wget -O TNBC_Pre_P012_t.zip "https://www.dropbox.com/s/7vyk08xrv81b1q2/TNBC_Pre_P012_t.zip?dl=1"

Donwload palette file

[4]:

# !wget -O TIL_reference_palette.pickle "https://www.dropbox.com/scl/fi/tswxdbzz47g97kduanbi6/TIL_reference_palette.pickle?rlkey=c2vibsg41exfxcqrhtiu32ikg&st=eyiyby6e&dl=0"
# !wget -O all_palette_dictionary.pkl "https://www.dropbox.com/scl/fi/oi2ykercl2ugguihogyzi/all_palette_dictionary.pkl?rlkey=lc7wbejo25p512hgw56wuo3zw&st=q8us6lpa&dl=0"

[5]:

os.chdir("scATAnno-main") #change to your own directory

output_name = "TNBC"
out_dir = os.path.join("scATAnno-main/case_study", output_name)
os.makedirs(out_dir, exist_ok=True)

[6]:

reference_data_path = "BCC_TIL_reference_atlas_final.h5ad"
reference_data = scATAnno_preprocess.load_reference_data(reference_data_path)
query_data = scATAnno_preprocess.import_query_data(path = 'TNBC_Pre_P012_t/',
                                    mtx_file = 'matrix.mtx.gz',
                                    cells_file = 'barcodes.tsv.gz',
                                    features_file = 'features.tsv.gz',
                                    variable_prefix = "TNBC_Pre_P012_t",
                                    celltype_col = "celltypes",
                                    add_metrics=False)

[7]:

print(reference_data)
print(query_data)
assert reference_data.var.shape[0] == query_data.var.shape[0]

AnnData object with n_obs × n_vars = 22008 × 344492
    obs: 'celltypes', 'tissue', 'dataset'
    var: 'selected'
    uns: 'Group_colors', 'celltypes_colors', 'spectral_eigenvalue', 'true_label_colors'
    obsm: 'X_spectral', 'X_umap'
View of AnnData object with n_obs × n_vars = 9935 × 344492
    obs: 'celltypes', 'tissue', 'dataset'

[8]:

# Integrate reference and query data
integrated_adata = scATAnno_assignment.scATAnno_integrate(reference_data, query_data, variable_prefix = "TNBC_Pre_P012_t", sample_size = 25000)

/Users/francis/Desktop/scATAnno-main/scATAnno/SnapATAC2_spectral.py:162: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual.
  adata.var["selected"] = selected_features

Compute similarity matrix
Normalization
Perform decomposition
Perform Nystrom extension

100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [03:11<00:00, 95.68s/it]

[9]:

# Apply harmony to remove batch effects
integrated_adata_harmony = scATAnno_assignment.scATAnno_harmony(integrated_adata, batch_col = "dataset")

2025-01-07 14:19:40,505 - harmonypy - INFO - Computing initial centroids with sklearn.KMeans...
2025-01-07 14:19:40 - INFO - Computing initial centroids with sklearn.KMeans...
2025-01-07 14:19:45,488 - harmonypy - INFO - sklearn.KMeans initialization complete.
2025-01-07 14:19:45 - INFO - sklearn.KMeans initialization complete.
2025-01-07 14:19:45,568 - harmonypy - INFO - Iteration 1 of 20
2025-01-07 14:19:45 - INFO - Iteration 1 of 20
2025-01-07 14:19:49,605 - harmonypy - INFO - Iteration 2 of 20
2025-01-07 14:19:49 - INFO - Iteration 2 of 20
2025-01-07 14:19:54,091 - harmonypy - INFO - Converged after 2 iterations
2025-01-07 14:19:54 - INFO - Converged after 2 iterations

[10]:

# Plot UMAP using spectral embeddings
integrated_adata = scATAnno_assignment.scATAnno_umap(integrated_adata_harmony, out_dir, use_rep = "X_spectral_harmony", save = True)
integrated_adata

/opt/anaconda3/envs/scATAnno/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
/opt/anaconda3/envs/scATAnno/lib/python3.11/site-packages/umap/umap_.py:1952: UserWarning: n_jobs value 1 overridden to 1 by setting random_state. Use no seed for parallelism.
  warn(
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
... storing 'celltypes' as categorical
... storing 'tissue' as categorical

[10]:

AnnData object with n_obs × n_vars = 31943 × 344492
    obs: 'celltypes', 'tissue', 'dataset'
    obsm: 'X_spectral', 'X_spectral_harmony', 'X_umap'

[11]:

scATAnno_plotting.defaultPlotting_umap()
sc.pl.umap(integrated_adata, color="dataset", palette = ['#228833', '#cccccc'], show=True, title = "Integration")

... storing 'celltypes' as categorical
... storing 'tissue' as categorical

[12]:

integrated_adata = sc.read_h5ad(os.path.join(out_dir, "1.Merged_query_reference.h5ad"))
integrated_adata = scATAnno_preprocess.add_variable(os.path.join(out_dir,"X_spectral_harmony.csv"), integrated_adata)

[13]:

reference = integrated_adata[integrated_adata.obs['dataset'] == "Atlas",:].copy()
reference.obs["celltypes"] =  scATAnno_assignment.curate_celltype_names(reference.obs["celltypes"], atlas = "TIL")

[14]:

import pickle
with open('TIL_reference_palette.pickle', 'rb') as handle:
    celltype_palette = pickle.load(handle)

scATAnno_plotting.defaultPlotting_umap()
sc.pl.umap(reference, color="celltypes", palette = celltype_palette, title = "TIL Reference Atlas")

[15]:

query = integrated_adata[integrated_adata.obs['dataset'] != "Atlas",:].copy()

[16]:

projecton = pd.read_csv(os.path.join("TNBC_Pre_P012_t/", "query_metadata.csv"), index_col=0)
query.obs = pd.merge(query.obs, projecton, left_index=True, right_index=True)

We first assign cell types based on KNN and compute KNN-based uncertainty score

[17]:

# Perform KNN assignment
reference_label_col = "celltypes"
use_rep = "X_spectral_harmony"
query_KNN = scATAnno_assignment.scATAnno_KNN_assign(reference, query, reference_label_col=reference_label_col, low_dim_col=use_rep)

Because we are annotating TIL subtypes, we adopt more loose threshold cutoffs to annotate query cells. Query cells are assigned high weighted distance-based uncertainty score if their distances to the assigned reference cell type centroid are greater than 95 percentile. Query cells with uncertainty score greater than 0.5 are annotated as unknown.

[18]:

# Perform weighted-distance based assignment
distance_threshold = 95
uncertainty_threshold = 0.5

atlas = "TIL"
query_distance = scATAnno_assignment.scATAnno_distance_assign(reference, query_KNN, reference_label_col=reference_label_col, distance_threshold=distance_threshold, atlas=atlas, uncertainty_threshold=uncertainty_threshold, use_rep = use_rep)

Finally, we perform cluster-level assignment of query cells.

[19]:

query_annotated = scATAnno_assignment.scATAnno_cluster_assign(query_distance, use_rep=use_rep, cluster_col = "Clusters")

[20]:

query_annotated

[20]:

AnnData object with n_obs × n_vars = 9935 × 344492
    obs: 'celltypes', 'tissue', 'dataset', 'paper_annotation', 'paper_annotation_major', 'paper_annotation_curated', 'Clusters', 'uncertainty_score_step1', 'pred_y', 'uncertainty_score_step2', 'Uncertainty_Score', '1.knn-based_celltype', '2.corrected_celltype', 'cluster_annotation'
    uns: 'neighbors', 'umap'
    obsm: 'X_umap', 'X_spectral_harmony', 'kernel_distance', 'distance', 'indices', 'neighbors_labels'
    obsp: 'distances', 'connectivities'

[21]:

import pickle
with open('all_palette_dictionary.pkl', 'rb') as handle:
    annotation_palette = pickle.load(handle)
sc.pl.umap(query_annotated, color = ["paper_annotation_curated", "cluster_annotation", 'Uncertainty_Score'], palette=annotation_palette,title = ["Annotation Published in Zhang et al.","scATAnno Annotation"],  legend_loc = "on data")

... storing 'paper_annotation' as categorical
... storing 'paper_annotation_major' as categorical
... storing 'paper_annotation_curated' as categorical
... storing 'Clusters' as categorical
... storing 'pred_y' as categorical
... storing '1.knn-based_celltype' as categorical
... storing '2.corrected_celltype' as categorical
... storing 'cluster_annotation' as categorical

WARNING: The title list is shorter than the number of panels. Using 'color' value instead for some plots.

Compare the annotation results between scATAnno and original publication.

[22]:

from pySankey.sankey import sankey

pred_col = 'cluster_annotation'
true_col = 'paper_annotation_curated'

sankey(
    query_annotated.obs[pred_col], query_annotated.obs[true_col], aspect=10, colorDict=annotation_palette,
    fontsize=10, figure_name=None
    )

# plt.gcf().set_size_inches(6, 10)
# plt.savefig(os.path.join(out_dir,"{}.png".format("Sankey_CellType_Annotation_paper_cluster")), bbox_inches='tight', dpi=300)