Welcome to Salmon’s documentation!

Contents:

Requirements

Binary Releases

Pre-compiled binaries of the latest release of Salmon for a number different platforms are available available under the Releases tab of Salmon’s GitHub repository. You should be able to get started quickly by finding a binary from the list that is compatible with your platform. Additionally, you can obtain a Docker image of the latest version from DockerHub using:

> docker pull combinelab/salmon

Requirements for Building from Source

• A C++11 conformant compiler (currently tested with GCC>=4.7 and Clang>=3.4)

• CMake. Salmon uses the CMake build system to check, fetch and install dependencies, and to compile and install Salmon. CMake is available for all major platforms (though Salmon is currently unsupported on Windows.)

Installation

After downloading the Salmon source distribution and unpacking it, change into the top-level directory:

> cd salmon

Then, create and out-of-source build directory and change into it:

> mkdir build
> cd build

Salmon makes extensive use of Boost. We recommend installing the most recent version (1.55) systemwide if possible. If Boost is not installed on your system, the build process will fetch, compile and install it locally. However, if you already have a recent version of Boost available on your system, it make sense to tell the build system to use that.

If you have Boost installed you can tell CMake where to look for it. Likewise, if you already have Intel’s Threading Building Blocks library installed, you can tell CMake where it is as well. The flags for CMake are as follows:

• -DFETCH_BOOST=TRUE – If you don’t have Boost installed (or have an older version of it), you can provide the FETCH_BOOST flag instead of the BOOST_ROOT variable, which will cause CMake to fetch and build Boost locally.

• -DBOOST_ROOT=<boostdir> – Tells CMake where an existing installtion of Boost resides, and looks for the appropritate version in <boostdir>. This is the top-level directory where Boost is installed (e.g. /opt/local).

• -DTBB_INSTALL_DIR=<tbbroot> – Tells CMake where an existing installation of Intel’s TBB is installed (<tbbroot>), and looks for the apropriate headers and libraries there. This is the top-level directory where TBB is installed (e.g. /opt/local).

• -DCMAKE_INSTALL_PREFIX=<install_dir> – <install_dir> is the directory to which you wish Salmon to be installed. If you don’t specify this option, it will be installed locally in the top-level directory (i.e. the directory directly above “build”).

There are a number of other libraries upon which Salmon depends, but CMake should fetch these for you automatically.

Setting the appropriate flags, you can then run the CMake configure step as follows:

> cmake [FLAGS] ..

The above command is the cmake configuration step, which should complain if anything goes wrong. Next, you have to run the build step. Depending on what libraries need to be fetched and installed, this could take a while (specifically if the installation needs to install Boost). To start the build, just run make.

> make

If the build is successful, the appropriate executables and libraries should be created. There are two points to note about the build process. First, if the build system is downloading and compiling boost, you may see a large number of warnings during compilation; these are normal. Second, note that CMake has colored output by default, and the steps which create or link libraries are printed in red. This is the color chosen by CMake for linking messages, and does not denote an error in the build process.

Finally, after everything is built, the libraries and executable can be installed with:

> make install

You can test the installation by running

> make test

This should run a simple test and tell you if it succeeded or not.

Salmon

Salmon is a tool for wicked-fast transcript quantification from RNA-seq data. It requires a set of target transcripts (either from a reference or de-novo assembly) to quantify. All you need to run Salmon is a FASTA file containing your reference transcripts and a (set of) FASTA/FASTQ file(s) containing your reads. Optionally, Salmon can make use of pre-computed alignments (in the form of a SAM/BAM file) to the transcripts rather than the raw reads.

The mapping-based mode of Salmon runs in two phases; indexing and quantification. The indexing step is independent of the reads, and only needs to be run once for a particular set of reference transcripts. The quantification step, obviously, is specific to the set of RNA-seq reads and is thus run more frequently. For a more complete description of all available options in Salmon, see below.

Note

Selective alignment

Selective alignment, first introduced by the --validateMappings flag in salmon, and now the default mapping strategy (in version 1.0.0 forward), is a major feature enhancement introduced in recent versions of salmon. When salmon is run with selective alignment, it adopts a considerably more sensitive scheme that we have developed for finding the potential mapping loci of a read, and score potential mapping loci using the chaining algorithm introduced in minimap2 5. It scores and validates these mappings using the score-only, SIMD, dynamic programming algorithm of ksw2 6. Finally, we recommend using selective alignment with a decoy-aware transcriptome, to mitigate potential spurious mapping of reads that actually arise from some unannotated genomic locus that is sequence-similar to an annotated transcriptome. The selective-alignment algorithm, the use of a decoy-aware transcriptome, and the influence of running salmon with different mapping and alignment strategies is covered in detail in the paper Alignment and mapping methodology influence transcript abundance estimation.

The use of selective alignment implies the use of range factorization, as mapping scores become very meaningful with this option. Selective alignment can improve the accuracy, sometimes considerably, over the faster, but less-precise mapping algorithm that was previously used. Also, there are a number of options and flags that allow the user to control details about how the scoring is carried out, including setting match, mismatch, and gap scores, and choosing the minimum score below which an alignment will be considered invalid, and therefore not used for the purposes of quantification.

The alignment-based mode of Salmon does not require indexing. Rather, you can simply provide Salmon with a FASTA file of the transcripts and a SAM/BAM file containing the alignments you wish to use for quantification.

Salmon is, and will continue to be, freely and actively supported on a best-effort basis. If you are in need of industrial-grade technical support, please consider the options at oceangenomics.com/support.

Using Salmon

As mentioned above, there are two “modes” of operation for Salmon. The first, requires you to build an index for the transcriptome, but then subsequently processes reads directly. The second mode simply requires you to provide a FASTA file of the transcriptome and a .sam or .bam file containing a set of alignments.

Note

Read / alignment order

Salmon, like eXpress 1, uses a streaming inference method to perform transcript-level quantification. One of the fundamental assumptions of such inference methods is that observations (i.e. reads or alignments) are made “at random”. This means, for example, that alignments should not be sorted by target or position. If your reads or alignments do not appear in a random order with respect to the target transcripts, please randomize / shuffle them before performing quantification with Salmon.

Note

Number of Threads

The number of threads that Salmon can effectively make use of depends upon the mode in which it is being run. In alignment-based mode, the main bottleneck is in parsing and decompressing the input BAM file. We make use of the Staden IO library for SAM/BAM/CRAM I/O (CRAM is, in theory, supported, but has not been thoroughly tested). This means that multiple threads can be effectively used to aid in BAM decompression. However, we find that throwing more than a few threads at file decompression does not result in increased processing speed. Thus, alignment-based Salmon will only ever allocate up to 4 threads to file decompression, with the rest being allocated to quantification. If these threads are starved, they will sleep (the quantification threads do not busy wait), but there is a point beyond which allocating more threads will not speed up alignment-based quantification. We find that allocating 8 — 12 threads results in the maximum speed, threads allocated above this limit will likely spend most of their time idle / sleeping.

For quasi-mapping-based Salmon, the story is somewhat different. Generally, performance continues to improve as more threads are made available. This is because the determination of the potential mapping locations of each read is, generally, the slowest step in quasi-mapping-based quantification. Since this process is trivially parallelizable (and well-parallelized within Salmon), more threads generally equates to faster quantification. However, there may still be a limit to the return on invested threads, when Salmon can begin to process fragments more quickly than they can be provided via the parser.

Preparing transcriptome indices (mapping-based mode)

One of the novel and innovative features of Salmon is its ability to accurately quantify transcripts without having previously aligned the reads using its fast, built-in selective-alignment mapping algorithm. Further details about the selective alignment algorithm can be found here.

If you want to use Salmon in mapping-based mode, then you first have to build a salmon index for your transcriptome. Assume that transcripts.fa contains the set of transcripts you wish to quantify. We generally recommend that you build a decoy-aware transcriptome file.

There are two options for generating a decoy-aware transcriptome:

• The first is to compute a set of decoy sequences by mapping the annotated transcripts you wish to index against a hard-masked version of the organism’s genome. This can be done with e.g. MashMap2, and we provide some simple scripts to greatly simplify this whole process. Specifically, you can use the generateDecoyTranscriptome.sh script, whose instructions you can find in this README.

