1
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ Why not store the read identifiers and metadata in a separate table file? $\endgroup$ – Chris_Rands Nov 15 '18 at 9:25
  • 1
    $\begingroup$ Of course, this is possible. But apart from not being an answer to my question, this is not a good solution because it is not standardized and invites inconsistencies. You have to write special code to read the table. The SAM format (which is also a table) is maybe not the ideal but it is the best we have. If we could all stick to the SAM format instead of introducing custom formats the interoperability of bioinformatics tools would be much better. $\endgroup$ – buenaonda Nov 16 '18 at 6:51
  • $\begingroup$ I’m a bit confused because the preprocessing you mention is pretty routine, and is generally done without the involvement of BAM files. $\endgroup$ – Konrad Rudolph Nov 20 '18 at 12:23
  • $\begingroup$ Please remember to confirm an answer once you've received one. $\endgroup$ – Stephopolis Dec 11 '18 at 18:40
4
$\begingroup$

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).

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.