载入包

import numpy as np
import pandas as pd
import scanpy as sc

sc.settings.verbosity = 3             # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=80, facecolor='white')

results_file = 'write/pbmc3k.h5ad'  # the file that will store the analysis results

两种方式读取数据

直接读取.h5ad文件

adata = sc.read('pancreas.h5ad')

读取10X测序文件

adata = sc.read_10x_mtx(
    'data/filtered_gene_bc_matrices/hg19/',  # the directory with the `.mtx` file
    var_names='gene_symbols',                # use gene symbols for the variable names (variables-axis index)
    cache=True)                              # write a cache file for faster subsequent reading

预处理、质控

adata.var_names_make_unique()  # this is unnecessary if using `var_names='gene_ids'` in `sc.read_10x_mtx`
sc.pl.highest_expr_genes(adata, n_top=20, )
# Basic filtering
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
# Filter mitochondrial genes
adata.var['mt'] = adata.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
# genes，counts分布
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True)
# Scatter plot
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')
# Filtering by 
adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :]

标准流程

# Normalize
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# Identify HVG
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.highly_variable_genes(adata)

# 储存counts到adata.raw
adata.raw = adata
adata = adata[:, adata.var.highly_variable]

# regress_out total_counts & pct_counts_mt
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])

# Scale
sc.pp.scale(adata, max_value=10)

# PCA
sc.pl.pca_variance_ratio(adata, log=True)
adata.write(results_file)

# Computing the neighborhood graph
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

# UMAP
sc.tl.umap(adata)

# Unsupervised clustering
sc.tl.leiden(adata)

# Save the result
adata.write(results_file)

寻找差异基因

# each one vs all
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

# Show the 10 top ranked genes per cluster 0, 1, …, 7 in a dataframe
pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(5)

# Get a table with the scores and groups
result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
pd.DataFrame(
    {group + '_' + key[:1]: result[key][group]
    for group in groups for key in ['names', 'pvals']}).head(5)

# 0 vs 1
sc.tl.rank_genes_groups(adata, 'leiden', groups=['0'], reference='1', method='wilcoxon')
sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=20)

作图函数

# FeaturePlot
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])

# DimPlot
# 当adata.obsm的keys为'umap'时：
sc.pl.umap(adata, color='leiden', legend_loc='on data')
# 当adata.obsm的keys为其他时：（这里为'UMAP'）
sc.pl.embedding(adata_atac, basis='UMAP', color='Clusters')


# VlnPlot
sc.pl.violin(adata, ['CST3', 'NKG7', 'PPBP'], groupby='leiden')