Single-Cell RNA-seq Workflow

Source: vignettes/single-cell-workflow.Rmd

This tutorial builds one cluster-level single-cell RNA-seq analysis through four cumulative configurations:

  1. short reads;
  2. short and long reads;
  3. short and long reads with positional-bias correction; and
  4. short and long reads with positional-bias correction and a GTEx prior.

Each configuration reuses the objects prepared in the previous sections, so the effect of each added source of information can be compared directly.

What This Tutorial Produces

By the end, you will have:

  • one mpaqt_index object for the reference transcriptome;
  • cluster-level short-read and long-read count lists;
  • four named lists of mpaqt_quant_result objects; and
  • a transcript-by-cluster TPM table for downstream analysis.

The four result lists are named sr_results, sr_lr_results, sr_lr_bias_results, and sr_lr_bias_prior_results.

Required Inputs

The short-read configuration requires:

  • a GTF annotation and matching transcriptome FASTA;
  • paired 10x Chromium FASTQ files, with barcode and UMI information in R1 and cDNA sequence in R2;
  • a cluster-assignment CSV; and
  • kallisto and bustools available on your system.

The cluster file must have barcodes in its first column and cluster labels in its second column. Additional columns are ignored. For example:

barcode,cluster
AAACCCAAGAAACACT-1,Astrocyte
AAACCCAAGAAACCAT-1,Neuron
AAACCCAAGAAACCCG-1,Neuron

The long-read configuration additionally requires a transcript-by-cell count matrix. Its first column contains transcript IDs, and its remaining column names are cell barcodes:

transcript_id,AAACCCAAGAAACACT-1,AAACCCAAGAAACCAT-1
ENST00000456328.2,2,0
ENST00000450305.2,0,3

The final configuration also uses the GTEx PCA table downloaded in Add the GTEx Prior.

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,
    threads = 8L
)

Set read_length to the length of the cDNA sequence represented by R2. Reuse the saved index in later runs:

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

See the index reference for other supported reference inputs.

Prepare Single-Cell Short Reads

Prepare cluster-level equivalence-class counts from the 10x FASTQ files and cluster assignments:

sr_list <- mpaqt_prepare_short_reads_sc(
    index = index,
    fastq_1 = "data/sample_R1.fastq.gz",
    fastq_2 = "data/sample_R2.fastq.gz",
    clusters_file = "data/clusters.csv",
    do_umi_dedup = TRUE,
    umi_dedup_mode = "fractional_proportional",
    output_dir = "results/short-read",
    technology = "10xv3",
    threads = 8L
)

sr_list contains one mpaqt_counts_sr object per cluster and an "all_clusters" object that pools their counts. With the settings above, each barcode-UMI molecule contributes a total count of one while its evidence can be distributed across equivalence classes.

UMI handling occurs at two separate stages. do_umi_dedup = TRUE collapses read records to molecule-level counts during preparation. Later, mpaqt_quant_sc() uses do_umi_correction = TRUE by default to normalize the short-read probability matrices during quantification.

Quantify with Short Reads

Run the first configuration for every entry in sr_list:

sr_results <- mpaqt_quant_sc(
    index = index,
    sr_counts_list = sr_list
)

The returned list uses the same cluster names as sr_list and also includes the pooled "all_clusters" result.

Prepare Single-Cell Long Reads

Aggregate the transcript-by-cell long-read matrix using the same cluster file:

lr_list <- mpaqt_prepare_long_reads_sc(
    index = index,
    count_file = "data/long_read_counts.csv",
    clusters_file = "data/clusters.csv",
    output_dir = "results/long-read"
)

lr_list contains one mpaqt_counts_lr object per cluster plus an "all_clusters" object. The short- and long-read lists must use the same cluster names. Cell barcodes in the count-matrix columns must match those in the cluster file exactly.

Quantify with Short and Long Reads

Add the long-read counts to obtain the second configuration:

sr_lr_results <- mpaqt_quant_sc(
    index = index,
    sr_counts_list = sr_list,
    lr_counts_list = lr_list
)

This call reuses the same cluster-level short-read evidence and adds transcript-level long-read evidence for each cluster.

Add Positional-Bias Correction

For a 3’-biased protocol or observed 3’ enrichment, add 3’ positional-bias correction to obtain the third configuration:

sr_lr_bias_results <- mpaqt_quant_sc(
    index = index,
    sr_counts_list = sr_list,
    lr_counts_list = lr_list,
    positional_bias = "3p"
)