• The second is to use the entire genome of the organism as the decoy sequence. This can be done by concatenating the genome to the end of the transcriptome you want to index and populating the decoys.txt file with the chromosome names. Detailed instructions on how to prepare this type of decoy sequence is available here. This scheme provides a more comprehensive set of decoys, but, obviously, requires considerably more memory to build the index.

Finally, pre-built versions of both the partial decoy and full decoy (i.e. using the whole genome) salmon indices for some common organisms are available via refgenie here.

If you are not using a pre-computed index, you run the salmon indexer as so:

> ./bin/salmon index -t transcripts.fa -i transcripts_index --decoys decoys.txt -k 31

This will build the mapping-based index, using an auxiliary k-mer hash over k-mers of length 31. While the mapping algorithms will make used of arbitrarily long matches between the query and reference, the k size selected here will act as the minimum acceptable length for a valid match. Thus, a smaller value of k may slightly improve sensitivity. We find that a k of 31 seems to work well for reads of 75bp or longer, but you might consider a smaller k if you plan to deal with shorter reads. Also, a shorter value of k may improve sensitivity even more when using selective alignment (enabled via the –validateMappings flag). So, if you are seeing a smaller mapping rate than you might expect, consider building the index with a slightly smaller k.

Quantifying in mapping-based mode

Then, you can quantify any set of reads (say, paired-end reads in files reads1.fq and reads2.fq) directly against this index using the Salmon quant command as follows:

> ./bin/salmon quant -i transcripts_index -l <LIBTYPE> -1 reads1.fq -2 reads2.fq --validateMappings -o transcripts_quant

If you are using single-end reads, then you pass them to Salmon with the -r flag like:

> ./bin/salmon quant -i transcripts_index -l <LIBTYPE> -r reads.fq --validateMappings -o transcripts_quant

Note

Order of command-line parameters

The library type -l should be specified on the command line before the read files (i.e. the parameters to -1 and -2, or -r). This is because the contents of the library type flag is used to determine how the reads should be interpreted.

You can, of course, pass a number of options to control things such as the number of threads used or the different cutoffs used for counting reads. Just as with the alignment-based mode, after Salmon has finished running, there will be a directory called salmon_quant, that contains a file called quant.sf containing the quantification results.

Providing multiple read files to Salmon

Often, a single library may be split into multiple FASTA/Q files. Also, sometimes one may wish to quantify multiple replicates or samples together, treating them as if they are one library. Salmon allows the user to provide a space-separated list of read files to all of it’s options that expect input files (i.e. -r, -1, -2). When the input is paired-end reads, the order of the files in the left and right lists must be the same. There are a number of ways to provide salmon with multiple read files, and treat these as a single library. For the examples below, assume we have two replicates lib_1 and lib_2. The left and right reads for lib_1 are lib_1_1.fq and lib_1_2.fq, respectively. The left and right reads for lib_2 are lib_2_1.fq and lib_2_2.fq, respectively. The following are both valid ways to input these reads to Salmon:

> salmon quant -i index -l IU -1 lib_1_1.fq lib_2_1.fq -2 lib_1_2.fq lib_2_2.fq --validateMappings -o out

> salmon quant -i index -l IU -1 <(cat lib_1_1.fq lib_2_1.fq) -2 <(cat lib_1_2.fq lib_2_2.fq) --validateMappings -o out

Similarly, both of these approaches can be adopted if the files are gzipped as well:

> salmon quant -i index -l IU -1 lib_1_1.fq.gz lib_2_1.fq.gz -2 lib_1_2.fq.gz lib_2_2.fq.gz --validateMappings -o out

> salmon quant -i index -l IU -1 <(gunzip -c lib_1_1.fq.gz lib_2_1.fq.gz) -2 <(gunzip -c lib_1_2.fq.gz lib_2_2.fq.gz) --validateMappings -o out

In each pair of commands, the first command lets Salmon natively parse the files, while the latter command creates, on-the-fly, an input stream that consists of the concatenation of both files. Both methods work, and are acceptable ways to merge the files. The latter method (i.e. process substitution) allows more complex processing to be done to the reads in the substituted process before they are passed to Salmon as input, and thus, in some situations, is more versatile.

Note

Interleaved FASTQ files

Salmon does not currently have built-in support for interleaved FASTQ files (i.e., paired-end files where both pairs are stored in the same file). We provide a script that can be used to run salmon with interleaved input. However, this script assumes that the input reads are perfectly synchronized. That is, the input cannot contain any un-paired reads.

Quantifying in alignment-based mode

Say that you’ve prepared your alignments using your favorite aligner and the results are in the file aln.bam, and assume that the sequence of the transcriptome you want to quantify is in the file transcripts.fa. You would run Salmon as follows:

> ./bin/salmon quant -t transcripts.fa -l <LIBTYPE> -a aln.bam -o salmon_quant

The <LIBTYPE> parameter is described below and is shared between both modes of Salmon. After Salmon has finished running, there will be a directory called salmon_quant, that contains a file called quant.sf. This contains the quantification results for the run, and the columns it contains are similar to those of Sailfish (and self-explanatory where they differ).

For the full set of options that can be passed to Salmon in its alignment-based mode, and a description of each, run salmon quant --help-alignment.

Note

Genomic vs. Transcriptomic alignments

Salmon expects that the alignment files provided are with respect to the transcripts given in the corresponding FASTA file. That is, Salmon expects that the reads have been aligned directly to the transcriptome (like RSEM, eXpress, etc.) rather than to the genome (as does, e.g. Cufflinks). If you have reads that have already been aligned to the genome, there are currently 3 options for converting them for use with Salmon. First, you could convert the SAM/BAM file to a FAST{A/Q} file and then use the lightweight-alignment-based mode of Salmon described below. Second, given the converted FASTA{A/Q} file, you could re-align these converted reads directly to the transcripts with your favorite aligner and run Salmon in alignment-based mode as described above. Third, you could use a tool like sam-xlate to try and convert the genome-coordinate BAM files directly into transcript coordinates. This avoids the necessity of having to re-map the reads. However, we have very limited experience with this tool so far.

Multiple alignment files

If your alignments for the sample you want to quantify appear in multiple .bam/.sam files, then you can simply provide the Salmon -a parameter with a (space-separated) list of these files. Salmon will automatically read through these one after the other quantifying transcripts using the alignments contained therein. However, it is currently the case that these separate files must (1) all be of the same library type and (2) all be aligned with respect to the same reference (i.e. the @SQ records in the header sections must be identical).

Description of some important options

Salmon exposes a number of useful optional command-line parameters to the user. The particularly important ones are explained here, but you can always run salmon quant -h to see them all.

--mimicBT2

This flag is a “meta-flag” that sets the parameters related to mapping and selective alignment to mimic alignment using Bowtie2 (with the flags --no-discordant and --no-mixed), but using the default scoring scheme and allowing both mismatches and indels in alignments.

--mimicStrictBT2

This flag is a “meta-flag” that sets the parameters related to mapping and selective alignment to mimic alignment using Bowtie2 (with the flags suggested by RSEM), but using the default scoring scheme and allowing both mismatches and indels in alignments. These setting essentially disallow indels in the resulting alignments.

--meta

As with the flags described above, this is a “meta-flag” that simply enables some options that may make more sense when quantifying metagenomic data. Specifically, the --meta flag sets the following options:

• The abundance optimization is initialized from the uniform distribution (compared to the default of using a weighted combination of the uniform intialization and the abundances learned during the online optimization)

• Rich equivalence classes are disabled. Using rich equivalence classes with metagenomic data should not be particularly problematic, but since they have been developed and tested most in the context of bulk RNA-seq quantification, they are currently disabled under this flag.

• The EM algorithm is used for abundance optimization instead of the default VBEM optimization. Neither is universally better than the other, but the parameters for the VBEM (e.g. the prior size and type) are set based on typical bulk RNA-seq transcriptome samples, and so may be less appropriate in the metagenomic context. Hence the --meta flags opts for the basic EM algorithm instead.

--recoverOrphans

