Skip to contents

Processing bulk RNA-seq inputs from FASTQ or BAM files

Source: vignettes/preparing-user-expression-data.Rmd
preparing-user-expression-data.Rmd

Bulk RNA-seq can be used with POWERUP after reads are converted into a gene-level expression matrix on a scale comparable to the reference expression data used for model training.

This article describes an optional workflow for paired-end bulk RNA-seq:

paired FASTQ -> STAR -> Salmon -> gene TPM -> log2(TPM + 1) -> POWERUP user matrix

If the starting point is a BAM file, paired reads are first recovered and then passed through the same workflow:

BAM -> paired FASTQ -> STAR -> Salmon -> gene TPM -> log2(TPM + 1) -> POWERUP user matrix

The resulting matrix has samples as rows and expression features as columns and can be supplied to prepare_powerup_data() as user_matrix. If you already have a compatible expression matrix, none of the upstream processing described here is required.

The workflow below uses STAR 2.7.11b, Salmon 1.10.0, human GRCh38, and GENCODE v38. The reference is human-only and does not include ERCC or viral sequences.

What you need

Input or software Requirement
FASTQ input Paired read 1 and read 2 FASTQ files, uncompressed or gzip-compressed
BAM input Paired-end RNA-seq BAM from which read pairs can be recovered
Sample identifiers One unique identifier for each sample
Reference expression table The expression table used for POWERUP model training, used here to define expected feature names
RNA reference Human GRCh38 genome, GENCODE v38 annotation, transcript FASTA, and STAR index
R package data.table
Command-line tools STAR 2.7.11b, Salmon 1.10.0, samtools, gffread, and standard Unix utilities (bash, awk, gzip/zcat)

Install the command-line tools using the package or environment manager appropriate for your system and make sure they are available on PATH. The examples assume a Unix-like environment, as is common for RNA-seq analysis, but do not assume cloud storage or a particular directory layout. The main bioinformatics tools are discovered with Sys.which().

The examples use eight threads. Increase or decrease this value to match the available compute resources.

Set up the working directories and software

Choose local directories for the reference and sample outputs. reference_expression_file should point to the same reference expression table that will be supplied to POWERUP for model training.

library(data.table)

threads <- 8
work_dir <- "powerup_rna"
reference_dir <- file.path(work_dir, "grch38_gencode38")
reference_expression_file <- "/path/to/reference_expression.csv"

