三代测序 - 数据质控

深圳华大基因科技有限公司（华大集团）旗下公司杭州华大序风科技有限公司（CycloneSEQ-华大序风）于2024年3月21日注册成立，致力于通过纳米孔单分子测序技术的创新研发及落地应用。2024 年 9 月 9 日，在“登峰探极·生命可测”华大集团生命科学全球新品发布会上，华大集团（华大智造）发布了名为 CycloneSEQ™ 的最新测序技术，并推出 CycloneSEQ-WT02（WT02，中文名“梧桐”）和 CycloneSEQ-WY01 (WY01，中文名“五岳”) 两款纳米孔测序仪（图1）。随后，在2025年3月1日的SEQ ALL联盟年度峰会上华大智造联合华大序风宣布G400-ER（更名后的CycloneSEQ-WY01）正式开售。

图1.CycloneSEQ-WT02和CycloneSEQ-WY01真机

作为一款纳米孔测序仪，现阶段测序所得的碱基质量会普遍偏低，根据官网性能参数的介绍，两款纳米孔测序仪的单次碱基准确率在97%左右，也就是Q15。因此，对下机数据进行质量查看和质控是数据分析前重要的一步。对于三代纳米孔测序平台，查看数据统计信息和质量最常用的就针对牛津纳米孔（ONT）数据开发的Nanopack分析套装，如NanoPlot，NanoComp和NanoQC，以及老牌质控软件fastp针对三代长度长数据优化的fastplong。这些软件的使用方法在我以往的推文中都有详细的教程，大家可以自行查阅。

