知识库/生信教程与工具/网页教程
全库搜索 Markdown 图册
Paper Analysis · 14 分钟阅读

原理+代码流程，一文学会华大C4平台scRNA-seq定量分析

站内编号 KB-000071
clippings生信教程单细胞scRNA-seq华大C4DNBelabDNBC4tools上游定量SeuratScanpy暂未试用
原文
板块与文章导航

摘要LLM 归库摘要

  • 核心内容： 华大 DNBelab C4 平台单细胞转录组定量分析完整教程：平台原理（双磁珠液滴技术）→ DNBseq 测序原理 → DNBC4tools 上游定量流程（下载/安装/参考基因组/建库/fastq 定量）→ 下游结果读取（R Seurat/自定义函数 + Python Scanpy/自定义函数）。
  • 适用场景： 手头有华大 DNBelab C4 / DNBC4tools 产出的单细胞数据（cDNA+oligo fastq 或 MEX 矩阵），需要从上游定量到下游分析。
  • 输入： ①上游：cDNA fastq + oligo fastq + 参考基因组 FASTA + GTF 注释；②下游：output/raw_matrix（MEX）或 filter_feature.h5ad。
  • 输出： 上游得到 01.data/02.count/03.analysis/04.report；下游得到 Seurat 对象或 AnnData 对象，可继续个性化分析。
  • 关键参数/步骤： tools mkgtf（GTF 校正/过滤）→ rna mkref（建索引，--threads）→ rna run（定量，--cDNAfastq/--oligofastq/--genomeDir）；读取用 Read10X/ReadMatrix_C4（R）/read_h5ad/read_anndata_C4（Python）；min.cells=3, min.features=200 质控。
  • 常见错误： GTF 需含 gene/transcript/exon 且带 gene_id/transcript_id，缺字段需先 mkgtf --action check 校正；ref.json 中 genomeDir 需指向参考路径；GTF 非唯一比对 reads 会被过滤。
  • 与口腔生信连接： 若课题组数据来自华大 DNBelab C4 平台（如牙周组织/口腔样本单细胞），可用此流程完成上游定量并接入 Seurat/Scanpy 下游；微生物组平台可类比。
  • 下一步处理： 暂未试用；需要实际华大 C4 数据/测试文件（教程提供百度网盘提取码 yckd）才能试跑。
  • 归库说明： 分类到 03 生信教程与工具库，置信度 0.92，理由：单细胞上游定量流程教程，含完整可复用代码。

BIOMAMBA Biomamba 生信基地 2025年12月28日 09:05

图片

一、写在前面

单细胞转录组技术 无疑是撬动 细胞异质性研究 、解析生命复杂调控网络的核心工具。单细胞转录组学已深度应用到 肿瘤微环境解析 、 发育生物学谱系追踪 、 免疫细胞亚群鉴定 、 神经细胞分型 等多个前沿领域，为科研工作者揭示生命活动的发生机制提供了前所未有的技术支撑。

作为国产单细胞平台，华大DNBelab C4平台具有捕获效率高、数据分析丰富、平台稳定性优秀等特征。本推送教程采用 “理论介绍 — 软件使用 — 上游定量分析” 的模块，给大家提供一套利用DNBC4tools软件读取上游定量文件分析流程，完成读取之后就可以正常进行下游的个性化分析，我们也给大家准备了更详细的教程： 这可能是全网最全的scRNA-seq分析教程！

此外，此前我们已陆续分享过多种主流平台的上游分析流程，包括 10X Genomics 平台基于 Cellranger 的上游文件处理方案、Drop-Seq 平台借助 dropseqRunner 的上游分析流程、 BD Rhapsody 平台的单细胞转录组定量流程 ，以及 10X Visium 平台通过 Space Ranger 软件实现的上游定量方法 。本次我们将聚焦 DNBC4tools 软件读取上游定量文件的核心过程，为大家提供针对性的学习指引。

学前Tips：本教程基于Linux环境演示，计算资源不足的同学可参考：

足够支持你完成硕博生涯的生信环境

忘记宣传了，独享用户连技术支持都是独享的

访问链接：https://biomamba.xiyoucloud.net/

欢迎联系客服微信[Biomamba_zhushou]获取帮助

如果你对下面的教程比较迷茫，那么你可以先行学习编程教程：

生信Linux及服务器使用技巧

如果需要 单细胞数据分析 教学 、生信热点 全文复现 、自测数据 个性化分析 辅导、 常态化实验学习 ，欢迎联系客服微信[Biomamba_zhushou]

图片

二、C4 单细胞测序平台

1. C4 单细胞测序平台介绍

