05单细胞分析2025-scanpy分析流程-02示例文件预处理和标准化处理及保存

scanpy分析流程-02示例文件分析

1.读取10x的数据到adata

通过 sc.read_10x_mtx() 函数从指定目录加载10x Genomics的单细胞RNA测序数据，读取的结果存储在adata 中

使用基因符号作为变量名 gene_symbols

并启用缓存以加快后续加载速度  cache=True

它是一个 AnnData 对象，方便在Scanpy的后续分析流程中使用

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
)

2.初步绘制整个单细胞数据集中表达量最高的基因

sc.pl.highest_expr_genes(adata, n_top=20)

sc.pl.highest_expr_genes 是一个用于绘制图表的函数，它展示了在整个单细胞数据集中表达量最高的基因 这有助于用户了解在数据集中哪些基因是最具代表性的、被高度表达的基因

Y轴：表示表达量最高的基因列表
X轴：表示这些基因的表达量占总表达量的百分比

#sc表示scanpy，是一个包的名字# pl表示plot作图

sc.pl.highest_expr_genes(adata,n_top=20)

每一种基因在所有细胞中的表达量占比的平均值=一个基因的表达量count值 3万种基因的表达量count 一个细胞约有3万种基因，把这一个细胞中的3万种基因的表达量count值加起来 然后用一个基因的表达量count值除以刚才的总表达量值比例值，得到这个基因在这个细胞中的表达量比例值计算 然后选取平均值最高的20种基因：做箱线图展示每一种基因在所有细胞 例如2798个细胞中个字的表达量的占比值 一个细胞测得2000个基因 一个细胞的表达矩阵有2万 

3.预处理 Preprocessing 对基因和细胞进行过滤 

sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

基础过滤 Basic filtering：

pp表示：scanpy.preprocessing

这段代码是对单细胞RNA测序数据中的细胞和基因进行过滤的。过滤步骤在单细胞分析中非常关键，用于去除质量差的细胞和低表达的基因，从而提高数据分析的可靠性

1. sc.pp.filter_cells(adata, min_genes=200)：

• 目的 ：过滤掉那些基因数（即检测到的基因数）少于200的 细胞

• 作用 ：

• 每个细胞中都包含了一定数量的基因表达数据，但是一些细胞可能只表达了很少的基因

• （例如由于实验错误或低质量的细胞）。通常这些低表达的细胞会被认为是低质量的

• min_genes=200 这个参数的含义是： 细胞中检测到的基因数必须至少有200个 ，

• 如果某个细胞表达的基因数量少于200个，则该细胞会被移除

• 这是为了过滤掉那些可能没有生物学意义或质量较差的细胞（例如，死亡细胞或双细胞等）

• 结果 ：adata 对象中的细胞会减少，只保留那些在至少200个基因中有表达的细胞

2. sc.pp.filter_genes(adata, min_cells=3)：

• 目的 ：过滤掉那些在少于3个细胞中表达的 基因 。

• 作用 ：

• 单细胞RNA测序数据中有很多基因只在少数细胞中表达，这些基因可能是噪音或不重要的基因

• 为了避免数据的稀疏性影响分析，通常会过滤掉这些基因

• min_cells=3 这个参数的含义是： 基因必须至少在3个细胞中有表达 ，如果某个基因只在少于

• 3个细胞中检测到，则该基因会被移除

• 这种过滤可以去掉极少表达的基因，通常它们对下游分析（如聚类、降维等）的影响较小

• 结果 ：adata 对象中的基因会减少，只保留那些在至少3个细胞中有表达的基因

总结：

sc.pp.filter_cells(adata, min_genes=200)

sc.pp.(adata, min_cells=3)

filter_cells 函数对细胞进行过滤操作 min_genes=200 一个细胞表达的基因至少有200

filter_genes 函数对基因进行过滤操作 min_cells=3 一种基因至少在3个细胞中表达

