BRB-seq: ultra-affordable high-throughput transcriptomics enabled by bulk RNA barcoding and sequencing

Abstract

Despite its widespread use, RNA-seq is still too laborious and expensive to replace RT-qPCR as the default gene expression analysis method. We present a novel approach, BRB-seq, which uses early multiplexing to produce 3' cDNA libraries for dozens of samples, requiring just 2 hours of hands-on time. BRB-seq has a comparable performance to the standard TruSeq approach while showing greater tolerance for lower RNA quality and being up to 25 times cheaper. We anticipate that BRB-seq will transform basic laboratory practice given its capacity to generate genome-wide transcriptomic data at a similar cost as profiling four genes using RT-qPCR.

Keywords: Barcoding; Gene expression; RNA-seq; Transcriptomics; qPCR.

Fig. 1

Global assessment of SCRB-seq’s performance for bulk RNA-seq. a Comparison of read alignment performances between TruSeq and five SCRB-seq datasets: one lymphoblastoid cell line (LCL; generated in-house), and four public datasets from [15, 18]. The no/multiple alignment values are derived from the STAR [35] alignment, and no gene/ambiguous and mapped to genes correspond to the annotation of the reads to the genes by Htseq [49]. b Total number of detected genes in the same LCL RNA samples by SCRB-seq and TruSeq at different detection thresholds (e.g., “Reads > 0” means that a gene is considered detected if it is covered by at least one read). c Evaluation of SCRB-seq’s performance relative to TruSeq using the data downsampled to 1M single-end reads and shown by the total number of identified DE genes and number of “true positive” DE genes. The latter represents a subset of DE genes identified using the full TruSeq 30M paired-end set; the error bars correspond to the variation produced by downsampled replicates (see the “Methods” section). d Assessment of the impact of the number of cycles during PCR pre-amplification of SCRB-seq libraries (downsampled to 1M single-end reads) prepared with BU3 primers. Performances were evaluated through variable quality measures: uniquely mapped reads, level of duplication, rate of MT-rRNA reads, and number of detected genes. e Assessment of the complexity of the libraries (downsampled to 100k single-end reads) obtained with different combinations of RT enzymes and DS cDNA generation procedures at various detection cutoffs (e.g., “Reads > 0” means that a gene is considered detected if it is covered by at least one read). f Read coverage across the gene body for different combinations of RT enzymes and DS cDNA generation procedures. Legend: DS cDNA, double-stranded cDNA; SE, single end; MMH, Maxima Fermentas Minus H Enzyme; SSII, Superscript II enzyme; SSS, second-strand synthesis using Nick translation; PCR, pre-amplification by polymerase chain reaction

Fig. 2

Schematic overview of the BRB-seq protocol. This schema highlights in details all steps of the final BRB-seq protocol. The bottom grayed window shows the final BRB-seq construct used for Illumina sequencing. The read Read1 and Read2 primers are used to sequence the barcode/UMI and cDNA fragment respectively. Index read (i7) is used to demultiplex Illumina libraries. Legend: DS cDNA, double-stranded cDNA

Fig. 3

BRB-seq’s overall performance relative to TruSeq. a Correlation of log2 read counts between technical replicates at t14 for the BRB-seq workflow (Pearson correlation r = 0.987). b Correlation of log2 read counts between BRB-seq and TruSeq (Pearson correlation r = 0.920). c Comparison of read alignment performances between BRB-seq and TruSeq. The no/multiple alignment values are derived from the STAR [35] alignment, and no gene/ambiguous and mapped to genes correspond to the annotation of the reads to the genes by Htseq [49]. d Comparison of library complexity between BRB-seq and TruSeq (e.g., “Reads > 0” means that a gene is considered detected if it is covered by at least one read). e Evaluation of BRB-seq’s performance relative to TruSeq using the data downsampled to 1M single-end reads and shown by the total number of identified DE genes and the number of “true positive” DE genes. The latter represents a subset of DE genes identified using the full TruSeq 30M paired-end set (see the “Methods” section). f The distribution of RPKM levels of expression of the DE genes detected (blue) or not detected (red) in the downsampled TruSeq (dotted) or BRB-seq (plain) that overlaps with the “gold standard” TruSeq ~ 30M paired-end reads. g The sequencing depth required for detecting genes with a given CPM expression level using TruSeq and BRB-seq libraries. A sequencing depth is considered sufficient if the gene is detected more than 95% of the time. h Power simulation analysis of public and in-house bulk SCRB-seq, BRB-seq, and TruSeq datasets (*p < 0.001; n.s. non-significant). i Correlation of expression values (normalized to HPRT1) determined by qPCR (in replicates, with 50 ng and 500 ng of total RNA used per RT), TruSeq and BRB-seq. Pearson’s r values are indicated. In all panels, for an unbiased comparison, all libraries were randomly downsampled to one million single-end reads (see the “Methods” section)

Fig. 4

BRB-seq multiplexing experiment and comparison with TruSeq. a Venn diagram showing the protein-coding genes detected (at least one read) across all 60 (TruSeq A) or 53 (TruSeq B) LCL samples after downsampling to 1M reads. b Distribution of counts per millions (CPM) of genes taken from every subset (corresponding color) of the Venn diagram shown in panel a. c Pearson’s correlations of log2 expressions, calculated sample by sample, i.e., of the same sample taken from two different dataset combinations (TruSeq A and B and BRB-seq). d Correlation heatmap showing in greater detail the individual LCL sample correlations between all three datasets (BRB-seq, TruSeq A, and TruSeq B). Highlighted in black are the three main clusters, showing, as expected, a clear separation by protocol (BRB-seq vs. TruSeq) or sequencing run (TruSeq A vs. B), overriding the relatively modest biological differences between 60 LCL samples, while maintaining an overall high correlation (Pearson’s r > 0.8). In all panels, all libraries were randomly downsampled to one million single-end reads for an unbiased comparison (see the “Methods” section)

Fig. 5

BRB-seq performance with fragmented RNA samples. a Pearson correlation between log2 read counts of intact (RNA quality number (RQN) = 8.9 and 9.8 for T0 and T14 respectively) versus fragmented samples (after 1 or 2 min of fragmentation). b Quality evaluation of BRB-seq libraries prepared with fragmented RNA samples (1 or 2 min fragmentation) compared with the intact RNA counterparts. For the analysis, the libraries were downsampled to 1M single-end reads (see the “Methods” section). “Max” threshold thus comes from the 1M downsampled intact RNA sample when compared to itself, without downsampling. Legend: RQN, RNA quality number (maximum is 10)

Fig. 6

The streamlined BRB-seq data analysis workflow and its low cost. a Schematic representation of the BRB-seq library post-sequencing data processing pipeline. It includes the BRB-seqTools module (available on github, see the “Methods” section) that can perform optional read trimming, alignment, sample demultiplexing, and generation of a count table. The count table can be further analyzed by standard algorithms or loaded into ASAP, a web-based analytical interface that facilitates data exploration and visualization. b The estimated per sample cost of library preparation for 96 samples for TruSeq and BRB-seq. Per sample cost of BRB-seq involving in-house made Tn5 or Nextera Tn5 is indicated