This flag (which should only be used in conjunction with selective alignment), performs orphan “rescue” for reads. That is, if mappings are discovered for only one end of a fragment, or if the mappings for the ends of the fragment don’t fall on the same transcript, then this flag will cause salmon to look upstream or downstream of the discovered mapping (anchor) for a match for the opposite end of the given fragment. This is done by performing “infix” alignment within the maximum fragment length upstream of downstream of the anchor mapping using edlib.

--hardFilter

This flag (which should only be used with selective alignment) turns off soft filtering and range-factorized equivalence classes, and removes all but the equally highest scoring mappings from the equivalence class label for each fragment. While we recommend using soft filtering (the default) for quantification, this flag can produce easier-to-understand equivalence classes if that is the primary object of study.

--skipQuant

Related to the above, this flag will stop execution before the actual quantification algorithm is run.

--allowDovetail

Dovetailing mappings and alignments are considered discordant and discarded by default — this is the same behavior that is adopted by default in Bowtie2. This is a change from the older behavior of salmon where dovetailing mappings were considered concordant and counted by default. If you wish to consider dovetailing mappings as concordant (the previous behavior), you can do so by passing the flag to salmon quant. Exotic library types (e.g. MU, MSF, MSR) are no longer supported. If you need support for such a library type, please submit a feature request describing the use-case.

-p / --threads

The number of threads that will be used for quasi-mapping, quantification, and bootstrapping / posterior sampling (if enabled). Salmon is designed to work well with many threads, so, if you have a sufficient number of processors, larger values here can speed up the run substantially.

Note

Default number of threads

The default behavior is for Salmon to probe the number of available hardware threads and to use this number. Thus, if you want to use fewer threads (e.g., if you are running multiple instances of Salmon simultaneously), you will likely want to set this option explicitly in accordance with the desired per-process resource usage.

--dumpEq

If Salmon is passed the --dumpEq option, it will write a file in the auxiliary directory, called eq_classes.txt that contains the equivalence classes and corresponding counts that were computed during quasi-mapping. The file has a format described in Equivalence class file.

--incompatPrior

This parameter governs the a priori probability that a fragment mapping or aligning to the reference in a manner incompatible with the prescribed library type is nonetheless the correct mapping. Note that Salmon sets this value, by default, to a small but non-zero probability. This means that if an incompatible mapping is the only mapping for a fragment, Salmon will still assign this fragment to the transcript. This default behavior is different than programs like RSEM, which assign incompatible fragments a 0 probability (i.e., incompatible mappings will be discarded). If you wish to obtain this behavior, so that only compatible mappings will be considered, you can set --incompatPrior 0.0. This will cause Salmon to only consider mappings (or alignments) that are compatible with the prescribed or inferred library type.

--fldMean

Note : This option is only important when running Salmon with single-end reads.

Since the empirical fragment length distribution cannot be estimated from the mappings of single-end reads, the --fldMean allows the user to set the expected mean fragment length of the sequencing library. This value will affect the effective length correction, and hence the estimated effective lengths of the transcripts and the TPMs. The value passed to --fldMean will be used as the mean of the assumed fragment length distribution (which is modeled as a truncated Gaussian with a standard deviation given by --fldSD).

--fldSD

Note : This option is only important when running Salmon with single-end reads.

Since the empirical fragment length distribution cannot be estimated from the mappings of single-end reads, the --fldSD allows the user to set the expected standard deviation of the fragment length distribution of the sequencing library. This value will affect the effective length correction, and hence the estimated effective lengths of the transcripts and the TPMs. The value passed to --fldSD will be used as the standard deviation of the assumed fragment length distribution (which is modeled as a truncated Gaussian with a mean given by --fldMean).

--minScoreFraction

This value controls the minimum allowed score for a mapping to be considered valid. It matters only when --validateMappings has been passed to Salmon. The maximum possible score for a fragment is ms = read_len * ma (or ms = (left_read_len + right_read_len) * ma for paired-end reads). The argument to --minScoreFraction determines what fraction of the maximum score s a mapping must achieve to be potentially retained. For a minimum score fraction of f, only mappings with a score > f * s will be kept. Mappings with lower scores will be considered as low-quality, and will be discarded.

It is worth noting that mapping validation uses extension alignment. This means that the read need not map end-to-end. Instead, the score of the mapping will be the position along the alignment with the highest score. This is the score which must reach the fraction threshold for the read to be considered as valid.

--bandwidth

This flag (which is only meaningful in conjunction with selective alignment), sets the bandwidth parameter of the relevant calls to ksw2’s alignment function. This determines how wide an area around the diagonal in the DP matrix should be calculated.

--maxMMPExtension

This flag (which should only be used with selective alignment) limits the length that a mappable prefix of a fragment may be extended before another search along the fragment is started. Smaller values for this flag can improve the sensitivity of mapping, but could increase run time.

--ma

This value should be a positive (typically small) integer. It controls the score given to a match in the alignment between the query (read) and the reference.

--mp

This value should be a negative (typically small) integer. It controls the score given to a mismatch in the alignment between the query (read) and the reference.

--go

This value should be a positive (typically small) integer. It controls the score penalty attributed to an alignment for each new gap that is opened. The alignment score computed uses an affine gap penalty, so the penalty of a gap is go + l * ge where l is the gap length. The value of go should typically be larger than that of ge.

--ge

This value should be a positive (typically small) integer. It controls the score penalty attributed to the extension of a gap in an alignment. The alignment score computed uses an affine gap penalty, so the penalty of a gap is go + l * ge where l is the gap length. The value of ge should typically be smaller than that of go.

--rangeFactorizationBins

The range-factorization feature allows using a data-driven likelihood factorization, which can improve quantification accuracy on certain classes of “difficult” transcripts. Currently, this feature interacts best (i.e., yields the most considerable improvements) when either (1) using alignment-based mode and simultaneously enabling error modeling with --useErrorModel or (2) when enabling --validateMappings in quasi-mapping-based mode. The argument to this option is a positive integer x, that determines fidelity of the factorization. The larger x, the closer the factorization to the un-factorized likelihood, but the larger the resulting number of equivalence classes. A value of 1 corresponds to salmon’s traditional rich equivalence classes. We recommend 4 as a reasonable parameter for this option (it is what was used in the range-factorization paper).

--useEM

Use the “standard” EM algorithm to optimize abundance estimates instead of the variational Bayesian EM algorithm. The details of the VBEM algorithm can be found in 3. While both the standard EM and the VBEM produce accurate abundance estimates, there are some trade-offs between the approaches. Specifically, the sparsity of the VBEM algorithm depends on the prior that is chosen. When the prior is small, the VBEM tends to produce a sparser solution than the EM algorithm, while when the prior is relatively larger, it tends to estimate more non-zero abundances than the EM algorithm. It is an active research effort to analyze and understand all the tradeoffs between these different optimization approaches. Also, the VBEM tends to converge after fewer iterations, so it may result in a shorter runtime; especially if you are computing many bootstrap samples.

The default prior used in the VB optimization is a per-nucleotide prior of 1e-5 reads per-nucleotide. This means that a transcript of length 100000 will have a prior count of 1 fragment, while a transcript of length 50000 will have a prior count of 0.5 fragments, etc. This behavior can be modified in two ways. First, the prior itself can be modified via Salmon’s --vbPrior option. The argument to this option is the value you wish to place as the per-nucleotide prior. Additionally, you can modify the behavior to use a per-transcript rather than a per-nucleotide prior by passing the flag --perTranscriptPrior to Salmon. In this case, whatever value is set by --vbPrior will be used as the transcript-level prior, so that the prior count is no longer dependent on the transcript length. However, the default behavior of a per-nucleotide prior is recommended when using VB optimization.

Note

Choosing between EM and VBEM algorithms

As mentioned above, a thorough comparison of all of the benefits and detriments of the different algorithms is an ongoing area of research. However, preliminary testing suggests that the sparsity-inducing effect of running the VBEM with a small prior may lead, in general, to more accurate estimates (the current testing was performed mostly through simulation). Hence, the VBEM is the default, and the standard EM algorithm is accessed via the –useEM flag.

--numBootstraps

