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

[shenwei356]

Hi @wwood , I'm very sorry for the late reply.

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

MAPQ can be computed as below (with the help of ChatGPT).

# bitscore can be replaced with alignment score with a simple calculation
S1 = best bitscore for query. 
S2 = second best bitscore for query (0 if none)

MAPQ = 40 * (S1 - S2) / S1
MAPQ *= (qcovHSP / 100)                       # cov_factor
MAPQ *= (pident / 100) * (1 - gaps / alenHSP) # qual_factor

MAPQ = int(min(60, max(0, MAPQ)))

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.

Edit distance = alenHSP * (1 - pident/100)

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.

Do you mean only include IDs of matched sequences, that would be smaller and easier.
But the alignment result file need to be read twice, one for these ids, and one for full records.

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?

Yes, genome id will not appear in SAM.

Sseqid is also assumed to be unique by Sam tools and co. Maybe encode as genome~sseqid as the reference name?

Can append the genome id to the sequence id as the prefix.

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

It will be lexicmap utils 2sam. If the header section is not required, it can accept STDIN.
If not, it has to read the alignment file twice.

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.

Can mark the one with the highest bitscore.
Well, if the alignments are the same to 2 refs, choose a random one.

Hope I can write the code this weekend.

[shenwei356]

@wwood , sorry for the late reply.
I've implemented an early version and tested it with a few examples.

(contents below are edited)

$ lexicmap utils 2sam -h
Convert the default search output to SAM format

Input:
   - Output file of 'lexicmap search' with the flag -a/--all.
   - Do not support STDIN.

Output:
   - Clipped regions in SEQ are represented as N's,
     as the input contains only the aligned portion of the sequences.
   - Different from SAM files produced by Minimap2,
     'X' (mismatch) in CIGAR is not converted to 'M' (match).
   - NM (Edit distance ) and AS (alignment score) fields are produced.

Usage:
  lexicmap utils 2sam [flags] 

Flags:
  -b, --buffer-size string      ► Size of buffer, supported unit: K, M, G. You need increase the value
                                when "bufio.Scanner: token too long" error reported (default "20M")
  -c, --concat-sgenome-sseqid   ► Concatenate sgenome and sseqid to make sure the reference sequence
                                names are distinct.
  -h, --help                    help for 2sam
  -o, --out-file string         ► Out file, supports and recommends a ".gz" suffix ("-" for stdout).
                                (default "-")
  -s, --separater string        ► Separater between sgenome and sseqid (default "~")

An example of searching a 16S rRNA gene (NC_000913.3:4166659-4168200, t.fasta.txt) in a bacterium (GCF_003697165.2, GCF_003697165.2.fa.gz). Note that: This genome has 7 copies of the gene on both strands, so the MAPQ is 0.

lexicmap

