前面，我们生信技能树的讲师小洁老师与萌老师新开了一个学习班：《掌握Python，解锁单细胞数据的无限可能》，身为技能树的一员，近水楼台先得月，学起！下面是我的学习笔记，希望可以给你带来一点参考

前面的python学习笔记：：

今天继续学习视频：python_day6 ！一口气学完吧！

touch day6.ipynb

课前准备操作到 23:52

本次课程需要用到的模块，提前安装好：

永久镜像设置：

#永久设置镜像
pip config set global.index-url https://pypi.mirrors.ustc.edu.cn/simple/
#升级pip
pip install pip --upgrade

或者临时的镜像使用，我更偏好这种好像：

# bash终端
conda activate sc
# 安装 单细胞分析需要的包
pip install scanpy python-igraph leidenalg -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple

# 安装jupter汉化的包：ipywidgets jupyterlab-language-pack-zh-CN
pip install ipywidgets jupyterlab-language-pack-zh-CN -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple

scanpy单细胞分析 01:40:04

代码参考自scanpy的标准流程：

https://scanpy.readthedocs.io/en/stable/tutorials/basics/clustering-2017.html
https://scanpy.readthedocs.io/en/stable/tutorials/basics/clustering.html#manual-cell-type-annotation

1.数据和包准备

数据来自著名的pbmc3k数据：

https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz

加载模块：

import pandas as pd
import scanpy as sc
import numpy as np
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")

2.读取数据

day6/01_data文件夹下面是标准的10x上游cellranger的输出结果：barcodes.tsv genes.tsv matrix.mtx

adata = sc.read_10x_mtx(
    "day6/01_data",  # 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
)
# this is unnecessary if using `var_names='gene_ids'` in `sc.read_10x_mtx`
adata.var_names_make_unique() 
adata

原始细胞数为2700个细胞：

简单看一下数据的指标：

# 表达矩阵里的数值范围
np.min(adata.X), np.max(adata.X)

# 基本过滤
# 过滤前 的细胞数与基因数
adata.X.shape

过滤：

# 过滤细胞：每个细胞至少表达200个基因
sc.pp.filter_cells(adata, min_genes=200)
# 过滤基因：每个基因至少在3个细胞中表达
sc.pp.filter_genes(adata, min_cells=3)

# 过滤后：可以看到过滤前后的细胞数量和基因数量
adata.X.shape

基因数过滤的比较多：

查看感兴趣的基因的表达矩阵

稀疏矩阵不支持直接查看，只能是转换成矩阵或者数据框才能查看。转换成矩阵就丢失了行名列名，转换成数据框更好。

# 转换成矩阵
adata[0:6, ['CD3D','TCL1A','MS4A1']].X.toarray()

# 转换成数据框
adata[0:6, ['CD3D','TCL1A'


    
,'MS4A1']].to_df()

3.质控

前面进行了简单初步的过滤，这里过滤低质量的细胞：包括双细胞，空包体，死细胞，红细胞等。

过滤低质量的细胞，常见的指标大家已经非常熟悉了，比如每个细胞中表达的基因数，count数，线粒体基因表达百分比，红细胞基因比例等

3.1 查看三个指标

在anndata对象中，基因的注释是在adata.var里，细胞的注释是在adata.obs里。

# 每个基因在多少个细胞中表达
adata.var.head()

# 每个细胞中表达多少个基因
adata.obs.head()

用calculate_qc_metrics来计算需要的过滤指标（"n_genes_by_counts", "total_counts", "pct_counts_mt"）

# 质控
# annotate the group of mitochondrial genes as "mt"
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)
adata.obs.head()

每个细胞中的指标如下：

绘制这些指标的小提琴图：

sc.pl.violin(
    adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4,
    # **{"color": "#f8766d"},  # 使用额外的参数来设置颜色
    multi_panel=True
)

小提琴图：免疫细胞的基因数以及count数都在一个比较低的水平，相对于其他上皮类，基质类细胞，癌细胞。

3.2 绘制指标之间的相关性

total_counts 与 pct_counts_mt的相关性：

sc.pl.scatter(adata, x="total_counts", y="pct_counts_mt")

