EMA uses a latent variable model to align barcoded short-reads (such as those produced by 10x Genomics' sequencing platform). More information is available in our paper. The full experimental setup is available here.
brew install brewsci/bio/ema
conda install -c bioconda ema
git clone --recursive https://github.com/arshajii/ema
cd ema
make
(The --recursive
flag is needed because EMA uses BWA's C API.)
usage: ema <count|preproc|align|help> [options]
count: perform preliminary barcode count (takes interleaved FASTQ via stdin)
-w <whitelist path>: specify barcode whitelist [required]
-o <output prefix>: specify output prefix [required]
-p: using haplotag barcodes
preproc: preprocess barcoded FASTQ files (takes interleaved FASTQ via stdin)
-w <whitelist path>: specify whitelist [required]
-n <num buckets>: number of barcode buckets to make [500]
-h: apply Hamming-2 correction [off]
-o: <output directory> specify output directory [required]
-b: output BX:Z-formatted FASTQs [off]
-p: using haplotag barcodes
-t <threads>: set number of threads [1]
all other arguments: list of all output prefixes generated by count stage
align: choose best alignments based on barcodes
-1 <FASTQ1 path>: first (preprocessed and sorted) FASTQ file [none]
-2 <FASTQ2 path>: second (preprocessed and sorted) FASTQ file [none]
-s <EMA-FASTQ path>: specify special FASTQ path [none]
-x: multi-input mode; takes input files after flags and spawns a thread for each [off]
-r <FASTA path>: indexed reference [required]
-o <SAM file>: output SAM file [stdout]
-R <RG string>: full read group string (e.g. '@RG\tID:foo\tSM:bar') [none]
-d: apply fragment read density optimization [off]
-p <platform>: sequencing platform (one of '10x', 'tru', 'cpt', 'haplotag') [10x]
-i <index>: index to follow 'BX' tag in SAM output [1]
-t <threads>: set number of threads [1]
all other arguments (only for -x): list of all preprocessed inputs
help: print this help message
If using haplotagging data, where barcodes are coded in the read headers as BX:Z:AxxCxxBxxDxx
, you do not need to provide
a barcode whitelist for the count
or proproc
steps.
EMA has several input modes:
-s <input>
: Input file is a single preprocessed "special" FASTQ generated by the preprocessing steps below.-x
: Input files are listed after flags (as inema align -a -b -c <input 1> <input 2> ... <input N>
). Each of these inputs are processed and all results are written to the SAM file specified with-o
.-1 <first mate>
/-2 <second mate>
: Input files are standard FASTQs. For interleaved FASTQs,-2
can be omitted. The only restrictions in this input mode are that read identifiers must end in:<barcode sequence>
and that the FASTQs must be sorted by barcode. For 10x data, the above two modes are preferred.
Multithreading can be enabled with -t <num threads>
. The actual threading mode is dependent on how the input is being read, however:
-s
,-1
/-2
: Multiple threads are spawned to work on the single input file (or pair of input files).-x
: Threads work on the input files individually.
(Note that, because of this, it never makes sense to spawn more threads than there are input files when using -x
.)
In this guide, we use the following additional tools:
We also use a 10x barcode whitelist, which can be found here.
Preprocessing 10x data entails several steps, the first of which is counting barcodes (-j
specifies the number of jobs to be spawned by parallel
):
cd /path/to/gzipped_fastqs/
parallel -j40 --bar 'pigz -c -d {} | \
ema count -w /path/to/whitelist.txt -o {/.} 2>{/.}.log' ::: *RA*.gz
Make sure that the FASTQs are interleaved and only contain the actual reads in the files above (as opposed to sample indices, typically with I1
in their filenames rather than RA
). This will produce *.ema-ncnt
and *.ema-fcnt
files, containing the count data.
If you do not have interleaved files, you can interleave them as follows:
parallel -j40 --bar 'paste <(pigz -c -d {} | paste - - - -) <(pigz -c -d {= s:_R1_:_R2_: =} | paste - - - -) | tr "\t" "\n" |\
ema count -w /path/to/whitelist.txt -o {/.} 2>{/.}.log' ::: *_R1_*.gz
where s:_R1_:_R2_:
is the regex that casts first-end filenames into the second-end filenames (make sure to adjust this if your naming scheme is different).
Now we can do the actual preprocessing, which splits the input into barcode bins (500 by default; specified with -n
). This preprocessing can be parallelized via -t
, which specifies how many threads to use:
pigz -c -d *RA*.gz | ema preproc -w /path/to/whitelist.txt -n 500 -t 40 -o output_dir *.ema-ncnt 2>&1 | tee preproc.log
or if you do not have interleaved files:
paste <(pigz -c -d *_R1_*.gz | paste - - - -) <(pigz -c -d *_R2_*.gz | paste - - - -) | tr "\t" "\n" |\
ema preproc -w /path/to/whitelist.txt -n 500 -t 40 -o output_dir *.ema-ncnt 2>&1 | tee preproc.log
First we map each barcode bin with EMA. Here, we'll do this using a combination of GNU Parallel and EMA's internal multithreading, which we found to be optimal due to the runtime/memory trade-off. In the following, for instance, we use 10 jobs each with 4 threads (for 40 total threads). We also pipe EMA's SAM output (stdout by default) to samtools sort
, which produces a sorted BAM:
parallel --bar -j10 "ema align -t 4 -d -r /path/to/ref.fa -s {} |\
samtools sort -@ 4 -O bam -l 0 -m 4G -o {}.bam -" ::: output_dir/ema-bin-???
Lastly, we map the no-barcode bin with BWA:
bwa mem -p -t 40 -M -R "@RG\tID:rg1\tSM:sample1" /path/to/ref.fa output_dir/ema-nobc |\
samtools sort -@ 4 -O bam -l 0 -m 4G -o output_dir/ema-nobc.bam
Note that @RG\tID:rg1\tSM:sample1
is EMA's default read group. If you specify another for EMA, be sure to specify the same for BWA as well (both tools take the full read group string via -R
).
EMA performs duplicate marking automatically. We mark duplicates on BWA's output with sambamba markdup
:
sambamba markdup -t 40 -p -l 0 output_dir/ema-nobc.bam output_dir/ema-nobc-dupsmarked.bam
rm output_dir/ema-nobc.bam
Now we merge all BAMs into a single BAM (might require modifying ulimit
s, as in ulimit -n 10000
):
sambamba merge -t 40 -p ema_final.bam output_dir/*.bam
Now you should have a single, sorted, duplicate-marked BAM ema_final.bam
.
Instructions for preprocessing and running EMA on data from other sequencing platforms can be found here.
EMA outputs a standard SAM file with several additional tags:
XG
: alignment probabilityMI
: cloud identifier (compatible with Long Ranger)XA
: alternate high-probability alignments