$ lexicmap utils 2sam t.lexicmap.tsv 
10:30:17.239 [INFO] round 1/2: extracting subject sequence IDs and lengths ...
10:30:17.240 [INFO]   elapsed time: 308.697µs
10:30:17.240 [INFO] round 2/2: converting to SAM format ...
10:30:17.240 [INFO]   elapsed time: 48.05µs
@HD	VN:1.6	SO:unsorted	GO:query
@SQ	SN:NZ_CP033092.2	LN:4903501
@PG	ID:lexicmap	PN:lexicmap	VN:0.9.0
NC_000913.3:4166659-4168200	0	NZ_CP033092.2	458559	0	79M1X8M1X120M1X1332M	*	0	1542	AAATTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAACGGTAACAGGAAGAAGCTTGCTTCTTTGCTGACGAGTGGCGGACGGGTGAGTAATGTCTGGGAAACTGCCTGATGGAGGGGGATAACTACTGGAAACGGTAGCTAATACCGCATAACGTCGCAAGACCAAAGAGGGGGACCTTCGGGCCTCTTGCCATCGGATGTGCCCAGATGGGATTAGCTAGTAGGTGGGGTAACGGCTCACCTAGGCGACGATCCCTAGCTGGTCTGAGAGGATGACCAGCCACACTGGAACTGAGACACGGTCCAGACTCCTACGGGAGGCAGCAGTGGGGAATATTGCACAATGGGCGCAAGCCTGATGCAGCCATGCCGCGTGTATGAAGAAGGCCTTCGGGTTGTAAAGTACTTTCAGCGGGGAGGAAGGGAGTAAAGTTAATACCTTTGCTCATTGACGTTACCCGCAGAAGAAGCACCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACGCCGTAAACGATGTCGACTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTAACGCGTTAAGTCGACCGCCTGGGGAGTACGGCCGCAAGGTTAAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGCAACGCGAAGAACCTTACCTGGTCTTGACATCCACGGAAGTTTTCAGAGATGAGAATGTGCCTTCGGGAACCGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGTGTTGTGAAATGTTGGGTTAAGTCCCGCAACGAGCGCAACCCTTATCCTTTGTTGCCAGCGGTCCGGCCGGGAACTCAAAGGAGACTGCCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACGACCAGGGCTACACACGTGCTACAATGGCGCATACAAAGAGAAGCGACCTCGCGAGAGCAAGCGGACCTCATAAAGTGCGTCGTAGTCCGGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTGGATCAGAATGCCACGGTGAATACGTTCCCGGGCCTTGTACACACCGCCCGTCACACCATGGGAGTGGGTTGCAAAAGAAGTAGGTAGCTTAACCTTCGGGAGGGCGCTTACCACTTTGTGATTCATGACTGGGGTGAAGTCGTAACAAGGTAACCGTAGGGGAACCTGCGGTTGGATCACCTCCTTA	*	NM:i:3	AS:i:3067
NC_000913.3:4166659-4168200	256	NZ_CP033092.2	1285123	0	79M1X8M1X120M1X1332M	*	0	1542	*	*	NM:i:3	AS:i:3067
NC_000913.3:4166659-4168200	272	NZ_CP033092.2	3780640	0	79M1X8M1X120M1X1332M	*	0	1542	*	*	NM:i:3	AS:i:3067
NC_000913.3:4166659-4168200	272	NZ_CP033092.2	4551515	0	79M1X8M1X120M1X1332M	*	0	1542	*	*	NM:i:3	AS:i:3067
NC_000913.3:4166659-4168200	272	NZ_CP033092.2	4591684	0	79M1X8M1X120M1X1332M	*	0	1542	*	*	NM:i:3	AS:i:3067
NC_000913.3:4166659-4168200	272	NZ_CP033092.2	4726193	0	75M1X3M1X8M1X1453M	*	0	1542	*	*	NM:i:3	AS:i:3067
NC_000913.3:4166659-4168200	272	NZ_CP033092.2	4844587	0	79M1X8M1X120M1X1332M	*	0	1542	*	*	NM:i:3	AS:i:3067

minimap2