Salmon has the ability to optionally compute bootstrapped abundance estimates. This is done by resampling (with replacement) from the counts assigned to the fragment equivalence classes, and then re-running the optimization procedure, either the EM or VBEM, for each such sample. The values of these different bootstraps allows us to assess technical variance in the main abundance estimates we produce. Such estimates can be useful for downstream (e.g. differential expression) tools that can make use of such uncertainty estimates. This option takes a positive integer that dictates the number of bootstrap samples to compute. The more samples computed, the better the estimates of variance, but the more computation (and time) required.

--numGibbsSamples

Just as with the bootstrap procedure above, this option produces samples that allow us to estimate the variance in abundance estimates. However, in this case the samples are generated using posterior Gibbs sampling over the fragment equivalence classes rather than bootstrapping. We are currently analyzing these different approaches to assess the potential trade-offs in time / accuracy. The --numBootstraps and --numGibbsSamples options are mutually exclusive (i.e. in a given run, you must set at most one of these options to a positive integer.)

--seqBias

Passing the --seqBias flag to Salmon will enable it to learn and correct for sequence-specific biases in the input data. Specifically, this model will attempt to correct for random hexamer priming bias, which results in the preferential sequencing of fragments starting with certain nucleotide motifs. By default, Salmon learns the sequence-specific bias parameters using 1,000,000 reads from the beginning of the input. If you wish to change the number of samples from which the model is learned, you can use the --numBiasSamples parameter. Salmon uses a variable-length Markov Model (VLMM) to model the sequence specific biases at both the 5’ and 3’ end of sequenced fragments. This methodology generally follows that of Roberts et al. 2, though some details of the VLMM differ.

Note: This sequence-specific bias model is substantially different from the bias-correction methodology that was used in Salmon versions prior to 0.6.0. This model specifically accounts for sequence-specific bias, and should not be prone to the over-fitting problem that was sometimes observed using the previous bias-correction methodology.

--gcBias

Passing the --gcBias flag to Salmon will enable it to learn and correct for fragment-level GC biases in the input data. Specifically, this model will attempt to correct for biases in how likely a sequence is to be observed based on its internal GC content.

You can use the FASTQC software followed by MultiQC with transcriptome GC distributions to check if your samples exhibit strong GC bias, i.e. under-representation of some sub-sequences of the transcriptome. If they do, we obviously recommend using the --gcBias flag. Or you can simply run Salmon with --gcBias in any case, as it does not impair quantification for samples without GC bias, it just takes a few more minutes per sample. For samples with moderate to high GC bias, correction for this bias at the fragment level has been shown to reduce isoform quantification errors 4 3.

This bias is distinct from the primer biases learned with the --seqBias option. Though these biases are distinct, they are not completely independent. When both --seqBias and --gcBias are enabled, Salmon will learn a conditional fragment-GC bias model. By default, Salmon will learn 3 different fragment-GC bias models based on the GC content of the fragment start and end contexts, though this number of conditional models can be changed with the (hidden) option --conditionalGCBins. Likewise, the number of distinct fragment GC bins used to model the GC bias can be changed with the (hidden) option --numGCBins.

Note : In order to speed up the evaluation of the GC content of arbitrary fragments, Salmon pre-computes and stores the cumulative GC count for each transcript. This requires an extra 4-bytes per nucleotide. While this extra memory usage should normally be minor, it can nonetheless be controlled with the --reduceGCMemory option. This option replaces the per-nucleotide GC count with a rank-select capable bit vector, reducing the memory overhead from 4-bytes per nucleotide to ~1.25 bits, while being only marginally slower).

--posBias

Passing the --posBias flag to Salmon will enable modeling of a position-specific fragment start distribution. This is meant to model non-uniform coverage biases that are sometimes present in RNA-seq data (e.g. 5’ or 3’ positional bias). Currently, a small and fixed number of models are learned for different length classes of transcripts, as is done in Roberts et al. 2. Note: The positional bias model is relatively new, and is still undergoing testing. It replaces the previous –useFSPD option, which is now deprecated. This feature should be considered as experimental in the current release.

--biasSpeedSamp

When evaluating the bias models (the GC-fragment model specifically), Salmon must consider the probability of generating a fragment of every possible length (with a non-trivial probability) from every position on every transcript. This results in a process that is quadratic in the length of the transcriptome — though each evaluation itself is efficient and the process is highly parallelized.

It is possible to speed this process up by a multiplicative factor by considering only every i^th fragment length, and interpolating the intermediate results. The --biasSpeedSamp option allows the user to set this sampling factor. Larger values speed up effective length correction, but may decrease the fidelity of bias modeling. However, reasonably small values (e.g. 10 or less) should have only a minor effect on the computed effective lengths, and can considerably speed up effective length correction on large transcriptomes. The default value for --biasSpeedSamp is 5.

--writeUnmappedNames

Passing the --writeUnmappedNames flag to Salmon will tell Salmon to write out the names of reads (or mates in paired-end reads) that do not map to the transcriptome. When mapping paired-end reads, the entire fragment (both ends of the pair) are identified by the name of the first read (i.e. the read appearing in the _1 file). Each line of the unmapped reads file contains the name of the unmapped read followed by a simple flag that designates how the read failed to map completely. If fragments are aligned against a decoy-aware index, then fragments that are confidently assigned as decoys are written in this file followed by the d (decoy) flag. Apart from the decoy flag, for single-end reads, the only valid flag is u (unmapped). However, for paired-end reads, there are a number of different possibilities, outlined below:

u   = The entire pair was unmapped. No mappings were found for either the left or right read.
m1  = Left orphan (mappings were found for the left (i.e. first) read, but not the right).
m2  = Right orphan (mappings were found for the right read, but not the left).
m12 = Left and right orphans. Both the left and right read mapped, but never to the same transcript.

By reading through the file of unmapped reads and selecting the appropriate sequences from the input FASTA/Q files, you can build an “unmapped” file that can then be used to investigate why these reads may not have mapped (e.g. poor quality, contamination, etc.). Currently, this process must be done independently, but future versions of Salmon may provide a script to generate this unmapped FASTA/Q file from the unmapped file and the original inputs.

--writeMappings

Passing the --writeMappings argument to Salmon will have an effect only in mapping-based mode and only when using a quasi-index. When executed with the --writeMappings argument, Salmon will write out the mapping information that it then processes to quantify transcript abundances. The mapping information will be written in a SAM compatible format. If no options are provided to this argument, then the output will be written to stdout (so that e.g. it can be piped to samtools and directly converted into BAM format). Otherwise, this argument can optionally be provided with a filename, and the mapping information will be written to that file. Note: Because of the way that the boost options parser that we use works, and the fact that --writeMappings has an implicit argument of stdout, if you provide an explicit argument to --writeMappings, you must do so with the syntax --writeMappings=<outfile> rather than the synatx --writeMappings <outfile>. This is due to a limitation of the parser in how the latter could be interpreted.

Note

Compatible mappings

The mapping information is computed and written before library type compatibility checks take place, thus the mapping file will contain information about all mappings of the reads considered by Salmon, even those that may later be filtered out due to incompatibility with the library type.

What’s this LIBTYPE?

Salmon, has the user provide a description of the type of sequencing library from which the reads come, and this contains information about e.g. the relative orientation of paired-end reads. As of version 0.7.0, Salmon also has the ability to automatically infer (i.e. guess) the library type based on how the first few thousand reads map to the transcriptome. To allow Salmon to automatically infer the library type, simply provide -l A or --libType A to Salmon. Even if you allow Salmon to infer the library type for you, you should still read the section below, so that you can interpret how Salmon reports the library type it discovers.

Note

Automatic library type detection in alignment-based mode

The implementation of this feature involves opening the BAM file, peaking at the first record, and then closing it to determine if the library should be treated as single-end or paired-end. Thus, in alignment-based mode automatic library type detection will not work with an input stream. If your input is a regular file, everything should work as expected; otherwise, you should provide the library type explicitly in alignment-based mode.

Also the automatic library type detection is performed on the basis of the alignments in the file. Thus, for example, if the upstream aligner has been told to perform strand-aware mapping (i.e. to ignore potential alignments that don’t map in the expected manner), but the actual library is unstranded, automatic library type detection cannot detect this. It will attempt to detect the library type that is most consistent with the alignment that are provided.