C4 平台是由华大智造自主研发的一站式单细胞测序数据分析平台，C系列平台基于液滴微流控策略，可以实现高效细胞捕获，完成高效的mRNA反转录和扩增后，结合华大DNBSEQ™高通量测序平台和Dr.Tom单细胞全交互式数据挖掘交付系统，进行高性价比单细胞转录组文库测序和数据分析挖掘，可广泛应用于研究免疫细胞、肿瘤微环境、神经系统发育及胚胎发育和干细胞分化等。

图片

优势 详细描述
稀有细胞高效捕获 高效捕获稀有类型细胞，实现研究新突破。
多种样本兼容性 除常规活细胞悬液外，珍稀液氮速冻样本也可实现细胞核水平研究。
测序数据高保真 搭配 DNBSEQ™ 测序系统，确保单个细胞数据不会混淆，数据质量更高。
自主数据分析 支持 Dr. Tom 交互式单细胞数据挖掘系统，助力用户轻松自主分析。
一站式服务模式 提供流式分选服务，一站式实现从组织样本到细胞数据获取及细胞群验证。

2. 测序原理#

图片

C4 单细胞测序平台基于微流控液滴技术，但是与传统的微流控液滴技术一个cell只与一个磁珠匹配不同，C4 平台会使用一大一小两个磁珠和捕获单个细胞进行配对，不同磁珠承担不同功能：

  • 大磁珠： 携带细胞条形码（Cell Barcode）和唯一分子标识符（UMI），用于标记每个细胞的转录本来源及去重。
  • 小磁珠：用于标记大磁珠（或者说用于标记液滴），不同的液滴会带有不同的小磁珠标记（Index）,同一个液滴内的所有大磁珠都带有相同的Index标记，数据分析时，可以通过大磁珠上所携带的Index信息将同一液滴内的所有大磁珠数据进行归类，从而实现对单个细胞的转录组数据进行准确归属。

  理想情况下微流控液滴技术希望1个液滴内含有1个磁珠1个细胞，实际情况(能产生数据的液滴)中情况可能是1个磁珠1个细胞;1个磁珠多个细胞;多个磁珠1个细胞;多个磁珠多个细胞，而C4平台通过使用两种不同大小的磁珠，可以有效区分这些情况，从而实现对单细胞的高效捕获和准确归属。

图片

3. 应用测序#

目前,基于华大智造 DNBelabC 系列高通量单细胞RNA文库制备试剂盒开展的研究或应用涵盖的样本类型丰富，包括人、鼠、拟南芥、水稻等超过50个物种的300多个组织类型。 根据已证实的实验经验表明,人的PBMC细胞、脑细胞核、脂肪细胞核、鼠肺癌细胞等样本类型的细胞捕获效率>50%,基因检出中位数2000左右;同时兼容复杂样本，如OCT包埋样本、鸟类、植物细胞核和植物原生质体均可在DNBelabC系列的单细胞平台实现。

  截至2025年7月，DNBelab C系列单细胞平台累计发表200+篇SCI文章，包含4篇Nature，4篇Cell，1篇Science，平均影响因子12+，其中IF>10文章数有100+篇。

图片

三、上游定量

1. DNBseq技术

DNBseq（DNA Nanoball Sequencing） 是华大基因自主研发的高通量测序技术，核心基于 DNA纳米球（DNA Nanoball，DNB）和高密度测序芯片。与传统NGS技术（如Illumina的桥式PCR扩增）不同，DNBseq避免了PCR扩增导致的重复误差，通过线性扩增生成单链DNA纳米球，结合联合探针锚定聚合（cPAS）技术进行测序。

  DNBseq测序文库构建分为常规建库和DNB制备两大阶段：

常规建库 （与传统NGS相似）

片段化：将DNA打断至200-500bp。 末端修复 & 加接头：修复片段末端，添加与测序芯片匹配的接头序列。 纯化：选择目标片段大小。

DNB制备 （关键技术）

  • 环化：将双链DNA连接成环状单链（ssDNA circle）。
  • 滚环扩增（RCA）：以环状DNA为模板，生成串联重复的线性单链结构。
  • DNB生成：通过机械或化学断裂，形成包含多拷贝重复序列的纳米球（直径~300nm）。
  • DNB芯片加载：将DNA纳米球固定到高密度芯片上，形成规则阵列。

优势 ：因无需PCR扩增，避免了GC偏好性和重复序列错误，显著提高数据均一性。

图片

对于采用DNBelab C系列高通量单细胞RNA文库制备试剂盒制备的文库是通常是这样式的：

图片

因此我们需要采用 DNBC4tools 进行上游的定量测序，用于分析高通量 DNBelab C 系列 TM 单细胞数据集的开源且灵活的流程。

2. DNBC4tools 下载与安装#