4.过滤线粒体高的细胞

# 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
)

# annotate the group of mitochondrial genes as "mt“将线粒体基因组注释为 “mt“ 

adata.var["mt"] = adata.var_names.str.startswith("MT-") #找到线粒体基因赋给mt 

sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, log1p=False,

 inplace=True) #percent_top=None 表示取百分之多少的表达量最高的基因做计算 这里是None，表示不做计算

让我们收集一些关于线粒体基因的信息，这些基因对质量控制非常重要Let’s assemble some information about mitochondrial genes, which are important for quality control.Citing from “Simple Single Cell” workflows

https://master.bioconductor.org/packages/release/workflows/html/simpleSingleCell.html#examining-gene-level-metrics

High proportions are indicative of poor-quality cells (Islam et al. 2014; Ilicic et al. 2016), possibly because of loss of cytoplasmic RNA from perforated cells. The reasoning is that mitochondria are larger than individual transcript molecules and less likely to escape through tears in the cell membrane.

线粒体高比例表明细胞质量较差（Islam 等人，2014 年；Ilicic 等人，2016 年），可能是因为穿孔细胞的细胞质 RNA 丢失。 其原因是线粒体比单个转录本分子大，不太可能通过细胞膜的裂口逃逸，线粒体基因含量的比例，会作为单细胞转录组测序数据质量好坏的一种表征（由于死细胞的细胞膜的通透性变大，细胞质中的分子，比如mRNA会到细胞外，造成胞内的mRNA减少。但是线粒体的分子仍旧会留在线粒体内，也就是还留在细胞内）线粒体基因含量高，细胞可能发生凋亡

5.小提琴图展示线粒体基因的相关参数

sc.pl.violin(
    adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4,
    multi_panel=True,
)

the number of genes expressed in the count matrix=n_genes_by_counts  #每个细胞检测到的基因数量 / 表达的基因数  

the total counts per cell=total_counts  #每个细胞的总基因计数/表达数 / 测序深度

the percentage of counts in mitochondrial genes=pct_counts_mt  #线粒体基因的计数百分比

6.散点图展示线粒体基因的占比：

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

sc.pl.scatter(adata, x="total_counts", y="pct_counts_mt") # 测序深度与线粒体比例 线粒体和测到的基因数 无相关性

sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts") # 测序深度与检测到的基因数  有相关性 正相关

7.设置过滤细胞和基因数标准

adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :].copy()

把adata理解为是一个dataframe，里面的每一行代表一个细胞（observation），每一列代表一种基因（variable)

从上面的作图来看，有表达的基因的种类数，在6000以上的，这样的细胞个数很少，属于极端情况。所以可以把这样的细胞的数据给剔除掉。在写条件代码做处理时，保留有表达的基因种类数，在2500以下的，这样的细胞的数据行

从上面的作图来看，线粒体基因的占比，在5%以上的，这样的细胞个数很少，属于极端情况。所以可以吧这样的数据给剔除掉。在写条件代码做处理时，是保留，线粒体基因的占比，在5%以下的，这样的细胞的数据行（数据中的线粒体基因的细胞数量几乎没有）

8.标准化的处理数据

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

Normalize_total()函数小结

sc.pp.normalize_total()标准化的处理，就是每一种基因的表达量，除以某一个细胞中的所有基因的总表达量，这算出来是一个比例值，也就是基因表达量的相对含量。细节层面的话，遇到0则加1的处理，避免出现0值的情况

sc.pp.normalize_total(adata, target_sum=1e4) # target_sum=1e4为扩大的倍数

那么对于一个细胞的各个基因的这些比例值，它们加起来，结果为1。对于所有的细胞，它们各自都是这样，由于这些比例值是小数，不够直观，或者不利于后面的运算处理。我们要把它们转化为大于1的数字，接近原来counts值的大小。我们会把这些比例值乘以一个固定的大的数字，例如10000，或1000000。这个固定的数字对所有细胞都是一样的。我们把这样的转化后的值称为normalized_counts值

