BCC Multiple Samples
[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')
import scATAnno
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")
Part I: First round annotation - Apply Healthy Adult reference atlas to BCC tumor samples
We first apply healthy adult reference atlas to two BCC samples. We load reference atlas data and BCC query data with healthy adult reference peak as feature input. Reference atlas is stored as h5ad AnnData in the data directory; query data is imported as MTX and TSV files from QuickATAC output.
Download reference data
[2]:
# !wget -O Healthy_Adult_reference_atlas.h5ad "https://www.dropbox.com/s/3ezp2t6gw6hw21v/Healthy_Adult_reference_atlas.h5ad?dl=1"
Download query data. Each BCC sample has two features, one aligned with Healthy adult reference peak, and the other aligned with BCC TIL reference peak.
[3]:
# !wget -O BCC_samples.zip "https://www.dropbox.com/s/eqkmt5trp3350cj/BCC_samples.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 HealthyAdult_reference_palette.pickle "https://www.dropbox.com/scl/fi/9nf09modjh9wc5bxhplr1/HealthyAdult_reference_palette.pickle?rlkey=huztqlg9axvj81oe9v6umso94&st=xk0srtmp&dl=0"
[5]:
os.chdir("scATAnno-main") #change to your own directory
output_name = "BCC_samples"
out_dir = os.path.join("case_study", output_name)
os.makedirs(out_dir, exist_ok=True)
reference_data_path = "Healthy_Adult_reference_atlas.h5ad"
reference_data = scATAnno_preprocess.load_reference_data(reference_data_path)
SU007_Total_Post = scATAnno_preprocess.import_query_data(path = 'BCC_samples/SU007_Total_Post_vHealthyAdult',
mtx_file = 'matrix.mtx.gz',
cells_file = 'barcodes.tsv.gz',
features_file = 'features.tsv.gz',
variable_prefix = "SU007_Total_Post",
celltype_col = "celltypes",
add_metrics=False)
SU010_Total_Pre = scATAnno_preprocess.import_query_data(path = 'BCC_samples/SU010_Total_Pre_vHealthyAdult',
mtx_file = 'matrix.mtx.gz',
cells_file = 'barcodes.tsv',
features_file = 'features.tsv.gz',
variable_prefix = "SU010_Total_Pre",
celltype_col = "celltypes",
add_metrics=False)
[6]:
print(SU010_Total_Pre)
print(SU007_Total_Post)
print(reference_data)
assert reference_data.var.shape[0] == SU010_Total_Pre.var.shape[0]
assert reference_data.var.shape[0] == SU007_Total_Post.var.shape[0]
View of AnnData object with n_obs × n_vars = 1071 × 890130
obs: 'celltypes', 'tissue', 'dataset'
View of AnnData object with n_obs × n_vars = 666 × 890130
obs: 'celltypes', 'tissue', 'dataset'
AnnData object with n_obs × n_vars = 100158 × 890130
obs: 'celltypes', 'tissue', 'dataset'
var: 'selected'
uns: 'cell type_colors', 'major cell type_colors', 'spectral_eigenvalue', 'tissue_colors'
obsm: 'X_umap'
We first create a dictionary of reference and query datasets
[7]:
data_list = [SU007_Total_Post, SU010_Total_Pre]
print("total number of query cells: {}".format(np.sum([i.obs.shape[0] for i in data_list])))
total number of query cells: 1737
Before integration, we need to select features of datasets and prepare for integration
[8]:
select_features(reference_data)
for anndata in data_list:
select_features(anndata)
datasets={}
datasets["Atlas"] = reference_data
for anndata in data_list:
key = (anndata.obs["dataset"])[0]
datasets[key] = anndata
datasets
/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
/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
[8]:
{'Atlas': AnnData object with n_obs × n_vars = 100158 × 890130
obs: 'celltypes', 'tissue', 'dataset'
var: 'selected'
uns: 'cell type_colors', 'major cell type_colors', 'spectral_eigenvalue', 'tissue_colors'
obsm: 'X_umap',
'SU007_Total_Post': AnnData object with n_obs × n_vars = 666 × 890130
obs: 'celltypes', 'tissue', 'dataset'
var: 'selected',
'SU010_Total_Pre': AnnData object with n_obs × n_vars = 1071 × 890130
obs: 'celltypes', 'tissue', 'dataset'
var: 'selected'}
We use SnapATAC2 to integrate the reference atlas and multiple query data
[9]:
# Integrate reference and query data
integrated_adata = scATAnno_assignment.scATAnno_integrate_multiple(datasets, sample_size = 25000)
Compute similarity matrix
Normalization
Perform decomposition
Perform Nystrom extension
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 6/6 [04:35<00:00, 45.98s/it]
Apply harmony to remove batch effects of datasets
[10]:
# Apply harmony to remove batch effects
integrated_adata_harmony = scATAnno_assignment.scATAnno_harmony(integrated_adata, batch_col = "dataset")
2025-01-07 14:44:11,303 - harmonypy - INFO - Computing initial centroids with sklearn.KMeans...
2025-01-07 14:44:11 - INFO - Computing initial centroids with sklearn.KMeans...
2025-01-07 14:44:18,261 - harmonypy - INFO - sklearn.KMeans initialization complete.
2025-01-07 14:44:18 - INFO - sklearn.KMeans initialization complete.
2025-01-07 14:44:18,580 - harmonypy - INFO - Iteration 1 of 20
2025-01-07 14:44:18 - INFO - Iteration 1 of 20
2025-01-07 14:44:42,622 - harmonypy - INFO - Iteration 2 of 20
2025-01-07 14:44:42 - INFO - Iteration 2 of 20
2025-01-07 14:45:06,559 - harmonypy - INFO - Converged after 2 iterations
2025-01-07 14:45:06 - INFO - Converged after 2 iterations
Plot UMAP using spectral embeddings
[11]:
# Plot UMAP using spectral embeddings
integrated_adata = scATAnno_assignment.scATAnno_umap(integrated_adata_harmony, out_dir, use_rep = "X_spectral_harmony", save = True)
/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
[12]:
# Set plotting parameters
scATAnno_plotting.defaultPlotting_umap()
sc.pl.umap(integrated_adata, color="dataset", palette = ['#cccccc', '#228833','#8c564b'], size=10, show=True, title = "Integration")
... storing 'celltypes' as categorical
... storing 'tissue' as categorical
[13]:
integrated_adata.obs
[13]:
| celltypes | tissue | dataset | |
|---|---|---|---|
| adipose_omentum_SM-CSSD4_1+AAGTCCTTAGGGCTGATATGCT_Atlas | Mesenchymal | adipose_omentum_SM-CSSD4_1+AAGTCCTTAGGGCTGATATGCT | Atlas |
| adipose_omentum_SM-CSSD4_1+AACGACCAAAGGCGTACTTCTA_Atlas | Mesenchymal | adipose_omentum_SM-CSSD4_1+AACGACCAAAGGCGTACTTCTA | Atlas |
| adipose_omentum_SM-IOBHJ_1+ATACTCGCTTAGCGGCTGAAAC_Atlas | Mesenchymal | adipose_omentum_SM-IOBHJ_1+ATACTCGCTTAGCGGCTGAAAC | Atlas |
| adipose_omentum_SM-CSSD4_1+CAGGAAAGCAAAAGGGATGCCA_Atlas | Mesenchymal | adipose_omentum_SM-CSSD4_1+CAGGAAAGCAAAAGGGATGCCA | Atlas |
| adipose_omentum_SM-IOBHJ_1+ATGCAGGTAAATAGCTGTGTCT_Atlas | Mesenchymal | adipose_omentum_SM-IOBHJ_1+ATGCAGGTAAATAGCTGTGTCT | Atlas |
| ... | ... | ... | ... |
| SU010_Total_Pre#TGGGTTATCCCAGCAG-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre |
| SU010_Total_Pre#GGGTGTCTCGGGACAA-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre |
| SU010_Total_Pre#ACCGGGTCATCCCAAA-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre |
| SU010_Total_Pre#CCCGTTATCTAGCAGT-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre |
| SU010_Total_Pre#GAACGTTGTTAAGTCC-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre |
101895 rows × 3 columns
[14]:
reference = integrated_adata[integrated_adata.obs['dataset'] == "Atlas",:].copy()
reference.obs["celltypes"] = scATAnno_assignment.curate_celltype_names(reference.obs["celltypes"], atlas = "HealthyAdult")
[15]:
import pickle
with open('HealthyAdult_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 = "Healthy Adult Reference Atlas")
... storing 'celltypes' as categorical
After integration, we want to assign cell types to the query data. First, query cells are assigned labels using K-Nearest neighbors along with KNN-based uncertainty score. Second, query cell labels are corrected using weighted-distance-based uncertainty score. Third, query cells are clustered and annotated at the cluster level.
[16]:
query = integrated_adata[integrated_adata.obs['dataset'] != "Atlas",:].copy()
Obtain UMAP coordinates and clusters provided by the original publication. Then Query cells are clustered and annotated at the cluster-level
[17]:
projection_allcells = pd.read_table("BCC_samples/BCC_study_metadata.csv", index_col = 0, sep = ",")
#projection_allcells.index = [projection_allcells.loc[i,"Group"] + "#" + i for i in projection_allcells.index]
sample_names = ["SU010_Total_Pre", "SU007_Total_Post"]
projection = projection_allcells[projection_allcells["Group"].isin(sample_names)]
query.obs = pd.merge(query.obs, projection, left_index=True, right_index=True)
projection.head()
[17]:
| UMAP_1 | UMAP_2 | Clusters | Group | depth | FRIP | Barcodes | Internal_Name | Group_Barcode | true_label | |
|---|---|---|---|---|---|---|---|---|---|---|
| SU007_Total_Post#GGACACCAGCGCACAA-1_SU007_Total_Post | 1.320694 | 12.943564 | Cluster14 | SU007_Total_Post | 50427 | 0.676443 | GGACACCAGCGCACAA-1 | SU007_Total_Post_57 | SU007_Total_Post#GGACACCAGCGCACAA-1 | Myeloid |
| SU007_Total_Post#ACTAGGTGTACCCATA-1_SU007_Total_Post | 5.631482 | 13.093065 | Cluster15 | SU007_Total_Post | 45889 | 0.580553 | ACTAGGTGTACCCATA-1 | SU007_Total_Post_59 | SU007_Total_Post#ACTAGGTGTACCCATA-1 | Endothelial |
| SU007_Total_Post#AGGACGATCGGTCTCT-1_SU007_Total_Post | 1.471412 | -7.473051 | Cluster4 | SU007_Total_Post | 21568 | 0.688590 | AGGACGATCGGTCTCT-1 | SU007_Total_Post_60 | SU007_Total_Post#AGGACGATCGGTCTCT-1 | Treg |
| SU007_Total_Post#GACCAATCATGGGTGA-1_SU007_Total_Post | -4.932662 | 2.046299 | Cluster7 | SU007_Total_Post | 25402 | 0.714097 | GACCAATCATGGGTGA-1 | SU007_Total_Post_61 | SU007_Total_Post#GACCAATCATGGGTGA-1 | Memory CD8 T |
| SU007_Total_Post#GCTCGAGAGTAATGTG-1_SU007_Total_Post | -4.171803 | 3.467478 | Cluster7 | SU007_Total_Post | 6544 | 0.714318 | GCTCGAGAGTAATGTG-1 | SU007_Total_Post_62 | SU007_Total_Post#GCTCGAGAGTAATGTG-1 | Memory CD8 T |
We first assign cell types based on KNN and compute KNN-based uncertainty score
Query cells are assigned high weighted distance-based uncertainty score if their distances to the assigned reference cell type centroid are greater than 90 percentile. Query cells with uncertainty score greater than 0.2 are annotated as unknown.
[18]:
# Perform weighted-distance based assignment
distance_threshold = 90
uncertainty_threshold = 0.2
reference_label_col = "celltypes"
use_rep = "X_spectral_harmony"
atlas = "HealthyAdult"
[19]:
# Perform KNN assignment
query_KNN = scATAnno_assignment.scATAnno_KNN_assign(reference, query, reference_label_col=reference_label_col, low_dim_col=use_rep)
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)
We then project query cells into the same UMAP coordinates as original publication for better comparisons, and perform cluster-level assignment.
[20]:
query_distance.obsm["X_umap"] = np.array(query_distance.obs[["UMAP_1", "UMAP_2"]])
query_annotated = scATAnno_assignment.scATAnno_cluster_assign(query_distance, use_rep=use_rep, cluster_col = "Clusters", UMAP=False)
query_annotated
[20]:
AnnData object with n_obs × n_vars = 1737 × 890130
obs: 'celltypes', 'tissue', 'dataset', 'UMAP_1', 'UMAP_2', 'Clusters', 'Group', 'depth', 'FRIP', 'Barcodes', 'Internal_Name', 'Group_Barcode', 'true_label', 'uncertainty_score_step1', 'pred_y', 'uncertainty_score_step2', 'Uncertainty_Score', '1.knn-based_celltype', '2.corrected_celltype', 'cluster_annotation'
uns: 'dataset_colors'
obsm: 'X_spectral', 'X_spectral_harmony', 'X_umap', 'kernel_distance', 'distance', 'indices', 'neighbors_labels'
Compare annotation results from scATAnno and the origional publication
[21]:
sc.pl.umap(query_annotated, color = ['true_label', "cluster_annotation", 'Uncertainty_Score'], title = ["Annotation Published in Satpathy et al.","scATAnno Annotation"], legend_loc = "on data")
... storing 'Clusters' as categorical
... storing 'Group' as categorical
... storing 'Barcodes' as categorical
... storing 'true_label' 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.
[22]:
immune_cells = pd.DataFrame(query_annotated[query_annotated.obs["cluster_annotation"] == "Immune Cells"].obs.index)
print(immune_cells.head())
0
0 SU007_Total_Post#TCGTTCGAGGTTCTCA-1_SU007_Tota...
1 SU007_Total_Post#GCAACCGAGAAGAGTG-1_SU007_Tota...
2 SU007_Total_Post#GTGTCAAGTTAGTAGA-1_SU007_Tota...
3 SU007_Total_Post#CTTGTCGAGCGAGCTA-1_SU007_Tota...
4 SU007_Total_Post#GCAGCTGCATTCGTCC-1_SU007_Tota...
Part II: Second round annotation - Apply BCC TIL reference atlas to BCC immune cells
Download TIL reference data
[23]:
# !wget "https://www.dropbox.com/s/ky4jezsj3pf2qwi/BCC_TIL_reference_atlas_final.h5ad?dl=1"
In the second round annotation, we apply BCC TIL reference atlas to the immune cells identified from the first round. Repeat the integration process. This time we use BCC TIL as reference, and use BCC query data with TIL reference peak as feature input. Reference atlas is stored as h5ad AnnData in the data directory; query data is imported as MTX and TSV files from QuickATAC output.
[24]:
reference_data_path = "BCC_TIL_reference_atlas_final.h5ad"
reference_data = scATAnno_preprocess.load_reference_data(reference_data_path)
SU007_Total_Post = scATAnno_preprocess.import_query_data(path = 'BCC_samples/SU007_Total_Post_vBCC',
mtx_file = 'matrix.mtx.gz',
cells_file = 'barcodes.tsv.gz',
features_file = 'features.tsv.gz',
variable_prefix = "SU007_Total_Post",
celltype_col = "celltypes",
add_metrics=False)
SU010_Total_Pre = scATAnno_preprocess.import_query_data(path = 'BCC_samples/SU010_Total_Pre_vBCC',
mtx_file = 'matrix.mtx.gz',
cells_file = 'barcodes.tsv',
features_file = 'features.tsv.gz',
variable_prefix = "SU010_Total_Pre",
celltype_col = "celltypes",
add_metrics=False)
[25]:
data_list = [SU007_Total_Post, SU010_Total_Pre]
print("Total number of query cells: {}".format(np.sum([i.obs.shape[0] for i in data_list])))
Total number of query cells: 1737
Select the immune cells from first round annotation
[26]:
immune_cells.index = immune_cells[0]
immune_cells.index = [i.split("_SU")[0] for i in immune_cells.index]
immune_cells.head()
[26]:
| 0 | |
|---|---|
| SU007_Total_Post#TCGTTCGAGGTTCTCA-1 | SU007_Total_Post#TCGTTCGAGGTTCTCA-1_SU007_Tota... |
| SU007_Total_Post#GCAACCGAGAAGAGTG-1 | SU007_Total_Post#GCAACCGAGAAGAGTG-1_SU007_Tota... |
| SU007_Total_Post#GTGTCAAGTTAGTAGA-1 | SU007_Total_Post#GTGTCAAGTTAGTAGA-1_SU007_Tota... |
| SU007_Total_Post#CTTGTCGAGCGAGCTA-1 | SU007_Total_Post#CTTGTCGAGCGAGCTA-1_SU007_Tota... |
| SU007_Total_Post#GCAGCTGCATTCGTCC-1 | SU007_Total_Post#GCAGCTGCATTCGTCC-1_SU007_Tota... |
[27]:
data_list_use = data_list.copy()
for i in range(len(data_list_use)):
print(data_list_use[i].obs["dataset"][0])
data_list_use[i] = data_list_use[i][data_list_use[i].obs.index.isin(immune_cells.index)]
print(data_list_use[i])
print("Total immune cells: {}".format(np.sum([i.obs.shape[0] for i in data_list_use])))
SU007_Total_Post
View of AnnData object with n_obs × n_vars = 618 × 344492
obs: 'celltypes', 'tissue', 'dataset'
SU010_Total_Pre
View of AnnData object with n_obs × n_vars = 188 × 344492
obs: 'celltypes', 'tissue', 'dataset'
Total immune cells: 806
We use SnapATAC2 to integrate the reference atlas and multiple query data
[28]:
for anndata in data_list_use:
select_features(anndata)
select_features(reference_data)
datasets={}
datasets["Atlas"] = reference_data
for anndata in data_list_use:
key = (anndata.obs["dataset"])[0]
datasets[key] = anndata
datasets
/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
/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
[28]:
{'Atlas': 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',
'SU007_Total_Post': AnnData object with n_obs × n_vars = 618 × 344492
obs: 'celltypes', 'tissue', 'dataset'
var: 'selected',
'SU010_Total_Pre': AnnData object with n_obs × n_vars = 188 × 344492
obs: 'celltypes', 'tissue', 'dataset'
var: 'selected'}
[29]:
# Integrate reference and query data
integrated_adata = scATAnno_assignment.scATAnno_integrate_multiple(datasets, sample_size = 25000)
Compute similarity matrix
Normalization
Perform decomposition
[30]:
# Apply harmony to remove batch effects
integrated_adata_harmony = scATAnno_assignment.scATAnno_harmony(integrated_adata, batch_col = "dataset")
2025-01-07 14:51:26,614 - harmonypy - INFO - Computing initial centroids with sklearn.KMeans...
2025-01-07 14:51:26 - INFO - Computing initial centroids with sklearn.KMeans...
2025-01-07 14:51:30,730 - harmonypy - INFO - sklearn.KMeans initialization complete.
2025-01-07 14:51:30 - INFO - sklearn.KMeans initialization complete.
2025-01-07 14:51:30,789 - harmonypy - INFO - Iteration 1 of 20
2025-01-07 14:51:30 - INFO - Iteration 1 of 20
2025-01-07 14:51:35,185 - harmonypy - INFO - Iteration 2 of 20
2025-01-07 14:51:35 - INFO - Iteration 2 of 20
2025-01-07 14:51:38,977 - harmonypy - INFO - Iteration 3 of 20
2025-01-07 14:51:38 - INFO - Iteration 3 of 20
2025-01-07 14:51:42,478 - harmonypy - INFO - Iteration 4 of 20
2025-01-07 14:51:42 - INFO - Iteration 4 of 20
2025-01-07 14:51:46,096 - harmonypy - INFO - Iteration 5 of 20
2025-01-07 14:51:46 - INFO - Iteration 5 of 20
2025-01-07 14:51:49,609 - harmonypy - INFO - Iteration 6 of 20
2025-01-07 14:51:49 - INFO - Iteration 6 of 20
2025-01-07 14:51:53,197 - harmonypy - INFO - Iteration 7 of 20
2025-01-07 14:51:53 - INFO - Iteration 7 of 20
2025-01-07 14:51:57,241 - harmonypy - INFO - Iteration 8 of 20
2025-01-07 14:51:57 - INFO - Iteration 8 of 20
2025-01-07 14:52:01,240 - harmonypy - INFO - Iteration 9 of 20
2025-01-07 14:52:01 - INFO - Iteration 9 of 20
2025-01-07 14:52:04,683 - harmonypy - INFO - Iteration 10 of 20
2025-01-07 14:52:04 - INFO - Iteration 10 of 20
2025-01-07 14:52:08,226 - harmonypy - INFO - Iteration 11 of 20
2025-01-07 14:52:08 - INFO - Iteration 11 of 20
2025-01-07 14:52:11,693 - harmonypy - INFO - Iteration 12 of 20
2025-01-07 14:52:11 - INFO - Iteration 12 of 20
2025-01-07 14:52:15,079 - harmonypy - INFO - Converged after 12 iterations
2025-01-07 14:52:15 - INFO - Converged after 12 iterations
[31]:
# Plot UMAP using spectral embeddings
integrated_adata = scATAnno_assignment.scATAnno_umap(integrated_adata_harmony, out_dir, use_rep = "X_spectral_harmony", save = True)
/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(
... storing 'celltypes' as categorical
... storing 'tissue' as categorical
[32]:
# Set plotting parameters
scATAnno_plotting.defaultPlotting_umap()
sc.pl.umap(integrated_adata, color="dataset", palette = ['#cccccc', '#228833','#8c564b'], size=10, show=True, title = "Integration")
... storing 'celltypes' as categorical
... storing 'tissue' as categorical
[33]:
reference = integrated_adata[integrated_adata.obs['dataset'] == "Atlas",:].copy()
reference.obs["celltypes"] = scATAnno_assignment.curate_celltype_names(reference.obs["celltypes"], atlas = "TIL")
[34]:
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")
[35]:
query = integrated_adata[integrated_adata.obs['dataset'] != "Atlas",:].copy()
[36]:
projection_allcells = pd.read_table("BCC_samples/BCC_study_metadata.csv", index_col = 0, sep = ",")
#projection_allcells.index = [projection_allcells.loc[i,"Group"] + "#" + i for i in projection_allcells.index]
projection_allcells["true_label"] = [i.strip() for i in projection_allcells["true_label"]]
sample_names = ["SU010_Total_Pre", "SU007_Total_Post"]
projection = projection_allcells[projection_allcells["Group"].isin(sample_names)]
query.obs = pd.merge(query.obs, projection, left_index=True, right_index=True)
query.obs
[36]:
| celltypes | tissue | dataset | UMAP_1 | UMAP_2 | Clusters | Group | depth | FRIP | Barcodes | Internal_Name | Group_Barcode | true_label | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| SU007_Total_Post#TCGTTCGAGGTTCTCA-1_SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | 2.118399 | -6.020474 | Cluster4 | SU007_Total_Post | 10159 | 0.722758 | TCGTTCGAGGTTCTCA-1 | SU007_Total_Post_312 | SU007_Total_Post#TCGTTCGAGGTTCTCA-1 | Treg |
| SU007_Total_Post#GCAACCGAGAAGAGTG-1_SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | 2.138090 | -8.092027 | Cluster4 | SU007_Total_Post | 12455 | 0.715857 | GCAACCGAGAAGAGTG-1 | SU007_Total_Post_1 | SU007_Total_Post#GCAACCGAGAAGAGTG-1 | Treg |
| SU007_Total_Post#GTGTCAAGTTAGTAGA-1_SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | 0.018025 | -1.757976 | Cluster1 | SU007_Total_Post | 9101 | 0.621910 | GTGTCAAGTTAGTAGA-1 | SU007_Total_Post_468 | SU007_Total_Post#GTGTCAAGTTAGTAGA-1 | Naive CD4 T |
| SU007_Total_Post#CTTGTCGAGCGAGCTA-1_SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | -1.315312 | 1.582670 | Cluster5 | SU007_Total_Post | 22287 | 0.669695 | CTTGTCGAGCGAGCTA-1 | SU007_Total_Post_3 | SU007_Total_Post#CTTGTCGAGCGAGCTA-1 | Naive CD8 T |
| SU007_Total_Post#GCAGCTGCATTCGTCC-1_SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | SU007_Total_Post | 0.688227 | -1.430185 | Cluster1 | SU007_Total_Post | 22482 | 0.693866 | GCAGCTGCATTCGTCC-1 | SU007_Total_Post_4 | SU007_Total_Post#GCAGCTGCATTCGTCC-1 | Naive CD4 T |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| SU010_Total_Pre#CCAGATAAGTACGCGA-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | -3.657552 | 2.379094 | Cluster7 | SU010_Total_Pre | 6428 | 0.583852 | CCAGATAAGTACGCGA-1 | SU010_Total_Pre_1057 | SU010_Total_Pre#CCAGATAAGTACGCGA-1 | Memory CD8 T |
| SU010_Total_Pre#ACAGGCCTCAAGAGAT-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | -3.117948 | 1.990111 | Cluster5 | SU010_Total_Pre | 3820 | 0.529188 | ACAGGCCTCAAGAGAT-1 | SU010_Total_Pre_1058 | SU010_Total_Pre#ACAGGCCTCAAGAGAT-1 | Naive CD8 T |
| SU010_Total_Pre#TGGGTTATCCCAGCAG-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | -1.875685 | 0.338385 | Cluster5 | SU010_Total_Pre | 4741 | 0.469099 | TGGGTTATCCCAGCAG-1 | SU010_Total_Pre_1061 | SU010_Total_Pre#TGGGTTATCCCAGCAG-1 | Naive CD8 T |
| SU010_Total_Pre#GGGTGTCTCGGGACAA-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | -4.008999 | 16.443800 | Cluster13 | SU010_Total_Pre | 6619 | 0.599940 | GGGTGTCTCGGGACAA-1 | SU010_Total_Pre_1064 | SU010_Total_Pre#GGGTGTCTCGGGACAA-1 | Plasma B |
| SU010_Total_Pre#CCCGTTATCTAGCAGT-1_SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | SU010_Total_Pre | -4.014687 | 1.344591 | Cluster10 | SU010_Total_Pre | 32230 | 0.209680 | CCCGTTATCTAGCAGT-1 | SU010_Total_Pre_1068 | SU010_Total_Pre#CCCGTTATCTAGCAGT-1 | NK1 |
806 rows × 13 columns
We first assign cell types based on KNN and compute KNN-based uncertainty score
[37]:
# 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.
[38]:
# Perform weighted-distance based assignment
distance_threshold = 95
uncertainty_threshold = 0.5
atlas = "HealthyAdult"
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)
We then project query cells into the same UMAP coordinates as original publication for better comparisons, and perform cluster-level assignment.
[39]:
query_distance.obsm["X_umap"] = np.array(query_distance.obs[["UMAP_1", "UMAP_2"]])
query_annotated = scATAnno_assignment.scATAnno_cluster_assign(query_distance, use_rep=use_rep, cluster_col = "Clusters", UMAP=False)
query_annotated
[39]:
AnnData object with n_obs × n_vars = 806 × 344492
obs: 'celltypes', 'tissue', 'dataset', 'UMAP_1', 'UMAP_2', 'Clusters', 'Group', 'depth', 'FRIP', 'Barcodes', 'Internal_Name', 'Group_Barcode', 'true_label', 'uncertainty_score_step1', 'pred_y', 'uncertainty_score_step2', 'Uncertainty_Score', '1.knn-based_celltype', '2.corrected_celltype', 'cluster_annotation'
uns: 'dataset_colors'
obsm: 'X_spectral', 'X_spectral_harmony', 'X_umap', 'kernel_distance', 'distance', 'indices', 'neighbors_labels'
[40]:
sc.pl.umap(query_annotated, color = ['true_label', "cluster_annotation", 'Uncertainty_Score'], title = ["Annotation Published in Satpathy et al.","scATAnno Annotation"], legend_loc = "on data")
... storing 'Clusters' as categorical
... storing 'Group' as categorical
... storing 'Barcodes' as categorical
... storing 'true_label' 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.
Finally, compare the annotation results between scATAnno and original publication.
[41]:
outcomes2 = []
for i in query_annotated.obs["cluster_annotation"].index:
if query_annotated.obs.loc[i,"cluster_annotation"] == query_annotated.obs.loc[i,"true_label"]:
outcomes2.append("Consistent")
elif query_annotated.obs.loc[i,"cluster_annotation"] == "unknown":
outcomes2.append("Unknown")
else:
outcomes2.append("Inconsistent")
query_annotated.obs["outcomes2"] = outcomes2
sc.pl.umap(query_annotated, color = "outcomes2", palette = ['#ff7f0e','#279e68'], title = "Comparison")
... storing 'outcomes2' as categorical
[42]:
adata_query = query_annotated.copy()
variable_col = "true_label"
outcome_col = "outcomes2"
listdict = []
for ct in adata_query.obs[variable_col].unique():
tmp_dict = {}
for pt in np.unique(adata_query.obs[outcome_col]):
tmp_dict[pt] = (
np.sum(
adata_query.obs[adata_query.obs[variable_col] == ct][outcome_col] == pt
)
)
l = len(adata_query.obs[adata_query.obs[variable_col] == ct])
tmp_dict['ct'] = f'{ct} ({l})'
tmp_dict['n_cells'] = len(adata_query.obs[adata_query.obs[variable_col] == ct])
listdict.append(tmp_dict)
df = pd.DataFrame(listdict).set_index('ct').sort_values(by='ct')
df = df.sort_values("n_cells", ascending=True)
fig, ax = plt.subplots(1, 1, figsize=(2, 6))
df[df.columns[0:len(df.columns)-1]].plot(kind='barh', stacked=True, ax=ax, width=1.0, color = ['#ff7f0e','#279e68'])
ax.set_xticklabels(ax.get_xticklabels(), fontsize=15)
ax.set_xlabel('cell numbers')
ax.set_ylabel('')
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0., frameon=False)
ax.grid(False)
sns.despine()
plt.show()
/var/folders/dq/79h3g1bj0rjb4wvj5r5c3w1h0000gn/T/ipykernel_2486/2520740504.py:24: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
ax.set_xticklabels(ax.get_xticklabels(), fontsize=15)
Compute accuracy and F1 scores
[43]:
query = query_annotated.copy()
paper_anno_col="true_label"
cluster_anno_col = "cluster_annotation"
from sklearn.metrics import f1_score, confusion_matrix, classification_report
report = classification_report(y_true=query.obs[paper_anno_col],
y_pred=query.obs[cluster_anno_col],
output_dict=True)
#Also generate a confusion matrix
cm = confusion_matrix(y_true=query.obs[paper_anno_col],
y_pred=query.obs[cluster_anno_col],)
FP = cm.sum(axis=0)-np.diag(cm)
FN = cm.sum(axis=1) - np.diag(cm)
TP = np.diag(cm)
TN = cm.sum() - (FP + FN + TP)
class_spec = TN/(TN+FP)
perfdf = pd.DataFrame(columns = ['class',
'precision',
'recall',
'f1_score',
'support',
'specificity'])
i=0
for key, val in report.items():
# we don't use accuracy so skip it
if key == 'accuracy':
continue
perflist = []
perflist.append(key)
for key2, val2 in val.items():
perflist.append(val2)
if i >= len(class_spec):
perflist.append(0)
else:
perflist.append(class_spec[i])
# Create a single-row DataFrame from perflist and concatenate
new_row = pd.DataFrame([perflist], columns=perfdf.columns)
perfdf = pd.concat([perfdf, new_row], ignore_index=True)
i += 1
/opt/anaconda3/envs/scATAnno/lib/python3.11/site-packages/sklearn/metrics/_classification.py:1565: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.
_warn_prf(average, modifier, f"{metric.capitalize()} is", len(result))
/opt/anaconda3/envs/scATAnno/lib/python3.11/site-packages/sklearn/metrics/_classification.py:1565: UndefinedMetricWarning: Recall is ill-defined and being set to 0.0 in labels with no true samples. Use `zero_division` parameter to control this behavior.
_warn_prf(average, modifier, f"{metric.capitalize()} is", len(result))
/opt/anaconda3/envs/scATAnno/lib/python3.11/site-packages/sklearn/metrics/_classification.py:1565: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.
_warn_prf(average, modifier, f"{metric.capitalize()} is", len(result))
/opt/anaconda3/envs/scATAnno/lib/python3.11/site-packages/sklearn/metrics/_classification.py:1565: UndefinedMetricWarning: Recall is ill-defined and being set to 0.0 in labels with no true samples. Use `zero_division` parameter to control this behavior.
_warn_prf(average, modifier, f"{metric.capitalize()} is", len(result))
/opt/anaconda3/envs/scATAnno/lib/python3.11/site-packages/sklearn/metrics/_classification.py:1565: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.
_warn_prf(average, modifier, f"{metric.capitalize()} is", len(result))
/opt/anaconda3/envs/scATAnno/lib/python3.11/site-packages/sklearn/metrics/_classification.py:1565: UndefinedMetricWarning: Recall is ill-defined and being set to 0.0 in labels with no true samples. Use `zero_division` parameter to control this behavior.
_warn_prf(average, modifier, f"{metric.capitalize()} is", len(result))
[44]:
perfdf
[44]:
| class | precision | recall | f1_score | support | specificity | |
|---|---|---|---|---|---|---|
| 0 | B | 1.000000 | 1.000000 | 1.000000 | 25.0 | 1.000000 |
| 1 | CD8 TEx | 1.000000 | 1.000000 | 1.000000 | 23.0 | 1.000000 |
| 2 | Effector CD8 T | 1.000000 | 1.000000 | 1.000000 | 48.0 | 1.000000 |
| 3 | Memory CD8 T | 1.000000 | 1.000000 | 1.000000 | 147.0 | 1.000000 |
| 4 | Myeloid | 1.000000 | 1.000000 | 1.000000 | 21.0 | 1.000000 |
| 5 | NK1 | 1.000000 | 1.000000 | 1.000000 | 29.0 | 1.000000 |
| 6 | NK2 | 1.000000 | 1.000000 | 1.000000 | 16.0 | 1.000000 |
| 7 | Naive CD4 T | 1.000000 | 1.000000 | 1.000000 | 140.0 | 1.000000 |
| 8 | Naive CD8 T | 1.000000 | 1.000000 | 1.000000 | 51.0 | 1.000000 |
| 9 | Plasma B | 1.000000 | 1.000000 | 1.000000 | 118.0 | 1.000000 |
| 10 | Tfh | 1.000000 | 1.000000 | 1.000000 | 44.0 | 1.000000 |
| 11 | Th1 | 0.000000 | 0.000000 | 0.000000 | 23.0 | 1.000000 |
| 12 | Th17 | 1.000000 | 1.000000 | 1.000000 | 19.0 | 1.000000 |
| 13 | Treg | 1.000000 | 1.000000 | 1.000000 | 102.0 | 1.000000 |
| 14 | unknown | 0.000000 | 0.000000 | 0.000000 | 0.0 | 0.971464 |
| 15 | macro avg | 0.866667 | 0.866667 | 0.866667 | 806.0 | 0.000000 |
| 16 | weighted avg | 0.971464 | 0.971464 | 0.971464 | 806.0 | 0.000000 |
[45]:
# Accuracy
from scATAnno import scATAnno_evaluation
df = query.obs.copy()
cell_accr = {}
for celltype_of_interest in np.unique(df[paper_anno_col]):
true_cells = df[df[paper_anno_col]==celltype_of_interest].index
accr = scATAnno_evaluation.compute_accuracy(df.loc[true_cells,paper_anno_col], df.loc[true_cells,cluster_anno_col])
cell_accr[celltype_of_interest] = accr
cell_level_df = pd.DataFrame([cell_accr], columns=cell_accr.keys())
cell_level_df = cell_level_df.T
cell_level_df.columns = ["accuracy"]
cell_level_df
[45]:
| accuracy | |
|---|---|
| B | 1.0 |
| CD8 TEx | 1.0 |
| Effector CD8 T | 1.0 |
| Memory CD8 T | 1.0 |
| Myeloid | 1.0 |
| NK1 | 1.0 |
| NK2 | 1.0 |
| Naive CD4 T | 1.0 |
| Naive CD8 T | 1.0 |
| Plasma B | 1.0 |
| Tfh | 1.0 |
| Th1 | 0.0 |
| Th17 | 1.0 |
| Treg | 1.0 |
[46]:
# Option 1: Using built-in float
data = np.array([float(cell_level_df.mean()),
float(perfdf[perfdf['class'] == "weighted avg"]["f1_score"])])
# Option 2: Using np.float64 (if you specifically need numpy float type)
data = np.array([np.float64(cell_level_df.mean()),
np.float64(perfdf[perfdf['class'] == "weighted avg"]["f1_score"])])
overallscore = pd.DataFrame(data.T,
index=['Accuracy', 'Weighted F1'],
columns=["value"])
[47]:
overallscore
[47]:
| value | |
|---|---|
| Accuracy | 0.928571 |
| Weighted F1 | 0.971464 |
[48]:
ax = sns.barplot(x = overallscore.index, y= overallscore["value"], capsize=.2, dodge=True, palette="Blues",
)
ax.set_title("Classification Performance")
ax.set_ylim(0,1)
ax.set_ylabel("")
fig = ax.get_figure()
fig.tight_layout()
plt.show()
