Introduction

Bambu is a method for transcript discovery and quantification from long read RNA-Seq data. Bambu uses aligned reads and genome reference annotations as input, and will return abundance estimates for all known transcripts and for newly discovered transcripts. Bambu uses the information from the reference annotations to correct misalignment at splice junctions, then reduces the aligned reads to read equivalent classes, and uses this information to identify novel transcripts across all samples of interest. Reads are then assigned to transcripts, and expression estimates are obtained using an expectation maximisation algorithm. Here, we present an example workflow for analysing Nanopore long read RNA-Sequencing data from two human cancer cell lines from the Singapore Nanopore Expression Project (SG-NEx).

Quick start: Transcript discovery and quantification with bambu

Installation

You can install bambu from github:

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("bambu")
BiocManager::install("NanoporeRNASeq")

General Usage

The default mode to run bambu is using a set of aligned reads (bam files), reference genome annotations (gtf file, TxDb object, or bambuAnnotation object), and reference genome sequence (fasta file or BSgenome). bambu will return a summarizedExperiment object with the genomic coordinates for annotated and new transcripts and transcript expression estimates. We highly recommend to use the same annotations that were used for genome alignment. If you have a gtf file and fasta file you can run bambu with the following options:

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)

bambu returns a SummarizedExperiment object which can be accessed as follows:

assays(se) returns the transcript abundance estimates as counts or CPM, as well the full length or unique read counts
rowRanges(se) returns a GRangesList with all annotated and newly discovered transcripts
rowData(se) returns additional information about each transcript such as the gene name and the class of newly discovered transcript

A complete workflow to identify and quantify transcript expression from Nanopore RNA-Seq data

To demonstrate the usage of Bambu, we used long-read RNA-Seq data generated using Oxford Nanopore Sequencing from the NanoporeRNASeq package, which consists of 6 samples from two human cell lines (K562 and MCF7) that were generated by the SG-NEx project. Each of these cell lines has three replicates, with 1 direct RNA sequencing run and 2 cDNA sequencing runs. Reads are aligned to chromosome 22 (Grch38) and stored as bam files. In this workflow, we will demonstrate how to apply bambu to these bam files to identify novel transcripts and estimate transcript expression levels, visualize the results, and identify differentially expressed genes and transcripts.

Input data

Aligned reads (bam files)

bambu takes genomic alignments saved in bam files. Here we use bam-files from the NanoporeRNASeq package, which contains reads aligned to the first half of the human chromosome 22 using minimap2.

library(bambu)
library(NanoporeRNASeq)
data("SGNexSamples")
SGNexSamples
##> DataFrame with 6 rows and 6 columns
##>                sample_id    Platform    cellLine    protocol cancer_type
##>              <character> <character> <character> <character> <character>
##> 1 K562_directcDNA_repl..      MinION        K562  directcDNA   Leukocyte
##> 2 K562_directcDNA_repl..     GridION        K562  directcDNA   Leukocyte
##> 3 K562_directRNA_repli..     GridION        K562   directRNA   Leukocyte
##> 4 MCF7_directcDNA_repl..      MinION        MCF7  directcDNA      Breast
##> 5 MCF7_directcDNA_repl..     GridION        MCF7  directcDNA      Breast
##> 6 MCF7_directRNA_repli..     GridION        MCF7   directRNA      Breast
##>                fileNames
##>              <character>
##> 1 NanoporeRNASeq/versi..
##> 2 NanoporeRNASeq/versi..
##> 3 NanoporeRNASeq/versi..
##> 4 NanoporeRNASeq/versi..
##> 5 NanoporeRNASeq/versi..
##> 6 NanoporeRNASeq/versi..
library(ExperimentHub)
NanoporeData <- query(ExperimentHub(), c("NanoporeRNA", "GRCh38", "BAM"))
bamFile <- Rsamtools::BamFileList(NanoporeData[["EH3808"]])

Genome sequence (fasta file/ BSGenome object)

bambu additionally requires a genome sequence, which is used to correct splicing junctions in read alignments. Ideally, we recommend to use the same genome seqeunce file that was used for alignment to be used for bambu.

# get path to fasta file
genomeSequenceData <- query(ExperimentHub(), c("NanoporeRNA", "GRCh38", "FASTA"))
genomeSequence <- genomeSequenceData[["EH7260"]]

As an option, users can also choose to use a BSgenome object:

library(BSgenome.Hsapiens.NCBI.GRCh38)
genomeSequenceBsgenome <- BSgenome.Hsapiens.NCBI.GRCh38

Genome annotations (bambu annotations object/ gtf file / TxDb object)

bambu also requires a reference transcript annotations object, which is used to correct read alignments, to train a model for transcript discovery, and for quantification. The annotation object can be created from a gtf file:

gtf.file <- system.file("extdata", "Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf", package = "bambu")
annotation <- prepareAnnotations(gtf.file)

The annotation object can also be created from a TxDb object:

txdb <- system.file("extdata", "Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf", package = "bambu")
annotation <- prepareAnnotations(txdb)

The annotation object can be stored and used again for re-running bambu. Here we will used the annotation object from the NanoporeRNASeq package that wasis prepared from a gtf file using the function in by function in bambu.

