Feature Request: output color set id for each positive kmer in query

[tbenavi1]

Hello,

Looking at other github issues I see that users have requested positional information in the reference where the k-mer maps. However, I am asking about something different. I would like an output file that gives information about each position in the query.

For every position in the query, I would like to know its corresponding ColorSet. And by "position in the query" I really mean the starting position of every k-mer in the query. A bed file format would be amazing:

seq1 0 5 0
seq1 5 7 1
seq1 7 10 2

This would indicate that the first 5 kmers of seq1 correspond to ColorSet 0, then the next two k-mers correspond to ColorSet 1 and so forth. Basically, I am envisioning an output format for just the lookup command of Fulgor. Does this make sense? I can give more explanation if necessary.

[jermp]

Hello,
yes, this makes sense and is also already implemented as it is a preliminary step for pseudoalignment.
However, this functionality is not explicitly exposed by the public API of Fulgor.

In other words, you want: for each kmer in the query that is found in the index, report its color set id. Right?

I'll implement this when I'm back from holidays; so this will happen in September.
Thank you for your interest in Fulgor! Do not hesitate to tell us more about your use case.

Best,
-Giulio

[jermp]

Hi @tbenavi1, in the meanwhile, could you please confirm if my understanding is correct?

In other words, you want: for each kmer in the query that is found in the index, report its color set id. Right?

Is this what you'd like?

[tbenavi1]

Hi, if possible I would like it in bed file format in query coordinates. It doesn't make sense to have the output contain the actual k-mers themselves, since this can be found in the query. And bed file format is nice because we can group consecutive k-mers in the query with the same color set id as a single line in the file. Does this make sense?

[tbenavi1]

Do color set ids start at 0 or 1? For a k-mer in the query that is not in the index, it would be nice if it was assigned a color set id of 0 (meaning not in the index) but I'm not sure how it is currently implemented.

[jermp]

There are many possible ways of compressing the output.
A sequence of triples, like (start pos in query, num kmers, color set id), makes sense to me. So if you have a triple (34,5,88) means that starting from pos 34 in the query the next 5 kmers all have color set 88.

Ids start at 0.

If a kmer is not found in the index, then its starting position in the query will not appear within the sequence of triples as said above.

[tbenavi1]

Sure, whatever is easiest to implement and we can parse the file for our needs. Thank you.

[jermp]

CC @vikshiv

[jermp]

Hi everyone,
I implemented the thing requested here: 052f1ba in the branch https://github.com/jermp/fulgor/tree/kmer-conservation. Can you check, @tbenavi1 and @vikshiv, if this is ok for you?

[jermp]

Merged into main.

[vikshiv]

This is really great, thank you @jermp! One question. My output looks like the expected output:
chr1_PATERNAL:2000001-2010000 267 (0 35 4361161) (35 15 4377737) (50 16 4361161) (66 18 4361161) (84 29 4361161) (113 2 4361161) (115 29 4361161)

I assume color set id 4361161 refers to "all" sequences (based on the kmer conservation results I've found). What would be the easiest way to retrieve colors for a given color set? For instance, I'm interested in the colors for sets smaller than some size (where kmer conservation dips).

I'm happy to implement in my own fork since I'm sure this is a niche request, but I can put in a PR. But if you could point me to the right spot in the codebase I can get a start on it.

Another unrelated question: for my large genome query, I am currently splitting into sliding windows and use multi-threading or to run the fulgor query. Is this the fastest way to go about this?

[jermp]

Yes, it means that all the first 35 kmers in your query all have color set id 4361161. You can then pass this id to this method to get an iterator to the actual color set: https://github.com/jermp/fulgor/blob/main/include/index.hpp#L31.

Then note that each iterator over a color set has a size() method returning its size, so it's very easy to filter out color sets based on their size.

For querying, yes, I would do as you said. Otherwise, you have too fewer queries and parallelism cannot be exploited.