A suite of tools for the pre-processing of BRB-seq data (bulk RNA-seq)
BRB-seq command-line tools are provided as a single executable jar file. The .jar file contains all required materials and can be run on any terminal.
For the tools to run properly, you must have Java 1.8 installed.
To check your java version, open your terminal application and r un the following command:
java -versionIf the output looks something like java version "1.8.x", you are good to go. If not, you may need to update your version; see the Oracle Java website to download the latest JRE (for users) or JDK (for developers).
The software relies on Picard Tools, but the Picard JAR is already embedded in the released JAR, so no need to install it yourself.
After sequencing your BRB-seq libraries, you should obtain two fastq files per library:
- R1 fastq file: This file should contain the barcodes and UMIs of the BRBseq construct. Usually, the barcode comes before the UMI sequence. And the UMI sequence is optional.
- R2 fastq file: This file should contain the exact same number of reads (and read names) than the R1 file. Except that the sequence of the reads are the sequences representing the RNA fragments.
BRB-seq tools is a suite dedicated to help you analyze these data, until the generation of the output count/UMI matrix. For further analyses (filtering, normalization, dimension reduction, clustering, differential expression), we recommend using ASAP web portal that you can freely access at asap.epfl.ch.
To check that BRB-seq Tools is working properly, run the following command:
java -jar BRBseqTools.1.1.jarAs shown in the output of this command, BRB-seq tool suite allows the user to run 3 possible tools.
Demultiplex Create demultiplexed fastq files from R1+R2 fastq (use this if you want to process them independently, thus doing the next steps without BRBseqTools)
CreateDGEMatrix Create the DGE Matrix (counts + UMI) from R2 aligned BAM and R1 fastq
Trim For trimming BRBseq construct in R2 fastq file or demultiplexed fastq files
AnnotateBAM For annotating the R2 BAM file using UMIs/Barcodes from the R1 fastq and/or GTF and/or Barcode files
This tool is used when you need to generate all the fastq files corresponding to all your multiplexed samples. You end up with one fastq file per sample.
Options:
-r1 %s [Required] Path of R1 FastQ files (containing barcode and optionally UMIs) [can be gzipped or raw].
-r2 %s [Required] Path of R2 FastQ files (containing read sequence) [can be gzipped or raw].
-c %s [Required] Path of Barcode/Samplename mapping file¹.
-n %i Number of allowed difference with the barcode [Default = 1]. Ambiguous barcodes (same distance from two or more existing barcodes) will be automatically discarded.
-o %s Output folder
-p %s Barcode pattern/order found in the reads of the R1 FastQ file. Barcode names should match the barcode file (default = 'BU' i.e. barcode followed by the UMI).
'B' [required] is used for specifying the barcode position.
'U' [optional] is used for specifying a UMI value position.
Character '?' [optional] can be used to ignore specific nucleotides.
-UMI %i If your barcode pattern contains UMI ('U'), you should specify this parameter as the length of the UMI.
¹You can download/edit this example of barcode/samplename mapping file
Example:
java -jar BRBseqTools.1.1.jar Demultiplex -r1 lib_example_R1.fastq.gz -r2 lib_example_R2.fastq.gz -c lib_example_barcodes.txt -p BU -UMI 14or, if no UMI:
java -jar BRBseqTools.1.1.jar Demultiplex -r1 lib_example_R1.fastq.gz -r2 lib_example_R2.fastq.gz -c lib_example_barcodes.txt -p BNote: When using this tool, UMIs are kept as indexes in all output .fastq files
This tool is used when you don't need the intermediary .fastq & .bam files from all your multiplexed samples. It greatly simplifies the demultiplexing and analysis, and directly generates a workable count/UMI matrix.
Options:
-f %s [Required] Path of R1 FastQ file [can be gzipped or raw].
-b %s [Required] Path of R2 aligned BAM file [do not need to be sorted or indexed].
-c %s [Required] Path of Barcode/Samplename mapping file¹.
-gtf %s [Required] Path of GTF file [can be gzipped or raw].
-s %s Do you want to count only reads falling on same strand than gene? [no, yes, reverse] (default = yes since BRB-seq is stranded protocol)
-n %i Number of allowed difference with the barcode [ambiguous reads will be automatically discarded].
-o %s Output folder
-t %s Path of existing folder for storing temporary files
-chunkSize %i Maximum number of reads to be stored in RAM (default = 10000000)
-p %s Barcode pattern/order found in the reads of the R1 FastQ file. Barcode names should match the barcode file (default = 'BU' i.e. barcode followed by the UMI).
'B' [required] is used for specifying the barcode position.
'U' [optional] is used for specifying a UMI value position.
Character '?' [optional] can be used to ignore specific nucleotides.
-UMI %i If your barcode pattern contains UMI ('U'), you should specify this parameter as the length of the UMI.
¹You can download/edit this example of barcode/samplename mapping file
Example:
java -jar BRBseqTools.1.1.jar CreateDGEMatrix -f lib_example_R1.fastq.gz -b lib_example_R2.bam -c lib_example_barcodes.txt -gtf Homo_sapiens.GRCh38.90.gtf.gz -p BU -UMI 14Note: The original BRB-seq protocol contains a UMI construct. But you can also generate a library without UMIs, or even not sequence the UMIs from R2 read. If UMIs are present, both UMI and read count matrices will be generated. If not, only the read count table will be generated.
This tool is used when you need specific downstream analyses of your BAM files that are not handled by BRBseqTools. It creates an annotated BAM file containing the UMI/Barcodes and can also have annotated genes information. For example, it allows you to use pipelines such as Picard MarkDuplicates, that can use the UMI/Barcode information for better quantifying the duplicated reads (Picard options BARCODE_TAG=BC READ_ONE_BARCODE_TAG=BX)
Options:
-f %s [Required] Path of R1 FastQ file [can be gzipped or raw].
-b %s [Required] Path of R2 aligned BAM file [do not need to be sorted or indexed].
-c %s Path of Barcode/Samplename mapping file¹.
-gtf %s Path of GTF file [can be gzipped or raw].
-n %i Number of allowed difference with the barcode [ambiguous reads will be automatically discarded].
-s %s Do you want to annotate genes only in case where reads are falling on same strand? [no, yes, reverse] (default = yes since BRB-seq is stranded protocol
-o %s Output folder
-t %s Path of existing folder for storing temporary files
-chunkSize %i Maximum number of reads to be stored in RAM (default = 10000000)
-p %s Barcode pattern/order found in the reads of the R1 FastQ file. Barcode names should match the barcode file (default = 'BU' i.e. barcode followed by the UMI).
'B' [required] is used for specifying the barcode position.
'U' [optional] is used for specifying a UMI value position.
Character '?' [optional] can be used to ignore specific nucleotides.
-UMI %i If your barcode pattern contains UMI ('U'), you should specify this parameter as the length of the UMI.
¹You can download/edit this example of barcode/samplename mapping file
Example:
java -jar BRBseqTools.1.1.jar AnnotateBAM -f lib_example_R1.fastq.gz -b lib_example_R2.bam -p BU -UMI 14Note: Most of the options are optionals here. If not put, then the corresponding tag will simply not be added to the annotated BAM.
This tool is used to trim the BRB-seq construct & the polyA sequences from the R2 fastq file. It can also be used to trim individually all the demultiplexed .fastq files if you used the "Demultiplex" tool (in this case, use the -uniqueBarcode option)
Options:
-f %s [Required] Path of FastQ file to trim (or containing folder for processing all fastq files at once).
-o %s Output folder
-uniqueBarcode If the fastq file(s) contain(s) only one barcode (for e.g. after demultiplexing), this option can be used for searching the specific barcode (most occuring) in the construct and trimming it when present.
-polyA %i Trim polyA strings that have more than this length (without mismatch), and all 3' string that follows [default=6]
-minLength %i If resulting trimmed reads are < this number, it is removed from the output fastq file [default=10]
Example:
java -jar BRBseqTools.1.1.jar Trim -f lib_example_R2.fastq.gzNote: If you use STAR for alignment, this step is optional, as it will not change much the results of the alignment (our tests have shown that the improvement is real but very minor)
First if you don't already have an indexed genome, you need to build the index (for STAR)
# Download last release of homo sapiens .fasta file on Ensembl
wget ftp://ftp.ensembl.org/pub/release-90/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
# Unzip
gzip -d Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
# Download last release of homo sapiens .gtf file on Ensembl
wget ftp://ftp.ensembl.org/pub/release-90/gtf/homo_sapiens/Homo_sapiens.GRCh38.90.gtf.gz
# Unzip
gzip -d Homo_sapiens.GRCh38.90.gtf.gz
# Generate the STAR genome index (in 'STAR_Index' folder) => This require ~30G RAM
mkdir STAR_Index/
STAR --runMode genomeGenerate --genomeDir STAR_Index/ --genomeFastaFiles Homo_sapiens.GRCh38.dna.primary_assembly.fa --sjdbGTFfile Homo_sapiens.GRCh38.90.gtf
# Note: you can add the argument '--runThreadN 4' for running 4 threads in parallel (or more if your computer has enough cores)
# Note: if genome is not human, you need to add/tune the argument '--genomeSAindexNbases xx' with xx = log2(nbBasesInGenome)/2 - 1
# For e.g. for drosophila melanogaster xx ~ 12.43, thus you should add the argument '--genomeSAindexNbases 12'When the index is built, you will never need to rebuild it. Then you can use only the following script:
# (Optional) Trim the read containing the sequence fragments (generates a 'lib_example_R2.trimmed.fastq.gz' file)
java -jar BRBseqTools.1.1.jar Trim -f lib_example_R2.fastq.gz
# Create output folder
mkdir BAM/
# Align only the R2 fastq file (using STAR, no sorting/indexing is needed)
STAR --runMode alignReads --genomeDir STAR_Index/ --outFilterMultimapNmax 1 --readFilesCommand zcat --outSAMtype BAM Unsorted --outFileNamePrefix BAM/ --readFilesIn lib_example_R2.trimmed.fastq.gz
# Note: you can add the argument '--runThreadN 4' for running 4 threads in parallel (or more if your computer has enough cores)
# Note: the '--outFilterMultimapNmax 1' option is recommended for removing multiple mapping reads from the output BAM
# (optional) Rename the output aligned BAM
mv BAM/Aligned.out.bam BAM/lib_example_R2.bam
# Demultiplex and generate output count/UMI matrix
java -jar BRBseqTools.1.0.jar CreateDGEMatrix -f lib_example_R1.fastq.gz -b BAM/lib_example_R2.bam -c lib_example_barcodes.txt -gtf Homo_sapiens.GRCh38.90.gtf -p BU -UMI 14
# Note: This example suppose that R1 has barcode followed by 14bp UMI (see above for other cases)
# Note: 'lib_example_barcodes.txt' should be created by the user and should contain the mapping between the barcode and the sample name¹¹You can download/edit this example of barcode/samplename mapping file
Then you can load the generated 'output.dge.reads.txt' count matrix in R and performs your analyses (or 'output.dge.umis.txt' if you prefer working with UMIs). Or you can upload this file to https://asap.epfl.ch and run the analysis pipeline online.
- src: all source files required for compilation
- lib: all JAR dependencies required for compilation / execution
- releases: latest release of BRB-seq Tools
Vincent Gardeux - [email protected]