data("HsChr22BambuAnnotation")
HsChr22BambuAnnotation
##> GRangesList object of length 1500:
##> $ENST00000043402
##> GRanges object with 2 ranges and 2 metadata columns:
##>       seqnames            ranges strand | exon_rank exon_endRank
##>          <Rle>         <IRanges>  <Rle> | <integer>    <integer>
##>   [1]       22 20241415-20243110      - |         2            1
##>   [2]       22 20268071-20268531      - |         1            2
##>   -------
##>   seqinfo: 1 sequence from an unspecified genome; no seqlengths
##> 
##> $ENST00000086933
##> GRanges object with 3 ranges and 2 metadata columns:
##>       seqnames            ranges strand | exon_rank exon_endRank
##>          <Rle>         <IRanges>  <Rle> | <integer>    <integer>
##>   [1]       22 19148576-19149095      - |         3            1
##>   [2]       22 19149663-19149916      - |         2            2
##>   [3]       22 19150025-19150283      - |         1            3
##>   -------
##>   seqinfo: 1 sequence from an unspecified genome; no seqlengths
##> 
##> $ENST00000155674
##> GRanges object with 8 ranges and 2 metadata columns:
##>       seqnames            ranges strand | exon_rank exon_endRank
##>          <Rle>         <IRanges>  <Rle> | <integer>    <integer>
##>   [1]       22 17137511-17138357      - |         8            1
##>   [2]       22 17138550-17138738      - |         7            2
##>   [3]       22 17141059-17141233      - |         6            3
##>   [4]       22 17143098-17143131      - |         5            4
##>   [5]       22 17145024-17145117      - |         4            5
##>   [6]       22 17148448-17148560      - |         3            6
##>   [7]       22 17149542-17149745      - |         2            7
##>   [8]       22 17165209-17165287      - |         1            8
##>   -------
##>   seqinfo: 1 sequence from an unspecified genome; no seqlengths
##> 
##> ...
##> <1497 more elements>

Transcript discovery and quantification

Running bambu

Next we apply bambu on the input data (bam files, annotations, genomeSequence). Bambu will perform isoform discovery to extend the provided annotations, and then quantify the transcript expression from these extended annotations using an Expectation-Maximisation algorithm. Here we will use 1 core, which can be changed to process multiple files in parallel.

se <- bambu(reads = bamFile, annotations = HsChr22BambuAnnotation, genome = genomeSequence,
    ncore = 1)

se
##> class: RangedSummarizedExperiment 
##> dim: 1501 1 
##> metadata(2): incompatibleCounts warnings
##> assays(4): counts CPM fullLengthCounts uniqueCounts
##> rownames(1501): BambuTx1 ENST00000043402 ... ENST00000641933
##>   ENST00000641967
##> rowData names(11): TXNAME GENEID ... txid eqClassById
##> colnames(1): bf442ce62cd_3844
##> colData names(1): name

Optionally, users can choose to apply bambu to do transcript discovery only

seNoquant <- bambu(reads = bamFile, annotations = HsChr22BambuAnnotation, genome = genomeSequence,
    quant = FALSE)

Optionally, users can choose to apply bambu to do quantification only (without isoform discovery)

seUnextended <- bambu(reads = bamFile, annotations = HsChr22BambuAnnotation, genome = genomeSequence,
    discovery = FALSE)

seUnextended
##> class: RangedSummarizedExperiment 
##> dim: 1500 1 
##> metadata(2): incompatibleCounts warnings
##> assays(4): counts CPM fullLengthCounts uniqueCounts
##> rownames(1500): ENST00000043402 ENST00000086933 ... ENST00000641933
##>   ENST00000641967
##> rowData names(4): TXNAME GENEID txid eqClassById
##> colnames(1): bf442ce62cd_3844
##> colData names(1): name

Running multiple samples

If you have multiple replicates for a sample, or plan to do comparative analysis between conditions, it may be beneficial to run all samples together instead of individually. This can be done by providing a vector of paths to all the bam files you want to analyze together.

bamFiles <- Rsamtools::BamFileList(NanoporeData[["EH3808"]], NanoporeData[["EH3809"]],
    NanoporeData[["EH3810"]], NanoporeData[["EH3811"]], NanoporeData[["EH3812"]],
    NanoporeData[["EH3813"]])
se.multiSample <- bambu(reads = bamFiles, annotations = HsChr22BambuAnnotation, genome = genomeSequence)

The advantage of running samples together include: novel transcripts that are identified in multiple samples are assigned unified IDs, enabling comparative analysis between different samples. This is especially important for downstream differential expression analysis when looking at novel transcripts. Running multiple samples can be multithreaded. While running multiple samples, by default, bambu will train a model separately on each sample and score novel transcripts in each sample separately.

If you need to combine samples in multiple configurations (for example different pairwise comparisons) we would recommend using the intermediate rcFiles to save processing time.

Modulating the sensitivity of discovery (pre and post analysis)

When doing transcript discovery there is a balance between sensitivity (the number of real novel transcripts that are detected) and the precision (how many of the novel transcripts are real). To control this balance, bambu uses the novel discovery rate (NDR) as the main parameter. The NDR threshold approximates the proportion of novel candidates output by bambu, relative to the number of known transcripts it found, i.e., an NDR of 0.1 would mean that 10% of all transcripts passing the threshold are classified as novel.

If you are using a genome where you expect a high amount of novel transcripts, a higher NDR is recommended so that these novel transcripts are not missed. Conversely if you are using a well annotated genome, we recommend a lower NDR threshold to reduce the presence of false positives. By default the NDR threshold is automatically chosen for the user based on predicted level of annotation completeness when compared to the default model trained on human reference annotations (Hg38).

