DNBC4tools—华大DNBelab系列单细胞分析pipeline

生信菜鸟团

发布于 2025-05-21 15:23:58

4.2K0

工欲善其事必先利其器

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偏好性和重复序列错误，显著提高数据均一性。

源自https://www.mgi-tech.com/products/resources

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

文库结构示意图

DNBC4tools

一个用于分析高通量 DNBelab C 系列 TM 单细胞数据集的开源且灵活的流程。

Github：

https://github.com/MGI-tech-bioinformatics/DNBelab_C_Series_HT_scRNA-analysis-software

如何安装

目前最新版本的软件，安装非常简单，下载解压即可直接使用。

##下载可能比较耗时，可以放到screen窗口（可选）
screen -R wget

cd ~/software
##wget 下载（必选）
wget -O dnbc4tools2.1.3.tar.gz "https://ftp.cngb.org/pub/CNSA/data5/CNP0006367/Single_Cell/CSE0000448/dnbc4tools2.1.3.tar.gz"

#解压
tar -xf dnbc4tools2.1.3.tar.gz

##添加到环境变量（可选）
ln -s ~/software/dnbc4tools2.1.3/dnbc4tools ~/software/bin/dnbc4tools

下载

解压直接使用

主要功能

单细胞RNA分析—- dnbc4tools rna run 使用单细胞 RNA 的 cDNA 和 oligo 文库测序数据，进行质量控制、比对和功能区域注释。随后，合并磁珠并识别细胞生成原始基因表达矩阵，识别细胞生成过滤后的基因表达矩阵。接着，对该矩阵进行细胞过滤、降维、聚类和注释分析，最终生成 HTML 格式的网页报告并输出分析结果。
单细胞ATAC分析— dnbc4tools atac run 使用单细胞 ATAC 文库测序数据，经过过滤和比对生成所有磁珠的 fragments 文件。合并磁珠并执行 peak 调用分析，利用 peaks 区域的片段信息进行细胞识别。随后进行细胞过滤、降维和聚类，最终整合各步骤结果生成 HTML 网页报告并输出分析结果。
单细胞VDJ分析— dnbc4tools vdj run 使用单细胞 VDJ 文库测序数据和对应样本的 5' 转录组分析结果。首先对数据进行过滤，利用 5' 转录组结果合并磁珠，接着比对 VDJ 基因区域并提取对应的 reads，进行从头组装并注释。根据组装注释结果和 5' 转录组的细胞获取情况进行细胞过滤，最后整合各步骤结果生成 HTML 网页报告并输出分析结果。

实例演示

至于如何使用，上述流程基本上都是需要

先构建参考数据库
然后准备输入文件
最后写脚本，提交运行

这里我们以单细胞转录组数据分析为例：

构建参考数据库

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

ensembl 参考基因组

https://www.ensembl.org/info/about/species.html
https://asia.ensembl.org/info/data/ftp/index.html

homo sapiens

#homo spaies比较新的是release-114，大家可以自己按需下载
https://ftp.ensembl.org/pub/release-114/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.gz
https://ftp.ensembl.org/pub/release-114/gtf/homo_sapiens/Homo_sapiens.GRCh38.114.gtf.gz

##这里就以服务器已有参考基因组release-112为例
cd ~/reference/human/ensembl 
1.5G 2月  20  2024 Homo_sapiens.GRCh38.112.gtf
3.2G 2月  14  2024 Homo_sapiens.GRCh38.dna.toplevel.fa

基因类型数量统计（可选）

##基因类型数量统计（可选）
mkdir ~/reference/human/homo_ensembl_112_dnbc4_index
cd ~/reference/human/ensembl 

dnbc4tools tools mkgtf --action stat --ingtf Homo_sapiens.GRCh38.112.gtf --output ../homo_ensembl_112_dnbc4_index/gtf_GRCh38_112_stat.txt --type gene_biotype

--action stat #计算每种基因类型的基因数量
--ingtf #gtf注释文件
--output  #输出文件
--type #GTF文件中基因类型标签设置，默认：gene_biotype

基因类型数量统计

这仅仅是一个统计信息，对后续分析没有影响。统计后可以对参考基因组有个大致的了解。【自己也可以用来对比一下不同版本的参考基因组信息的区别】

校正gtf文件（可选）

dnbc4tools tools mkgtf --action check --ingtf ./Homo_sapiens.GRCh38.112.gtf --output ../homo_ensembl_112_dnbc4_index/Homo_sapiens.GRCh38.112_corrected.gtf

--action check #检查GTF文件是否符合分析要求，并生成一个新的GTF文件

校正GTF文件

