add StarDist+RCTD and 6 benchmarking datasets

README.md

@@ -29,7 +29,12 @@ We provide source codes for reproducing the SpatialScope analysis in the main te
 
 All relevent materials involved in the reproducing codes are availabel from [here](https://drive.google.com/drive/folders/1PXv_brtr-tXshBVEd_HSPIagjX9oF7Kg?usp=sharing)
 
-+ [Benchmarking](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Benchmarking-MERFISH.ipynb)
++ [Benchmarking Dataset 1](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Benchmarking-Dataset_1.ipynb)
++ [Benchmarking Dataset 2](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Benchmarking-Dataset_2.ipynb)
++ [Benchmarking Dataset 3](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Benchmarking-Dataset_3.ipynb)
++ [Benchmarking Dataset 4](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Benchmarking-Dataset_4.ipynb)
++ [Benchmarking Dataset 5](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Benchmarking-Dataset_5.ipynb)
++ [Benchmarking Dataset 6](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Benchmarking-Dataset_6.ipynb)
 + [Human Heart (Visium, a single slice)](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Human-Heart.ipynb)
 + [Mouse Brain (Visium, 3D alignment of multiple slices)](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Mouse-Brain.ipynb)
 + [Mouse Cerebellum (Slideseq-V2)](https://github.com/YangLabHKUST/SpatialScope/blob/master/demos/Mouse-Cerebellum-Slideseq.ipynb)

compared_methods/SDRCTD.py

@@ -0,0 +1,205 @@
+import anndata
+import numpy as np
+import pandas as pd
+import sys
+import pickle
+import os
+import copy
+import argparse
+from sklearn.model_selection import KFold
+from sklearn.metrics import mean_squared_error
+import pandas as pd
+import matplotlib.pyplot as plt
+import scanpy as sc
+import warnings
+warnings.filterwarnings('ignore')
+import seaborn as sns
+
+from SDRCTD_utils import *
+
+
+
+
+class SDRCTD:
+ def __init__(self,tissue,out_dir,RCTD_results_dir,RCTD_results_name, ST_Data, SC_Data, cell_class_column = 'cell_type', hs_ST = True):
+ self.tissue = tissue
+ self.out_dir = out_dir 
+ self.RCTD_results_dir = RCTD_results_dir 
+ self.RCTD_results_name = RCTD_results_name 
+ self.ST_Data = ST_Data
+ self.SC_Data = SC_Data
+ self.cell_class_column = cell_class_column
+ self.hs_ST = hs_ST
+
+ if not os.path.exists(out_dir):
+ os.mkdir(out_dir)
+ if not os.path.exists(os.path.join(out_dir,tissue)):
+ os.mkdir(os.path.join(out_dir,tissue))
+
+ self.out_dir = os.path.join(out_dir,tissue)
+ loggings = configure_logging(os.path.join(self.out_dir,'logs'))
+ self.loggings = loggings 
+
+ self.LoadRCTDresults()
+ if SC_Data is not None:
+ self.LoadSCData()
+
+
+ def LoadRCTDresults(self):
+ with open(os.path.join(self.RCTD_results_dir, self.RCTD_results_name + '.pickle'), 'rb') as handle:
+ RCTD_results = pickle.load(handle)
+
+ if self.hs_ST:
+ try:
+ weights = RCTD_results['results']['weights']
+ except:
+ weights = RCTD_results['results']
+ else:
+ weights = RCTD_results['results']
+
+
+ self.weights = (weights / np.array(weights.sum(1))[:, None])
+
+ def single_cell_type_assignment(self, cell_num_column = 'cell_count', VisiumCellsPlot = True):
+ seged_sp_adata = sc.read(self.ST_Data) #ST_Data already complete nuclei segmentation with StarDist. 'cell_locations' already in uns
+
+ mat = self.weights.values
+ cell_nums = np.array(seged_sp_adata.obs[cell_num_column])
+
+ cell_counts = distribute_cells(mat, cell_nums)
+ cell_types = self.weights.columns
+ cell_type_list = assign_cell_type(cell_counts, cell_types)
+ seged_sp_adata.uns['cell_locations']['SDRCTD_cell_type'] = cell_type_list
+
+ self.cell_type_list = cell_type_list
+ self.seged_sp_adata = seged_sp_adata
+
+ seged_sp_adata.uns['RCTD_weights'] = self.weights
+ seged_sp_adata.write(os.path.join(self.out_dir, 'single_cell_type_label_bySDRCTD.h5ad'))
+
+ # plot results
+ if self.hs_ST or not VisiumCellsPlot:
+ fig, ax = plt.subplots(figsize=(10,8.5),dpi=100)
+ sns.scatterplot(data=seged_sp_adata.uns['cell_locations'], x="x",y="y",s=10,hue='SDRCTD_cell_type',palette='tab20',legend=True)
+ plt.axis('off')
+ plt.legend(bbox_to_anchor=(0.97, .98),framealpha=0)
+ plt.savefig(os.path.join(self.out_dir, 'SDRCTD_estemated_ct_label.png'))
+ plt.close()
+
+ elif VisiumCellsPlot:
+ if seged_sp_adata.obsm['spatial'].shape[1] == 2:
+ fig, ax = plt.subplots(1,1,figsize=(14, 8),dpi=200)
+ PlotVisiumCells(seged_sp_adata,"SDRCTD_cell_type",size=0.4,alpha_img=0.4,lw=0.4,palette='tab20',ax=ax)
+ plt.savefig(os.path.join(self.out_dir, 'SDRCTD_estemated_ct_label.png'))
+ plt.close()
+
+ def cell_type_mean_assignment(self):
+ # cell type mean as decomposed cell gene expression
+ ref_df = pd.DataFrame([[ct, i]for i, ct in enumerate(self.sc_data_process_marker.obs[self.cell_class_column].astype('category').cat.categories)], columns = ['cell_type', 'cell_type_code'])
+ ref_df.index = ref_df.cell_type
+ ref_df = ref_df.iloc[:,1:]
+
+ x_decom = self.mu[ref_df.loc[np.array(self.cell_type_list)].cell_type_code.tolist()]
+ x_decom_adata = anndata.AnnData(X = x_decom.copy(), obs = self.seged_sp_adata.uns['cell_locations'].copy(), var = self.sc_data_process_marker.var)
+ x_decom_adata.write(os.path.join(self.out_dir, 'cell_type_mean_bySDRCTD.h5ad'))
+
+
+ def LoadSCData(self):
+ # load sc data
+ sc_data_process = anndata.read_h5ad(self.SC_Data)
+ if 'Marker' in sc_data_process.var.columns:
+ sc_data_process_marker = sc_data_process[:,sc_data_process.var['Marker']]
+ else:
+ sc_data_process_marker = sc_data_process
+
+ if sc_data_process_marker.X.max() <= 30:
+ self.loggings.info(f'Maximum value: {sc_data_process_marker.X.max()}, need to run exp')
+ try:
+ sc_data_process_marker.X = np.exp(sc_data_process_marker.X) - 1
+ except:
+ sc_data_process_marker.X = np.exp(sc_data_process_marker.X.toarray()) - 1
+
+
+ cell_type_array = np.array(sc_data_process_marker.obs[self.cell_class_column])
+ cell_type_class = np.unique(cell_type_array)
+ df_category = sc_data_process_marker.obs[[self.cell_class_column]].astype('category').apply(lambda x: x.cat.codes)
+
+ # parameters: mean and cell type index
+ cell_type_array_code = np.array(df_category[self.cell_class_column]) 
+ try:
+ data = sc_data_process_marker.X.toarray()
+ except:
+ data = sc_data_process_marker.X
+
+ n, d = data.shape
+ q = cell_type_class.shape[0]
+ self.loggings.info(f'scRNA-seq data shape: {data.shape}')
+ self.loggings.info(f'scRNA-seq cell class number: {q}')
+
+ mu = np.zeros((q, d))
+ for k in range(q):
+ mu[k] = data[cell_type_array_code == k].mean(0).squeeze()
+ self.mu = mu
+ self.sc_data_process_marker = sc_data_process_marker
+
+
+
+
+
+
+
+if __name__ == "__main__":
+ HEADER = """
+ <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
+ <> 
+ <> StarDist + RCTD
+ <> Version: %s
+ <> MIT License
+ <>
+ <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
+ <> Software-related correspondence: %s or %s
+ <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
+ <> Visium data example:
+ python <install path>/src/Cell_Type_Identification.py \\
+ --cell_class_column cell_type \\
+ --tissue heart \\
+ --out_dir ./output \\
+ --ST_Data ./output/heart/sp_adata_ns.h5ad \\
+ --SC_Data ./Ckpts_scRefs/Heart_D2/Ref_Heart_sanger_D2.h5ad 
+ <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> 
+ """ 
+ def str2bool(v):
+ if isinstance(v, bool):
+ return v
+ if v.lower() in ("yes", "true", "t", "y", "1"):
+ return True
+ elif v.lower() in ("no", "false", "f", "n", "0"):
+ return False
+ else:
+ raise argparse.ArgumentTypeError("Boolean value expected.")
+
+ parser = argparse.ArgumentParser(description='simulation sour_sep')
+ parser.add_argument('--out_dir', type=str, help='output path', default=None)
+ parser.add_argument('--RCTD_results_dir', type=str, help='RCTD results path', default=None)
+ parser.add_argument('--RCTD_results_name', type=str, help='RCTD results file\'s name', default='InitProp')
+ parser.add_argument('--ST_Data', type=str, help='ST data path', default=None)
+ parser.add_argument('--SC_Data', type=str, help='single cell reference data path', default=None)
+ parser.add_argument('--cell_class_column', type=str, help='input cell class label column in scRef file', default = 'cell_type')
+ parser.add_argument('--cell_num_column', type=str, help='cell number column in spatial file', default = 'cell_count')
+ parser.add_argument('--hs_ST', action="store_true", help='high resolution ST data such as Slideseq, DBiT-seq, and HDST, MERFISH etc.')
+ parser.add_argument("--VisiumCellsPlot", type=str2bool, const=True, default=True, nargs="?", help="whether to plot in VisiumCells mode or just scatter plot")
+ args = parser.parse_args()
+
+ args.tissue = 'SDRCTD_results'
+ if not os.path.exists(args.out_dir):
+ os.mkdir(args.out_dir)
+ if not os.path.exists(os.path.join(args.out_dir,args.tissue)):
+ os.mkdir(os.path.join(args.out_dir,args.tissue))
+ if args.RCTD_results_dir is None:
+ args.RCTD_results_dir = args.out_dir
+
+
+ sdr = SDRCTD(args.tissue,args.out_dir, args.RCTD_results_dir, args.RCTD_results_name, args.ST_Data, args.SC_Data, cell_class_column = args.cell_class_column, hs_ST = args.hs_ST)
+ sdr.single_cell_type_assignment(cell_num_column = args.cell_num_column, VisiumCellsPlot = args.VisiumCellsPlot)
+ if args.SC_Data is not None:
+ sdr.cell_type_mean_assignment()

compared_methods/SDRCTD_utils.py

@@ -0,0 +1,139 @@
+import os
+import scanpy as sc
+# import squidpy as sq
+import numpy as np
+import pandas as pd
+import pathlib
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+# import skimage
+import seaborn as sns
+from itertools import chain
+# from stardist.models import StarDist2D
+from csbdeep.utils import normalize
+from anndata import AnnData
+from scipy.spatial.distance import pdist
+import logging
+import sys
+from sklearn.metrics.pairwise import cosine_similarity
+
+def PlotVisiumCells(adata,annotation_list,size=0.8,alpha_img=0.3,lw=1,subset=None,palette='tab20',show_circle = True, legend = True, ax=None,**kwargs):
+ merged_df = adata.uns['cell_locations'].copy()
+ test = sc.AnnData(np.zeros(merged_df.shape), obs=merged_df)
+ test.obsm['spatial'] = merged_df[["x", "y"]].to_numpy().copy()
+ test.uns = adata.uns
+
+ if subset is not None:
+ #test = test[test.obs[annotation_list].isin(subset)]
+ test.obs.loc[~test.obs[annotation_list].isin(subset),annotation_list] = None
+
+ sc.pl.spatial(
+ test,
+ color=annotation_list,
+ size=size,
+ frameon=False,
+ alpha_img=alpha_img,
+ show=False,
+ palette=palette,
+ na_in_legend=False,
+ ax=ax,title='',sort_order=True,**kwargs
+ )
+ if show_circle:
+ sf = adata.uns['spatial'][list(adata.uns['spatial'].keys())[0]]['scalefactors']['tissue_hires_scalef']
+ spot_radius = adata.uns['spatial'][list(adata.uns['spatial'].keys())[0]]['scalefactors']['spot_diameter_fullres']/2
+ for sloc in adata.obsm['spatial']:
+ rect = mpl.patches.Circle(
+ (sloc[0] * sf, sloc[1] * sf),
+ spot_radius * sf,
+ ec="grey",
+ lw=lw,
+ fill=False
+ )
+ ax.add_patch(rect)
+ ax.axes.xaxis.label.set_visible(False)
+ ax.axes.yaxis.label.set_visible(False)
+
+ if not legend:
+ ax.get_legend().remove()
+
+ # make frame visible
+ for _, spine in ax.spines.items():
+ spine.set_visible(True)
+
+
+
+def assign_cell_type(cell_counts, cell_types):
+ cell_type_list = []
+ for i in range(cell_counts.shape[0]):
+ cell_count = cell_counts[i]
+ idx = np.where(cell_count > 0)[0]
+ cell_type_list_row = [[cell_types[idx][_]] * cell_count[idx[_]] for _ in range(idx.shape[0])]
+ cell_type_list_row = np.array([item for sublist in cell_type_list_row for item in sublist])
+ np.random.shuffle(cell_type_list_row)
+ cell_type_list = cell_type_list + list(cell_type_list_row)
+
+ return cell_type_list
+
+def distribute_cells(mat, cell_nums): # mat: spots * cell_type; cell_nums: spots * 1
+ cell_nums = cell_nums.astype(int)
+ cell_nums_original = cell_nums.copy()
+
+ mat[np.absolute(mat) < 1e-3] = 0
+ cell_counts = np.zeros(mat.shape).astype(int)
+ assert not np.any(cell_nums < 0)
+ assert not np.any(mat < 0)
+
+ mat = mat * cell_nums[:, None]
+ cell_num_dist = np.floor(mat).astype(int)
+ cell_counts = cell_counts + cell_num_dist
+ cell_nums_remain = cell_nums - cell_num_dist.sum(1)
+ mat_remain = mat - cell_num_dist
+
+ assert not np.any(cell_nums_remain < 0)
+ assert not np.any(mat_remain < 0)
+
+ mat = mat_remain
+ cell_nums = cell_nums_remain
+
+ while(np.any(cell_nums_remain > 0)):
+ mat[mat.argsort()[:, ::-1].argsort() >= cell_nums[:, None]] = 0
+ mat = np.divide(mat, mat.sum(1)[:,None], out=np.zeros_like(mat), where=mat.sum(1)[:,None]!=0)
+ mat = mat * cell_nums[:, None]
+ cell_num_dist = np.floor(mat).astype(int)
+ cell_counts = cell_counts + cell_num_dist
+ cell_nums_remain = cell_nums - cell_num_dist.sum(1)
+ mat_remain = mat - cell_num_dist
+
+ assert not np.any(cell_nums_remain < 0)
+ assert not np.any(mat_remain < 0)
+
+ mat = mat_remain
+ cell_nums = cell_nums_remain
+
+ assert np.array_equal(cell_counts.sum(1), cell_nums_original)
+
+ return cell_counts
+
+
+
+
+def configure_logging(logger_name):
+ LOG_LEVEL = logging.DEBUG
+ log_filename = logger_name+'.log'
+ importer_logger = logging.getLogger('importer_logger')
+ importer_logger.setLevel(LOG_LEVEL)
+ formatter = logging.Formatter('%(asctime)s : %(levelname)s : %(message)s')
+
+ fh = logging.FileHandler(filename=log_filename)
+ fh.setLevel(LOG_LEVEL)
+ fh.setFormatter(formatter)
+ importer_logger.addHandler(fh)
+
+ sh = logging.StreamHandler(sys.stdout)
+ sh.setLevel(LOG_LEVEL)
+ sh.setFormatter(formatter)
+ importer_logger.addHandler(sh)
+ return importer_logger
+
+
+