From 8a09a58f4bf66597ddf97e32c46720c4481c4fc0 Mon Sep 17 00:00:00 2001
From: Hasindu Gamaarachchi <hasindu2008@gmail.com>
Date: Sat, 27 Jul 2024 17:45:55 +1000
Subject: [PATCH] try the implementation

---
 src/eventalign.c | 39 ++++++++++++++++++++++++++++++++++++++-
 src/f5c.c        |  2 +-
 src/f5cmisc.h    |  2 +-
 3 files changed, 40 insertions(+), 3 deletions(-)

diff --git a/src/eventalign.c b/src/eventalign.c
index efc764ac..7d3a2ced 100644
--- a/src/eventalign.c
+++ b/src/eventalign.c
@@ -1656,6 +1656,7 @@ void emit_event_alignment_tsv_header(FILE* fp, int8_t print_read_names, int8_t w
     if(write_samples) {
         fprintf(fp, "\t%s", "samples");
     }
+    fprintf(fp, "\tread_pos\tread_kmer");
     fprintf(fp, "\n");
 }
 
@@ -2020,18 +2021,52 @@ std::vector<float> get_scaled_samples_for_event(const event_table* events,scalin
     return out;
 }
 
+static inline int32_t *get_event_to_base_map(const event_table* et, index_pair_t *base_to_event_map, int32_t read_len, uint32_t kmer_size){
+    int32_t *event_to_base_map = (int32_t*)malloc(sizeof(int32_t)*et->n);
+    MALLOC_CHK(event_to_base_map);
+    for(size_t i=0; i<et->n; i++){
+        event_to_base_map[i] = -1;
+    }
+    int n_kmer = read_len - kmer_size + 1;
+    for(int i=0; i<n_kmer; i++){
+        if(base_to_event_map[i].start != -1){
+            for(int j=base_to_event_map[i].start; j<=base_to_event_map[i].stop; j++){
+                event_to_base_map[j] = i;
+            }
+        }
+    }
+
+    return event_to_base_map;
+
+}
+
+static inline void sprintf_read_kmer(kstring_t *sp, int32_t *event_to_base_map, uint64_t event_idx, char *read, int32_t read_len, uint32_t kmer_size){
+
+    int32_t read_pos = event_to_base_map[event_idx];
+    if(read_pos == -1){
+        sprintf_append(sp, "\t.\t.");
+        return;
+    } else {
+        char kmer[kmer_size+1];
+        kmer_cpy(kmer, read+read_pos, kmer_size);
+        sprintf_append(sp, "\t%d\t%s", read_pos, kmer);
+    }
+
+}
 
 char *emit_event_alignment_tsv(uint32_t strand_idx,
                               const event_table* et, model_t* model, uint32_t kmer_size, scalings_t scalings,
                               const std::vector<event_alignment_t>& alignments,
                               int8_t print_read_names, int8_t scale_events, int8_t write_samples, int8_t write_signal_index, int8_t collapse,
-                              int64_t read_index, char* read_name, char *ref_name,float sample_rate, float *rawptr)
+                              int64_t read_index, char* read_name, char *ref_name,float sample_rate, float *rawptr, int32_t read_len, char *read, index_pair_t *base_to_event_map)
 {
 
     kstring_t str;
     kstring_t *sp = &str;
     str_init(sp, sizeof(char)*alignments.size()*120);
 
+    int32_t *event_to_base_map = get_event_to_base_map(et, base_to_event_map, read_len, kmer_size);
+
     size_t n_collapse = 1;
     for(size_t i = 0; i < alignments.size(); i+=n_collapse) {
 
@@ -2153,9 +2188,11 @@ char *emit_event_alignment_tsv(uint32_t strand_idx,
             sample_str.resize(sample_str.size() - 1);
             sprintf_append(sp, "\t%s", sample_str.c_str());
         }
+        sprintf_read_kmer(sp, event_to_base_map, ea.event_idx, read, read_len, kmer_size);
         sprintf_append(sp, "\n");
     }
 
+    free(event_to_base_map);
 
     //str_free(sp); //freeing is later done in free_db_tmp()
     return sp->s;
diff --git a/src/f5c.c b/src/f5c.c
index 3ced84ee..bd94fe64 100644
--- a/src/f5c.c
+++ b/src/f5c.c
@@ -874,7 +874,7 @@ void eventalign_single(core_t* core, db_t* db, int32_t i){
         db->event_alignment_result_str[i] = emit_event_alignment_sam(qname, core->m_hdr, db->bam_rec[i], *event_alignment_result, sam_out_version, &(db->et[i]), db->sig[i]->nsample, ref_len, rna, db->scalings[i]);
     } else {
         db->event_alignment_result_str[i] = emit_event_alignment_tsv(0,&(db->et[i]),core->model,core->kmer_size, db->scalings[i],*event_alignment_result, print_read_names, scale_events, write_samples, write_signal_index, collapse_events,
-                   db->read_idx[i], qname, contig, db->sig[i]->sample_rate, db->sig[i]->rawptr);
+                   db->read_idx[i], qname, contig, db->sig[i]->sample_rate, db->sig[i]->rawptr, db->read_len[i], db->read[i], db->base_to_event_map[i]);
     }
 }
 
diff --git a/src/f5cmisc.h b/src/f5cmisc.h
index 821857c4..f32bb8ba 100644
--- a/src/f5cmisc.h
+++ b/src/f5cmisc.h
@@ -79,7 +79,7 @@ char *emit_event_alignment_tsv(uint32_t strand_idx,
                               const event_table* et, model_t* model, uint32_t kmer_size,  scalings_t scalings,
                               const std::vector<event_alignment_t>& alignments,
                               int8_t print_read_names, int8_t scale_events, int8_t write_samples, int8_t write_signal_index, int8_t collapse,
-                              int64_t read_index, char* read_name, char *ref_name, float sample_rate, float *samples);
+                              int64_t read_index, char* read_name, char *ref_name, float sample_rate, float *samples, int32_t read_len, char *read, index_pair_t *base_to_event_map);
 void emit_event_alignment_tsv_header(FILE* fp, int8_t print_read_names, int8_t write_samples, int8_t write_signal_index);
 void emit_sam_header(samFile* fp, const bam_hdr_t* hdr);
 char *emit_event_alignment_sam(char* read_name, bam_hdr_t* base_hdr, bam1_t* base_record,