SAM format output?

[wwood]

Hi,

Great tool. Curious if would be striaghtforward to implement SAM format output from LexicMap? Then it could be slotted in as another mapper in CoverM.

Thanks, ben

[shenwei356]

It should be easy. I can add a new command lexicmap utils 2sam, just like the 2blast.

[shenwei356]

Hi @wwood Ben, since LexicMap is an aligner, not a mapper. So not all information can be provided in a SAM format. Current format is here.

I don't know what SAM fields you need. Below is the initial plan. Please give some suggestions.

header:

Tag	Tag2	Value	Comment
@hd	VN	1.5
	SO	unsorted
@sq	SN		There would be a massive number of records
	LN

data:

Col	Field	data in LexicMap	How
1	QNAME	query	query
2	FLAG	N/A	16 for reverse complement
3	RNAME	sseqid	sseqid
4	POS	sstart, but for aligned region	for strand +: sstart - (qstart - 1)
5	MAPQ	N/A	255?
6	CIGAR	cigar, but for aligned region	add "H" records in 5' and 3' end
7	RNEXT	N/A	*
8	PNEXT	N/A	*
9	TLEN	alenHSP	ssend - sstart
10	SEQ	qseq, but for aligned region	unaligned region is not available, can be * if it's not important.
11	QUAL	N/A	*

[wwood]

Thanks @shenwei356 I think that mostly does it. Couple comments, keeping in mind I could be forgetting some parts of what is needed.

1. The next coverm version has mapq based filtering, but we can just throw an error when the user asks for that and lexicmap.

2. I do need the NM tag, because that is used for pident filtering. I think that should be doable, just mismatches plus insert bases plus deleted bases.

3. I think it might be possible to only include sequences in the ref in the header, but would need to check this. Would help since then header would be way smaller. Alternatively could think about cram.

4. It probably makes sense for a user to have the option of prespecifying that the database is already formatted. But if only sseqid is given then the genome would not be known by the Sam format? Sseqid is also assumed to be unique by Sam tools and co. Maybe encode as genome~sseqid as the reference name?

5. Is the idea to stream through the util that does conversion? Would be ideal.

6. Sam format requires 1 alignment of a query sequence to be primary, so some choice will have to made even if the alignments are the same to 2 refs.

Happy to test out something concrete. Imagine I'm missing something.

Hth. ben