作为国产纳米孔测序仪，后续数据分析最理想的工具软件，是针对自家数据开发的算法，但是这需要时间和科研圈的集体贡献。华大序风专门建立了一个Github账号，CycloneSEQ Bioinformatics ，来发布他们针对CycloneSEQ下机数据分析的算法工具（图2）。现在已经有三款软件 Bamboo (数据质控），Glycine（鉴定全长cDNA）和Citrus（拆分barcode）。今天我就一起来看看Bamboo的安装和使用方法。

图2.华大序风github账号

一、软件介绍

Bamboo 是由华大序风开发的一款针对长读长测序数据的质量控制工具，支持无参（reference-free, or pre-alignment）或有参（reference-based, or post-alignment）QC，支持 fastq 或 bam 格式文件作为输入。现在最新版本为v0.2.2（2024.10.08）。作者为程瀚森（Hansen Cheng），张嘉远（Jia-Yuan Zhang），连明（Ming Lian）。

Github: https://github.com/CycloneSEQ-Bioinformatics/Bamboo

二、软件安装

Bamboo暂时还不支持conda一键安装部署。如果你的conda环境下有samtools和minimap2，这可以直接下载解压二进制文件。如果没有，则创建Bamboo的conda分析环境，然后安装samtools和minimap2。

# 创建名bamboo_env 的conda环境，安装samtools和Minimap2
$ conda create -n bamboo_env samtools minimap2

#下载软件
$ wget -c -t 0 https://github.com/CycloneSEQ-Bioinformatics/Bamboo/releases/download/v0.2.2/bamboo-0.2.2.tar.gz

#解压缩
$ tar -zxvf bamboo-0.2.2.tar.gz

#添加可执行权限,如果已经是可执行程序则可以跳过这一步
$ cd bamboo-0.2.2
$ chmod +x bamboo

#添加bamboo到环境变量，可以直接调用
$ echo 'PATH=/mnt/data/home/mli/Desktop/Software/bamboo-0.2.2:$PATH' >> ~/.bashrc && source ~/.bashrc

# 测试bamboo安装，查看帮助文档
$ bamboo --help

三、软件使用

根据不同的输入组合Bamboo可以进行不同的分析模式（图3），即：

Fastq only：无参QC
Bam & ref.fa：有参QC
Fastq & ref.fa：无参&有参QC

图3.Bamboo不同分析模式流程图

1. 基本使用

无参QC模式

若只提供fastq文件，则执行无参QC，命令示例如下：

$ bamboo -i input1.fastq input2.fastq -o output_dir --sample_size 50000 --seed 123

其中，

-i/--sequence_path：指定fastq文件（可以是纯文本或gz文件）路径，可以给定一个或多个文件，文件之间用空格分隔开。注意：无论给定的是一个还是多个fastq文件，最终都会被合并成一个文件进行处理；多样本的话就分开跑。
-o/--output_dir：指定输出目录，会在该目录下创建一个名叫sequence的子目录来保存分析结果；
--sample_size：指定reads抽样量，默认为1e5；也可以设定为-1以不执行抽样，分析所有reads；
--seed：指定随机抽样过程的随机数种子；

注意：无参QC无并行模式，故设置 -t 参数无效。

有参QC模式

若同时提供bam文件和ref.fa，则执行有参QC，命令示例如下：

$ bamboo -b input.bam -r reference.fasta -o output_dir --realign

其中，

-b/--bam_path：指定bam文件路径，要求输入的bam文件满足两个条件：1）CIGAR列中以=或X区分展示M状态；2）不包含secondary alignment。若不满足其中的任意一个条件，请添加--realign参数以调用minimap2进行重新比对；
-r/--reference_path：指定参考序列的fasta文件（可以是纯文本或gz文件）路径；
-o/--output_dir：指定输出目录，会在该目录下创建一个名叫alignment的子目录来保存分析结果；
--realign：调用minimap2进行重新比对；
-t/--threads：指定有参QC的并行线程数。注意，该参数仅在执行有参QC时有效。

无参&有参QC模式

若同时提供fastq文件和ref.fa，则执行无参&有参QC，命令示例如下：

$ bamboo \
-i input.fastq.gz -r reference.fasta.gz \
-o test_bamboo_fastqtobam2 \
--sample_size 10000 \
-t 4 \
--align_all \
--minimap2_args "-ax map-ont --eqx --secondary=no"

其中，

-o/--output_dir：指定输出目录，会在该目录下创建两个分别名叫sequence和alignment的子目录来保存分析结果；
-t/--threads：指定有参QC的并行线程数。注意，该参数仅在执行有参QC时有效。该参数的设定会对minimap2运行产生影响，即若在--minimap2_args未设定线程数，则会将其设为与-t的值一致；
--align_all：在执行有参QC之前，会调用minimap2执行参考基因组比对，默认情况下使用抽样后的fastq文件进行比对，若设置该参数，则会使用原始fastq文件进行比对；
--minimap2_args：设置调用minimap2执行参考基因组比对时的参数组合，默认为“-ax map-ont --eqx --secondary=no。

2.结果说明/解读

分析结束后可以获得一个html网页报告，你可以下载到本地后，用浏览器打开

分别会得到以下结果但不限以下结果：

1）无参分析Reference-free QC：

Sequencing summary (统计总结)：以表格形式展示从fastq文件统计得到的全局测序数据信息，包括Total reads（总序列条数）、Total bases（总碱基数）、Average read quality （平均碱基质量）Read N50 （序列N50）和Average read length（平均序列长度）等（图4）。

图4.Sequencing summary (统计总结)

Sequence length (序列长度)：以图形可视化展示从fastq文件统计得到的read长度分布信息（图5）；这张图展示了read长度的分布情况，横坐标表示read长度，纵坐标表示具有该长度的read数量。

图5.Sequence length (序列长度)

Sequence quality：以图形可视化展示从fastq文件统计得到的read/base测序质量分布信息（图6）；此图展示了reads质量值的分布情况。横坐标表示reads的质量值，纵坐标表示具有该质量值的reads数量。

图6.Sequence quality (序列质量)

Sequence content：以图形可视化展示从fastq文件统计得到的read的碱基组成信息，包括A/T/C/G碱基组成（图7）、GC含量（图8）等。图7展示了reads首尾或5’/3’端1000 bp内，4种碱基（ATCG）的占比。x轴表示从5'端或3’端起始的相对位置，y轴表示各碱基的比例。图8展示了每个reads的GC含量分布情况。横坐标表示GC含量百分比，纵坐标表示具有该GC含量的reads数量。

图7. A/T/C/G碱基组成

图8.GC含量

2）有参分析Reference-based QC：

Mapping summary(比对总结)：以表格形式展示从bam文件得到的全局比对信息(图9），包括Total reads(总序列数）、Unmapped reads(未比对上的序列)等。

图9.Mapping summary(比对总结)

Coverage and bias(覆盖度和偏好)：以图形可视化展示从bam文件统计得到的基因组覆盖信息，包括全基因组覆盖度（图10）、覆盖度分布等；图10展示了基因组各染色体的覆盖情况，每条染色体占一行，横坐标是染色体位置（以Mb为单位），纵坐标代表该位置的测序深度。图中使用颜色区分不同染色体，以帮助识别高覆盖和低覆盖区域，评估数据的完整性和均匀性。

图10.全基因组覆盖度

Sequence accuracy(序列准确度)：以图形可视化展示从bam文件统计得到的read准确率/一致性信息，包括read准确率分布（图11）、总体错误水平等；图11展示了read的准确率分布情况。read的准确率基于参考基因组的比对结果计算得到，即match bases/read align length。横坐标是read比对到参考基因组的一致性百分比，纵坐标是具有该一致性的reads比例。此分析用以帮助用户评估测序数据的准确性。

图11.Sequence accuracy(序列准确度)

Low complexity regions(低复杂度区域)：以图形可视化展示从bam文件统计得到的低复杂度区域，即homo/heteropolymer区域的测序错误模式（图12），以及它们对测序错误的贡献度。图12展示了不同长度的homopolymer区域在测序结果中发生的长度膨胀（expansion）或收缩（contraction）情况的相对误差分布。对于基因组中(X)n的homopolymer区域，若在测序结果中仍然保持(X)m的形式，如果长度一致，即m=n，则测序正确，若长度发生改变，即m≠n，则我们将这种情况成为发生expansion（m>n）或contraction（m<n）错误，除此之外的情况，都归为其他复杂情况（complex）。横坐标是参考序列中homopolymer区域的长度，纵坐标不同色带的宽度为不同错误情况的占比。左图是不区homopolymer区域碱基组成的总体错误情况，右图则是细分各个homopolymer区域碱基组成下，各自错误情况。

图12.homo/heteropolymer区域的测序错误模式

四、帮助文档

usage: bamboo [-h] [-b BAM_PATH] [-r REFERENCE_PATH] [--realign] [--minimap2_path MINIMAP2_PATH] [--minimap2_args MINIMAP2_ARGS] [--samtools_path SAMTOOLS_PATH]
              [--align_all] [-i SEQUENCE_PATH [SEQUENCE_PATH ...]] [-o OUTPUT_DIR] [-t THREADS] [--sample_size SAMPLE_SIZE] [--seed SEED] [--keep-intermediates]

Bamboo v0.2.0: a tool for quality control and error profiling of long-read sequencing data.

optional arguments:
  -h, --help            show this help message and exit

Sequence analyses:
  Arguments for sequence analyses.

  -i SEQUENCE_PATH [SEQUENCE_PATH ...], --sequence_path SEQUENCE_PATH [SEQUENCE_PATH ...]
                        Path to the input FASTQ file. If multiple input files are supplied, they will be concatenated before analyses. (default: None)

Alignment analyses:
  Arguments for alignment analyses.

  -b BAM_PATH, --bam_path BAM_PATH
                        Path to the input BAM file. (default: None)
  -r REFERENCE_PATH, --reference_path REFERENCE_PATH
                        Path to the reference FASTA file. (default: None)
  --realign             Re-align sampled reads using Minimap2. Use this option if the input BAM file does not contain x/= CIGAR operations. (default: False)
  --minimap2_path MINIMAP2_PATH
                        Path to Minimap2 executable. (default: minimap2)
  --minimap2_args MINIMAP2_ARGS
                        Command line arguments for Minimap2 (default: -ax map-ont --eqx --secondary=no -t 8)
  --samtools_path SAMTOOLS_PATH
                        Path to samtools executable. (default: samtools)
  --align_all           When `--bam_path` is not specified, perform alignment for all input reads (rather than just the sampled reads) to the reference genome.
                        Aligning all reads will improve accuracy of sequencing coverage analyses, but can take a considerable amount of time. (default: False)

General arguments:
  General input/output arguments.

  -o OUTPUT_DIR, --output_dir OUTPUT_DIR
                        Directory to save output figures and reports. (default: bamboo_report)
  -t THREADS, --threads THREADS
  --sample_size SAMPLE_SIZE
                        The number of reads to be analyzed. Use --sample_size -1 to disable random sampling and analyze all reads in the input data. (default:
                        100000)
  --seed SEED           Random seed for sampling. (default: 42)
  --keep-intermediates  Do not remove intermediate data files generated in the analyses. (default: False)

参考文献

CycloneSEQ-Bioinformatics/Bamboo/wiki
华大序风官网：https://cycloneseq.genomics.cn

三代测序 - 数据质控 | Bamboo