total_counts 与 n_genes_by_counts的相关性：

sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts")

3.3 过滤

小的知识点：

链式赋值：如果你在一段代码中连续对数据进行多次筛选或其他操作，可以在最后一步使用 .copy() 来确保最终结果是一个独立的副本。这样可以避免在中间步骤中不小心修改原始数据。

adata = adata[adata.obs.n_genes_by_counts 2500, :]
adata = adata[adata.obs.pct_counts_mt 5, :].copy()
adata.shape
# (2638, 13714)

过滤后，细胞数变成了2638个。

4.找高变基因(HVG)

挑选高变基因的指标：离散度-是基因表达方差与基因表达平均水平的比值，用于评估基因表达的变异程度。

# 首先将数据矩阵归一化（校正文库大小）：
sc.pp.normalize_total(adata, target_sum=1e4)

# 对数据进行log
sc.pp.log1p(adata)

# 高变基因
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

# 查看每个基因的指标
adata.var.head()

查看前10个高变化基因(因为scanpy和seurat的高变化基因默认参数不同，所以找出的基因也不相同)

sc.pl.highly_variable_genes(adata)

# adata.var.loc[adata.var.highly_variable == True]
he = adata.var.sort_values(by='dispersions_norm', ascending=False)

print(he[he.highly_variable == True].index.tolist()[0:10])
# ['DOK3', 'ARVCF', 'YPEL2', 'UBE2D4', 'FAM210B', 'CTB-113I20.2', 'GBGT1', 'LRRIQ3', 'MTIF2', 'TTC8']

存储一下adata的原始数据

# 存储一下adata的原始数据
adata.raw = adata
adata = adata[:, adata.var.highly_variable].copy() #非必须

5.降维

5.1 线性降维PCA

# 回归
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"]) 

# scale
sc.pp.scale(adata, max_value=10) 
sc.tl.pca(adata, svd_solver="arpack")

# 绘制 pca 聚类结果
sc.pl.pca(adata)

主成分贡献度：

sc.pl.pca_variance_ratio(adata, log=True)

选择高变基因以及pca分析之后，数据的变化：

5.2 计算邻接矩阵图

sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

5.3 非线性降维可视化

UMAP：

sc.tl.umap(adata)

# Using the igraph implementation and a fixed number of iterations can be significantly faster, especially for larger datasets
sc.tl.leiden(adata,flavor="igraph",n_iterations=2,resolution=0.9)

# umap图
sc.pl.umap(adata, color=["leiden"], legend_loc="on data", size=5)

6.差异分析：筛选亚群高表达基因

# Obtain cluster-specific differentially expressed genes
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon",pts=True)

sc.pl.rank_genes_groups_dotplot(
    adata, groupby="leiden", standard_scale="var", n_genes=5
)

每个亚群高表达基因top5展示：

marker基因表格：




    
pd.DataFrame(adata.uns["rank_genes_groups"]["names"]).head(10)

将差异结果整理成一个表格：

result = adata.uns["rank_genes_groups"]
result.keys() #看看结果都包含哪些部分

keys_to_get = ("pvals","logfoldchanges","pvals_adj","names")
subset = {k: pd.DataFrame(result[k]) for k in keys_to_get}
{key: value.shape for key, value in subset.items()}
{key: value.iloc[0:4,0:4,] for key, value in subset.items()}
type(result['names'])

groups = result["names"].dtype.names #记录数组提取列名
groups

# pbmc_markers = pd.DataFrame(
#     {
#         group + "_" + key: result[key][group]
#         for group in groups
#         for key in [ "pvals","logfoldchanges","pvals_adj","names"]
#     }
# )

# pbmc_markers.shape
# pbmc_markers.head(5)

import pandas as pd
import numpy as np
n = 4
split_df = []
for i in [int(g) for g in groups]:
    a = pbmc_markers.iloc[:, i:(i + n)].copy()
    a['cluster'] = str(i)
    a.columns =   a.columns.str.replace(r'\d+_', '', regex=True)
    split_df.append(a)
pbmc_markers = pd.concat(split_df, ignore_index=True)
pbmc_markers.head()