GTF 文件格式要求为：对于单细胞 RNA 分析，GTF 文件至少需包含“gene”或“transcript”类型以及“exon”类型的注释，并且属性中必须包含“gene_id”或“gene_name”以及“transcript_id”或“transcript_name”。存在内容缺失的 GTF 文件会导致主分析流程无法注释而报错。
软件会主动填补 gene 行和 transcript 行缺失的信息，会根据 gene_id 和 gene_name 以及 transcript_id 和 transcript_name 互相填补。
警告信息提示该位置存在多个基因信息，在分析时会导致过滤。

基因类型过滤

cd ~/reference/human/homo_ensembl_112_dnbc4_index/

dnbc4tools tools mkgtf --ingtf ./Homo_sapiens.GRCh38.112_corrected.gtf --output ./Homo_sapiens.GRCh38.112_corrected_filter.gtf --type gene_biotype

基因类型过滤

构建索引文件

这一步，比较耗时，建议放到后台处理。【人类参考基因构建索引，6个线程，耗时1个多小时】

screen -R mkref

dnbc4tools rna mkref --ingtf ./Homo_sapiens.GRCh38.112_corrected_filter.gtf --fasta ../ensembl/Homo_sapiens.GRCh38.dna.toplevel.fa --threads 6 --species Homo_sapiens_GRCh38.112

--ingtf #gtf基因组注释文件
--fasta #FASTA格式的基因组文件
--threads #调用线程，默认：4
--species #指定构建参考数据库的物种名称

构建参考基因组索引

单细胞转录组定量

基本用法：

##单样本
dnbc4tools rna run \
  --name sample \
  --cDNAfastq1 /data/sample_cDNA_R1.fastq.gz \
  --cDNAfastq2 /data/sample_cDNA_R2.fastq.gz \
  --oligofastq1 /data/sample_oligo1_1.fq.gz,/data/sample_oligo2_1.fq.gz \
  --oligofastq2 /data/sample_oligo1_2.fq.gz,/data/sample_oligo2_2.fq.gz \
  --genomeDir /opt/database/Homo_sapiens \
  --threads 10

##多样本批量生成运行脚本
dnbc4tools rna multi --list sample.tsv \
  --genomeDir /opt/database/Homo_sapiens \
  --threads 10


--name  #样本ID
--cDNAfastq1 #cDNA文库 R1 FASTQ文件及路径
--cDNAfastq2 #cDNA文库 R2 FASTQ文件及路径
--oligofastq1 #oligo文库 R1 FASTQ文件及路径
--oligofastq2 #oligo文库 R2 FASTQ文件及路径
--genomeDir #基因组文件所在目录的路径
--threads #调用线程

准备样本数据

如果是单样本的话，直接按照上述单样本的程序代码，指定输入输出，提交后台运行即可。
如果是多样本，并且一个样本对应多个fastq文件的话，处理起来相对就复杂一点。可能需要先对样本数据重命名，再批量生成运行脚本。

多样本处理首先需要制作一个自己的样本信息对应列表sample.tsv ：

比如这样的信息表：

第一列是样本名称
第二列是 cDNA 文库测序数据，多个 fastq 文件以逗号分隔，R1 和 R2 文件以分号分隔。
第三列是寡核苷酸文库测序数据。多个 fastq 文件以逗号分隔，R1 和 R2 文件以分号分隔。

$cat sample.tsv 
D2-1    D2-1_L001_f1.fq.gz,D2-1_L002_f1.fq.gz;D2-1_L001_r2.fq.gz,D2-1_L002_r2.fq.gz     D2-1_oligo_f1.fq.gz;D2-1_oligo_r2.fq.gz
D2-2    D2-2_L001_f1.fq.gz,D2-2_L002_f1.fq.gz;D2-2_L001_r2.fq.gz,D2-2_L002_r2.fq.gz     D2-2_oligo_f1.fq.gz;D2-2_oligo_r2.fq.gz
D5-1    D5-1_L001_f1.fq.gz,D5-1_L002_f1.fq.gz;D5-1_L001_r2.fq.gz,D5-1_L002_r2.fq.gz     D5-1_oligo_f1.fq.gz;D5-1_oligo_r2.fq.gz
D5-2    D5-2_L001_f1.fq.gz,D5-2_L002_f1.fq.gz;D5-2_L001_r2.fq.gz,D5-2_L002_r2.fq.gz     D5-2_oligo_f1.fq.gz;D5-2_oligo_r2.fq.gz
D8-1    D8-1_L001_f1.fq.gz,D8-1_L002_f1.fq.gz;D8-1_L001_r2.fq.gz,D8-1_L002_r2.fq.gz     D8-1_oligo_f1.fq.gz;D8-1_oligo_r2.fq.gz
D8-2    D8-2_L001_f1.fq.gz,D8-2_L002_f1.fq.gz;D8-2_L001_r2.fq.gz,D8-2_L002_r2.fq.gz     D8-2_oligo_f1.fq.gz;D8-2_oligo_r2.fq.gz