这样，每一个细胞的所有基因的normalized_counts值之和都相等，也就是细胞001的normalizedcounts.sum()它们都是相等的一个值

关于这个名字“normalized_total",可以理解为是用total总的基因表达量来作为分母，来做的标准化处理target_sum=1e4为扩大的倍数

sc.pp.log1p(adata)

对以上数值进行取LOG值 sc.pp.log1p(adata) np.log1p()是numpy的函数，含义为log1 plus x,即log(1+x)，默认底数为e=2.718.可以写参数out来指定生成的结果保存到哪一个变量中，一般可以用来表示对array做原位的改变操作，对标准化后的表达量数值，做1n（1+x）的对数化处理 sc.pp.log1p(adata)为什么要做对数化处理呢？是为了让数值的范围缩小一些吗？做ln(1+x)处理，这一步很重要，在scanpy体系中必须做。

因为下面的标注表达量高度变化的基因的函数，期望接受的数值，就是In（1+×）的数值而不是原始的数值，在下面标注表达量高度变化的基因的函数的代码内部，会对数值做指数运算，然后减1，求出count值。如果在这里不做In(1+x)处理的话，那么在调用标注表达量高度变化的基因的函数时，得到的数值就不对了

9.计算高变基因

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

高变基因：

这里的表达量高度变化的基因，这种变异，是指相对于基因自身的表达量的平均值的变异程度，就是dispersion=variance/mean Expects logarithmized data,except when flavor='seurat_v3’ in which count data is expected,默认情况下，需要输入的数值是对数化处理后的数值。如果是 flavor=seuratv3'则要求输入的是count值，而不是对数化的值

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

只有满足条件，表达量平均值，在0.0125到3之间的，dispersion值不能小于0.5，这些基因才有可能被标注为表达量高度变化的基因，才能做下游的降维分析

highly-variablegenes(）函数的作用是把表达量高度变化的基因给标注出来

adata，里面含对数化的数值。但是当flavor='seurat_v3"时，应该要用 coun数值，而不是对数化的值

对于Satijia15和zheng17的基于dispersion值的方法，标准化的dispersion值，是通过缩放dispersion的平均值和标准差，对那些落入一个给定箱子的基因，对于基因的平均表达量。这意味看，对于平均表达量的每个箱子，表送量高度变化的基因会被选出来

而另一种方法是Stuart19的方法，对每一种基因计算标准化的方差值。首先，用规范化的标准差对数据做标准化（即对每一个基因做z-score标准化）。然后，计算上述标准化后的每一个基因的方差值，这个方差值就是标准化的方差值。对基因按照方差值的大小来排序

10.散点图展示高变基因

sc.pl.highly_variable_genes(adata)

散点图第1个图，把adata中的每一个基因在所有细胞中的平均值，以及标准化的dispersion值取出来，以及哪些基因是highly variable gene的信息取出来，做散点图，每一个点代表一个基因，横坐标是基因的表达量的平均值，纵坐标是标准化的dispersion值，然后给是highly variable gene的点涂黑色，给不是 highly variable gene的点涂灰色

散点图第2个图，只是把上面的标准化的dispersion值换成了没做标准化的dispersion 值，其它的数据信息以及画图做法都一样

11.备份数据到adata.raw

adata.raw = adata.copy()

Set the.raw attribute of the AnnData object to the normalized and logarithmized raw gene expression for later use indifferential testing and visualizations of gene expression.This simply freezes the state of the AnnData object

差异检验，以及可视化基因表达量时，要用标准化和对数化的基因表达量数据。这里我们把 adata赋值给adata.raw 后面可以拿到现在的原始的基因表达量的数值

记录 简要的学习 成长 体会 心得 总结 思路 心流 链接 资料  懂得都懂