前面，我们生信技能树的讲师小洁老师与萌老师新开了一个学习班：《掌握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()
整理的表格如下：
今天学习到这里~
