Scanpy 分析 3k PBMCs:数据预处理

简介: Scanpy 分析 3k PBMCs:数据预处理

引言

本系列讲解 使用Scanpy分析单细胞(scRNA-seq)数据教程,持续更新,欢迎关注,转发!

数据集

本次使用的数据集包含一位健康供体的3k PBMCs,这些数据可以从10x Genomics的官方网站免费获取。如果你使用的是unix系统,只需取消以下代码的注释并运行,就可以完成数据的下载和解压操作。另外,代码的最后一行会自动创建一个目录,用于存放处理后的数据。

# !mkdir -p data
# !curl https://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz -o data/pbmc3k_filtered_gene_bc_matrices.tar.gz
# !cd data; tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz
# !mkdir -p write

数据导入

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

把计数矩阵导入到一个 AnnData 对象里,这个对象能存储很多注释信息以及数据的各种不同表现形式。它还自带一种基于 HDF5 的文件格式,即 .h5ad

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)

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

我们来收集一些有关线粒体基因的信息,这些基因对于质量控制来说很重要。

如果细胞里线粒体基因的比例很高,那就意味着这个细胞的质量不太好。出现这种情况可能是因为细胞膜有破损,细胞质 RNA 从破洞里流失了。因为线粒体比单个转录本分子大,所以相对来说更难从细胞膜的裂缝里跑出去。

通过使用 pp.calculate_qc_metrics 这个工具,我们可以非常高效地计算出很多质量控制指标。

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

绘制了一些质量控制指标的小提琴图,具体包括:

  • 计数矩阵中表达的基因数量

  • 每个细胞的总计数

  • 线粒体基因的计数百分比

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

细胞过滤

接下来,要移除那些表达过多线粒体基因或者总计数过多的细胞:

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

具体操作是通过切片 AnnData 对象来进行过滤。

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

对数据矩阵进行总计数归一化(也就是库大小校正) 每个细胞的读段数调整到10,000个,这样就可以让不同细胞之间的计数具有可比性了。

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

对数转换

sc.pp.log1p(adata)

识别highly-variable基因

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

sc.pl.highly_variable_genes(adata)

AnnData 对象的 .raw 属性设置为经过归一化和对数化处理的原始基因表达数据,以便后续用于差异测试和基因表达的可视化。这一步相当于把 AnnData 对象当前的状态给固定下来。

adata.raw = adata.copy()

进行过滤

adata = adata[:, adata.var.highly_variable]

对每个细胞的总计数和线粒体基因表达百分比的影响进行回归分析。然后将数据缩放到单位方差。

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

对每个基因也进行缩放,使其方差为单位方差。并且,将超过标准差10的值进行裁剪处理。

sc.pp.scale(adata, max_value=10)
相关文章
|
4月前
|
安全 网络协议 Linux
Linux网络应用层协议展示:HTTP与HTTPS
此外,必须注意,从HTTP迁移到HTTPS是一项重要且必要的任务,因为这不仅关乎用户信息的安全,也有利于你的网站评级和粉丝的信心。在网络世界中,信息的安全就是一切,选择HTTPS,让您的网站更加安全,使您的用户满意,也使您感到满意。
137 18
|
4月前
|
安全 搜索推荐 数据可视化
提升游戏APP留存与付费:高性价比工具与实践
2024年中国移动游戏市场收入达2843.3亿元,头部产品占据超七成份额。在存量市场竞争中,留住用户成为关键。Xinstall推出三大工具:Deeplink实现无缝跳转提升转化;安全加速优化启动体验;分享统计助力社交裂变,以低成本助力开发者提升留存与付费转化率。
125 3
|
4月前
|
编解码 算法 数据可视化
Scanpy 分析 scRNA-seq:细胞类型注释
Scanpy 分析 scRNA-seq:细胞类型注释
Scanpy 分析 scRNA-seq:细胞类型注释
|
4月前
|
监控 Kubernetes Go
日志采集效能跃迁:iLogtail 到 LoongCollector 的全面升级
LoongCollector 在日志场景中实现了全面的重磅升级,从功能、性能、稳定性等各个方面均进行了深度优化和提升,本文我们将对 LoongCollector 的升级进行详细介绍。
414 86
|
4月前
|
机器学习/深度学习 存储 人工智能
浅入浅出——生成式 AI
团队做 AI 助理,而我之前除了使用一些 AI 类产品,并没有大模型相关的积累。故先补齐一些基本概念,避免和团队同学沟通起来一头雾水。这篇文章是学习李宏毅老师《生成式 AI 导论》的学习笔记。
476 27
浅入浅出——生成式 AI
|
4月前
|
存储 SQL 大数据
从 o11y 2.0 说起,大数据 Pipeline 的「多快好省」之道
SLS 是阿里云可观测家族的核心产品之一,提供全托管的可观测数据服务。本文以 o11y 2.0 为引子,整理了可观测数据 Pipeline 的演进和一些思考。
318 35
|
4月前
|
数据采集 Web App开发 JavaScript
基于Selenium的Python爬虫抓取动态App图片
基于Selenium的Python爬虫抓取动态App图片
337 68
|
4月前
|
开发工具 开发者
HarmonyOS实战:腾讯IM之聊天详情页面搭建(二)
本文讲解了在鸿蒙系统中实现腾讯IM聊天功能的完整流程,涵盖对话列表展示、历史消息获取、实时消息更新及文本消息发送等核心功能。通过实际代码示例,详细说明了如何利用IM SDK实现聊天业务逻辑。适合开发者逐步学习并实践,建议点赞收藏以便参考。
157 9
HarmonyOS实战:腾讯IM之聊天详情页面搭建(二)
|
4月前
|
人工智能 前端开发 Java
“最近我给有代码洁癖的同事墙裂安利了通义灵码”
通义灵码2.5.0版本现已全面支持Qwen3,采用混合专家架构,参数量仅为DeepSeek-R1的1/3,是国内首个“混合推理模型”。它在性能评测中超越了DeepSeek-R1、OpenAI-o1等国际主流模型,并全面支持MCP能力,集成国内最大MCP中文社区。作为程序员体验后发现,通义灵码可通过简单指令生成完整项目代码,包括前后端、接口调用等,大幅降低开发门槛。文中通过两个Demo展示了其强大功能:一是聚合多平台热榜数据并推送微信通知;二是基于高德和12306 MCP生成旅游攻略HTML页面。整个过程无需手动编写代码,推荐开发者尝试。
257 47