Visium HD 空间转录组分析探索之--基础分析

import osimport numpy as npimport pandas as pdimport scanpy as scimport seaborn as snsimport omicverse as ovimport matplotlib.pyplot as pltfrom matplotlib.pyplot import rc_context
import warningswarnings.filterwarnings('ignore')
plt.style.use('default')plt.rcParams['figure.facecolor'] = 'white'sc.settings.set_figure_params(dpi=100, dpi_save=200, figsize=(5, 5), facecolor='white')

data_dir = './CRC/P1_CRC/'tissue_position_df = pd.read_parquet(f'{data_dir}/binned_outputs/square_008um/spatial/tissue_positions.parquet')tissue_position_df.to_csv(f'{data_dir}/binned_outputs/square_008um/spatial/tissue_positions_list.csv', index=False, header=None)adata = sc.read_visium(f'{data_dir}/binned_outputs/square_008um', library_id='P1')adata.obs['sample'] = 'P1'adata.var_names_make_unique()

# Data QCadata.var['mt'] = adata.var_names.str.upper().str.startswith('MT-') adata.var['ribo'] = adata.var_names.str.upper().str.startswith(("RPS","RPL"))sc.pp.calculate_qc_metrics(adata, qc_vars=['mt', 'ribo'], percent_top=None, log1p=False, inplace=True)adata.obs['log10GenesPerUMI'] = np.log10(adata.obs['n_genes_by_counts']) / np.log10(adata.obs['total_counts'])sc.pl.violin(    adata,    ['n_genes_by_counts', 'total_counts', 'pct_counts_mt', 'log10GenesPerUMI'],     groupby='sample',     jitter=False,     rotation=90,     multi_panel=True,     show=False)

## Cell filtersc.pp.filter_cells(adata, min_counts=50)adata = adata[(adata.obs['pct_counts_mt'] < 30) & (adata.obs['log10GenesPerUMI'] > 0.8) & (adata.obs['log10GenesPerUMI'] < 0.99)]sc.pl.violin(    adata,    ['n_genes_by_counts', 'total_counts', 'pct_counts_mt', 'log10GenesPerUMI'],     groupby='sample',     jitter=False,     rotation=90,     multi_panel=True,     show=False)
with rc_context({'figure.figsize': (8, 8)}):    sc.pl.spatial(        adata,         color=['n_genes_by_counts', 'total_counts', 'log10GenesPerUMI'],        library_id='P1',        alpha_img=0.6    )

## data normalization and log10 transformadata.layers["counts"] = adata.X.copy()sc.pp.normalize_total(adata)sc.pp.log1p(adata)adata.raw = adata ## select highly variable gene and data scale sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5, n_top_genes=2000)sc.pl.highly_variable_genes(adata, show=False)adata = adata[:, adata.var['highly_variable']]sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])sc.pp.scale(adata, max_value=10) ## PCAsc.tl.pca(adata, svd_solver='arpack')sc.pl.pca(adata, color='sample', show=False)sc.pl.pca_variance_ratio(adata, log=True, show=False) ## UMAPsc.pp.neighbors(adata, n_neighbors=15, n_pcs=15)sc.tl.umap(adata)sc.tl.leiden(adata, resolution=0.6, key_added='res_0.6')sc.tl.leiden(adata, resolution=0.8, key_added='res_0.8')

from matplotlib import patheffectsimport matplotlib.pyplot as pltfig, ax = plt.subplots(figsize=(6,6)) ov.pl.embedding(adata,                basis='X_umap',                color=['res_0.6'],                size=4,                show=False,                 legend_loc=None,                 add_outline=False,                 frameon='small',                legend_fontoutline=2,                ax=ax                ) ov.utils.gen_mpl_labels(    adata,    'res_0.6',    exclude=("None",),      basis='X_umap',    ax=ax,    adjust_kwargs=dict(arrowprops=dict(arrowstyle='-', color='black')),    text_kwargs=dict(fontsize= 12 ,weight='bold', path_effects=[patheffects.withStroke(linewidth=2, foreground='w')] ),)

sc.tl.rank_genes_groups(    adata,     groupby='res_0.6',     use_raw=True,     method="wilcoxon",    pts=True) sc.pl.rank_genes_groups_dotplot(    adata,     groupby='res_0.6',    standard_scale="var",     n_genes=5,    cmap='bwr',    show=False,) result = adata.uns['rank_genes_groups']groups = result['names'].dtype.namesall_diff = []for c in groups:    tmp = sc.get.rank_genes_groups_df(adata, group=c)    tmp = tmp[(tmp['logfoldchanges']>0) & (tmp['pvals_adj']<0.05)]    tmp['cluster'] = c    all_diff.append(tmp) all_diff = pd.concat(all_diff, axis=0)