## 切换工作目录
cd /home/chenao/workspace/protocol/C4_sc
#### 使用 \`wget\` 下载
wget-O dnbc4tools2.1.3.tar.gz "ftp://ftp2.cngb.org/pub/CNSA/data7/CNP0006367/Single_Cell/CSE0000448/dnbc4tools2.1.3.tar.gz"
#### 使用 \`curl\` 下载
# curl -o dnbc4tools2.1.3.tar.gz "ftp://ftp2.cngb.org/pub/CNSA/data7/CNP0006367/Single_Cell/CSE0000448/dnbc4tools2.1.3.tar.gz"

解压下载好的安装包

## 执行解压
tar-xzvf dnbc4tools2.1.3.tar.gz

查看生成的文件目录，我们会看到类似的结构

ls dnbc4tools2.1.3

## dnbc4tools
## external
## lib
## misc
## sourceC4.bash

这就完成了安装。dnbc4tools 现在可以使用了。

3. DNBC4tools 使用示例#

3.1 下载物种参考基因组文件#

下载参考基因组和注释文件，根据自己的物种按需选择：

  • 基因组文件：基因组文件应以 FASTA 格式提供，包含特定物种的完整基因组序列，包括染色体、线 粒体及其他遗传信息，通常为主装配版本。这些文件为基因组分析和比对提供基础数据。
  • 注释文件：基因组注释文件应 以 GTF 格式提供，包含基因组中基因、转录本、外显子及其他功能区域的详细信息。该文件标识基因的位置、类型（如“gene”、“transcript”、“exon”）及其相关属性（如“gene_id”、“gene_name”、“transcript_id”、“transcript_name”）。这些信息对理解基因组的功能和结构至关重要。

  对于可从 Ensembl 数据库 获取的物种，建议使用该处提供的文件。Ensembl 的 GTF 文件包含可选标签，便于过滤（通过 dnbc4tools tools mkgtf）。如果 Ensembl 无法提供所需物种的文件，则可以使用其他来源的 GTF 和 FASTA 文件。请注意，GTF 文件为必需，不支持 GFF 文件。基因组文件与注释文件需对应，GTF 文件格式要求为：对于单细胞 RNA 分析，GTF 文件至少需包含“gene”或“transcript”类型以及“exon”类型的注释，并且属性中必须包含“gene_id”或“gene_name”以及“transcript_id”或“transcript_name”。

## 下载参考基因组,这里以小鼠为例
curl-O https://ftp.ensembl.org/pub/release-115/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.toplevel.fa.gz
## 下载注释文件
curl-O https://ftp.ensembl.org/pub/release-115/gtf/mus_musculus/Mus_musculus.GRCm39.115.gtf.gz

## 解压
gunzip  Mus_musculus.GRCm39.dna.toplevel.fa.gz
gunzip  Mus_musculus.GRCm39.115.gtf.gz

3.2 过滤GTF文件（可选）#

从 ENSEMBL 和 UCSC 等网站下载的 GTF 文件通常包含多种基因类型的基因。选择您研究中比较感兴趣的基因类型，过滤部分基因类型可以减少重叠的基因注释。与多个基因非唯一比对的 reads 会被过滤。

dnbc4tools

提供了基因类型数量统计、基因类型过滤和校正 GTF 文件的功能。

  • 基因类型数量统计（可选）
dnbc4tools2.1.3/dnbc4tools tools mkgtf --action stat --ingtf Mus_musculus.GRCm39.115.gtf --output gtfstat.txt --type gene_biotype

在上面的命令中，需要查看 GTF 文件中的 tag 确定 type 的类型。

cat gtfstat.txt

## Type Count
## lncRNA   35034
## protein_coding   20121
## processed_pseudogene 9489
## misc_RNA 2216
## unprocessed_pseudogene   1949
## snRNA    1910
## miRNA    1879
## transcribed_unprocessed_pseudogene   1587
## transcribed_processed_pseudogene 1149
## TEC  1019
## snoRNA   942
## rRNA_pseudogene  497
## transcribed_unitary_pseudogene   201
## IG_V_pseudogene  187
## IG_V_gene    146
## TR_V_gene    107
## unitary_pseudogene   89
## TR_J_gene    79
## rRNA 53
## scaRNA   49
## IG_D_gene    37
## TR_V_pseudogene  33
## Mt_tRNA  22
## artifact 19
## IG_J_gene    18
## IG_C_gene    14
## IG_C_pseudogene  9
## ribozyme 8
## TR_C_gene    6
## sRNA 5
## TR_D_gene    5
## vault_RNA    4
## TR_J_pseudogene  4
## pseudogene   4
## IG_J_pseudogene  3
## translated_processed_pseudogene  2
## Mt_rRNA  2
## IG_pseudogene    1
  • 校正 GTF 文件（可选）