dir.create(work_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(reference_dir, recursive = TRUE, showWarnings = FALSE)

star <- Sys.which("STAR")
salmon <- Sys.which("salmon")
samtools <- Sys.which("samtools")
gffread <- Sys.which("gffread")

required_tools <- c(STAR = star, salmon = salmon, samtools = samtools, gffread = gffread)
missing_tools <- names(required_tools)[!nzchar(required_tools)]

if (length(missing_tools) > 0) {
  message <- paste(
    "Install these command-line tools and add them to PATH:",
    paste(missing_tools, collapse = ", ")
  )
  stop(message)
}

cat("STAR:", system2(star, "--version", stdout = TRUE), "\n")
cat("Salmon:", system2(salmon, "--version", stdout = TRUE), "\n")

For the closest match to this workflow, use STAR 2.7.11b and Salmon 1.10.0.

Build the RNA reference

The RNA reference only needs to be built once and can then be reused for all samples. The example below downloads the same GRCh38 genome FASTA and GENCODE v38 annotation used for this workflow, generates a transcript FASTA with gffread, and builds the STAR index.

genome_url <- paste0(
  "https://storage.googleapis.com/gcp-public-data--broad-references/",
  "hg38/v0/Homo_sapiens_assembly38.fasta"
)

gtf_url <- paste0(
  "https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/",
  "release_38/gencode.v38.annotation.gtf.gz"
)

genome <- file.path(reference_dir, "Homo_sapiens_assembly38.fasta")
gtf_gz <- file.path(reference_dir, "gencode.v38.annotation.gtf.gz")
gtf <- file.path(reference_dir, "gencode.v38.annotation.gtf")
transcript_fasta <- file.path(reference_dir, "gencode.v38.transcripts.fa")
star_index <- file.path(reference_dir, "star_index")

if (!file.exists(genome)) download.file(genome_url, genome, mode = "wb")
if (!file.exists(gtf_gz)) download.file(gtf_url, gtf_gz, mode = "wb")

if (!file.exists(gtf)) {
  status <- system2("gzip", c("-dc", gtf_gz), stdout = gtf)
  if (status != 0) stop("Failed to decompress the GENCODE annotation.")
}

if (!file.exists(transcript_fasta)) {
  status <- system2(gffread, c(gtf, "-g", genome, "-w", transcript_fasta))
  if (status != 0) stop("Failed to generate the transcript FASTA.")
}

if (!file.exists(file.path(star_index, "Genome"))) {
  dir.create(star_index, recursive = TRUE, showWarnings = FALSE)

  status <- system2(
    star,
    c(
      "--runMode", "genomeGenerate",
      "--runThreadN", threads,
      "--genomeDir", star_index,
      "--genomeFastaFiles", genome,
      "--sjdbGTFfile", gtf,
      "--sjdbOverhang", 100
    )
  )

  if (status != 0) stop("STAR index generation failed.")
}

cat("Reference build complete\n")
cat("Genome:", genome, "\n")
cat("GTF:", gtf, "\n")
cat("Transcript FASTA:", transcript_fasta, "\n")
cat("STAR index:", star_index, "\n")

Building the STAR index is the most resource-intensive setup step. Keep the complete reference directory after it has been built so the same files can be reused for future samples.

Define the reference feature map

POWERUP ultimately needs the new expression matrix to use the same feature names as the reference expression data. Read the reference header and map each feature to its gene symbol.

The code below works with plain gene-symbol columns and with columns formatted as GENE (EntrezID). If an Entrez ID suffix is present, it is removed only for gene matching; the original reference feature name is preserved in the final POWERUP matrix.

reference_names <- names(fread(reference_expression_file, nrows = 0, check.names = FALSE))
reference_features <- reference_names[-1]

feature_map <- data.table(
  feature = reference_features,
  gene_name = sub(" \\([0-9]+\\)$", "", reference_features)
)

cat("Reference expression features:", nrow(feature_map), "\n")

The expression values in the reference file are not used during this preprocessing step. The file is used here only to define the expected gene features and their exact column names.

Build the transcript-to-gene map

Salmon reports transcript-level TPM values. A transcript-to-gene map from the same GENCODE v38 annotation is used to aggregate those transcript TPM values to genes.

tx2gene_file <- file.path(reference_dir, "gencode.v38.tx2gene.tsv")

if (!file.exists(tx2gene_file)) {
  cmd <- paste0("awk -F '\\t' '$3==\"transcript\" {print $9}' ", shQuote(gtf))
  attrs <- system2("bash", c("-c", shQuote(cmd)), stdout = TRUE)

  has_tx <- grepl('transcript_id "', attrs)
  has_gene <- grepl('gene_name "', attrs)
  attrs <- attrs[has_tx & has_gene]

  tx2gene <- data.table(
    transcript_id = sub('.*transcript_id "([^"]+)".*', '\\1', attrs),
    gene_name = sub('.*gene_name "([^"]+)".*', '\\1', attrs)
  )

  tx2gene <- unique(tx2gene)
  fwrite(tx2gene, tx2gene_file, sep = "\t")
} else {
  tx2gene <- fread(tx2gene_file)
}

tx2gene[, transcript_key := sub("\\.[0-9]+$", "", transcript_id)]

cat("Transcript-to-gene mappings:", nrow(tx2gene), "\n")
cat("Unique genes:", uniqueN(tx2gene$gene_name), "\n")

Transcript version suffixes are removed before matching because transcript identifiers in Salmon output and the annotation may differ only by their version suffix.

Define the RNA-seq processing function

The following function accepts either a BAM file or paired FASTQ files. All inputs are local file paths. For BAM input, the BAM is name-sorted and converted back to paired FASTQ before alignment. Both input routes then use the same STAR and Salmon steps.

run_cmd <- function(command, args) {
  status <- system2(command, args)
  if (status != 0) stop("Command failed: ", command)
}

process_powerup_rna <- function(
  sample_id,
  bam = NULL,
  fastq1 = NULL,
  fastq2 = NULL,
  threads = 8,
  keep_intermediates = TRUE
) {
  sample_dir <- file.path(work_dir, sample_id)
  dir.create(sample_dir, recursive = TRUE, showWarnings = FALSE)

  cat("\n============================================================\n")
  cat("Processing sample:", sample_id, "\n")
  cat("Sample directory:", sample_dir, "\n")
  cat("Threads:", threads, "\n")
  cat("============================================================\n")

  bam_mode <- !is.null(bam)

  if (bam_mode) {
    cat("\nStage 1: preparing BAM input\n")
    cat("Input BAM:", bam, "\n")

    fastq1 <- file.path(sample_dir, paste0(sample_id, "_R1.fastq"))
    fastq2 <- file.path(sample_dir, paste0(sample_id, "_R2.fastq"))
    name_bam <- file.path(sample_dir, paste0(sample_id, ".namesorted.bam"))

    cat("Name-sorting BAM with samtools.\n")
    run_cmd(samtools, c("sort", "-n", "-@", threads, "-o", name_bam, bam))

    cat("Converting BAM to paired FASTQ.\n")
    run_cmd(
      samtools,
      c(
        "fastq", "-@", threads,
        "-1", fastq1,
        "-2", fastq2,
        "-0", "/dev/null",
        "-s", "/dev/null",
        "-n", name_bam
      )
    )

    unlink(name_bam)

    cat("FASTQ 1:", fastq1, "\n")
    cat("FASTQ 2:", fastq2, "\n")
  } else {
    cat("\nStage 1: preparing paired FASTQ input\n")

    if (is.null(fastq1) || is.null(fastq2)) stop("Provide either BAM or paired FASTQs.")

    cat("FASTQ 1:", fastq1, "\n")
    cat("FASTQ 2:", fastq2, "\n")
  }

  cat("\nStage 2: STAR alignment\n")

  star_prefix <- file.path(sample_dir, paste0(sample_id, "."))

  star_args <- c(
    "--runThreadN", threads,
    "--genomeDir", star_index,
    "--readFilesIn", fastq1, fastq2,
    "--outFileNamePrefix", star_prefix,
    "--outSAMtype", "BAM", "Unsorted",
    "--quantMode", "TranscriptomeSAM"
  )

  compressed <- grepl("\\.gz$", fastq1) && grepl("\\.gz$", fastq2)
  if (compressed) star_args <- c(star_args, "--readFilesCommand", "zcat")

  cat("Starting STAR.\n")
  run_cmd(star, star_args)
  cat("STAR alignment complete.\n")

  cat("\nStage 3: Salmon transcript quantification\n")

  transcript_bam <- paste0(star_prefix, "Aligned.toTranscriptome.out.bam")
  salmon_dir <- file.path(sample_dir, "salmon")

  run_cmd(
    salmon,
    c(
      "quant",
      "-t", transcript_fasta,
      "-l", "IU",
      "-a", transcript_bam,
      "-p", threads,
      "-o", salmon_dir
    )
  )

  cat("Salmon quantification complete.\n")

  cat("\nStage 4: aggregating transcript TPM to genes\n")

  quant <- fread(file.path(salmon_dir, "quant.sf"))
  quant[, transcript_key := sub("\\.[0-9]+$", "", Name)]

  gene_tpm <- merge(
    quant[, .(transcript_key, TPM)],
    tx2gene[, .(transcript_key, gene_name)],
    by = "transcript_key"
  )

  gene_tpm <- gene_tpm[, .(TPM = sum(TPM)), by = gene_name]

  cat("Transcripts quantified:", nrow(quant), "\n")
  cat("Genes quantified:", nrow(gene_tpm), "\n")

  cat("\nStage 5: creating POWERUP expression matrix\n")

  matched <- feature_map[gene_name %in% gene_tpm$gene_name]
  unmatched <- feature_map[!gene_name %in% gene_tpm$gene_name]

  tpm <- gene_tpm$TPM[match(matched$gene_name, gene_tpm$gene_name)]
  expression <- log2(tpm + 1)

  user_matrix <- as.data.table(as.list(expression))
  setnames(user_matrix, matched$feature)
  user_matrix[, cell_line := sample_id]
  setcolorder(user_matrix, c("cell_line", matched$feature))

  gene_tpm_file <- file.path(sample_dir, paste0(sample_id, "_gene_tpm.tsv"))
  output_file <- file.path(sample_dir, paste0(sample_id, "_powerup_expression.csv"))
  unmatched_file <- file.path(sample_dir, "unmatched_powerup_features.tsv")

  fwrite(gene_tpm, gene_tpm_file, sep = "\t")
  fwrite(user_matrix, output_file)
  fwrite(unmatched, unmatched_file, sep = "\t")

  cat("POWERUP features matched:", nrow(matched), "of", nrow(feature_map), "\n")
  cat("POWERUP features unmatched:", nrow(unmatched), "\n")
  cat("Gene TPM file:", gene_tpm_file, "\n")
  cat("POWERUP matrix:", output_file, "\n")

  if (!keep_intermediates) {
    cat("\nStage 6: removing intermediate files\n")

    unlink(paste0(star_prefix, "Aligned.out.bam"))
    unlink(transcript_bam)
    if (bam_mode) unlink(c(fastq1, fastq2))

    cat("Intermediate cleanup complete.\n")
  }

  cat("\n============================================================\n")
  cat("Sample complete:", sample_id, "\n")
  cat("POWERUP features:", ncol(user_matrix) - 1, "\n")
  cat("STAR log:", paste0(star_prefix, "Log.final.out"), "\n")
  cat("POWERUP matrix:", output_file, "\n")
  cat("============================================================\n")

  user_matrix
}

The function writes three useful final files for each sample:

  • <sample>_gene_tpm.tsv: gene-level TPM before transformation and feature filtering;
  • <sample>_powerup_expression.csv: the final one-row POWERUP expression matrix;
  • unmatched_powerup_features.tsv: reference features that were not represented in the quantified genes.

STAR and Salmon also write their standard logs and quantification files within the sample directory. Set keep_intermediates = FALSE if the large alignment files and BAM-derived FASTQs do not need to be retained after successful processing.

Starting from a BAM file

For BAM input, provide a local BAM path and a sample identifier:

sample_a <- process_powerup_rna(
  sample_id = "sample_A",
  bam = "/path/to/sample_A.bam"
)

The BAM route performs:

  1. name-sort the BAM with samtools;
  2. recover paired FASTQ reads;
  3. align the recovered reads to GRCh38 with STAR;
  4. quantify transcripts from the STAR transcriptome BAM with Salmon;
  5. sum transcript TPM values to genes;
  6. transform gene TPM as log2(TPM + 1);
  7. match and order genes using the reference expression features.

The existing alignments in the input BAM are not used directly for quantification. Realignment puts BAM-derived reads and FASTQ-derived reads through the same STAR and Salmon workflow.

Starting from paired FASTQ files

For paired FASTQ input, provide read 1 and read 2 directly:

sample_b <- process_powerup_rna(
  sample_id = "sample_B",
  fastq1 = "/path/to/sample_B_R1.fastq.gz",
  fastq2 = "/path/to/sample_B_R2.fastq.gz"
)

Both uncompressed FASTQs and paired .fastq.gz files are accepted. Gzip-compressed inputs are passed to STAR with zcat.

Combine multiple samples

Each call returns a one-row data.table. Samples processed against the same reference feature map can be combined into one POWERUP user matrix.

if (!identical(names(sample_a), names(sample_b))) {
  stop("The POWERUP expression matrices do not have identical columns.")
}

user_expression <- rbindlist(list(sample_a, sample_b))
fwrite(user_expression, file.path(work_dir, "powerup_expression.csv"))

user_expression

The resulting matrix has the expected structure:

cell_line,GENE1 (1234),GENE2 (5678),...
sample_A,4.21,0.83,...
sample_B,3.77,1.05,...

Expression values are log2(TPM + 1).

Use the expression matrix in POWERUP

Pass the combined matrix to prepare_powerup_data() as user_matrix. Here, expression and response are the reference datasets used for model training and desired_targets contains the responses to model.

prepared <- prepare_powerup_data(
  gene_expression = expression,
  response = response,
  targets = desired_targets,
  user_matrix = user_expression
)

The prepared object can then proceed through the standard POWERUP workflow:

models <- fit_powerup_models(prepared)
models <- add_powerup_predictions(models, prepared)
predictions <- summarize_predictions(models)

Important considerations

RNA-seq expression values depend on the genome and annotation, alignment and quantification settings, gene identifiers, and expression transformation. New samples should therefore be processed as comparably as possible to the reference expression data used to train the POWERUP models.

This workflow is designed for paired-end bulk RNA-seq and produces log2(TPM + 1) expression. If the reference models were trained using materially different RNA-seq preprocessing, use a preprocessing strategy appropriate to that reference rather than assuming this workflow will make the datasets equivalent.

The reference expression table determines the final feature names. Genes that cannot be matched are recorded in unmatched_powerup_features.tsv and are not added to the generated user matrix. prepare_powerup_data() subsequently intersects the available reference and user features before model fitting and prediction.

Summary

For paired-end bulk RNA-seq, this workflow prepares POWERUP input as:

FASTQ -> STAR -> Salmon -> transcript TPM -> gene TPM -> log2(TPM + 1) -> reference feature matching -> POWERUP user matrix

or, when starting from an existing BAM:

BAM -> recover paired FASTQ -> STAR -> Salmon -> transcript TPM -> gene TPM -> log2(TPM + 1) -> reference feature matching -> POWERUP user matrix

Continue with Preparing data for POWERUP for feature selection, user samples, held-out reference samples, and the structure of the prepared object used by the rest of the POWERUP workflow.

Session information 

## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.7.3
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] powerup_1.0.99
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.37     desc_1.4.3        R6_2.6.1          fastmap_1.2.0    
##  [5] xfun_0.51         cachem_1.1.0      knitr_1.50        htmltools_0.5.8.1
##  [9] rmarkdown_2.29    lifecycle_1.0.4   cli_3.6.5         sass_0.4.9       
## [13] pkgdown_2.2.0     textshaping_1.0.0 jquerylib_0.1.4   systemfonts_1.2.2
## [17] compiler_4.4.2    tools_4.4.2       ragg_1.5.1        bslib_0.9.0      
## [21] evaluate_1.0.3    yaml_2.3.10       jsonlite_2.0.0    rlang_1.1.6      
## [25] fs_1.6.5          htmlwidgets_1.6.4