@HD	VN:1.6	SO:unsorted	GO:query
@SQ	SN:NZ_CP033092.2	LN:4903501
@SQ	SN:NZ_CP033091.2	LN:131333
@PG	ID:minimap2	PN:minimap2	VN:2.29-r1283	CL:minimap2 -N 100 -a -o t.sam GCF_003697165.2.fa.gz t.16S.fa
NC_000913.3:4166659-4168200	16	NZ_CP033092.2	4844587	0	1542M	*	0	0	TAAGGAGGTGATCCAACCGCAGGTTCCCCTACGGTTACCTTGTTACGACTTCACCCCAGTCATGAATCACAAAGTGGTAAGCGCCCTCCCGAAGGTTAAGCTACCTACTTCTTTTGCAACCCACTCCCATGGTGTGACGGGCGGTGTGTACAAGGCCCGGGAACGTATTCACCGTGGCATTCTGATCCACGATTACTAGCGATTCCGACTTCATGGAGTCGAGTTGCAGACTCCAATCCGGACTACGACGCACTTTATGAGGTCCGCTTGCTCTCGCGAGGTCGCTTCTCTTTGTATGCGCCATTGTAGCACGTGTGTAGCCCTGGTCGTAAGGGCCATGATGACTTGACGTCATCCCCACCTTCCTCCAGTTTATCACTGGCAGTCTCCTTTGAGTTCCCGGCCGGACCGCTGGCAACAAAGGATAAGGGTTGCGCTCGTTGCGGGACTTAACCCAACATTTCACAACACGAGCTGACGACAGCCATGCAGCACCTGTCTCACGGTTCCCGAAGGCACATTCTCATCTCTGAAAACTTCCGTGGATGTCAAGACCAGGTAAGGTTCTTCGCGTTGCATCGAATTAAACCACATGCTCCACCGCTTGTGCGGGCCCCCGTCAATTCATTTGAGTTTTAACCTTGCGGCCGTACTCCCCAGGCGGTCGACTTAACGCGTTAGCTCCGGAAGCCACGCCTCAAGGGCACAACCTCCAAGTCGACATCGTTTACGGCGTGGACTACCAGGGTATCTAATCCTGTTTGCTCCCCACGCTTTCGCACCTGAGCGTCAGTCTTCGTCCAGGGGGCCGCCTTCGCCACCGGTATTCCTCCAGATCTCTACGCATTTCACCGCTACACCTGGAATTCTACCCCCCTCTACGAGACTCAAGCTTGCCAGTATCAGATGCAGTTCCCAGGTTGAGCCCGGGGATTTCACATCTGACTTAACAAACCGCCTGCGTGCGCTTTACGCCCAGTAATTCCGATTAACGCTTGCACCCTCCGTATTACCGCGGCTGCTGGCACGGAGTTAGCCGGTGCTTCTTCTGCGGGTAACGTCAATGAGCAAAGGTATTAACTTTACTCCCTTCCTCCCCGCTGAAAGTACTTTACAACCCGAAGGCCTTCTTCATACACGCGGCATGGCTGCATCAGGCTTGCGCCCATTGTGCAATATTCCCCACTGCTGCCTCCCGTAGGAGTCTGGACCGTGTCTCAGTTCCAGTGTGGCTGGTCATCCTCTCAGACCAGCTAGGGATCGTCGCCTAGGTGAGCCGTTACCCCACCTACTAGCTAATCCCATCTGGGCACATCCGATGGCAAGAGGCCCGAAGGTCCCCCTCTTTGGTCTTGCGACGTTATGCGGTATTAGCTACCGTTTCCAGTAGTTATCCCCCTCCATCAGGCAGTTTCCCAGACATTACTCACCCGTCCGCCACTCGTCAGCAAAGAAGCAAGCTTCTTCCTGTTACCGTTCGACTTGCATGTGTTAGGCCTGCCGCCAGCGTTCAATCTGAGCCATGATCAAACTCTTCAATTT	*	NM:i:3	ms:i:3066	AS:i:3066	nn:i:0	tp:A:P	cm:i:267	s1:i:1508	s2:i:1513	de:f:0.0019	rl:i:0
NC_000913.3:4166659-4168200	272	NZ_CP033092.2	3780640	0	1542M	*	0	0	*	*	NM:i:3	ms:i:3066	AS:i:3066	nn:i:0	tp:A:S	cm:i:267	s1:i:1508	de:f:0.0019	rl:i:0
NC_000913.3:4166659-4168200	256	NZ_CP033092.2	458559	0	1542M	*	0	0	*	*	NM:i:3	ms:i:3066	AS:i:3066	nn:i:0	tp:A:S	cm:i:267	s1:i:1508	de:f:0.0019	rl:i:0
NC_000913.3:4166659-4168200	272	NZ_CP033092.2	4591684	0	1542M	*	0	0	*	*	NM:i:3	ms:i:3066	AS:i:3066	nn:i:0	tp:A:S	cm:i:267	s1:i:1508	de:f:0.0019	rl:i:0
NC_000913.3:4166659-4168200	272	NZ_CP033092.2	4726193	0	1542M	*	0	0	*	*	NM:i:3	ms:i:3066	AS:i:3066	nn:i:0	tp:A:S	cm:i:269	s1:i:1513	de:f:0.0019	rl:i:0
NC_000913.3:4166659-4168200	256	NZ_CP033092.2	1285123	0	1542M	*	0	0	*	*	NM:i:3	ms:i:3066	AS:i:3066	nn:i:0	tp:A:S	cm:i:267	s1:i:1508	de:f:0.0019	rl:i:0
NC_000913.3:4166659-4168200	272	NZ_CP033092.2	4551515	0	1542M	*	0	0	*	*	NM:i:3	ms:i:3066	AS:i:3066	nn:i:0	tp:A:S	cm:i:267	s1:i:1508	de:f:0.0019	rl:i:0

Binaries are here: #10