GTF 文件格式要求为：对于单细胞 RNA 分析，GTF 文件至少需包含“gene”或“transcript”类型以及“exon”类型的注释，并且属性中必须包“gene_id”或“gene_name”以及“transcript_id”或“transcript_name”。存在内容缺失的 GTF 文件会导致主分析流程无法注释而报错。

dnbc4tools2.1.3/dnbc4tools tools mkgtf --action check --ingtf Mus_musculus.GRCm39.115.gtf --output Mus_musculus.GRCm39_corrected.gtf

软件会主动填补 gene 行和 transcript 行缺失的信息，会根据 gene_id 和 gene_name 以及 transcript_id 和 transcript_name 互相填补。

警告信息提示该位置存在多个基因信息，在分析时会导致过滤去重。

# 基因类型过滤
# dnbc4tools tools mkgtf --ingtf Homo_sapiens.GRCh38_corrected.gtf --output Homo_sapiens.GRCh38_filter.gtf --type gene_biotype

命令中我们可以添加参数 include 来获取需要的 gene 类型，默认是下面列出的基因型： - protein_coding - lncRNA/lincRNA - antisense - IG_V_gene - IG_LV_gene - IG_D_gene - IG_J_gene - IG_C_gene - IG_V_pseudogene - IG_J_pseudogene - IG_C_pseudogene - TR_V_gene - TR_D_gene - TR_J_gene - TR_C_gene

例如应用以上筛选来过滤，可以使用默认选项，也可以使用如下命令，类似 Cellranger 过滤步骤：

dnbc4tools2.1.3/dnbc4tools tools mkgtf --ingtf Mus_musculus.GRCm39_corrected.gtf --output Mus_musculus.GRCm39_filtered.gtf --type gene_biotype \
--include protein_coding,lncRNA,lincRNA,antisense,IG_V_gene,IG_LV_gene,IG_J_gene,IG_C_gene,IG_V_pseudogene,IG_J_pseudogene,IG_C_pseudogene,\
TR_V_gene,TR_D_gene,TR_J_gene,TR_C_gene

这将从原始未过滤的 GTF 文件生成一个过滤的 GTF 文件。在输出文件中，其他基因类型被排除在 GTF 注释之外。

3.3 构建参考数据库#

在运行单细胞RNA定量之前，需要注释文件 GTF 和参考基因组 FASTA 来构建索引文件，用于测序 reads 的比对和注释。

dnbc4tools2.1.3/dnbc4tools rna mkref --fasta Mus_musculus.GRCm39.dna.toplevel.fa --ingtf Mus_musculus.GRCm39_filtered.gtf --species Mus_musculus --threads 30

看一下输出的参考基因组文件结构

ls GRCm39

## _bt_log
## chrLength.txt
## chrNameLength.txt
## chrName.txt
## chrStart.txt
## exonGeTrInfo.tab
## exonInfo.tab
## geneInfo.tab
## Genome
## genomeParameters.txt
## Log.out
## mtgene.list
## ref.json
## SA
## SAindex
## sjdbInfo.txt
## sjdbList.fromGTF.out.tab
## sjdbList.out.tab
## transcriptInfo.tab

3.4 准备Fastq文件#

  • 单样本：直接准备好对应的 cDNAfastq 文件和 oligofastq 文件即可
  • 多样本（可选）：为了简化每个样本单独生成主分析流程，可以使用配置文件来生成一个包含多个样本的主流程 shell 脚本
# dnbc4tools rna multi --list sample.tsv --genomeDir /opt/database/Homo_sapiens --threads 10

