基于 Trae 的单细胞 RNA 测序分析与可视化
资料来源：火山引擎-开发者社区
1 本实践材料
AI 原生 IDE Trae（http://trae.com.cn）
2 实践目标在 AI 辅助下学习代码并完成对 scRNA-seq 数据进行质量控制、过滤、标准化、降维以及可视化。具体如下：1.质量控制：
使用 Scanpy 计算细胞和基因的质量控制指标（如 UMI 数、检测到的基因数）。
筛选低质量细胞（如 UMI 数过低的细胞）和低表达基因。
可视化：绘制小提琴图和散点图。
2.标准化与对数变换：
使用 Scanpy 进行标准化处理，使每个细胞的总 UMI 数一致。
对数据进行对数变换，减少偏态分布。
3.降维和聚类： 通过 PCA 和 UMAP 进行降维，并进行聚类分析。操作：
PCA 降维 ：使用 Scanpy 进行 PCA 降维，提取主要变异信息。
UMAP 降维 ：基于 PCA 结果，使用 Scanpy 进行 UMAP 降维。
t-SNE 降维 ：基于 PCA 结果想，使用 Scanpy 进行 t-SNE 降维。
4.可视化： 绘制 UMAP 图和 t-SNE 图，展示细胞聚类结果。
3 本实践需要的基础
Python 基础与环境配置
Python 语法基础：变量、数据类型、函数、类定义、事件处理。
包管理工具：需通过 pip 安装库。
开发环境：代码适合在 IDE（如 PyCharm, VS Code）中运行。
R、Rtools 安装与 R 环境配置（R 版本 4.4.2，Rtools44，版本过低无法操作)
技术栈
Scanpy：用于单细胞数据分析，包括数据预处理、降维和可视化。
Pandas：用于数据读取、筛选和转换。
参考：基于单细胞数据分析的实践指南（https://www.sc-best-practices.org/preamble.html）
数据集来源：figshare
4 实践步骤
4.1 导入库
import numpy as np
import scanpy as sc
import seaborn as sns
#画图的
from matplotlib import pyplot as plt
from scipy.stats import median\_abs\_deviation
#从 Scipy 统计模块中导入 median\_absolute\_deviation （MAD）函数。在单细胞数据分析中，它通常用于质量控制中的异常值检测。#设置图像参数
sc.settings.verbosity = 0
sc.settings.set\_figure\_params(
dpi=80,
facecolor=
"white"
,
frameon=False,
)
4.2 数据准备
1.使用 scanpy 读取数据文件（如 H5AD 格式）
使用 Scanpy 库（一个用于单细胞数据分析的 Python 工具）从一个 10x Genomics 格式的 .h5 文件中读取单细胞测序数据
adata = sc.read\_10x\_h5(
filename=
"filtered\_feature\_bc\_matrix.h5"
,
backup\_url=
"https://figshare.com/ndownloader/files/39546196"
,
)
由于网络问题，该链接可能有时不能正常访问，所以可以先将数据集下载至本地，命名为 filtered\_feature\_bc\_matrix.h5 ，再读取本地文件。
Trae 提示词
改为读取本地文件 filtered\_feature\_bc\_matrix.h5
将代码改为
filename =
"filtered\_feature\_bc\_matrix.h5"
# 使用 sc.read\_10x\_h5 读取本地文件
adata = sc.read\_10x\_h5(filename)
# 查看 AnnData 对象的概览
(adata)
2.检查数据结构，确认细胞和基因的分布。
4.3 数据预处理和可视化
1.质量控制：
使用 Scanpy 计算细胞和基因的质量控制指标（如 UMI 数、检测到的基因数）。
筛选低质量细胞（如 UMI 数过低的细胞）和低表达基因。
adata.var\_names\_make\_unique()
# 确保基因名唯一
Trae 提示词
如何筛选低质量细胞（如 UMI 数过低的细胞）和低表达基因。有：线粒体基因，核糖体基因，血红蛋白基因，并进行可视化。
==== 质量控制 ====# 标记特殊基因
adata.var[
'mt'
] = adata.var\_names.str.contains(r
'^MT-|^mt-|^MT[ABCDEFGH]|^Mt\.[A-Za-z]'
, regex=True)
adata.var[
'ribo'
] = adata.var\_names.str.startswith((
'RPS'
,
'RPL'
,
'rps'
,
'rpl'
))
adata.var[
'hb'
] = adata.var\_names.str.contains(
'^HB[AB]?[1-9]$'
, regex=True)
# 计算质控指标
sc.pp.calculate\_qc\_metrics(
adata,
qc\_vars=[
'mt'
,
'ribo'
,
'hb'
],
percent\_top=None,
log1p=False,
inplace=True
)
作用 ：识别线粒体基因（mt）、核糖体基因（ribo）和血红蛋白基因（hb）
原理 ：使用正则表达式匹配基因命名模式（如 MT-开头的线粒体基因）计算每个细胞的这些基因表达百分比（pct_counts_mt 等）
==== 质控前可视化 ====
(
"=== 质控前数据分布 ==="
)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
metrics = [
'n\_genes\_by\_counts'
,
'total\_counts'
,
'pct\_counts\_mt'
]
titles = [
'Genes per Cell'
,
'UMI Counts'
,
'Mitochondrial %'
]
for
i, metric
in
enumerate(metrics):
sns.violinplot(y=adata.obs[metric], ax=axs[i])
axs[i].set\_title(titles[i])
if
metric ==
'pct\_counts\_mt'
:
axs[i].set\_ylim(0, adata.obs[metric].max() * 1.2)
plt.tight\_layout()
plt.show()
展示未过滤数据的分布
关键指标 ：
1.每个细胞检测到的基因数
2.每个细胞的总 UMI 数
3.线粒体基因百分比
# ==== 细胞过滤 ====
(
"\n=== 应用过滤条件 ==="
)
min\_genes, max\_genes = 200, 5000
max\_mt, min\_counts = 20, 1000
cell\_filter = (
adata.obs.n\_genes\_by\_counts.between(min\_genes, max\_genes) &
(adata.obs.pct\_counts\_mt < max\_mt) &
(adata.obs.total\_counts > min\_counts)
)
adata = adata[cell\_filter, :].copy()
# ==== 基因过滤 ====
sc.pp.filter\_genes(adata, min\_cells=10)
细胞过滤标准 ： 保留基因数在 200-5000 之间的细胞 线粒体基因百分比 <20% 总 UMI 数 >1000
目的 ：移除低质量细胞（基因太少/线粒体基因过高）和多重捕获（基因过多）
基因过滤作用 ：去除在少于 10 个细胞中表达的基因
目的 ：减少噪声基因对后续分析的影响
# ==== 质控后可视化 ====
(
"\n=== 质控后数据分布 ==="
)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
for
i, metric
in
enumerate(metrics):
sns.violinplot(y=adata.obs[metric], ax=axs[i])
axs[i].set\_title(f
"{titles[i]} (Filtered)"
)
if
metric ==
'pct\_counts\_mt'
:
axs[i].set\_ylim(0, adata.obs[metric].max() * 1.2)
plt.tight\_layout()
plt.show()
作用 ：验证过滤效果
对比 ：与质控前分布对比，确认过滤阈值合理性
==== 最终结果输出 ====print
(
"\n=== 质控结果 ==="
)
(f
"剩余细胞数: {adata.n\_obs}"
)
(f
"剩余基因数: {adata.n\_vars}"
)
(f
"线粒体百分比范围: {adata.obs.pct\_counts\_mt.min():.2f}% - {adata.obs.pct\_counts\_mt.max():.2f}%" )
输出指标 ：
1.保留的细胞总数
2.保留的基因总数
3.线粒体基因百分比的范围
整体流程 ：质量控制 → 数据过滤 → 结果验证 → 输出报告
# 改用散点图验证数据
plt.scatter(adata.obs.total\_counts, adata.obs.pct\_counts\_mt, alpha=0.3)
plt.xlabel(
'Total UMI'
)
plt.ylabel(
'MT%'
)
plt.show()
4.4 标准化
Trae 提示词
如何进行标准化处理
(
"\n=== 数据标准化 ==="
)
# 保存原始计数到 layers
adata.layers[
"counts"
] = adata.X.copy()
# 标准化处理
sc.pp.normalize\_total(adata, target\_sum=1e4)
sc.pp.log1p(adata)
数据标准化： 将每个细胞的总基因表达量标准化为 10,000。标准化总基因表达量可以消除细胞测序深度的差异，使得不同细胞之间的基因表达量可以更公平地比较。例如，测序深度高的细胞可能会有更高的总表达量，通过标准化可以消除这种偏差。 对数据进行对数变换，具体操作：
对每个基因表达量取自然对数后加 1（ log1p 是 log(1 + x) 的操作，加 1 的目的是避免对 0 取对数导致的数学错误）。
作用：对数变换可以压缩数据的动态范围，使得数据更接近正态分布，便于后续的分析。例如，基因表达量通常具有长尾分布，对数变换可以使其分布更加均匀。
4.5 降维分析
选择表达量变异最大的 2000 个基因。
对数据进行缩放。
进行主成分分析（PCA）。
构建邻接图。
进行 UMAP 和 t-SNE 降维。
Trae 提示词
如何进行降维分析
保存原始数据
adata\_raw = adata.copy()
# ==== 降维分析 ====
sc.pp.highly\_variable\_genes(adata, n\_top\_genes=2000)
sc.pp.scale(adata, max\_value=10)
sc.tl.pca(adata, svd\_solver=
'arpack'
)
sc.pp.neighbors(adata, n\_pcs=15)
sc.tl.umap(adata)
sc.tl.tsne(adata, n\_pcs=15, random\_state=42)
# ==== 可视化 ====
plt.figure(figsize=(18, 6))
sc.pl.pca\_scatter(adata, color=
"total\_counts"
)
# 质控对比
plt.subplot(131)
sns.violinplot(data=adata\_raw.obs[[
'n\_genes\_by\_counts'
,
'total\_counts'
,
'pct\_counts\_mt'
]],
inner=
"quartile"
, palette=
"Pastel1"
).set\_title(
'Pre-QC'
)
# UMAP 可视化
sc.pl.umap(adata, color=
'pct\_counts\_mt'
, legend\_loc=
'right margin'
,
title=
'UMAP (MT%)'
, edges=False)
# t-SNE 可视化
sc.pl.tsne(adata, color=
'n\_genes\_by\_counts'
, show=False, legend\_loc=
'right margin'
,
title=
't-SNE (Genes/Cell)'
, edges=False)
plt.tight\_layout()
plt.show()
# ==== 结果输出 ====
(f
"\n 最终细胞数: {adata.n\_obs}"
)
(f
"最终基因数: {adata.n\_vars}"
)
(f
"MT% 范围: {adata.obs.pct\_counts\_mt.min():.1f}% - {adata.obs.pct\_counts\_mt.max():.1f}%"
)
<<< 左右滑动见更多 >>>
降维分析代码详解：
sc.pp.highly\_variable\_genes(adata, n\_top\_genes=2000)
功能：选择表达量变异最大的基因。
具体操作：计算每个基因的表达量变异，并选择变异最大的 2000 个基因。
目的：减少数据维度，保留最具信息量的基因。高变基因通常更能反映细胞之间的差异，提高降维分析的效率和效果。
对数据进行缩放
sc.pp.scale(adata, max\_value=10)
功能：对数据进行缩放。
具体操作：将数据缩放到均值为 0，标准差为 1 的范围内，最大值限制为 10。
目的：便于后续的主成分分析（PCA）。缩放后的数据更容易进行线性分析，同时限制最大值可以避免极端值的影响。
进行主成分分析（PCA）
sc.tl.pca(adata, svd\_solver=
'arpack'
)
功能：进行主成分分析（PCA）。
具体操作：将高维数据投影到低维空间，提取主要的变异方向。
目的：PCA 是一种常用的降维方法，可以去除噪声并保留数据的主要结构，为后续的 UMAP 和 t-SNE 分析提供基础。
构建邻接图
sc.pp.neighbors(adata, n\_pcs=15)
功能：基于 PCA 结果构建邻接图。
具体操作：使用前 15 个主成分构建邻接图，表示细胞之间的相似性。
目的：邻接图用于后续的 UMAP 和 t-SNE 分析，帮助算法更好地理解细胞之间的局部和全局结构。
进行 UMAP 降维
sc.tl.umap(adata)
功能：进行 UMAP 降维。
具体操作：将数据投影到二维空间，保留局部结构的同时尽量保留全局结构。
目的：UMAP 是一种非线性降维方法，适合高维数据的可视化。它能够清晰地展示细胞之间的聚类和连续性。
进行 t-SNE 降维
sc.tl.tsne(adata, n\_pcs=15, random\_state=42)
功能：进行 t-SNE 降维。
具体操作：使用前 15 个主成分进行 t-SNE 降维，随机种子为 42，确保结果可复现。
目的：t-SNE 是另一种非线性降维方法，特别擅长展示局部结构，但可能会牺牲一些全局结构。它也适合高维数据的可视化。
小结降维分析：
选择高变基因：减少数据维度，保留最具信息量的基因。
数据缩放：便于后续的主成分分析。
主成分分析（PCA）：提取主要的变异方向，为后续分析提供基础。
构建邻接图：表示细胞之间的相似性，用于 UMAP 和 t-SNE 分析。
UMAP 降维：将数据投影到二维空间，保留局部和全局结构。
t-SNE 降维：将数据投影到二维空间，特别擅长展示局部结构。
通过这些步骤，可以有效地处理单细胞测序数据，为后续的细胞类型鉴定、差异表达分析等提供基础。
总结 1.数据质量评估 ：
如果大部分细胞的检测到的基因数、总 UMI 数、线粒体基因比例、核糖体基因比例和血红蛋白基因比例都符合预期范围，说明数据质量较好。
如果存在异常值，需要根据设定的阈值过滤掉低质量细胞。
2.筛选阈值的确定 ： 根据小提琴图和散点图的分布情况，可以确定合理的筛选阈值。
例如：
检测到的基因数：>200。
总 UMI 数：>200。
线粒体基因比例：<5%。
核糖体基因比例：<20%。
血红蛋白基因比例：<10%。
3.数据分布的直观展示 ： 小提琴图和散点图提供了数据分布的直观展示，帮助研究人员更好地理解数据的特征和潜在问题。 通过这些分析，可以有效地去除低质量细胞，提高后后续分析的可靠性和准确性。4.保留局部和全局结构 ： UMAP 在降维过程中既保留了数据的局部结构（即相似的数据点在低维空间中仍然保持相近），也尽可能保留了全局结构（即不同簇之间的相对位置）。 这使得 UMAP 图能够更好地展示数据的整体分布和细胞类型的层次关系。UMAP 图特别擅长展示数据中的连续性和过渡性。例如，在细胞分化过程中，UMAP 图可以清晰地展示细胞从一种状态到另一种状态的过渡路径。 这对于研究细胞分化、发育过程以及细胞状态的动态变化非常有帮助。 t-SNE 特别擅长展示数据的局部结构，能够将相似的数据点聚集在一起，形成清晰的簇。 这使得 t-SNE 图非常适合用于识别不同的细胞类型或状态，尤其是在单细胞测序数据中。 t-SNE 能够揭示数据中隐藏的模式和结构，即使这些模式在高维空间中不明显。 通过调整参数（如 perplexity ），可以更好地控制局部和全局结构的平衡，从而发现不同的细胞亚群。
5.在实际应用中 ，研究人员通常会同时使用 UMAP 和 t-SNE 图，以获得更全面的数据视图。 UMAP 图可以帮助理解数据的整体结构和连续性，而 t-SNE 图可以帮助识别局部结构和细胞亚群。通过结合这两种方法，可以更深入地理解单细胞测序数据中的生物学信息。
还未添加个人签名 2022-01-25 加入
还未添加个人简介
评论