The library type string consists of three parts: the relative orientation of the reads, the strandedness of the library, and the directionality of the reads.

The first part of the library string (relative orientation) is only provided if the library is paired-end. The possible options are:

I = inward
O = outward
M = matching

The second part of the read library string specifies whether the protocol is stranded or unstranded; the options are:

S = stranded
U = unstranded

If the protocol is unstranded, then we’re done. The final part of the library string specifies the strand from which the read originates in a strand-specific protocol — it is only provided if the library is stranded (i.e. if the library format string is of the form S). The possible values are:

F = read 1 (or single-end read) comes from the forward strand
R = read 1 (or single-end read) comes from the reverse strand

An example of some library format strings and their interpretations are:

IU (an unstranded paired-end library where the reads face each other)

SF (a stranded single-end protocol where the reads come from the forward strand)

OSR (a stranded paired-end protocol where the reads face away from each other,
     read1 comes from reverse strand and read2 comes from the forward strand)

Note

Strand Matching

Above, when it is said that the read “comes from” a strand, we mean that the read should align with / map to that strand. For example, for libraries having the OSR protocol as described above, we expect that read1 maps to the reverse strand, and read2 maps to the forward strand.

For more details on the library type, see Fragment Library Types.

Output

For details of Salmon’s different output files and their formats see Salmon Output File Formats.

Misc

Salmon, in quasi-mapping-based mode, can accept reads from FASTA/Q format files, or directly from gzipped FASTA/Q files (the ability to accept compressed files directly is a feature of Salmon 0.7.0 and higher). If your reads are compressed in a different format, you can still stream them directly to Salmon by using process substitution. Say in the quasi-mapping-based Salmon example above, the reads were actually in the files reads1.fa.bz2 and reads2.fa.bz2, then you’d run the following command to decompress the reads “on-the-fly”:

> ./bin/salmon quant -i transcripts_index -l <LIBTYPE> -1 <(bunzip2 -c reads1.fa.gz) -2 <(bunzip2 -c reads2.fa.bz2) -o transcripts_quant

and the bzipped files will be decompressed via separate processes and the raw reads will be fed into Salmon. Actually, you can use this same process even with gzip compressed reads (replacing bunzip2 with gunzip or pigz -d). Depending on the number of threads and the exact configuration, this may actually improve Salmon’s running time, since the reads are decompressed concurrently in a separate process when you use process substitution.

Finally, the purpose of making this software available is for people to use it and provide feedback. The paper describing this method is published in Nature Methods. If you have something useful to report or just some interesting ideas or suggestions, please contact us (rob.patro@cs.stonybrook.edu and/or carlk@cs.cmu.edu). If you encounter any bugs, please file a detailed bug report at the Salmon GitHub repository.

References

1

Roberts, Adam, and Lior Pachter. “Streaming fragment assignment for real-time analysis of sequencing experiments.” Nature Methods 10.1 (2013): 71-73.

2(1,2)

Roberts, Adam, et al. “Improving RNA-Seq expression estimates by correcting for fragment bias.” Genome Biology 12.3 (2011): 1.

3(1,2)

Patro, Rob, et al. “Salmon provides fast and bias-aware quantification of transcript expression.” Nature Methods (2017). Advanced Online Publication. doi: 10.1038/nmeth.4197..

4

Love, Michael I., Hogenesch, John B., Irizarry, Rafael A. “Modeling of RNA-seq fragment sequence bias reduces systematic errors in transcript abundance estimation.” Nature Biotechnology 34.12 (2016). doi: 10.1038/nbt.368.2..

5

Li, Heng. “Minimap2: pairwise alignment for nucleotide sequences.” Bioinformatics 34.18 (2018): 3094-3100.

6

Global alignment and alignment extension.

Alevin

Alevin is a tool — integrated with the salmon software — that introduces a family of algorithms for quantification and analysis of 3’ tagged-end single-cell sequencing data. Currently alevin supports the following single-cell protocols:

1. Drop-seq

2. 10x-Chromium v1/2/3

3. inDropV2

4. CELSeq 1/2

5. Quartz-Seq2

6. sci-RNA-seq3

Alevin works under the same indexing scheme (as salmon) for the reference, and consumes the set of FASTA/Q files(s) containing the Cellular Barcode(CB) + Unique Molecule identifier (UMI) in one read file and the read sequence in the other. Given just the transcriptome and the raw read files, alevin generates a cell-by-gene count matrix (in a fraction of the time compared to other tools).

Alevin works in two phases. In the first phase it quickly parses the read file containing the CB and UMI information to generate the frequency distribution of all the observed CBs, and creates a lightweight data-structure for fast-look up and correction of the CB. In the second round, alevin utilizes the read-sequences contained in the files to map the reads to the transcriptome, identify potential PCR/sequencing errors in the UMIs, and performs hybrid de-duplication while accounting for UMI collisions. Finally, a post-abundance estimation CB whitelisting procedure is done and a cell-by-gene count matrix is generated.

Using Alevin

Alevin requires the following minimal set of necessary input parameters (generally providing the flags in that order is recommended):

• -l: library type (same as salmon), we recommend using ISR for both Drop-seq and 10x-v2 chemistry.

• -1: CB+UMI file(s), alevin requires the path to the FASTQ file containing CB+UMI raw sequences to be given under this command line flag. Alevin also supports parsing of data from multiple files as long as the order is the same as in -2 flag.

• -2: Read-sequence file(s), alevin requires the path to the FASTQ file containing raw read-sequences to be given under this command line flag. Alevin also supports parsing of data from multiple files as long as the order is the same as in -1 flag.

• --dropseq / --chromium / --chromiumV3: the protocol, this flag tells the type of single-cell protocol of the input sequencing-library.

• -i: index, file containing the salmon index of the reference transcriptome, as generated by salmon index command.

• -p: number of threads, the number of threads which can be used by alevin to perform the quantification, by default alevin utilizes all the available threads in the system, although we recommend using ~10 threads which in our testing gave the best memory-time trade-off.

• -o: output, path to folder where the output gene-count matrix (along with other meta-data) would be dumped.

• --tgMap: transcript to gene map file, a tsv (tab-separated) file — with no header, containing two columns mapping of each transcript present in the reference to the corresponding gene (the first column is a transcript and the second is the corresponding gene).

Once all the above requirement are satisfied, alevin can be run using the following command:

> salmon alevin -l ISR -1 cb.fastq.gz -2 reads.fastq.gz --chromium  -i salmon_index_directory -p 10 -o alevin_output --tgMap txp2gene.tsv

Providing multiple read files to Alevin

Often, a single library may be split into multiple FASTA/Q files. Also, sometimes one may wish to quantify multiple replicates or samples together, treating them as if they are one library. Alevin allows the user to provide a space-separated list of files to all of it’s options that expect input files (i.e. -1, -2). The order of the files in the left and right lists must be the same. There are a number of ways to provide alevin with multiple CB and read files, and treat these as a single library. For the examples below, assume we have two replicates lib_A and lib_B. The left and right reads for lib_A are lib_A_cb.fq and lib_A_reads.fq, respectively. The left and right reads for lib_B are lib_B_cb.fq and lib_B_read.fq, respectively. The following are both valid ways to input these reads to alevin:

> salmon alevin -l ISR -1 lib_A_cb.fq lib_B_cb.fq -2 lib_A_read.fq lib_B_read.fq

Similarly, both of these approaches can be adopted if the files are gzipped as well:

> salmon alevin -l ISR -1 lib_A_cb.fq.gz lib_B_cb.fq.gz -2 lib_A_read.fq.gz lib_B_read.fq.gz

Note

Don’t provide data through input stream

To keep the time-memory trade-off within acceptable bounds, alevin performs multiple passes over the Cellular Barcode file. Alevin goes through the barcode file once by itself, and then goes through both the barcode and read files in unison to assign reads to cells using the initial barcode mapping. Since the pipe or the input stream can’t be reset to read from the beginning again, alevin can’t read in the barcodes, and might crash.

Description of important options