然后执行批量脚本生成，得到多样本数据，每个样本的运行脚本

##在fq文件目录执行该脚本
dnbc4tools rna multi --list sample.tsv --genomeDir ~/reference/human/homo_ensembl_112_dnbc4_index --threads 10

批量生成运行脚本

后台提交运行

所有样本的运行脚本都生成后，我们可以在后台批量提交，这时候需要考虑样本数量和服务器的资源使用情况，决定是否分批提交运行。关于提交后台以及服务器资源查看，详见

screen -R BGI

##这个分批运行代码要自己视情况修改
ls ../st2_rawdata/D*.sh|sort -V|cut -d "/" -f 3| cut -d "." -f 1 |xargs -iF -P 5 sh -c 'bash ../st2_rawdata/F.sh 1>log_F.txt 2>&1'

结果文件

正常运行的进度输出：

运行输出日志

整个运行下来可以得到下列文件：

$tree 
.
├── 01.data
│   ├── alignment_report.csv
│   ├── anno_report.csv
│   ├── beads_stat.txt
│   ├── CB_UB_count.txt
│   ├── cDNA.barcode.counts.txt
│   ├── cDNAin1
│   ├── cDNAin2
│   ├── cDNA_para
│   ├── cDNA.sequencing.report.csv
│   ├── final_sorted.bam
│   ├── Log.final.out
│   ├── Log.out
│   ├── Log.progress.out
│   ├── oligo.barcode.counts.txt
│   ├── oligo_para
│   ├── oligo.reads.fq.gz
│   ├── oligo.sequencing.report.csv
│   └── SJ.out.tab
├── 02.count
│   ├── barcodeTranslate.hex.txt
│   ├── barcodeTranslate.txt
│   ├── beads_barcodes.txt
│   ├── cellNumber.merge.png
│   ├── filter_matrix
│   │   ├── barcodes.tsv.gz
│   │   ├── features.tsv.gz
│   │   └── matrix.mtx.gz
│   ├── raw_matrix
│   │   ├── barcodes.tsv.gz
│   │   ├── features.tsv.gz
│   │   └── matrix.mtx.gz
│   ├── saturation_cDNA.png
│   ├── saturation_cDNA.xls
│   ├── similarity.all.csv
│   ├── similarity.droplet.csv
│   └── singlecell.csv
├── 03.analysis
│   ├── cluster.csv
│   ├── cluster.png
│   ├── filter_feature.h5ad
│   ├── filter_QCplot.png
│   ├── marker.csv
│   ├── QC_Clutser.h5ad
│   ├── raw_QCplot.png
│   └── raw_qc.xls
├── 04.report
│   ├── D2-1_scRNA_report.html
│   ├── div
│   │   ├── barcoderanks.div
│   │   ├── barcoderanks.html
│   │   ├── cluster.div
│   │   ├── cluster.html
│   │   ├── gene.saturation.div
│   │   ├── gene.saturation.html
│   │   ├── merge.div
│   │   ├── merge.html
│   │   ├── sequence.saturation.div
│   │   ├── sequence.saturation.html
│   │   ├── umi.cluster.div
│   │   ├── umi.cluster.html
│   │   ├── violin.summary.base64
│   │   ├── violin.summary.div
│   │   └── violin.summary.html
│   ├── metrics_summary.xls
│   └── table
│       └── marker.table.txt
├── log
│   ├── 20250520.txt
│   └── _bt_log
└── output
    ├── anno_decon_sorted.bam
    ├── anno_decon_sorted.bam.bai
    ├── attachment
    │   ├── RNAvelocity_matrix
    │   │   ├── barcodes.tsv.gz
    │   │   ├── features.tsv.gz
    │   │   ├── spanning.mtx.gz
    │   │   ├── spliced.mtx.gz
    │   │   └── unspliced.mtx.gz
    │   └── splice_matrix
    │       ├── barcodes.tsv.gz
    │       ├── features.tsv.gz
    │       └── matrix.mtx.gz
    ├── D2-1_scRNA_report.html
    ├── filter_feature.h5ad
    ├── filter_matrix
    │   ├── barcodes.tsv.gz
    │   ├── features.tsv.gz
    │   └── matrix.mtx.gz
    ├── metrics_summary.xls
    ├── raw_matrix
    │   ├── barcodes.tsv.gz
    │   ├── features.tsv.gz
    │   └── matrix.mtx.gz
    └── singlecell.csv

15 directories, 81 files

output 目录就是我们需要的结果文件。