其中sample.tsv文件使用制表符 ( 分隔符。第一列包含样本名称，第二列包含 cDNA 文库测序数据，第三列包含寡核苷酸文库测序数据。多个 fastq 文件应以逗号分隔，R1 和 R2 文件应以分号分隔。

# sample1 /data/cDNA1_R1.fq.gz;/data/cDNA1_R2.fq.gz /data/oligo1_R1.fq.gz,/data/oligo4_R1.fq.gz;/data/oligo1_R2.fq.gz,/data/oligo4_R2.fq.gz 
# sample2 /data/cDNA2_R1.fq.gz;/data/cDNA2_R2.fq.gz /data/oligo2_R1.fq.gz;/data/oligo2_R2.fq.gz 
# sample3 /data/cDNA3_R1.fq.gz;/data/cDNA3_R2.fq.gz /data/oligo3_R1.fq.gz;/data/oligo3_R2.fq.gz

运行完成后输出：

# sample1.sh
# sample2.sh
# sample3.sh

其中文件 sample1.sh 如下：

# cat sample1.sh

# dnbc4tools2.1.3/dnbc4tools rna run --name sample1 --cDNAfastq1 /data/cDNA1_R1.fq.gz --cDNAfastq2 /data/cDNA1_R2.fq.gz --oligofastq1 /data/oligo1_R1.fq.gz,/data/oligo4_R1.fq.gz --oligofastq2 /data/oligo1_R2.fq.gz,/data/oligo4_R2.fq.gz --genomeDir /database/scRNA/Mus_musculus/mm10 --threads 10

3.5 定量分析流程#

RNA 主分析流程:处理单个样本单细胞 RNA 的 cDNA 和 oligo 文库测序数据。该流程包括质量控制、比对和功能区域注释。随后，系统将合并磁珠以识别细胞，并生成原始基因表达矩阵及过滤后的基因表达矩阵。接下来，分析将对该矩阵进行细胞过滤、降维、聚类和注释，最终生成 HTML 格式的报告并输出分析结果。

PS:这一步要注意参考基因组中ref.json文件，其中 genomeDir 指定了参考基因组路径

dnbc4tools2.1.3/dnbc4tools rna run \
--name test \
--cDNAfastq1 data/E150004608_L01_21_1.fq.gz \
--cDNAfastq2 data/E150004608_L01_21_2.fq.gz \
--oligofastq1 data/E150004608_L01_12_1.fq.gz \
--oligofastq2 data/E150004608_L01_12_2.fq.gz \
--genomeDir GRCm39 \
--threads 30

成功的运行会以Analysis Finished结束。

# 结果文件目录如示
ls test

## 01.data
## 02.count
## 03.analysis
## 04.report
## log
## output

4. 下游文件的读取#

4.1 R使用Seurat包读取#

## 直接使用 Read10X 函数读取
library(Seurat)

## Loading required package: SeuratObject

## Loading required package: sp

## 'SeuratObject' was built with package 'Matrix' 1.7.1 but the current
## version is 1.7.2; it is recomended that you reinstall 'SeuratObject' as
## the ABI for 'Matrix' may have changed

## 
## Attaching package: 'SeuratObject'

## The following objects are masked from 'package:base':
## 
##     intersect, t

counts.data <-Read10X(data.dir ="./test/output/filter_matrix",gene.column =1)

counts.data[1:5, 1:5]

## 5 x 5 sparse Matrix of class "dgCMatrix"
##         CELL1_N3 CELL4_N4 CELL9_N2 CELL11_N2 CELL25_N3
## Xkr4       1        .        1         5         .
## Gm1992      .        .        .         .         .
## Gm71132     .        .        .         .         .
## Gm19938     .        .        .         1         .
## Gm37381     .        .        .         .         .

# 使用计数矩阵创建 Seurat 对象
# 'counts' 参数接收您的稀疏矩阵 (dgCMatrix)
seurat_object <-CreateSeuratObject(
counts = counts.data, 
project ="SparseMatrixProject", 
min.cells =3,           # 质量控制：保留至少在 3 个细胞中表达的基因
min.features =200# 质量控制：保留表达至少 200 个基因的细胞
)

# 查看创建的对象摘要
print(seurat_object)

## An object of class Seurat 
## 40706 features across 8350 samples within 1 assay 
## Active assay: RNA (40706 features, 0 variable features)
##  1 layer present: counts

这样我们就生成了对应的Seurat单细胞转录组对象。

4.2 R自定义函数读取多种类型结果文件#

## 使用函数进行读取。(适用于 raw_matrix、filter_matrix 和 RNAvelocity_matrix 等)

ReadMatrix_C4 <-function(mex_dir=NULL,
barcode.path =NULL,
feature.path =NULL,
matrix.path=NULL,
use_10X=FALSE) {
# ------------------------------------------------------------------
# 1. 参数检查和路径设置
# ------------------------------------------------------------------

# 检查是否提供了任何矩阵路径或目录
if (is.null(mex_dir) &&is.null(barcode.path)  &&is.null(feature.path) &&
is.null(matrix.path)) {
stop("No matrix set.") # 如果都没有，停止执行并报错
  }

# 如果提供了 MEX 目录，检查它是否存在
if (!is.null(mex_dir) &&!file.exists(mex_dir) ) {
stop(paste0(mex_dir, " does not exist.")) # 如果目录不存在，停止执行并报错
  }

# 如果没有提供独立的路径（barcode, feature, matrix），则假定它们在 mex_dir 中
# 这是标准的10x Genomics/Cell Ranger MEX 格式的默认文件名
if (is.null(barcode.path)  &&is.null(feature.path) &&is.null(matrix.path)) {
    barcode.path <-paste0(mex_dir, "/barcodes.tsv.gz")
    feature.path <-paste0(mex_dir, "/features.tsv.gz") # 10x v3 格式使用 features.tsv.gz
    matrix.path <-paste0(mex_dir, "/matrix.mtx.gz")
  }

# 定义包含 Spliced/Unspliced/Spanning 计数的文件的路径
# 这通常用于 RNA velocity 或更精细的转录本分析
  spliced.path <-paste0(mex_dir, "/spliced.mtx.gz")
  unspliced.path <-paste0(mex_dir, "/unspliced.mtx.gz")
  spanning.path <-paste0(mex_dir, "/spanning.mtx.gz")

# 检查必须存在的条形码和特征文件
if (!file.exists(barcode.path) ||!file.exists(feature.path)) {
stop(paste0("No expression file found at ", mex_dir)) # 如果缺失，停止执行并报错
  }

# ------------------------------------------------------------------
# 2. 内部函数：读取单个矩阵（标准 MEX 格式）
# ------------------------------------------------------------------

# 定义一个内部函数来读取标准的 .mtx, features.tsv, barcodes.tsv 格式
  .ReadPISA0 <-function(barcode.path, feature.path, matrix.path, use_10X) {
# 1. 读取稀疏矩阵数据
    mat <- Matrix::readMM(file = matrix.path)

# 2. 读取特征（基因/转录本）名称
    feature.names <-read.delim(feature.path,
header =FALSE,
stringsAsFactors =FALSE
    )
# 3. 读取条形码（细胞）名称
    barcode.names <-read.delim(barcode.path,
header =FALSE,
stringsAsFactors =FALSE
    )

# 4. 设置矩阵的列名（Cell Barcodes）
colnames(mat) <- barcode.names$V1

# 5. 设置矩阵的行名（Features/Genes）
if (use_10X ==TRUE) {
# 10X Genomics V3 格式的特征文件通常 V2 列包含基因符号/名称
rownames(mat) <-make.unique(feature.names$V2)
    } else {
# 其他格式（如 V1 格式或一些预处理后的数据）可能 V1 列包含基因 ID/名称
rownames(mat) <-make.unique(feature.names$V1)
    }

# 返回构建好的矩阵
    mat
  }

# ------------------------------------------------------------------
# 3. 矩阵读取逻辑：单矩阵 vs. Spliced/Unspliced/Spanning 矩阵
# ------------------------------------------------------------------

# 优先检查是否为标准单矩阵文件
# 如果 'spliced.path' 不存在，但 'matrix.path' 存在，则读取并返回单个表达矩阵
if (!file.exists(spliced.path) &&file.exists(matrix.path)) {
return(.ReadPISA0(barcode.path, feature.path, matrix.path, use_10X))
  }

# 否则，假设存在 Spliced/Unspliced/Spanning 矩阵，并准备读取它们
  mat <-list() # 初始化一个列表来存储多个矩阵

# 读取 Spliced 矩阵
cat("Load spliced matrix ...\n")
  mat$spliced <- Matrix::readMM(file = spliced.path)

# 读取 Unspliced 矩阵
cat("Load unspliced matrix ...\n")
  mat$unspliced <- Matrix::readMM(file = unspliced.path)

# 读取 Spanning 矩阵
cat("Load spanning matrix ...\n")
  mat$spanning <- Matrix::readMM(file = spanning.path)
}

# 测试读取文件
test <-ReadMatrix_C4(mex_dir='./test/output/raw_matrix')

test[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         CELL1_N3 CELL2_N1 CELL3_N2 CELL4_N4 CELL5_N2
## Xkr4           1        .        1        .        .
## Gm1992         .        .        .        .        .
## Gm71132        .        .        .        .        .
## Gm19938        .        .        .        .        .
## Gm37381        .        .        .        .        .

# 使用计数矩阵创建 Seurat 对象
# 'counts' 参数接收您的稀疏矩阵 (dgCMatrix)
seurat_object <-CreateSeuratObject(
counts = test, 
project ="SparseMatrixProject", 
min.cells =3,           # 质量控制：保留至少在 3 个细胞中表达的基因
min.features =200# 质量控制：保留表达至少 200 个基因的细胞
)

## Warning: Data is of class dgTMatrix. Coercing to dgCMatrix.

# 查看创建的对象摘要
print(seurat_object)

## An object of class Seurat 
## 41274 features across 47999 samples within 1 assay 
## Active assay: RNA (41274 features, 0 variable features)
##  1 layer present: counts

4.3 Python使用Scanpy读取#

## 直接使用scanpy提供read_h5ad函数读取
import scanpy as sc
adata = sc.read_h5ad('./test/output/filter_feature.h5ad')
print(adata)

## AnnData object with n_obs × n_vars = 8350 × 45577
##     var: 'gene_symbols'

这样我们就生成了对应的

Anndata

格式的单细胞转录组对象。

4.4 Python自定义函数读取#

import pandas as pd
import scipy.io
import anndata
from scipy.sparse import csr_matrix
import os

def read_anndata_C4(path):
"""
    读取 10x Genomics MEX (Matrix Exchange) 格式的数据，并将其转换为 AnnData 对象。

    该函数假定 'path' 目录下包含三个标准文件：
    - matrix.mtx.gz (稀疏矩阵数据)
    - features.tsv.gz (特征/基因信息)
    - barcodes.tsv.gz (细胞条形码信息)

    Args:
        path (str): 包含 MEX 文件的目录路径。

    Returns:
        anndata.AnnData: 包含加载数据的 AnnData 对象。
    """

# 1. 读取稀疏矩阵数据 (matrix.mtx.gz)
# scipy.io.mmread 用于读取 Matrix Market 格式的稀疏矩阵文件 (.mtx)
# 默认 dtype 是 float64，将其转换为 float32 以节省内存
    mat = scipy.io.mmread(os.path.join(path, "matrix.mtx.gz")).astype("float32")

# 2. 矩阵转置 (Transpose)
# Matrix Market 格式的约定是 (Feature x Cell)，但 AnnData 要求 (Cell x Feature)
    mat = mat.transpose()

# 3. 转换为 CSR 稀疏矩阵格式
# CSR (Compressed Sparse Row) 格式是 AnnData 用于高效矩阵操作的首选格式
    mat = csr_matrix(mat)

# 4. 创建 AnnData 对象
# 使用转置并转换后的稀疏矩阵来初始化 AnnData 对象
# dtype="float32" 再次确保数据类型
    adata = anndata.AnnData(mat, dtype="float32")

# 5. 读取并设置特征/基因信息 (var)
# features.tsv.gz 文件通常包含基因ID和/或基因名称
    genes = pd.read_csv(os.path.join(path, "features.tsv.gz"), header=None, sep='\t')

# 将第一列（V1/0）的值作为 AnnData 的 var_names (特征名称)
    var_names = genes[0].values

# 确保特征名称是唯一的，这是 AnnData 的最佳实践
    var_names = anndata.utils.make_index_unique(pd.Index(var_names))
    adata.var_names = var_names

# 将基因 ID/符号作为 AnnData.var 的一个新列 'gene_symbols' 存储起来
# 即使 var_names 被唯一化了，原始的 ID 仍然有用
    adata.var['gene_symbols'] = genes[0].values

# 6. 读取并设置细胞条形码信息 (obs)
# barcodes.tsv.gz 文件包含细胞的唯一条形码
# 将第一列（V1/0）的值作为 AnnData 的 obs_names (观测名称，即细胞名称)
    adata.obs_names = pd.read_csv(os.path.join(path, "barcodes.tsv.gz"), header=None)[0].values

# 7. 最终检查和确保 var_names 唯一性 (尽管前面已经做过，这是一个安全步骤)
    adata.var_names_make_unique()

# 8. 返回构建好的 AnnData 对象
return adata

## 测试一下
test = read_anndata_C4('./test/output/filter_matrix')

## /home/chenao/miniconda3/envs/cpdb/lib/python3.9/site-packages/anndata/_core/anndata.py:402: FutureWarning: The dtype argument is deprecated and will be removed in late 2024.
##   warnings.warn(

test

## AnnData object with n_obs × n_vars = 8350 × 45577
##     var: 'gene_symbols'

完成读取之后就可以正常进行下游的个性化分析，我们也给大家准备了详细的教程： 小白级单细胞分析全教程

图片

四、参考

华大官方文档（https://bgitechsolutions.com/sequencing/244）
华大B站视频介绍 dnbc4tools_github（https://www.bilibili.com/video/av958683573/?vd\_source=038c940ddc24688e3c4a0eb678552cc1）

图片

五、演示环境

一切不给测试文件和分析环境版本的教程都是耍流氓，本推送的代码和测试文件可以在以下链接中下载： 图片

通过网盘分享的文件:

https://pan.baidu.com/s/1bkf-YB1WCjm6gCRYREKr0w

提取码: yckd

图片

R语言演示环境(趁开学前| 入手足够完成硕博生涯的生信环境)：

sessionInfo()

## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3;  LAPACK version 3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Asia/Shanghai
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] sceasy_0.0.7           reticulate_1.43.0.9001 Seurat_5.2.1          
## [4] SeuratObject_5.0.2     sp_2.2-0              
## 
## loaded via a namespace (and not attached):
##   [1] deldir_2.0-4           pbapply_1.7-2          gridExtra_2.3 
##   [4] rlang_1.1.5            magrittr_2.0.3       RcppAnnoy_0.0.22 
##   [7] spatstat.geom_3.4-1    matrixStats_1.5.0      ggridges_0.5.6   
##  [10] compiler_4.4.2         png_0.1-8              vctrs_0.6.5 
##  [13] reshape2_1.4.4         stringr_1.5.1          pkgconfig_2.0.3 
##  [16] fastmap_1.2.0          labeling_0.4.3         promises_1.3.2 
##  [19] rmarkdown_2.29         purrr_1.0.2            xfun_0.50   
##  [22] cachem_1.1.0           jsonlite_1.8.9         goftest_1.2-3 
##  [25] later_1.4.1         spatstat.utils_3.1-4   irlba_2.3.5.1  
##  [28] parallel_4.4.2         cluster_2.1.8          R6_2.5.1      
##  [31] ica_1.0-3           bslib_0.9.0            stringi_1.8.4  
##  [34] RColorBrewer_1.1-3     spatstat.data_3.1-6    parallelly_1.42.0   
##  [37] spatstat.univar_3.1-3  lmtest_0.9-40          jquerylib_0.1.4   
##  [40] scattermore_1.2        Rcpp_1.0.14            knitr_1.49   
##  [43] tensor_1.5             future.apply_1.11.3    zoo_1.8-12     
##  [46] sctransform_0.4.1      httpuv_1.6.15          Matrix_1.7-2   
##  [49] splines_4.4.2          igraph_2.1.4           tidyselect_1.2.1      
##  [52] abind_1.4-8            rstudioapi_0.17.1      yaml_2.3.10   
##  [55] spatstat.random_3.4-1  codetools_0.2-20       miniUI_0.1.1.1     
##  [58] spatstat.explore_3.4-3 listenv_0.9.1          lattice_0.22-6   
##  [61] tibble_3.2.1           plyr_1.8.9             withr_3.0.2 
##  [64] shiny_1.10.0           ROCR_1.0-11            evaluate_1.0.3        
##  [67] Rtsne_0.17             future_1.34.0          fastDummies_1.7.5     
##  [70] survival_3.8-3         polyclip_1.10-7        fitdistrplus_1.2-2    
##  [73] pillar_1.10.1          KernSmooth_2.23-26     plotly_4.10.4   
##  [76] generics_0.1.3         RcppHNSW_0.6.0         ggplot2_3.5.1   
##  [79] munsell_0.5.1          scales_1.3.0           globals_0.16.3 
##  [82] xtable_1.8-4           glue_1.8.0             lazyeval_0.2.2        
##  [85] tools_4.4.2            data.table_1.16.4      RSpectra_0.16-2 
##  [88] RANN_2.6.2          dotCall64_1.2          cowplot_1.1.3   
##  [91] grid_4.4.2             tidyr_1.3.1            colorspace_2.1-1      
##  [94] nlme_3.1-168           patchwork_1.3.0        cli_3.6.3             
##  [97] spatstat.sparse_3.1-0  spam_2.11-1            viridisLite_0.4.2     
## [100] dplyr_1.1.4            uwot_0.2.2             gtable_0.3.6          
## [103] sass_0.4.9             digest_0.6.37          progressr_0.15.1      
## [106] ggrepel_0.9.6          htmlwidgets_1.6.4      farver_2.1.2          
## [109] htmltools_0.5.8.1      lifecycle_1.0.4        httr_1.4.7            
## [112] mime_0.12              MASS_7.3-64

图片

欢迎致谢

如果以上内容对你有帮助，欢迎在文章的Acknowledgement中加上这一段，联系客服微信可以发放奖励：

Since Biomamba and his wechat public account team produce bioinformatics tutorials and share code with annotation, we thank Biomamba for their guidance in bioinformatics and data analysis for the current study.

欢迎在发文/毕业时向我们分享你的喜悦~

已致谢文章：

KIF18A在肝细胞癌转移中的双重角色

《Advanced Science》新型Arf1抑制剂促进癌症干细胞衰老并增强抗肿瘤免疫

13分+文章利用scRNA-Seq揭示地铁细颗粒物引起肺部炎症的分子机制
IF14.3| scRNA-seq+脂质组多组学分析揭示宫内生长受限导致肝损伤的性别差异
鼻咽癌的Bulk RNA-Seq与scRNA-Seq联合分析
除了铁死亡，还有铜死亡？！
银屑病和脂肪肝病中共同病理和免疫特征

scRNA-seq揭示脓毒症预后水平预测的关键靶点！

机器学习+生信多组学联合构建牙周炎"线粒体功能障碍与免疫微环境"关联网络

bulk+scRNA-seq挖掘BCL2-MAPK14-TXN氧化应激诊断模型，鉴定脓毒症中氧化应激关键基因

如何联系我们

图片

开箱即用的单细胞分析镜像 的 微信号[Biomamba_zhushou] ，方便各位随时交流。同时我们也构建了 交流群矩阵 ，欢迎大家入群讨论。

大家可以阅读完这几篇之后添加

给生信入门初学者的小贴士

没有检索，就没有发言权

已有生信基地联系方式的同学 无需重复添加

图片

您点的每个赞和在看，我都认真当成了喜欢

微信扫一扫赞赏作者

上游定量 · 目录

闪记

复制 LaTeX 公式

kimi