Alevin exposes a number of useful optional command-line parameters to the user. The particularly important ones are explained here, but you can always run salmon alevin -h to see them all.

-p / --numThreads

The number of threads that will be used for quantification. Alevin is designed to work well with many threads, so, if you have a sufficient number of processors, larger values here can speed up the run substantially. In our testing we found that usually 10 threads gives the best time-memory trade-off.

Note

Default number of threads

The default behavior is for Alevin to probe the number of available hardware threads and to use this number.

Thus, if you want to use fewer threads (e.g., if you are running multiple instances of Salmon simultaneously), you will likely want to set this option explicitly in accordance with the desired per-process resource usage.

--whitelist

This is an optional argument, where user can explicitly specify the whitelist CB to use for cell detection and CB sequence correction. If not given, alevin generates its own set of putative CBs.

Note

Not 10x 737k whitelist

This flag does not use the technologically defined whitelisted cellular barcodes provided by 10x, instead it’s a per experiment level list of subsampled cellular barcodes that need to quantified for consistency with other tools for example an input would be a file generated by cellranger with the name barcodes.tsv (uncompressed).

--noQuant

Generally used in parallel with --dumpfq. If Alevin is passed the --noQuant option, the pipeline will stop before starting the mapping. The general use-case is when we only need to concatenate the CB on the read-id of the second file and break the execution afterwards.

--noDedup

If Alevin is passed the --noDedup option, the pipeline only performs CB correction, maps the read-sequences to the transcriptome generating the interim data-structure of CB-EqClass-UMI-count. Used in parallel with --dumpBarcodeEq or --dumpBfh for the purposes of obtaining raw information or debugging.

--mrna

The list of mitochondrial genes which are to be used as a feature for CB whitelising naive Bayes classification.

Note

It is generally advisable to not use nuclear mitrochondrial genes in this as they can be both up and/or down regulated which might cancel out the usefulness of this feature. Please check issue #367 in salmon repo to know more about it.

--rrna

The list of ribosomal genes which are to be used as a feature for CB whitelising naive Bayes classification.

--dumpfq

Generally used along with --noQuant. If activated, alevin will sequence correct the CB and attach the corrected CB sequence to the read-id in the second file and dumps the result to standard-out (stdout).

--dumpBfh

Alevin internally uses a potentially big data-structure to concisely maintain all the required information for quantification. This flags dumps the full CB-EqClass-UMI-count data-structure for the purposed of allowing raw data analysis and debugging.

--dumpFeatures

If activated, alevin dumps all the features used by the CB classification and their counts at each cell level. It’s generally used in pair with other command line flags.

--dumpMtx

This flags is used to internally convert the default binary format of alevin for gene-count matrix into a human readable mtx (matrix market exchange) sparse format.

--forceCells

Alevin performs a heuristic based initial CB white-listing by finding the knee in the distribution of the CB frequency. Although knee finding algorithm works pretty well in most of the case, it sometimes over shoot and results in very less number of CB. With this flag, by looking at the CB frequency distribution, a user can explicitly specify the number of CB to consider for initial white-listing.

--expectCells

Just like forceCells flag, it’s yet another way of skipping the knee calculation heuristics, if it’s failing. This command line flag uses the cellranger type white-listing procedure. As specified in their algorithm overview page, “All barcodes whose total UMI counts exceed m/10 are called as cells”, where m is the frequency of the top 1% cells as specified by the parameter of this command line flag.

--numCellBootstraps

Alevin provides an estimate of the inferential uncertainty in the estimation of per cell level gene count matrix by performing bootstrapping of the reads in per-cell level equivalence classes. This command line flag informs Alevin to perform certain number of bootstrap and generate the mean and variance of the count matrix. This option generates three additional file, namely, quants_mean_mat.gz, quants_var_mat.gz and quants_boot_rows.txt. The format of the files stay the same as quants_mat.gz while the row order is saved in quants_boot_rows.txt and the column order is stays the same as in file quants_mat_cols.txt.

Note

Alevin can also dump the full bootstrap cell-gene count matrix of a experiment. To generate inferential replicates of the experiemnt, –numCellBootstraps has to be paired with –dumpFeatures which generates a file with name quants_boot_mat.gz. The output format is the same as quants_mat.gz and we fit the 3D cube of the cell-inference-gene counts in 2D as follows: if an experiment has C cells, G genes and N inferential replicates; alevin output file quants_boot_mat.gz would contain C*N rows and G columns while, starting from the top, the first N rows would represent first cell and it’s N inferential replicate. For more information on importing and using inferential replicates for single-cell data in generating accurate differential expression analysis, check out tximport and our Swish paper.

--debug

Alevin peforms intelligent white-listing downstream of the quantification pipeline and has to make some assumptions like chosing a fraction of reads to learn low confidence CB and in turn might erroneously exit – if the data results in no mapped or deduplicated reads to a CB in low confidence region. The problem doesn’t happen when provided with external whitelist but if there is an error and the user is aware of this being just a warning, the error can be skipped by running Alevin with this flag.

--minScoreFraction

This value controls the minimum allowed score for a mapping to be considered valid. It matters only when --validateMappings has been passed to Salmon. The maximum possible score for a fragment is ms = read_len * ma (or ms = (left_read_len + right_read_len) * ma for paired-end reads). The argument to --minScoreFraction determines what fraction of the maximum score s a mapping must achieve to be potentially retained. For a minimum score fraction of f, only mappings with a score > f * s will be kept. Mappings with lower scores will be considered as low-quality, and will be discarded.

It is worth noting that mapping validation uses extension alignment. This means that the read need not map end-to-end. Instead, the score of the mapping will be the position along the alignment with the highest score. This is the score which must reach the fraction threshold for the read to be considered as valid.

Single-cell protocol specific notes

In cases where single-cell protocol supports variable length cellbarcodes, alevin adds nucleotide padding to make the lengths uniform. Furthermore, the padding scheme ensures that there are no collisions added in the process. The padding scheme is as follows:

1. sci-RNA-seq3: The barcode is composed of 9-10 bp hairpin adaptor and 10 bp reverse transcription index making it 19-20 bp long. If the bacode is 20 bp long, alevin adds A and it adds AC if it is 19 bp long. Thus, the length of barcode in the output is 21 bp.

2. inDropV2: 8-11 bp barcode1 along with 8 bp barcode2 makes up the barcode. For barcode lengths of 16, 17, 18, and 19 bp, alevin adds AAAC, AAG, AT, and A respectively. Thus, the length of barcode in the output is 20 bp. Furthermore, the position of barcode1 is dependent on finding exact match of sequence w1. If exact match is not found, a search for w1 is performed allowing a maximum hamming distance 2 b/w w1 and read2 substring of w1 length within the required bounds; the first match is returned.

Output

Typical 10x experiment can range form hundreds to tens of thousand of cells – resulting in huge size of the count-matrices. Traditionally single-cell tools dumps the Cell-v-Gene count matrix in various formats. Although, this itself is an open area of research but by default alevin dumps a per-cell level gene-count matrix in a binary-compressed format with the row and column indexes in a separate file.

A typical run of alevin will generate 4 files:

• quants_mat.gz – Compressed count matrix.

• quants_mat_cols.txt – Column Header (Gene-ids) of the matrix.

• quants_mat_rows.txt – Row Index (CB-ids) of the matrix.

• quants_tier_mat.gz – Tier categorization of the matrix.

Note

Working with R packages

Alevin generates multiple metadata files like the hash codes of the reference transcriptome and it’s crucial for working with downstream R package like tximeta . Hence along with the above files, it’s advisable to keep the complete output folder generated by alevin.

Along with the Cell-v-Gene count matrix, alevin dumps a 3-fold categorization of each estimated count value of a gene(each cell disjointly) in the form of tiers. Tier 1 is the set of genes where all the reads are uniquely mapping. Tier 2 is genes that have ambiguously mapping reads, but connected to unique read evidence as well, that can be used by the EM to resolve the multimapping reads. Tier 3 is the genes that have no unique evidence and the read counts are, therefore, distributed between these genes according to an uninformative prior.

