Is there a read mapper that allows input in SAM format instead of FASTA/FASTQ?

Question

When we get raw sequenced reads in FASTQ format, we usually need to do some pre-processing (QC, demultiplexing, etc.) prior to sequence alignment. The only way to register the results of this preprocessing is to add metadata to each read. This is best done by using the SAM format which has a well-defined system of metadata tags.

However, it seems that not a single read mapper (except maybe for bwa aln with option -b) accepts SAM (or one of the binary equivalents BAM/CRAM) as input to make use of the metadata that has been added during preprocessing. This implies that one has to convert the unaligned reads back to FASTQ format.

In any case, whether one uses FASTQ throughout or converts back to FASTQ, a common annoyance with FASTQ is that metadata is registered as part of the query name - this is not standardized or even well-defined. After alignment, we again have to do some text processing to get the metadata out of the query name and into proper metadata tags defined by SAM. This back and forth hacking is highly unfortunate and error prone. I have seen this, for example, in some single-cell sequencing projects where we have highly multiplexed sequencing runs with reads that need to be QC'd and classified. The new SAM standard from May 2018 does indeed define proper tags for handling cell-multiplexed runs. SAM is the de-facto standard to store both unaligned and aligned reads for years to come so it should be standard for read mappers to accept it as input.

I would be grateful for some hints and if we could collect here some read mappers that do accept SAM/BAM/CRAM as input format.

Answer 1

The common solution for scRNA-seq is to put cell barcodes and such in read headers and then post-process things with UMItools.

But regarding your actual question, STAR can accept SAM/BAM as input with the --readFilesType SAM PE option (it's SAM for both that format and BAM). You can swap SE for PE if your data is single-end (or is effectively that way, as is the case in many scRNA-seq experiments). Note that in the case of PE experiments that mates must be next to each other. I really don't know what this does with auxiliary tags though, so you would need to see if they're lost before you start processing a lot of data).

Answer 2

With bwa-mem or minimap2, the recommended way is

samtools fastq -T BC,RX name-grouped.sam | bwa mem -C -p -

This passes the BC and RX tags through bwa-mem and copies them to the SAM output. If you are running paired-end alignment and your BAM is coordinate sorted, you have to run samtools collate to group reads by name:

samtools collate -Ol0 in.bam tmp | samtools fastq -T BC,RX - | bwa mem -pC -

Here bwa mem -C or minimap2 -y copies FASTQ comments to SAM output. This strategy is much easier to deal with than SAM/BAM. A caveat is that it doesn't work well for read groups because each read group requires a separate header line. You have to manually modify the header.

I am not sure which short-read RNA-seq mappers have implemented a similar feature. You might consider to send a feature request.

Answer 3

At the Broad Institute, we do something similar to what is described by user172818. We use uBAM (unmapped BAM) to store unmapped reads, because as you say the SAM format has a more appropriate way to encode metadata. However, we use a Picard tool called SamToFastq to do the initial conversion then pipe the output of that to BWA mem. Then we use an additional step to reattach metadata, read group info and so on, using a Picard tool called MergeBamAlignments.

You can see the full implementation here (starts as line 71):

https://github.com/gatk-workflows/gatk4-genome-processing-pipeline/blob/master/tasks/Alignment.wdl

Answer 4

SNAP is a hash-index based genomic aligner that can take SAM and BAM as input. It has been used as an aligner for RNA-seq (see SNAPR) and BS-seq (see BS-Seeker3). SNAPR seems to accept SAM/BAM but it is not clear whether BS-Seeker3 wraps this functionality as well.

I tried SNAP and it seems to do the job well, i.e. it only adds tags during alignment and does not remove metadata that has been deposited prior to alignment, neither from the SAM header nor from individual reads.

NOTE I post this answer to my own question after I became aware of SNAP.