Tutorial 3: Use ADEPT on customized ST dataset

packages loading in

import warnings
import pandas as pd
import numpy as np
import scanpy as sc
from scipy import sparse
import os
from imputation.impute import impute_
import GAAE
import argparse
import matplotlib.pyplot as plt
plt.rcParams['pdf.fonttype'] = 42
plt.rcParams['ps.fonttype'] = 42
import time
import seaborn as sns
from GAAE.utils import impute, DE_num_calc

build up your own loader function

Below is an example for loading a customized dataset. You could provide cell annotation to it optionally.

def load_customized(data_path):
    """this is an example for reading .h5 file, if you have different saving format, please refer to scanpy or implement your own function to read in"""
    ad = sc.read_visium(path=data_path, count_file='_filtered_feature_bc_matrix.h5', load_images=True)
    ad.var_names_make_unique()

    """comment out this if you dont have gt/cell annotation. This is an example for reading gt which is stored in .txt file"""
    gt_dir = os.path.join("./", "./", 'gt')
    gt_df = pd.read_csv(os.path.join(gt_dir, 'tissue_positions_list_GTs.txt'), sep=',', header=None, index_col=0)
    ad.obs['original_clusters'] = gt_df.loc[:, 6]
    keep_bcs = ad.obs.dropna().index
    ad = ad[keep_bcs].copy()
    ad.obs['original_clusters'] = ad.obs['original_clusters'].astype(int).astype(str)

    return ad

To make the tutorial runnable, I will still use the DLPFC dataset here

def load_customized(data_path):
    # 151507, ..., 151676 12 in total
    ad = sc.read_visium(path="/home/yunfei/spatial_benchmarking/benchmarking_data/DLPFC12/151673", count_file='151673_filtered_feature_bc_matrix.h5', load_images=True)
    ad.var_names_make_unique()

    return ad

def filter_num_calc(args_, comp_):
    if comp_ is not None:
        return comp_
    print("optimizing minimum filter number")
    # input_dir = os.path.join(args_.data_dir, args_.input_data)
    data_name = args_.input_data
    data_path = args_.data_dir
    if data_name:
        adata = load_customized(data_path)
    else:
        print("exception in dataname")
        exit(-1)
    adata.var_names_make_unique()

    for temp_count in range(5, 150):
        sc.pp.filter_genes(adata, min_counts=temp_count)
        try:
            X = adata.X.todense()
        except:
            X = adata.X
        if np.count_nonzero(X)/(X.shape[0]*X.shape[1]) > args_.filter_nzr:
            return temp_count
    return 150

def initialize(args_, gene_min_count):
    print("initializing spatial transcriptomic data")
    data_path = args_.data_dir
    data_name = args_.input_data
    if data_name:
        adata_ = load_customized(data_path)
        adata_ori_ = adata_
        if args_.use_preprocessing:
            # Normalization
            sc.pp.filter_genes(adata_, min_counts=gene_min_count)
            if args_.use_hvgs != 0:
                sc.pp.highly_variable_genes(adata_, flavor="seurat_v3", n_top_genes=args_.use_hvgs)
            sc.pp.normalize_total(adata_, target_sum=1e4)
            sc.pp.log1p(adata_)
        else:
            sc.pp.filter_genes(adata_, min_counts=gene_min_count)
    else:
        print("exception in data name")
        exit(-1)
    return adata_, adata_ori_

argument parser: change the necessary args to your customized value by setting args.data_dir=’/path/to/your/input/data/’

parser = argparse.ArgumentParser()
parser.add_argument('--data_dir', type=str, default="./", help="root dir for input data")
parser.add_argument('--gt_dir', type=str, default="./", help="root dir for data ground truth")
parser.add_argument('--input_data', type=str, default="151673", help="input data section id")
parser.add_argument('--impute_cluster_num', type=str, default="7", help="diff cluster numbers for imputation")
parser.add_argument('--cluster_num', type=int, default=7, help="input data cluster number")
parser.add_argument('--radius', type=int, default=150, help="input data radius")
parser.add_argument('--no_de', type=int, default=0, help='switch on/off DEG selection module')
parser.add_argument("--use_mean", type=int, default=0, help="use mean value in de list or not")
parser.add_argument("--impute_runs", type=int, default=1, help="time of runs for imputation")
parser.add_argument("--runs", type=int, default=2, help="total runs for the data")
parser.add_argument('--use_hvgs', type=int, default=3000, help="select highly variable genes before training")
parser.add_argument('--use_preprocessing', type=int, default=1, help='use preprocessed input or raw input')
parser.add_argument('--save_fig', type=int, default=1, help='saving output visualization')
parser.add_argument('--de_nzr_min', type=float, default=0.299, help='suggested min nzr threshold after de selection')
parser.add_argument('--de_nzr_max', type=float, default=0.399, help='suggested max nzr threshold after de selection')
parser.add_argument('--use_gpu_id', type=str, default='0', help='use which GPU, only applies when you have multiple gpu')
args, unknown = parser.parse_known_args()
args.impute_cluster_num = args.impute_cluster_num.split(",")  # ["5", "6", "7"]
args.filter_nzr = 0.15
args.filter_num = 5
root_d = args.data_dir