Alevin can also dump the count-matrix in a human readable – matrix-market-exchange (_mtx_) format, if given flag –dumpMtx which generates a new output file called quants_mat.mtx.

Output Quality Check

Alevin generated gene-count matrix can be visualized for various quality checks using alevinQC , a shiny based R package and it is actively supported by Charlotte Soneson.

Tutorial & Parsers

We have compiled a step-by-step resource to help get started with aleivn. We have tutorials on how to get input, run and generate output using alevin’s framework which can be found here at Alevin Tutorials. The tutorial also covers the topic of integrating alevin with downstream analysis tools like Seurat and Monocle. If you are interested in parsing various output binary formats like quants_mat.gz, quants_tier_mat.gz, cell_umigraph.gz etc. of alevin in python, checkout our companion repo for python parsing. This repo is also available on pip and can be installed through pip install vpolo. We cover how to use this library on our alevin-tutorial website too.

Alevin Logs

Alevin generates alevin_meta_info.json file with the following json entries. Please note based on the command line flags provided during the time alevin was run, some of the below json entries may not be present.

• total_reads – Total number of reads in the experiment as observed by alevin.

• reads_with_N – Total number of reads with at least one nucleotide N in their cellular barcode sequence (and are not used for quantification).

• noisy_cb_reads – Total number of reads from noisy cellular barcodes (and are not used for quantification). A cellular barcode can be marked noisy based on many different conditions, for example all the barcodes below “knee” threshold or all the barcodes below provided threshold on –expectCells / –forceCells.

• noisy_umi_reads – Total number of reads with at least one nucleotide N in their UMI sequence (and are not used for quantification).

• used_reads – Total reads used for the quantification: total_reads - reads_with_N - noisy_cb_reads - noisy_umi_reads.

• mapping_rate – Fraction of reads mapping to the reference i.e. #mapped reads / total_reads.

• reads_in_eqclasses - Total number of reads present in the bfh (cell level equivalence classes).

• total_cbs – Total number of cellular barcodes observed by alevin in the experiment.

• used_cbs – Total number of cellular barcodes used by alevin for the quantification.

• initial_whitelist – Total number of whitelisted cellular barcodes by “knee” based thresholding.

• low_conf_cbs – Total number of low confidence cellular barcodes quantified for intelligent whitelisting.

• num_features – Total number of features used intelligent whitelisting of the cellular barcodes.

• final_num_cbs – Total number of cellular barcodes present in the output quant matrix.

• deduplicated_umis – Total number of UMIs present in the experiment post UMI deduplication across all cells.

• mean_umis_per_cell – Mean of the number of UMIs (post deduplication) present in each cell.

• mean_genes_per_cell – Mean of the number of genes expressed (>0 counts) in each cell.

• no_read_mapping_cbs – Total number of cellular barcodes with no reads mapped to them.

• num_bootstraps – Total number of bootstrap inferential replicates generated for each cell.

Misc

Finally, the purpose of making this software available is because we believe it may be useful for people dealing with single-cell RNA-seq data. We want the software to be as useful, robust, and accurate as possible. So, if you have any feedback — something useful to report, or just some interesting ideas or suggestions — please contact us (asrivastava@cs.stonybrook.edu and/or rob.patro@cs.stonybrook.edu). If you encounter any bugs, please file a detailed bug report at the Salmon GitHub repository.

BibTex

@article{srivastava2019alevin,

title={Alevin efficiently estimates accurate gene abundances from dscRNA-seq data},

author={Srivastava, Avi and Malik, Laraib and Smith, Tom and Sudbery, Ian and Patro, Rob},

journal={Genome biology},

volume={20},

number={1},

pages={65},

year={2019},

publisher={BioMed Central}

}