Full usage and another example: https://github.com/shenwei356/LexicMap/blob/main/docs/content/usage/utils/2sam.md

[shenwei356]

Dear @wwood , is there any chance for you to try it? If there are no further issues, I would like to tag a new release.

[wwood]

Sorry for the slow follow-up here. I integrated lexicmap into coverm, but ran into this issue. I managed to reproduce it outside CoverM (but run from a git checkout to provide the test data). I used the most recent commit, not the binary you provided, fwiw.

$ lexicmap index -I tests/data/rhys_bug/genomes/ -O /tmp/tmp.lmi
17:28:44.198 [INFO] LexicMap v0.9.0
17:28:44.198 [INFO]   https://github.com/shenwei356/LexicMap
..

$ lexicmap search -d 'tests/data/rhys_bug/genomes.lmi/' -a -j 8  'tests/data/rhys_bug/20120700_S3D.head100000.1.fq.gz' 'tests/data/rhys_bug/20120700_S3D.head100000.2.fq.gz' -o /tmp/lexicmap.tsv
17:32:16.016 [INFO] LexicMap v0.9.0
17:32:16.016 [INFO]   https://github.com/shenwei356/LexicMap
17:32:16.016 [INFO]
17:32:16.016 [INFO] checking input files ...
17:32:16.016 [INFO]   2 input file(s) given
17:32:16.017 [INFO]
17:32:16.017 [INFO] loading index: tests/data/rhys_bug/genomes.lmi/
17:32:16.017 [INFO]   reading masks...
17:32:16.036 [INFO]   reading indexes of seeds (k-mer-value) data...
17:32:16.195 [INFO]   creating reader pools for 1 genome batches, each with 8 readers...
17:32:16.196 [INFO] index loaded in 180.577374ms
17:32:16.196 [INFO]
17:32:16.196 [INFO] searching with 8 threads...
17:32:16.196 [INFO]   maximum number of concurrent queries: 8, force garbage collection for every 64 queries
processed queries: 49664, speed: 76476.943 queries per minute
17:32:55.419 [INFO]
17:32:55.419 [INFO] processed queries: 50000, speed: 76486.438 queries per minute
17:32:55.419 [INFO] 34.9480% (17474/50000) queries matched
17:32:55.419 [INFO] done searching
17:32:55.419 [INFO] search results saved to: /tmp/lexicmap.tsv
17:32:55.420 [INFO]
17:32:55.420 [INFO] elapsed time: 39.404128899s
17:32:55.420 [INFO]

$ lexicmap utils 2sam /tmp/lexicmap.tsv -c | samtools sort -l0 -@ 0 >/tmp/lexicmap.sam
17:34:28.459 [INFO] round 1/2: extracting subject sequence IDs and lengths ...
17:34:28.469 [INFO]   elapsed time: 10.721355ms
17:34:28.469 [INFO] round 2/2: converting to SAM format ...
[E::sam_parse1] CIGAR and query sequence are of different length
[W::sam_read1_sam] Parse error at line 1376
samtools sort: truncated file. Aborting

Am I using it incorrectly? Or something amiss in 2sam?
Thanks.

[shenwei356]

Could you please share the index or genomes, and some reads that could reproduce the error? I can't reproduce this with my current data.

Previously, I have had similar errors. Because only the aligned parts of the query and target sequences are saved in the output of lexicmap search. So, for clipping in 5' or 3' end , I use N's to replace the missing parts.

And please make sure you're using the latest commit.

[wwood]

All of the data is available in the coverm repo, I believe, just checkout the main and then run the above from there. I'm afk but I did a checkout of the code about 24 hours ago. If things have changed I can try rerunning? -\------------- Ben Woodcroft Group leader, Centre for Microbiome Research, QUT

…