preprocessing to get the parameters for initialization

filter_num = filter_num_calc(args, args.filter_num)

adata, adata_ori = initialize(args, filter_num)

if os.path.exists('./cache/' + args.data_dir.split('/')[-2] + args.input_data + '.txt'):
    with open('./cache/' + args.data_dir.split('/')[-2] + '.txt', 'r') as fp:
        line = fp.readlines()[0]
        split_ = line.strip().split(",")
        de_top_k_list = []
        for e in split_:
            de_top_k_list.append(int(e))
    print("previously cached de list = ", de_top_k_list)
else:
    de_top_k_list = DE_num_calc(args, adata)
    print("optimized de list = ", de_top_k_list)
    with open('./cache/' + args.data_dir.split('/')[-2] + '.txt', 'a+') as fp:
        # fp.write('de list: ')
        fp.write(','.join([str(i) for i in de_top_k_list]))

initializing spatial transcriptomic data

/home/yunfei/anaconda3/envs/test/lib/python3.7/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
  utils.warn_names_duplicates("var")

optimizing top DEs before imputation
original non-zero =  0.12213016680918756
(3639, 18032)
DE topk =  50
section id =  151673
------Calculating spatial graph...
The graph contains 21124 edges, 3639 cells.
5.8049 neighbors per cell on average.
Size of Input:  (3639, 3000)
GAAE(
  (conv1): GATConv(3000, 512, heads=1)
  (conv2): GATConv(512, 30, heads=1)
  (conv3): GATConv(30, 512, heads=1)
  (conv4): GATConv(512, 3000, heads=1)
)

100%|██████████| 1000/1000 [00:11<00:00, 85.23it/s]
R[write to console]:                    __           __
   ____ ___  _____/ /_  _______/ /_
  / __ `__ \/ ___/ / / / / ___/ __/
 / / / / / / /__/ / /_/ (__  ) /_
/_/ /_/ /_/\___/_/\__,_/____/\__/   version 5.4.10
Type 'citation("mclust")' for citing this R package in publications.

fitting ...
  |======================================================================| 100%

From cffi callback <function _processevents at 0x7fd6f558a830>:
Traceback (most recent call last):
  File "/home/yunfei/anaconda3/envs/test/lib/python3.7/site-packages/rpy2/rinterface_lib/callbacks.py", line 275, in _processevents
    @ffi_proxy.callback(ffi_proxy._processevents_def,
  File "/home/yunfei/anaconda3/envs/test/lib/python3.7/site-packages/rpy2/rinterface.py", line 92, in _sigint_handler
    raise KeyboardInterrupt()
KeyboardInterrupt

selecting DEs after first round of clustering ...

/home/yunfei/Spatial/ADEPT/GAAE/Train_ADEPT_modified.py:508: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
  gene_list = adata.var.index[:, np.newaxis]

(18032, 1182) (18032, 2457)
label =  1

initial clustering to impute data

de_list_epoch = []
if de_top_k_list != []:
    print("performing DEGs selection")
    adata_list = []
    for de_ in de_top_k_list:
        for cluster_n in args.impute_cluster_num:
            print("cluster_n = ", cluster_n)
            GAAE.get_kNN(adata, rad_cutoff=args.radius)

            ari_ini, ari_final, de_list, adata_out = GAAE.train_ADEPT_use_DE(adata, n_epochs=1000,
                                                                        num_cluster=int(cluster_n),
                                                                        dif_k=de_, device_id=args.use_gpu_id)
            de_list_epoch.append(de_list)
            adata_list.append(adata_out)
    g_union = set.union(*de_list_epoch)
    imputed_ad = impute(args, adata_list, g_union, de_top_k_list)
else:
    print("skip performing DEGs selection")
    imputed_ad = adata

final clustering

"""result of imputed data"""
GAAE.get_kNN(imputed_ad, rad_cutoff=args.radius)
ari_ini, ARI, de_list, adata_out = GAAE.train_ADEPT_use_DE(imputed_ad, n_epochs=1000, num_cluster=args.cluster_num, device_id=args.use_gpu_id)

aris = []
print('Dataset:', args.input_data)
print('ARI:', ARI)
aris.append(ARI)
print(aris)