Skip to content

Feature request: add query length to kmer_conservation output #84

Description

@cschu

Hi,

I'd be interested in having fulgor report the length of the query in kmer_conservation in order to facilitate assessment of how well the query matches. (E.g, 70 shared 31mers are a great match for 100bp, but less so for a 150bp or 250bp read.) My input are preprocessed WGS/RNAseq reads of varying lengths, so I cannot work with a uniform query length and parsing the input fastqs to count bases for all reads ist not really efficient, especially if fulgor already has that information.

I have a customised fulgor version that already does this, but I thought it could be of use for others and the code changes are trivial. The output for the non-matching reads is theoretically not necessary, but leaving it unchanged would screw with the output format.

Thanks,
Christian

diff --git a/tools/kmer_conservation.cpp b/tools/kmer_conservation.cpp
index 1851bdf..6eb561f 100644
--- a/tools/kmer_conservation.cpp
+++ b/tools/kmer_conservation.cpp
@@ -24,14 +24,14 @@ void kmer_conservation(FulgorIndex const& index,
             index.kmer_conservation(record.seq, kmer_conservation_info);
             buff_size += 1;
             if (!kmer_conservation_info.empty()) {
-                ss << record.name << '\t' << kmer_conservation_info.size();
+                ss << record.name << '\t' << record.seq.length() << '\t' << kmer_conservation_info.size();
                 for (auto kct : kmer_conservation_info) {
                     ss << "\t(" << kct.start_pos_in_query << ' ' << kct.num_kmers << ' '
                        << kct.color_set_id << ')';
                 }
                 ss << '\n';
             } else {
-                ss << record.name << "\t0\n";
+                ss << record.name << '\t' << record.seq.length() << "\t0\n";
             }
             kmer_conservation_info.clear();
             options.increment_processed_reads();

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions