Bulk RNA-seq Workflow

Source: vignettes/bulk-workflow.Rmd

This tutorial shows the two bulk RNA-seq configurations covered here in one continuous R workflow:

  1. quantify a sample from short reads;
  2. add transcript-level long-read counts and quantify the same sample again.

The second configuration reuses the index and short-read counts created for the first, making it easy to compare the two results.

What This Tutorial Produces

By the end, you will have:

  • an mpaqt_index object for the reference transcriptome;
  • an mpaqt_counts_sr object prepared from paired short-read FASTQ files;
  • an mpaqt_counts_lr object prepared from a Bambu transcript-count table;
  • sr_result, based on short reads only; and
  • sr_lr_result, based on short and long reads together.

The workflow is:

reference + short reads  ->  sr_result
                     + long-read counts  ->  sr_lr_result

Required Inputs

For the short-read configuration, you need:

  • a GTF annotation and matching transcriptome FASTA;
  • paired short-read FASTQ files; and
  • kallisto and bustools available on your system.

For the short-plus-long-read configuration, you also need a transcript-count table produced by Bambu. Its first column contains transcript IDs, and its second column contains counts for the sample used here.

Transcript IDs must use the same convention in the GTF, transcriptome FASTA, and Bambu table, including version suffixes when present.

Create or Load the Index

Load MPAQT and create the index once for a reference:

library(mpaqt)

index <- mpaqt_index(
    annotation = "reference/gencode.annotation.gtf",
    transcriptome = "reference/gencode.transcripts.fa",
    output_file = "reference/mpaqt.index.rds",
    read_length = 91L,
    store_gtf = TRUE,
    threads = 8L
)

Set read_length to the length of the transcript sequence represented by one short-read mate. store_gtf = TRUE lets the same index support the optional raw-FLNC route shown later. Reuse the saved index for later samples:

index <- mpaqt_read_index("reference/mpaqt.index.rds")

The index reference describes additional ways to create an index, including starting from a genome or an existing kallisto index.

Prepare Short Reads

Prepare equivalence-class counts from paired FASTQ files:

sr_counts <- mpaqt_prepare_short_reads(
    index = index,
    fastq_1 = "data/sample_R1.fastq.gz",
    fastq_2 = "data/sample_R2.fastq.gz",
    output_dir = "results/short-read",
    threads = 8L
)

This step runs kallisto and bustools, returns an mpaqt_counts_sr object, and saves results/short-read/mpaqt.short_read.rds for reuse.

Quantify with Short Reads

Run the short-read-only configuration:

sr_result <- mpaqt_quant(
    index = index,
    sr_counts = sr_counts
)

sr_result is an mpaqt_quant_result object containing transcript-level TPM values, abundances, convergence information, and fitted values.

head(tpm(sr_result))

Prepare Long Reads

For the second configuration, prepare long-read counts from a Bambu table such as:

transcript_id,sample1
ENST00000456328.2,18
ENST00000450305.2,6
ENST00000488147.1,0

Pass that table to mpaqt_prepare_long_reads():

lr_counts <- mpaqt_prepare_long_reads(
    index = index,
    bambu_counts = "data/bambu_counts.csv",
    output_dir = "results/long-read"
)

MPAQT matches the table’s transcript IDs to the index order and saves results/long-read/mpaqt.long_read.rds. Transcripts absent from the table receive zero counts; transcript IDs that are not in the index are ignored, and substantial mismatches trigger a warning.

Quantify with Short and Long Reads

Reuse index and sr_counts, then add lr_counts:

sr_lr_result <- mpaqt_quant(
    index = index,
    sr_counts = sr_counts,
    lr_counts = lr_counts
)

You now have both configurations in separate objects: sr_result for short reads alone and sr_lr_result for short and long reads together.

Inspect and Save Results

Use the result accessors to inspect the second configuration:

head(tpm(sr_lr_result))
head(abundances(sr_lr_result))
head(gene_abundances(sr_lr_result, index))

gene_abundances() needs the same index used for quantification so it can map transcripts to genes.

Save both results with distinct prefixes:

mpaqt_save_result(
    sr_result,
    output_dir = "results/quantification",
    prefix = "sample.short"
)

mpaqt_save_result(
    sr_lr_result,
    output_dir = "results/quantification",
    prefix = "sample.short_long"
)

Each call writes a full .quant.rds result and a .quant.csv table of TPM values. The files can be loaded later with readRDS() and data.table::fread(), respectively.

Alternative Input Formats

The examples above use the most direct path from paired FASTQ files and a Bambu table. If those inputs have already been processed, reuse the existing files instead.

Load prepared short-read counts from RDS:

sr_counts <- mpaqt_prepare_short_reads(
    index = index,
    rds_file = "results/short-read/mpaqt.short_read.rds"
)

Or prepare them from matching BUS and equivalence-class files:

sr_counts <- mpaqt_prepare_short_reads(
    index = index,
    bus_file = "data/output.bus",
    ec_file = "data/matrix.ec",
    output_dir = "results/short-read"
)

Load prepared long-read counts from RDS:

lr_counts <- mpaqt_prepare_long_reads(
    index = index,
    rds_file = "results/long-read/mpaqt.long_read.rds"
)

Raw FLNC reads can also be aligned and quantified before MPAQT prepares the long-read count object:

lr_counts <- mpaqt_prepare_long_reads(
    index = index,
    flnc_fastq = "data/sample.flnc.fastq.gz",
    genome = "reference/genome.fa",
    output_dir = "results/long-read",
    threads = 8L
)

This route additionally requires minimap2, samtools, and the Bioconductor package Bambu. See the short-read and long-read references for input details.

Troubleshooting

Transcript IDs do not match. Confirm that the GTF, transcriptome FASTA, and count table come from the same reference release and either all retain or all remove transcript version suffixes.

kallisto or bustools is not found. Check the external tools described in the installation article, then restart R if your PATH changed.

The Bambu table has multiple samples. The preparation function uses the second column. Supply a table containing the transcript-ID column and the one sample column you want to quantify.

A saved object cannot be loaded. Make sure the RDS contains the matching MPAQT object type and was prepared against the same index.

Next Steps