To manually select an NDR value, use the NDR argument in bambu:

se.NDR_0.3 <- bambu(reads = bamFiles, annotations = HsChr22BambuAnnotation, genome = genomeSequence,
    NDR = 0.3)

Alternatively transcript discovery can be run without thresholds, producing a GRangesList annotation object with all transcripts scored with its NDR score. Note that this means turning quant = FALSE in running bambu. The annotations can be filtered by their NDR score (see example below), read count and gene read proportion between the discovery and quantification steps or used for other types of analysis.

newAnnotations <- bambu(reads = bamFiles, annotations = HsChr22BambuAnnotation, genome = genomeSequence,
    NDR = 1, quant = FALSE)
annotations.filtered <- newAnnotations[(!is.na(mcols(newAnnotations)$NDR) & mcols(newAnnotations)$NDR <
    0.1) | is.na(mcols(newAnnotations)$NDR)]
se.NDR_1 <- bambu(reads = bamFiles, annotations = annotations.filtered, genome = genomeSequence,
    NDR = 1, discovery = FALSE)

Additionally there are other thresholds that advanced users can access through opt.discovery when running bambu (see arguments).

Visualise results

bambu provides functions to visualise and explore the results. When multiple samples are used, we can visualise the correlation and clustering of all samples with a heatmap:

library(ggplot2)
plotBambu(se.multiSample, type = "heatmap")

Additionally, we can also visualise the correlation with a 2-dimmensional PCA plot.

plotBambu(se.multiSample, type = "pca")

In addition to visualising the correlation between samples, bambu also provide a function to visualise the extended annotation and expression estimation for individual genes. Here we look at gene ENSG00000099968 and visualise the transcript coordinates for annotated and novel isoforms and expression levels for these isoforms across all samples.

plotBambu(se.multiSample, type = "annotation", gene_id = "ENSG00000099968")

##> [[1]]
##> TableGrob (3 x 1) "arrange": 3 grobs
##>   z     cells    name                grob
##> 1 1 (2-2,1-1) arrange      gtable[layout]
##> 2 2 (3-3,1-1) arrange      gtable[layout]
##> 3 3 (1-1,1-1) arrange text[GRID.text.237]

Obtain gene expression estimates from transcript expression

For the downstream analysis, we will add the condition of interest to the object that describes the samples. Here we are interested in a comparison of the 2 cell lines:

colData(se.multiSample)$condition <- as.factor(SGNexSamples$cellLine)

To obtain the accurate gene expression estimates which uses all reads that can be assigned to each gene (including reads that are incompatible with all existing annotations) you can run the following command. Looking at the output, we can see there are novel genes identified as well

seGene.multiSample <- transcriptToGeneExpression(se.multiSample)
seGene.multiSample
##> class: RangedSummarizedExperiment 
##> dim: 579 6 
##> metadata(0):
##> assays(2): counts CPM
##> rownames(579): BambuGene1 BambuGene10 ... ENSG00000284654
##>   ENSG00000284665
##> rowData names(2): GENEID newGeneClass
##> colnames(6): bf442ce62cd_3844 bf44279494ce8_3846 ... bf4423c9a8d16_3852
##>   bf44267751387_3854
##> colData names(2): name condition

We can again use the function to visualise the gene expression data across the 6 samples with a heatmap or PCA plot. As expected, samples from the same cell line showed higher correlation than across the cell lines.

colData(seGene.multiSample)$groupVar <- SGNexSamples$cellLine
plotBambu(seGene.multiSample, type = "heatmap")

Save data (gtf/text)

bambu includes a function to write the extended annotations, the transcript and the gene expression estimates that include any newly discovered genes and transcripts to text files.

save.dir <- tempdir()
writeBambuOutput(se.multiSample, path = save.dir, prefix = "NanoporeRNASeq_")

bambu also includes a function that only exports the extended annotations to gtf file:

save.file <- tempfile(fileext = ".gtf")
writeToGTF(rowRanges(se.multiSample), file = save.file)

Advanced usages for different use cases

Using a pre-trained model

By default, bambu requires at least 1000 transcripts from the annotations to be detected in a sample in order to train a sample specific model. In use cases where this is not possible bambu will instead use a default pre-trained model to calculate the transcript probability score (TPS) for each read class. Users can force this behavior if they believe their annotations are not sufficient for sample specific training (for example if they suspect a high proportion of real novel transcripts are present in their sample). This is advantageous when you want NDR calibration without the impacts of a model trained using low quality annotations.

se <- bambu(reads = bamFiles, annotations = HsChr22BambuAnnotation, genome = genomeSequence,
    opt.discovery = list(fitReadClassModel = FALSE))

The default pre-trained model was trained on SGNex_HepG2_directRNA_replicate5_run1 and has the following characteristics:

Genome: Homo_sapiens.GRCh38.dna_sm.primary_assembly
Annotations: Homo_sapiens.GRCh38.91
Read count: 7,861,846
Technology: Nanopore (ONT)
Library preparation: directRNA
Base Calling Accuracy: 79%
Average Read Length: 1093

We have found the pre-trained model works successfully across species borders (on Arabidopsis thaliana) and on different technologies (PacBio), with only small decreases in performance compared to using a sample specific model. The pre-trained model is not always effective in samples with large differences in sequencing quality or if the library preparation results in biases in the overall structure of the transcriptome. In this case, we would recommend training a new model using similar data from a different sample that has quality reference annotations (See Training a model on another species/dataset and applying it).