@article{Srivastava2020,

doi = {10.1093/bioinformatics/btaa450},

url = {https://doi.org/10.1093/bioinformatics/btaa450},

year = {2020},

month = jul,

publisher = {Oxford University Press ({OUP})},

volume = {36},

number = {Supplement{_}1},

pages = {i292–i299},

author = {Avi Srivastava and Laraib Malik and Hirak Sarkar and Rob Patro},

title = {A Bayesian framework for inter-cellular information sharing improves {dscRNA}-seq quantification},

journal = {Bioinformatics}

}

DOI

• https://doi.org/10.1186/s13059-019-1670-y

• https://doi.org/10.1093/bioinformatics/btaa450

References

1

Zhu, Anqi, et al. “Nonparametric expression analysis using inferential replicate counts.” BioRxiv (2019): 561084.

2

Qiu, Xiaojie, et al. “Reversed graph embedding resolves complex single-cell trajectories.” Nature methods 14.10 (2017): 979.

3

Butler, Andrew, et al. “Integrating single-cell transcriptomic data across different conditions, technologies, and species.” Nature biotechnology 36.5 (2018): 411.

4

Macosko, Evan Z., et al. “Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets.” Cell 161.5 (2015): 1202-1214.

5

Zheng, Grace XY, et al. “Massively parallel digital transcriptional profiling of single cells.” Nature communications 8 (2017): 14049.

6

Patro, Rob, et al. “Salmon provides fast and bias-aware quantification of transcript expression.” Nature Methods (2017). Advanced Online Publication. doi: 10.1038/nmeth.4197.

7

Petukhov, Viktor, et al. “Accurate estimation of molecular counts in droplet-based single-cell RNA-seq experiments.” bioRxiv (2017): 171496.

8

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/algorithms/overview

Salmon Output File Formats

Quantification File

Salmon’s main output is its quantification file. This file is a plain-text, tab-separated file with a single header line (which names all of the columns). This file is named quant.sf and appears at the top-level of Salmon’s output directory. The columns appear in the following order:

Name

Length

EffectiveLength

TPM

NumReads

Each subsequent row describes a single quantification record. The columns have the following interpretation.

• Name — This is the name of the target transcript provided in the input transcript database (FASTA file).

• Length — This is the length of the target transcript in nucleotides.

• EffectiveLength — This is the computed effective length of the target transcript. It takes into account all factors being modeled that will effect the probability of sampling fragments from this transcript, including the fragment length distribution and sequence-specific and gc-fragment bias (if they are being modeled).

• TPM — This is salmon’s estimate of the relative abundance of this transcript in units of Transcripts Per Million (TPM). TPM is the recommended relative abundance measure to use for downstream analysis.

• NumReads — This is salmon’s estimate of the number of reads mapping to each transcript that was quantified. It is an “estimate” insofar as it is the expected number of reads that have originated from each transcript given the structure of the uniquely mapping and multi-mapping reads and the relative abundance estimates for each transcript.

Command Information File

In the top-level quantification directory, there will be a file called cmd_info.json. This is a JSON format file that records the main command line parameters with which Salmon was invoked for the run that produced the output in this directory.

Auxiliary Files

The top-level quantification directory will contain an auxiliary directory called aux_info (unless the auxiliary directory name was overridden via the command line). This directory will have a number of files (and subfolders) depending on how salmon was invoked.

Meta information

The auxiliary directory will contain a JSON format file called meta_info.json which contains meta information about the run, including stats such as the number of observed and mapped fragments, details of the bias modeling etc. If Salmon was run with automatic inference of the library type (i.e. --libType A), then one particularly important piece of information contained in this file is the inferred library type. Most of the information recorded in this file should be self-descriptive.

Unique and ambiguous count file

The auxiliary directory also contains 2-column tab-separated file called ambig_info.tsv. This file contains information about the number of uniquely-mapping reads as well as the total number of ambiguously-mapping reads for each transcript. This file is provided mostly for exploratory analysis of the results; it gives some idea of the fraction of each transcript’s estimated abundance that derives from ambiguously-mappable reads.

Observed library format counts

When run in mapping-based mode, the quantification directory will contain a file called lib_format_counts.json. This JSON file reports the number of fragments that had at least one mapping compatible with the designated library format, as well as the number that didn’t. It also records the strand-bias that provides some information about how strand-specific the computed mappings were.

Finally, this file contains a count of the number of mappings that were computed that matched each possible library type. These are counts of mappings, and so a single fragment that maps to the transcriptome in more than one way may contribute to multiple library type counts. Note: This file is currently not generated when Salmon is run in alignment-based mode.

Fragment length distribution

The auxiliary directory will contain a file called fld.gz. This file contains an approximation of the observed fragment length distribution. It is a gzipped, binary file containing integer counts. The number of (signed, 32-bit) integers (with machine-native endianness) is equal to the number of bins in the fragment length distribution (1,001 by default — for fragments ranging in length from 0 to 1,000 nucleotides).

Sequence-specific bias files

If sequence-specific bias modeling was enabled, there will be 4 files in the auxiliary directory named obs5_seq.gz, obs3_seq.gz, exp5_seq.gz, exp5_seq.gz. These encode the parameters of the VLMM that were learned for the 5’ and 3’ fragment ends. Each file is a gzipped, binary file with the same format.

It begins with 3 32-bit signed integers which record the length of the context (window around the read start / end) that is modeled, follwed by the length of the context that is to the left of the read and the length of the context that is to the right of the read.

Next, the file contains 3 arrays of 32-bit signed integers (each of which have a length of equal to the context length recorded above). The first records the order of the VLMM used at each position, the second records the shifts and the widths required to extract each sub-context — these are implementation details.

Next, the file contains a matrix that encodes all VLMM probabilities. This starts with two signed integers of type std::ptrdiff_t. This is a platform-specific type, but on most 64-bit systems should correspond to a 64-bit signed integer. These numbers denote the number of rows (nrow) and columns (ncol) in the array to follow.

Next, the file contains an array of (nrow * ncol) doubles which represent a dense matrix encoding the probabilities of the VLMM. Each row corresponds to a possible preceeding sub-context, and each column corresponds to a position in the sequence context. Unused values (values where the length of the sub-context exceed the order of the model at that position) contain a 0. This array can be re-shaped into a matrix of the appropriate size.

Finally, the file contains the marginalized 0:sup:th-order probabilities (i.e. the probability of each nucleotide at each position in the context). This is stored as a 4-by-context length matrix. As before, this entry begins with two signed integers that give the number of rows and columns, followed by an array of doubles giving the marginal probabilities. The rows are in lexicographic order.

Fragment-GC bias files

If Salmon was run with fragment-GC bias correction enabled, the auxiliary directory will contain two files named expected_gc.gz and observed_gc.gz. These are gzipped binary files containing, respectively, the expected and observed fragment-GC content curves. These files both have the same form. They consist of a 32-bit signed int, dtype which specifies if the values to follow are in logarithmic space or not. Then, the file contains two signed integers of type std::ptrdiff which give the number of rows and columns of the matrix to follow. Finally, there is an array of nrow by ncol doubles. Each row corresponds to a conditional fragment GC distribution, and the number of columns is the number of bins in the learned (or expected) fragment-GC distribution.

Equivalence class file

If salmon was run with the --dumpEq option, then a file called eq_classes.txt will exist in the auxiliary directory. The format of that file is as follows:

N (num transcripts)
M (num equiv classes)
tn_1
tn_2
...
tn_N
eq_1_size t_11 t_12 ... count
eq_2_size t_21 t_22 ... count

That is, the file begins with a line that contains the number of transcripts (say N) then a line that contains the number of equivalence classes (say M). It is then followed by N lines that list the transcript names — the order here is important, because the labels of the equivalence classes are given in terms of the ID’s of the transcripts. The rank of a transcript in this list is the ID with which it will be labeled when it appears in the label of an equivalence class. Finally, the file contains M lines, each of which describes an equivalence class of fragments. The first entry in this line is the number of transcripts in the label of this equivalence class (the number of different transcripts to which fragments in this class map — call this k). The line then contains the k transcript IDs. Finally, the line contains the count of fragments in this equivalence class (how many fragments mapped to these transcripts). The values in each such line are tab separated. Note: The indices for transcripts referenced in this file start at 0.

If salmon was run with the --dumpEqWeights or -d option, then the eq_classes.txt file will include a textual representation of the range-factorized equivalence classes will exist in the auxiliary directory. The format of that file is specified as follows:

N (num transcripts)
M (num equiv classes)
tn_1
tn_2
...
tn_N
eq_1_size t_11 t_12 ... p_11 p_12 ... count
eq_2_size t_21 t_22 ... p_21 p_22 ... count

That is, the file begins with a line that contains the number of transcripts (say N) then a line that contains the number of equivalence classes (say M). It is then followed by N lines that list the transcript names — the order here is important, because the labels of the equivalence classes are given in terms of the ID’s of the transcripts. The rank of a transcript in this list is the ID with which it will be labeled when it appears in the label of an equivalence class. Finally, the file contains M lines, each of which describes a range-factorized equivalence class of fragments. The first entry in this line is the number of transcripts in the label of this equivalence class (the number of different transcripts to which fragments in this class map — call this k). The line then contains the k transcript IDs that partially define the label of this range-factorized equivalence class followed by k floating point values which correspond to the conditional probabilities of drawing a fragment from each of these k transcripts within this range-factorized equivalence class. Finally, the line contains the count of fragments in this equivalence class (how many fragments mapped to these transcripts with approximately this conditional probability distribution). The values in each such line are tab separated. Note: The indices for transcripts referenced in this file start at 0. Note: Unlike the simple equivalence classes, the same transcript set can appear more than once in the set of range-factorized equivalence classes. This is because different sets of fragments can induce quite different conditional probability distributions among these transcripts. For more details on this representation, please check the paper describing range-factorized equivalence classes.

Fragment Library Types

There are numerous library preparation protocols for RNA-seq that result in sequencing reads with different characteristics. For example, reads can be single end (only one side of a fragment is recorded as a read) or paired-end (reads are generated from both ends of a fragment). Further, the sequencing reads themselves may be unstranded or strand-specific. Finally, paired-end protocols will have a specified relative orientation. To characterize the various different typs of sequencing libraries, we’ve created a miniature “language” that allows for the succinct description of the many different types of possible fragment libraries. For paired-end reads, the possible orientations, along with a graphical description of what they mean, are illustrated below:

The library type string consists of three parts: the relative orientation of the reads, the strandedness of the library, and the directionality of the reads.

The first part of the library string (relative orientation) is only provided if the library is paired-end. The possible options are:

I = inward
O = outward
M = matching

The second part of the read library string specifies whether the protocol is stranded or unstranded; the options are:

S = stranded
U = unstranded

If the protocol is unstranded, then we’re done. The final part of the library string specifies the strand from which the read originates in a strand-specific protocol — it is only provided if the library is stranded (i.e. if the library format string is of the form S). The possible values are:

F = read 1 (or single-end read) comes from the forward strand
R = read 1 (or single-end read) comes from the reverse strand

So, for example, if you wanted to specify a fragment library of strand-specific paired-end reads, oriented toward each other, where read 1 comes from the forward strand and read 2 comes from the reverse strand, you would specify -l ISF on the command line. This designates that the library being processed has the type “ISF” meaning, Inward (the relative orientation), Stranded (the protocol is strand-specific), Forward (read 1 comes from the forward strand).

The single end library strings are a bit simpler than their pair-end counter parts, since there is no relative orientation of which to speak. Thus, the only possible library format types for single-end reads are U (for unstranded), SF (for strand-specific reads coming from the forward strand) and SR (for strand-specific reads coming from the reverse strand).

A few more examples of some library format strings and their interpretations are:

IU (an unstranded paired-end library where the reads face each other)

SF (a stranded single-end protocol where the reads come from the forward strand)

OSR (a stranded paired-end protocol where the reads face away from each other,
     read1 comes from reverse strand and read2 comes from the forward strand)

Note

Correspondence to TopHat library types

The popular TopHat RNA-seq read aligner has a different convention for specifying the format of the library. Below is a table that provides the corresponding sailfish/salmon library format string for each of the potential TopHat library types:

TopHat	Salmon (and Sailfish)
	Paired-end	Single-end
-fr-unstranded	-l IU	-l U
-fr-firststrand	-l ISR	-l SR
-fr-secondstrand	-l ISF	-l SF

The remaining salmon library format strings are not directly expressible in terms of the TopHat library types, and so there is no direct mapping for them.

Indices and tables

• Index

• Module Index

• Search Page