MPAQT estimates one set of positional weights from the pooled sr_list[["all_clusters"]] counts and applies those weights to each cluster. For a 5’-biased protocol or observed 5’ enrichment, use positional_bias = "5p" instead.

Add the GTEx Prior

The GTEx prior table contains transcript expression embeddings derived from principal component analysis of GTEx v8 multi-tissue bulk RNA-seq data. Each row represents one transcript; the transcript_id column is followed by 200 features named PC1 through PC200. MPAQT uses these scores as predictive features when fitting its abundance prior.

Download gtex_PCA_output.csv.gz from the Zenodo record and decompress it:

wget -O gtex_PCA_output.csv.gz "https://zenodo.org/records/21478011/files/gtex_PCA_output.csv.gz?download=1"
gunzip gtex_PCA_output.csv.gz

Read the decompressed table, then pass it directly to prior_model in the fourth configuration:

gtex_prior_path <- "gtex_PCA_output.csv"
prior_data <- data.table::fread(gtex_prior_path)

sr_lr_bias_prior_results <- mpaqt_quant_sc(
    index = index,
    sr_counts_list = sr_list,
    lr_counts_list = lr_list,
    positional_bias = "3p",
    prior_model = prior_data,
    prior_start = 25L,
    convergence_start = 50L
)

MPAQT matches prior_data$transcript_id to the index after removing version suffixes. With these iteration settings, prior fitting begins at iteration 25 and convergence is not checked until iteration 50, allowing the fitted prior to affect the abundance estimates before the run can stop.

These four configurations are cumulative to make their comparison clear. The R API does not make long reads, positional-bias correction, and the prior intrinsically dependent: arguments can be omitted when a study calls for a different combination.

Inspect, Combine, and Save Results

Choose the result list for the configuration you want to take forward. Here we use the fourth configuration and exclude the pooled "all_clusters" result from the cluster-by-cluster output:

final_results <- sr_lr_bias_prior_results
cluster_ids <- setdiff(names(final_results), "all_clusters")
cluster_results <- final_results[cluster_ids]

head(tpm(cluster_results[[1]]))

combined_tpm <- do.call(
    combine_quant_results,
    c(
        unname(cluster_results),
        list(sample_names = cluster_ids)
    )
)

head(combined_tpm)

combined_tpm has one transcript_id column followed by one TPM column per cluster. Save the per-cluster result objects, per-cluster TPM tables, and the combined table together:

mpaqt_save_result_sc(
    cluster_results,
    output_dir = "results/quantification",
    prefix = "sample"
)

Use sr_results, sr_lr_results, or sr_lr_bias_results in the same way when you want to inspect or save an earlier configuration.

Alternative Input Formats

If kallisto has already produced matching BUS and equivalence-class files, prepare the single-cell short-read list without rerunning pseudoalignment:

sr_list <- mpaqt_prepare_short_reads_sc(
    index = index,
    clusters_file = "data/clusters.csv",
    bus_file = "data/output.bus",
    ec_file = "data/matrix.ec",
    do_umi_dedup = TRUE,
    umi_dedup_mode = "fractional_proportional",
    output_dir = "results/short-read",
    technology = "10xv3",
    threads = 8L
)

Previously prepared cluster objects can be reconstructed into named lists from their stored sample_id values:

sr_files <- list.files(
    "results/short-read",
    pattern = "\\.short_read\\.rds$",
    full.names = TRUE
)
sr_list <- lapply(sr_files, mpaqt_read_short_read_counts)
names(sr_list) <- vapply(sr_list, function(x) x$sample_id, character(1))

lr_files <- list.files(
    "results/long-read",
    pattern = "\\.long_read\\.rds$",
    full.names = TRUE
)
lr_list <- lapply(lr_files, mpaqt_read_long_read_counts)
names(lr_list) <- vapply(lr_list, function(x) x$sample_id, character(1))

See the single-cell short-read and single-cell long-read references for complete input details.

Troubleshooting

Cluster names differ between lists. Prepare short and long reads from the same cluster file, or rename list entries so every short-read cluster has the corresponding long-read entry.

No long-read barcodes match a cluster. Compare the matrix column names with the first column of clusters.csv, including any -1 suffixes.

Transcript IDs do not match. Use the same reference release for the index and long-read counts. The GTEx table can be matched with or without transcript version suffixes.

kallisto, bustools, or GPBoost is unavailable. Follow the dependency checks in the installation article, then restart R if your environment changed.

The final result contains all_clusters. This pooled entry is useful for shared weight estimation. Remove it with setdiff() as shown before combining or saving per-cluster results.

Next Steps