For the example in this vignette, the data used is rather small, hence the pre-trained model is used by default. So you won’t see differences in the output, but you may try it out with more data (see SG-NEx data)

De-novo transcript discovery

In cases where the organism does not yet have reference annotations, or unreliable annotations, bambu can be run in de-novo mode. In de-novo mode, bambu does not train a model, and instead uses the pre-trained model to classify novel transcripts (see Using a pre-trained model. To learn how to train a new model for a more closely related organism/sample see Training a model on another species/dataset and applying it. Without annotations bambu is unable to calibrate the NDR output, nor be able to recommend a threshold and will instead use the TPS as the thresholded value. Therefore you should supply a manual NDR threshold (Modulating the sensitivity of discovery (pre and post analysis)) and note that the precision of the output is unlikely to linearly match an applied threshold. The TPS threshold used is (> 1-NDR). If an NDR is not provided, a default NDR threshold of <0.1 is used (an effective TPS threshold of > 0.9). As in Modulating the sensitivity of discovery (pre and post analysis) an NDR of 1 can be provided to output all possible read classes with their TPS scores

novelAnnotations <- bambu(reads = bamFiles, annotations = NULL, genome = genomeSequence,
    NDR = 1, quant = FALSE)

Storing and using preprocessed files (rcFiles)

The first step of bambu involves the construction of read classes which is a large fraction of the running time. This could be time-consuming if we want to perform transcript discovery & quantification multiple times on the same dataset using different configurations (eg. NDR, or combinations of samples), especially when the sample is large. To mitigate this, we can store the read class information as read class files (rcFiles) during a bambu run. Then they can be used as an input argument in the bambu main function for the subsequent bambu run.

# se <- bambu(rcFile = rcFiles, annotations = annotations, genome = fa.file)

rcFiles can be generated in two ways, either as a direct output of the bambu() function when quant and discovery are FALSE, or as written outputs when a path is provided to the rcOutdir argument. When rcFiles are output using rcOutdir this is done using BiocFileCache. For more details on how to access, use, and identify these files see here. A short example is shown below.

Example using rcOutDir to produce preprocessed files

se <- bambu(reads = bamFiles, rcOutDir = save.dir, annotations = HsChr22BambuAnnotation,
    genome = genomeSequence)

This will store a preprocessed rcFile in the provided directory for each sample file provided to reads. To access these files for future use, we recommend using the BioCFileCache package which provides the metadata needed to identify the sample.

library(BiocFileCache)
bfc <- BiocFileCache(save.dir, ask = FALSE)
info <- bfcinfo(bfc)

The info object is a tibble which associates the filename (fpath) with the sample (rname) to help you identify which .rds file you need.

info
##> # A tibble: 6 × 10
##>   rid   rname create_time access_time rpath rtype fpath last_modified_time etag 
##>   <chr> <chr> <chr>       <chr>       <chr> <chr> <chr>              <dbl> <chr>
##> 1 BFC1  bf44… 2024-05-01… 2024-05-01… /tmp… rela… 22bc…                 NA <NA> 
##> 2 BFC2  bf44… 2024-05-01… 2024-05-01… /tmp… rela… 22bc…                 NA <NA> 
##> 3 BFC3  bf44… 2024-05-01… 2024-05-01… /tmp… rela… 22bc…                 NA <NA> 
##> 4 BFC4  bf44… 2024-05-01… 2024-05-01… /tmp… rela… 22bc…                 NA <NA> 
##> 5 BFC5  bf44… 2024-05-01… 2024-05-01… /tmp… rela… 22bc…                 NA <NA> 
##> 6 BFC6  bf44… 2024-05-01… 2024-05-01… /tmp… rela… 22bc…                 NA <NA> 
##> # ℹ 1 more variable: expires <dbl>
# running bambu using the first file
se <- bambu(reads = info$rpath[1:3], annotations = HsChr22BambuAnnotation, genome = genomeSequence)

This output is also generated when both quant and discovery are set to false in a list form indexed by sample.

se <- bambu(reads = bamFiles, annotations = HsChr22BambuAnnotation, genome = genomeSequence,
    discovery = FALSE, quant = FALSE)
##>   |                                                                              |                                                                      |   0%  |                                                                              |============                                                          |  17%  |                                                                              |=======================                                               |  33%  |                                                                              |===================================                                   |  50%  |                                                                              |===============================================                       |  67%  |                                                                              |==========================================================            |  83%  |                                                                              |======================================================================| 100%

As this is an intermediate object it will not be suitable to use for general use cases. We will document the object below for any potential advanced use cases that may arise.

rowData(se[[1]])
##> DataFrame with 5601 rows and 23 columns
##>                     chr.rc strand.rc    startSD     endSD readCount.posStrand
##>                   <factor>  <factor>  <numeric> <numeric>           <integer>
##> rc.1                    22         +         NA        NA                   1
##> rc.2                    22         *         NA        NA                   1
##> rc.3                    22         +         NA        NA                   1
##> rc.4                    22         +         NA        NA                   1
##> rc.5                    22         *         NA        NA                   1
##> ...                    ...       ...        ...       ...                 ...
##> rcunsplicedNew.95       22         *         NA        NA                   0
##> rcunsplicedNew.96       22         *         NA        NA                   1
##> rcunsplicedNew.97       22         *    9.89949      0.00                   0
##> rcunsplicedNew.98       22         * 1275.61568   1474.49                  74
##> rcunsplicedNew.99       22         *         NA        NA                   0
##>                             intronStarts             intronEnds
##>                              <character>            <character>
##> rc.1              10716989,10718464,10.. 10717993,10719774,10..
##> rc.2                            10717385               10721323
##> rc.3                            10717714               10720835
##> rc.4                   10718330,10718568      10718490,10723770
##> rc.5                   10718577,10723906      10723868,10729611
##> ...                                  ...                    ...
##> rcunsplicedNew.95                     NA                     NA
##> rcunsplicedNew.96                     NA                     NA
##> rcunsplicedNew.97                     NA                     NA
##> rcunsplicedNew.98                     NA                     NA
##> rcunsplicedNew.99                     NA                     NA
##>                               confidenceType readCount            readIds
##>                                     <factor> <integer>             <list>
##> rc.1              lowConfidenceJunctionReads         1                  1
##> rc.2              lowConfidenceJunctionReads         1                  2
##> rc.3              lowConfidenceJunctionReads         1                  3
##> rc.4              lowConfidenceJunctionReads         1                  4
##> rc.5              lowConfidenceJunctionReads         1                  5
##> ...                                      ...       ...                ...
##> rcunsplicedNew.95               unsplicedNew         1               1020
##> rcunsplicedNew.96               unsplicedNew         1               1021
##> rcunsplicedNew.97               unsplicedNew         2          1022,1023
##> rcunsplicedNew.98               unsplicedNew        97 1024,1025,1031,...
##> rcunsplicedNew.99               unsplicedNew         1               1176
##>                              GENEID novelGene  numExons geneReadProp
##>                         <character> <logical> <integer>    <numeric>
##> rc.1                     BambuGene1      TRUE         8     0.142857
##> rc.2                   BambuGene328      TRUE         2     0.333333
##> rc.3                     BambuGene2      TRUE         2     1.000000
##> rc.4                     BambuGene1      TRUE         3     0.142857
##> rc.5                     BambuGene1      TRUE         3     0.142857
##> ...                             ...       ...       ...          ...
##> rcunsplicedNew.95 unstranded.Gene69      TRUE         1    1.0000000
##> rcunsplicedNew.96 unstranded.Gene70      TRUE         1    1.0000000
##> rcunsplicedNew.97   ENSG00000172967     FALSE         1    0.0124224
##> rcunsplicedNew.98   ENSG00000172967     FALSE         1    0.6024845
##> rcunsplicedNew.99      BambuGene445      TRUE         1    0.5000000
##>                   geneReadCount     equal compatible numAstart   numAend
##>                       <integer> <logical>  <numeric> <integer> <integer>
##> rc.1                          7        NA         NA        NA        NA
##> rc.2                          3        NA         NA        NA        NA
##> rc.3                          1        NA         NA        NA        NA
##> rc.4                          7        NA         NA        NA        NA
##> rc.5                          7        NA         NA        NA        NA
##> ...                         ...       ...        ...       ...       ...
##> rcunsplicedNew.95             1        NA         NA        NA        NA
##> rcunsplicedNew.96             1        NA         NA        NA        NA
##> rcunsplicedNew.97           161     FALSE          0         3         5
##> rcunsplicedNew.98           161     FALSE          0         1         6
##> rcunsplicedNew.99             2        NA         NA        NA        NA
##>                   numTstart   numTend txScore.noFit     txScore
##>                   <integer> <integer>     <numeric>   <numeric>
##> rc.1                     NA        NA            NA          NA
##> rc.2                     NA        NA            NA          NA
##> rc.3                     NA        NA            NA          NA
##> rc.4                     NA        NA            NA          NA
##> rc.5                     NA        NA            NA          NA
##> ...                     ...       ...           ...         ...
##> rcunsplicedNew.95        NA        NA            NA          NA
##> rcunsplicedNew.96        NA        NA            NA          NA
##> rcunsplicedNew.97         2         4     0.0256271 0.000549702
##> rcunsplicedNew.98         3         1     0.1397222 0.006596871
##> rcunsplicedNew.99        NA        NA            NA          NA

Tracking read-to-transcript assignment

Some use cases require knowing which individual reads support specific transcripts (novel and annotated). By default this feature is off due to the memory overhead it introduces but can be turned on using the trackReads argument. The output has three columns: read_id, a list of indices of equal matches, a list of indices of compatible matches. These indices match the annotations found in rowRanges(se)

se <- bambu(reads = bamFiles, annotations = HsChr22BambuAnnotation, genome = genomeSequence,
    trackReads = TRUE)
##>   |                                                                              |                                                                      |   0%  |                                                                              |============                                                          |  17%  |                                                                              |=======================                                               |  33%  |                                                                              |===================================                                   |  50%  |                                                                              |===============================================                       |  67%  |                                                                              |==========================================================            |  83%  |                                                                              |======================================================================| 100%
metadata(se)$readToTranscriptMaps[[1]]
##> # A tibble: 22,387 × 3
##>    readId                               equalMatches compatibleMatches
##>    <chr>                                <list>       <list>           
##>  1 9db30f53-2e3a-45bc-a01a-586782f00797 <NULL>       <NULL>           
##>  2 dec47cfc-becc-4ff1-9344-0275c2518a41 <NULL>       <NULL>           
##>  3 d544f2e7-a8cc-4a4c-94da-161213ba767b <NULL>       <NULL>           
##>  4 3a49d2c1-92e8-409c-bbec-0901cf6efa6c <NULL>       <NULL>           
##>  5 881bcf9d-9d04-44cb-8d69-a2c03eb9f417 <NULL>       <NULL>           
##>  6 4f41b1d0-69bf-4eeb-8eab-fd552ec1e991 <NULL>       <NULL>           
##>  7 a6d7ed99-683c-4f23-a138-d387059fc65e <NULL>       <NULL>           
##>  8 af2ffb1f-98e3-4e0b-befc-5571fd6af49e <NULL>       <NULL>           
##>  9 9a0ed6a4-6d9b-4264-ba97-858f93eca3a7 <NULL>       <NULL>           
##> 10 1e6750b9-ed24-4626-96f3-c25625ffbfa3 <NULL>       <NULL>           
##> # ℹ 22,377 more rows

Training a model on another species/dataset and applying it

In situations where training is not or cannot be performed, and the default model is also not suitable for the sample (the sample is generated from a different technology, species, condition, etc), bambu provides the option to train a new model, if well annotated similar data is available. For example one might train a model on arabidopsis to apply to an unannotated plant sample.

   # first train the model using a related annotated dataset from .bam

se = bambu(reads = sample1.bam, annotations = annotations, genome = fa.file, discovery = FALSE, quant = FALSE, opt.discovery = list(returnModel = TRUE)) # note that discovery and quant need to be set to FALSE, alternatively you can have them set to TRUE and retrieve the model from the rcFile as long as returnModel = TRUE.
# newDefaultModel = metadata(se[[1]])$model # [[1]] will select the model trained on the first sample

  #alternatively train the model using an rcFile
rcFile <- readRDS(pathToRcFile)
newDefaultModel = trainBambu(rcFile)

# use the trained model on another sample
# sample2.bam and fa.file2 represent the aligned reads and genome for the poorly annotated sample
se <- bambu(reads = sample2.bam, annotations = NULL, genome = fa.file2, opt.discovery = list(defaultModels = newDefaultModel, fitReadClassModel = FALSE))

#trainBambu Arguments

rcFile <- NULL, min.readCount = 2, nrounds = 50, NDR.threshold = 0.1, verbose = TRUE

Including single exons

By default bambu does not report single exon transcripts because they are known to have a high frequency of false positives and do not have splice junctions that are used by bambu to distinguish read classes. Nevertheless bambu trains a separate model on single-exon transcripts, and these predictions can be accessed and included in the annotations.

se <- bambu(reads = bamFiles, annotations = HsChr22BambuAnnotation, genome = genomeSequence,
    opt.discovery = list(min.txScore.singleExon = 0))
##>   |                                                                              |                                                                      |   0%  |                                                                              |============                                                          |  17%  |                                                                              |=======================                                               |  33%  |                                                                              |===================================                                   |  50%  |                                                                              |===============================================                       |  67%  |                                                                              |==========================================================            |  83%  |                                                                              |======================================================================| 100%

Downstream analysis

Identifying differentially expressed genes

One of the most common tasks when analysing RNA-Seq data is the analysis of differential gene expression across a condition of intertest. Here we use DESeq2 to find the differentially expressed genes between MCF7 and K562 cell lines. Similar to using results from Salmon, estimates from bambu will first be rounded.

library(DESeq2)
dds <- DESeqDataSetFromMatrix(round(assays(seGene.multiSample)$counts), colData = colData(se.multiSample),
    design = ~condition)
dds.deseq <- DESeq(dds)
deGeneRes <- DESeq2::results(dds.deseq, independentFiltering = FALSE)
head(deGeneRes[order(deGeneRes$padj), ])
##> log2 fold change (MLE): condition MCF7 vs K562 
##> Wald test p-value: condition MCF7 vs K562 
##> DataFrame with 6 rows and 6 columns
##>                  baseMean log2FoldChange     lfcSE      stat
##>                 <numeric>      <numeric> <numeric> <numeric>
##> ENSG00000185686  509.5330       -7.18914  0.491885 -14.61550
##> ENSG00000283633   97.7268       -9.15686  1.241648  -7.37476
##> ENSG00000240972 2422.0142        2.42269  0.337762   7.17276
##> ENSG00000197077   26.6353        9.13674  1.316819   6.93849
##> ENSG00000099977  233.5728        1.88222  0.295669   6.36596
##> ENSG00000169635   44.2941       -3.41640  0.568819  -6.00613
##>                                                                  pvalue
##>                                                               <numeric>
##> ENSG00000185686 0.00000000000000000000000000000000000000000000000223691
##> ENSG00000283633 0.00000000000016463564442875887204802358248569125144718
##> ENSG00000240972 0.00000000000073501963738434064723469394279351604856113
##> ENSG00000197077 0.00000000000396301028105484197572231012093810231539798
##> ENSG00000099977 0.00000000019407034315011366727597473872107131098729838
##> ENSG00000169635 0.00000000190005590148128080419213970528400509119038020
##>                                                                  padj
##>                                                             <numeric>
##> ENSG00000185686 0.000000000000000000000000000000000000000000000639756
##> ENSG00000283633 0.000000000023542897153312519636878684077131826564433
##> ENSG00000240972 0.000000000070071872097307140962563386178858781949685
##> ENSG00000197077 0.000000000283355235095421190358932019871807117539220
##> ENSG00000099977 0.000000011100823628186501602749632544239744014191729
##> ENSG00000169635 0.000000090569331303941049047086719533949938920613931

A quick summary of differentially expressed genes

summary(deGeneRes)
##> 
##> out of 286 with nonzero total read count
##> adjusted p-value < 0.1
##> LFC > 0 (up)       : 20, 7%
##> LFC < 0 (down)     : 23, 8%
##> outliers [1]       : 0, 0%
##> low counts [2]     : 0, 0%
##> (mean count < 0)
##> [1] see 'cooksCutoff' argument of ?results
##> [2] see 'independentFiltering' argument of ?results

We can also visualise the MA-plot for differentially used isoforms using . However, visualizing the MA-plots using the original log-fold change results will be affected by the noise associated with log2 fold changes from low count genes without requiring arbitrary filtering thresholds. As recommended in the DESeq2 tutorial. we applied the same shrinkage to effect sizes to improve the visualization.

library(apeglm)
resLFC <- lfcShrink(dds.deseq, coef = "condition_MCF7_vs_K562", type = "apeglm")
plotMA(resLFC, ylim = c(-3, 3))

Identifying differential transcript usage

We used DEXSeq to detect alternative used isoforms.

library(DEXSeq)
dxd <- DEXSeqDataSet(countData = round(assays(se.multiSample)$counts), sampleData = as.data.frame(colData(se.multiSample)),
    design = ~sample + exon + condition:exon, featureID = rowData(se.multiSample)$TXNAME,
    groupID = rowData(se.multiSample)$GENEID)
dxr <- DEXSeq(dxd)
head(dxr)
##> 
##> LRT p-value: full vs reduced
##> 
##> DataFrame with 6 rows and 12 columns
##>                                  groupID   featureID exonBaseMean dispersion
##>                              <character> <character>    <numeric>  <numeric>
##> BambuGene1:BambuTx1           BambuGene1    BambuTx1     0.547211         NA
##> BambuGene28:BambuTx2         BambuGene28    BambuTx2     1.410226         NA
##> BambuGene33:BambuTx3         BambuGene33    BambuTx3     1.820219         NA
##> BambuGene10:BambuTx4         BambuGene10    BambuTx4     0.387598         NA
##> BambuGene12:BambuTx5         BambuGene12    BambuTx5     6.931283         NA
##> ENSG00000283023:BambuTx6 ENSG00000283023    BambuTx6     0.812250    2.31969
##>                               stat    pvalue      padj      K562      MCF7
##>                          <numeric> <numeric> <numeric> <numeric> <numeric>
##> BambuGene1:BambuTx1             NA        NA        NA        NA        NA
##> BambuGene28:BambuTx2            NA        NA        NA        NA        NA
##> BambuGene33:BambuTx3            NA        NA        NA        NA        NA
##> BambuGene10:BambuTx4            NA        NA        NA        NA        NA
##> BambuGene12:BambuTx5            NA        NA        NA        NA        NA
##> ENSG00000283023:BambuTx6  0.611639  0.434171  0.942744  0.260969  0.431158
##>                          log2fold_MCF7_K562   genomicData    countData
##>                                   <numeric> <GRangesList>     <matrix>
##> BambuGene1:BambuTx1                      NA                  4:3:0:...
##> BambuGene28:BambuTx2                     NA                  4:9:3:...
##> BambuGene33:BambuTx3                     NA                 14:6:3:...
##> BambuGene10:BambuTx4                     NA                  3:2:0:...
##> BambuGene12:BambuTx5                     NA               32:37:13:...
##> ENSG00000283023:BambuTx6            1.45859                  1:7:1:...

We can visualize the MA-plot

plotMA(dxr, cex = 0.8)

Getting help

Questions and issues can be raised at the Bioconductor support site (once bambu is available through bioconductor): https://support.bioconductor.org. Please tag your your posts with bambu.

Alternatively, issues can be raised at the bambu Github repository:https://github.com/GoekeLab/bambu.

Citing bambu

A manuscript describing bambu is currently in preparation. If you use bambu for your research, please cite using the following doi: 10.5281/zenodo.3900025.

Session Information

sessionInfo()
##> R version 4.4.0 beta (2024-04-15 r86425)
##> Platform: x86_64-pc-linux-gnu
##> Running under: Ubuntu 22.04.4 LTS
##> 
##> Matrix products: default
##> BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
##> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.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: America/New_York
##> tzcode source: system (glibc)
##> 
##> attached base packages:
##> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
##> [8] base     
##> 
##> other attached packages:
##>  [1] DEXSeq_1.50.0                         
##>  [2] RColorBrewer_1.1-3                    
##>  [3] AnnotationDbi_1.66.0                  
##>  [4] BiocParallel_1.38.0                   
##>  [5] apeglm_1.26.0                         
##>  [6] DESeq2_1.44.0                         
##>  [7] ggplot2_3.5.1                         
##>  [8] BSgenome.Hsapiens.NCBI.GRCh38_1.3.1000
##>  [9] Rsamtools_2.20.0                      
##> [10] NanoporeRNASeq_1.13.0                 
##> [11] ExperimentHub_2.12.0                  
##> [12] AnnotationHub_3.12.0                  
##> [13] BiocFileCache_2.12.0                  
##> [14] dbplyr_2.5.0                          
##> [15] bambu_3.6.0                           
##> [16] BSgenome_1.72.0                       
##> [17] rtracklayer_1.64.0                    
##> [18] BiocIO_1.14.0                         
##> [19] Biostrings_2.72.0                     
##> [20] XVector_0.44.0                        
##> [21] SummarizedExperiment_1.34.0           
##> [22] Biobase_2.64.0                        
##> [23] GenomicRanges_1.56.0                  
##> [24] GenomeInfoDb_1.40.0                   
##> [25] IRanges_2.38.0                        
##> [26] S4Vectors_0.42.0                      
##> [27] BiocGenerics_0.50.0                   
##> [28] MatrixGenerics_1.16.0                 
##> [29] matrixStats_1.3.0                     
##> 
##> loaded via a namespace (and not attached):
##>   [1] splines_4.4.0            bitops_1.0-7             filelock_1.0.3          
##>   [4] tibble_3.2.1             graph_1.82.0             XML_3.99-0.16.1         
##>   [7] rpart_4.1.23             lifecycle_1.0.4          httr2_1.0.1             
##>  [10] doParallel_1.0.17        lattice_0.22-6           ensembldb_2.28.0        
##>  [13] MASS_7.3-60.2            OrganismDbi_1.46.0       backports_1.4.1         
##>  [16] magrittr_2.0.3           Hmisc_5.1-2              sass_0.4.9              
##>  [19] rmarkdown_2.26           jquerylib_0.1.4          yaml_2.3.8              
##>  [22] ggbio_1.52.0             DBI_1.2.2                abind_1.4-5             
##>  [25] zlibbioc_1.50.0          purrr_1.0.2              AnnotationFilter_1.28.0 
##>  [28] biovizBase_1.52.0        RCurl_1.98-1.14          nnet_7.3-19             
##>  [31] VariantAnnotation_1.50.0 rappdirs_0.3.3           circlize_0.4.16         
##>  [34] GenomeInfoDbData_1.2.12  genefilter_1.86.0        annotate_1.82.0         
##>  [37] codetools_0.2-20         DelayedArray_0.30.0      xml2_1.3.6              
##>  [40] tidyselect_1.2.1         shape_1.4.6.1            UCSC.utils_1.0.0        
##>  [43] farver_2.1.1             base64enc_0.1-3          GenomicAlignments_1.40.0
##>  [46] jsonlite_1.8.8           GetoptLong_1.0.5         Formula_1.2-5           
##>  [49] survival_3.6-4           iterators_1.0.14         bbmle_1.0.25.1          
##>  [52] foreach_1.5.2            tools_4.4.0              progress_1.2.3          
##>  [55] Rcpp_1.0.12              glue_1.7.0               gridExtra_2.3           
##>  [58] SparseArray_1.4.0        xfun_0.43                dplyr_1.1.4             
##>  [61] numDeriv_2016.8-1.1      withr_3.0.0              formatR_1.14            
##>  [64] BiocManager_1.30.22      fastmap_1.1.1            GGally_2.2.1            
##>  [67] fansi_1.0.6              digest_0.6.35            R6_2.5.1                
##>  [70] mime_0.12                colorspace_2.1-0         Cairo_1.6-2             
##>  [73] dichromat_2.0-0.1        biomaRt_2.60.0           RSQLite_2.3.6           
##>  [76] utf8_1.2.4               tidyr_1.3.1              generics_0.1.3          
##>  [79] data.table_1.15.4        prettyunits_1.2.0        httr_1.4.7              
##>  [82] htmlwidgets_1.6.4        S4Arrays_1.4.0           ggstats_0.6.0           
##>  [85] pkgconfig_2.0.3          gtable_0.3.5             blob_1.2.4              
##>  [88] hwriter_1.3.2.1          ComplexHeatmap_2.20.0    htmltools_0.5.8.1       
##>  [91] geneplotter_1.82.0       RBGL_1.80.0              ProtGenerics_1.36.0     
##>  [94] clue_0.3-65              scales_1.3.0             png_0.1-8               
##>  [97] knitr_1.46               rstudioapi_0.16.0        reshape2_1.4.4          
##> [100] rjson_0.2.21             coda_0.19-4.1            checkmate_2.3.1         
##> [103] curl_5.2.1               bdsmatrix_1.3-7          cachem_1.0.8            
##> [106] GlobalOptions_0.1.2      stringr_1.5.1            BiocVersion_3.19.1      
##> [109] parallel_4.4.0           foreign_0.8-86           restfulr_0.0.15         
##> [112] pillar_1.9.0             grid_4.4.0               vctrs_0.6.5             
##> [115] xtable_1.8-4             cluster_2.1.6            htmlTable_2.4.2         
##> [118] evaluate_0.23            GenomicFeatures_1.56.0   magick_2.8.3            
##> [121] mvtnorm_1.2-4            cli_3.6.2                locfit_1.5-9.9          
##> [124] compiler_4.4.0           rlang_1.1.3              crayon_1.5.2            
##> [127] labeling_0.4.3           emdbook_1.3.13           plyr_1.8.9              
##> [130] stringi_1.8.3            txdbmaker_1.0.0          munsell_0.5.1           
##> [133] lazyeval_0.2.2           Matrix_1.7-0             hms_1.1.3               
##> [136] bit64_4.0.5              statmod_1.5.0            KEGGREST_1.44.0         
##> [139] highr_0.10               memoise_2.0.1            bslib_0.7.0             
##> [142] bit_4.0.5                xgboost_1.7.7.1

Analysing Long Read RNA-Seq data with bambu

Ying Chen, Andre Sim, Yuk Kei Wan, Jonathan Göke