________________________________ From: Wei Shen ***@***.***> Sent: Saturday, March 7, 2026 8:56:01 PM To: shenwei356/LexicMap ***@***.***> Cc: Ben J Woodcroft ***@***.***>; Mention ***@***.***> Subject: Re: [shenwei356/LexicMap] SAM format output? (Issue #26) [https://avatars.githubusercontent.com/u/2655946?s=20&v=4]shenwei356 left a comment (shenwei356/LexicMap#26)<#26 (comment)> Could you please share the index or genomes, and some reads that could reproduce the error? I can't reproduce this with my current data. And please make sure you're using the latest commit. ― Reply to this email directly, view it on GitHub<#26 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAADX5BUAWS52YYKDYMDLC34PP54DAVCNFSM6AAAAACRE24SYGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHM2DAMJWGI2DSMBZGQ>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[shenwei356]

There is a bug, and it's fixed now. It works with data above. binaries are also updated: #10

[wwood]

Hi @shenwei356 sorry, I think I might have come across something futher. If you create the bam file with lexicmap 2sam, on tid 1221 there are 4 reads mapped. i.e.

(coverm) (conda)cpu1n001:20260311:~/git/coverm$ samtools view /tmp/bam_cache/tmp.lmi.20120700_S3D.head100000.1.fq.gz.bam |grep 73.20120700_S3D.34~73.20120700_S3D.34_contig_2590
HISEQ:191:C5DWMANXX:5:1101:10623:8541   0       73.20120700_S3D.34~73.20120700_S3D.34_contig_2590       4138    60      126M    *       0       126     CGGCTTCAGCGAGGCCAGCCCAGAGCAGCGAGCAGCAGCCCTGCCCTTCGTCCACACTATGTTCGATTGGGTCAAGGCCGAGATGACATGGATGTTCGCAGCGATCACCTGGTCAACGGTGGCTGA      *       NM:i:0  AS:i:251
HISEQ:191:C5DWMANXX:5:1101:10623:8541   16      73.20120700_S3D.34~73.20120700_S3D.34_contig_2590       4471    60      126M    *       0       126     GAGCTTCTGCACGGCGTGCTACACGGGGGTGTATCCAACGCAGTGGGTGGATGTGGAGGAGATACTGCCGGCGGGTGCGCTGGTCGGATAGCGCGCAGAGCGACCGGCCTGAGCAAGATCGAGATC      *       NM:i:0  AS:i:251
HISEQ:191:C5DWMANXX:5:1101:18319:9129   0       73.20120700_S3D.34~73.20120700_S3D.34_contig_2590       6410    60      115M1D11M       *       0       127TGAAACATCCTTCGCCGAGTGGATCAGACCGCGCTCGGCCAAGGCAATGAGCGCTCTGTTCAGTGCCGCCTCAGCCTCCAGATTGATCCACGGCGGCGTACCCCAGAGCGAGCTACAGTTCGTCCT   *       NM:i:0  AS:i:243
HISEQ:191:C5DWMANXX:5:1101:18319:9129   16      73.20120700_S3D.34~73.20120700_S3D.34_contig_2590       6731    60      126M    *       0       126     CGGCGACGAATTGTTTGAACTTTGGCAATCCGGAGAAGCCGGAGATCATGGCGCAGTTTTCGGCGACGATCGATGGGATCGCGGAGGCGTGTACGGCGCTGGGGACGCCGATTACGGGAGGCAATG      *       NM:i:0  AS:i:251

On the 3rd read, the NM is set to 0, but there is a D in the cigar. AFAICS, this seems to be a bug in 2sam?

[shenwei356]

Sorry, I shouldn't have calculated it like this:

Edit distance = alenHSP * (1 - pident/100)

Now, it's counted from the CIGAR string.

HISEQ:191:C5DWMANXX:5:1101:18319:9129   16      73.20120700_S3D.34~73.20120700_S3D.34_contig_2590       6731    60      126M    *       0       126     CGGCGACGAATTGTTTGAACTTTGGCAATCCGGAGAAGCCGGAGATCATGGCGCAGTTTTCGGCGACGATCGATGGGATCGCGGAGGCGTGTACGGCGCTGGGGACGCCGATTACGGGAGGCAATG     *       NM:i:0  AS:i:251
HISEQ:191:C5DWMANXX:5:1101:18319:9129   0       73.20120700_S3D.34~73.20120700_S3D.34_contig_2590       6410    60      115M1D11M       *       0       127     TGAAACATCCTTCGCCGAGTGGATCAGACCGCGCTCGGCCAAGGCAATGAGCGCTCTGTTCAGTGCCGCCTCAGCCTCCAGATTGATCCACGGCGGCGTACCCCAGAGCGAGCTACAGTTCGTCCT     *       NM:i:1  AS:i:243