全长转录组 | 三代全长转录组分析流程（PacBio & ONT ）-- Bambu

今天我们继续介绍一款使用三代全长转录本数据进行转录本注释和定量的工具 - Bambu。来自新加坡科技研究局 (A-STAR) 的Jonathan Göke（图1）开发的长度长RNA-seq转录组分析工具Bambu，于2023年6月12日发表在《Nature Methods》杂志上，题目为Context-aware transcript quantification from long-read RNA-seq data with Bambu。该工具基于机器学习来识别和表征新转录本，从而能够对不同物种和样本进行适应性分析。

图1. Jonathan Göke团队

大多数转录本定量方法依赖于固定的基因参考注释文件；然而，实际情况下的转录组是动态的，需要根据当时的情境情况来判断，静态转录组注释文件包含了一些基因的非活跃转录本（isoforms)，而对于另一些基因存在注释不完整的情况。Bambu，一种利用长读长RNA-seq测序数据，通过基于机器学习的转录本鉴定方法，实现对实际测序数据里转录本的定量。为了识别新的转录本，Bambu 估计了新发现率（novel discovery rate），用可解释的、精度校准的单一参数替代了任意的单样本阈值。Bambu 保留了全长和独特的转录本序列，使其在存在非活跃转录本（isoforms) 的情况下能够进行准确的定量。与现有的转录本鉴定方法相比，Bambu 在不损失灵敏性的情况下实现了更高的精度。依据实际数据情况的动态注释改善了新的和已知的转录本的定量。利用长读长RNA测序数据和机器学习，Bambu 促进了准确的转录本鉴定和定量 （图2）。

图2. Bambu文章摘要

一、软件介绍

Bambu 是一个利用长读长RNA-Seq数据进行多样本转录本鉴定和定量的R包。在进行序列比对后，可以使用Bambu 获取已知和新转录本以及基因的表达量。Bambu 的输出可以直接用于可视化和下游分析，例如差异基因表达或转录本使用情况等。

Bambu 分析流程可以分成四步（图3. a-d）：首先，采用概率模型，利用参考注释、基因组序列以及从数据中提取的特征（详细特征可见方法部分）来校正RNA连接（junction）区域或位点的比对。使用相同剪切位点的校正序列被合并成一类序列 (Read class, RCs）。第二步，整合来自所有样本的合并序列（Read classes），并对跨样本间的新发现率（novel discovery rate，NDR）进行计算。小于特定NDR阈值的Read class序列被归为新的转录本，并将新注释的转录本加入参考转录本注释中，形成基于实际样本的注释库（扩展了原本的参考注释）。第三步，利用扩展参考基因注释，将每一个Read class序列重新归类标注。在这个过程中，由于比对误差造成的不精确的序列匹配可以被修正。第四步，使用扩展参考基因注释，对每一样本中的转录本进行最终的定量，获得基因/转录本的表达矩阵。

图3. Bambu分析流程

二、软件安装

Bambu: https://github.com/GoekeLab/bambu

版本：v3.2.4

1. 使用Bioconductor安装Bambu

在R中，或Rstudio，或Rstudio-server环境中运行以下命令进行安装，安装完成如图4所示。

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")  #如果已经安装BiocManager可忽略此行命令
BiocManager::install("bambu")

图4. Bambu的R包安装

建议：因为全长转录组数据量一般较大以及计算量大，一般台式机或者自己的笔记本电脑部署的本地R软件，可能导致占用大量内存和计算资源。建议使用服务器，部署 Rstudio-server，这样就可以在服务器上进行数据存储和计算了。具体安装部署步骤可参考官网：https://posit.co/download/rstudio-server/。

2. 测试Bambu的安装

用自带的测试数据确认Bambu安装成功。

#加载bambu软件包
library(bambu)

#运行以下命令进行测序数据，参考基因组，参考基因组注释文件的导入，并进行全长转录组分析
test.bam <- system.file("extdata", "SGNex_A549_directRNA_replicate5_run1_chr9_1_1000000.bam", package = "bambu") #读入测序数据

fa.file <- system.file("extdata", "Homo_sapiens.GRCh38.dna_sm.primary_assembly_chr9_1_1000000.fa", package = "bambu") #读入参考基因组

gtf.file <- system.file("extdata", "Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf", package = "bambu") #读入参考基因组注释文件

bambuAnnotations <- prepareAnnotations(gtf.file) #参考基因组注释预处理

se <- bambu(reads = test.bam, annotations = bambuAnnotations, genome = fa.file) #全长转录组分析

运行完最后一行，得到 --- Start isoform quantification --- 和--- Finished running Bambu ---，则说明Bambu运行和安装成功。

三、软件使用

默认模式下输入文件:

• 经过和参考基因组比对的序列文件，.bam文件格式。
• 参考基因组注释文件，.gtf文件。
• 参考基因组文件，.fa文件。

1. 准备参考基因组比对序列文件

因为Bambu说明文档没有提供比对方法，所以这里采用常用的长度长（long-read）比对软件minimap2。

#根据三代测序平台和建库方法选择合适的运行命令，一步法
$ minimap2 -ax splice:hq -uf ref.fa iso-seq.fq | samtools sort -@ 12 -o align.bam --write-index -        # PacBio Iso-seq/traditional cDNA
$ minimap2 -ax splice ref.fa nanopore-cdna.fa | samtools sort -@ 12 -o align.bam --write-index -      # Nanopore 2D cDNA-seq
$ minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq | samtools sort -@ 12 -o align.bam --write-index -  # Nanopore Direct RNA-seq

#分步进行示例
$ minimap2 -ax splice:hq -uf ref.fa iso-seq.fq | samtools view -@ 12 -bS | samtools sort -@ 12 -o align.bam        # PacBio Iso-seq/traditional cDNA
$ minimap2 -ax splice ref.fa nanopore-cdna.fa | samtools view -@ 12 -bS | samtools sort -@ 12 -o align.bam      # Nanopore 2D cDNA-seq
$ minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq | samtools view -@ 12 -bS | samtools sort -@ 12 -o align.bam  # Nanopore Direct RNA-seq

#对bam文件进行索引
$ samtools index align.bam

注意:

• PacBio平台和ONT平台产生的数据，具体参数有所不同，详细请参考minimap2使用文档。
• 这里使用的参考基因组和Bambu使用的要统一起来。
• 需要提前安装samtools，留意安装版本。
• 因为minimap2输出的是.sam文件格式，所以使用samtools将.sam转换成.bam，并且使用samtools sort 对 .bam 进行排序，然后用samtools index进行索引。上面展示了一步转.bam和分步转.bam的示例。
• samtools v1.18版本有--write-index - 选项; samtools v1.9版本则没有。

2. 读取样本序列.bam文件

#可同时读入多个样本文件
samples <- c(S1.bam, S1.bam, S1.bam)

3. 准备参考基因组文件

如果人和小鼠可以在Gencode上下载，如遇到其它物种可以从Ensembl上下载。

4. 准备参考基因组注释

annotations <- prepareAnnotations( *.gft ）

#保存注释文件为rds文件
saveRDS(annotations, ”/path/to/annotations.rds” )

#下次读取注释文件
annotations <- readRDS("/path/to/annotations.rds")

5. 运行Bambu

se <- bambu(reads = samples, annotations = annotations, genome = genome.fa)

6. 只鉴定新转录本，不定量

se.discoveryOnly <- bambu(reads = samples, annotations = annotations, genome = genome.fa, quant = FALSE)

7. 只对已知的转绿本和基因进行定量，不做新的转录本鉴定

se.quantOnly <- bambu(reads = samples, annotations = annotations, genome = genome.fa, discovery = FALSE)

8. 不依赖于参考注释的转录本组装

novelAnnotations <- bambu(reads = test.bam, annotations = NULL, genome = fa.file, NDR = 0.5, quant = FALSE)

具体参数设置和调试参考官方github。

四、输出文件

Bambu 运行完后会返回一个SummarizedExperiment对象，可以通过以下方式访问：

• assays(se) 返回转录本表达丰度矩阵，以count或CPM的形式。
• rowRanges(se) 返回GRangesList，包含所有已注释和新发现的转录本。
• rowData(se) 返回每个转录本的额外信息。

通过使用assays()中的变量（如counts或CPM）提取转录本表达量：

• assays(se)$counts - count表达量。
• assays(se)$CPM - 测序深度归一化后的表达量。
• assays(se)$fullLengthCounts - 每个转录本的全长序列count表达量。
• assays(se)$uniqueCounts - 每个转录本唯一回贴序列的count表达量。

可以使用writeBambuOutput()输出完整的结果文件。这个命令可以产生三个文件：1. 扩展的.gtf文件。2.转录本的count表达矩阵。3.基因的count表达矩阵。

writeBambuOutput(se, path = "./bambu/")

如果只对新的转录本感兴趣，可以过滤掉参考注释的转录本。

se.novel = se[mcols(se)$novelTranscript,]
writeBambuOutput(se.novel, path = "./bambu/")

五、可视化

通过plotBambu可以对基因/转录本进行可视化（图5）。

plotBambu(se, type = "annotation", gene_id = "ENSG00000107104")

图5. 转录本tx.9的外显子组成和在样本中的表达量

转录本tx.9以及其对应基因的其它转录本的结构和表达量（图6）。

plotBambu(se, type = "annotation", transcript_id = "tx.9")

图6. 转录本tx.9以及其对应基因的其它转录本的结构和表达量

通过plotBambu展示样本分组PCA聚类（图7）。

图7. 基于基因/转录本表达量的PCA聚类

参考文献

1. Chen, Ying, et al. "Context-aware transcript quantification from long-read RNA-seq data with Bambu." Nature methods. (2023)
2. Bambu github： https://github.com/GoekeLab/bambu