Read kmer

src/eventalign.c

@@ -1642,7 +1642,7 @@ void emit_sam_header(samFile* fp, const bam_hdr_t* hdr)
 }
 
 
-void emit_event_alignment_tsv_header(FILE* fp, int8_t print_read_names, int8_t write_samples, int8_t write_signal_index)
+void emit_event_alignment_tsv_header(FILE* fp, int8_t print_read_names, int8_t write_samples, int8_t write_signal_index, int8_t write_read_kmer)
 {
     fprintf(fp, "%s\t%s\t%s\t%s\t%s\t", "contig", "position", "reference_kmer",
             (print_read_names? "read_name" : "read_index"), "strand");
@@ -1656,6 +1656,11 @@ void emit_event_alignment_tsv_header(FILE* fp, int8_t print_read_names, int8_t w
     if(write_samples) {
         fprintf(fp, "\t%s", "samples");
     }
+
+    if(write_read_kmer){
+        fprintf(fp, "\t%s\t%s", "read_pos", "read_kmer");
+    }
+
     fprintf(fp, "\n");
 }
 
@@ -2034,18 +2039,142 @@ std::vector<float> get_scaled_samples_for_event(const event_table* events,scalin
     return out;
 }
 
+typedef struct {
+    int32_t capacity;
+    int32_t size;
+    int32_t *map;
+    int32_t ref_pos_start;
+} ref2read_t;
+
+static inline ref2read_t get_ref2read_map(bam1_t *record, int32_t read_len){
+
+    ref2read_t ref2read;
+    ref2read.capacity = read_len;
+    ref2read.size = 0;
+
+    ref2read.map = (int32_t*)malloc(sizeof(int32_t)*ref2read.capacity);
+    MALLOC_CHK(ref2read.map);
+    for(int i = 0; i < ref2read.capacity; i++){
+        ref2read.map[i] = -1;
+    }
+
+    // This code is derived from bam_fillmd1_core
+    uint32_t *cigar = bam_get_cigar(record);
+    const bam1_core_t *c = &record->core;
+
+    // read pos is an index into the original sequence that is present in the FASTQ
+    // on the strand matching the reference
+    int read_pos = 0;
+
+    int ref_pos = c->pos;
+    ref2read.ref_pos_start = ref_pos;
+
+    for (uint32_t ci = 0; ci < c->n_cigar; ++ci) {
+
+        int cigar_len = cigar[ci] >> 4;
+        int cigar_op = cigar[ci] & 0xf;
+
+        // Set the amount that the ref/read positions should be incremented
+        // based on the cigar operation
+        int read_inc = 0;
+        int ref_inc = 0;
+
+        // Process match between the read and the reference
+        bool is_aligned = false;
+        if(cigar_op == BAM_CMATCH || cigar_op == BAM_CEQUAL || cigar_op == BAM_CDIFF) {
+            is_aligned = true;
+            read_inc = 1;
+            ref_inc = 1;
+        } else if(cigar_op == BAM_CDEL) {
+            ref_inc = 1;
+        } else if(cigar_op == BAM_CREF_SKIP) {
+            ref_inc = 1;
+        } else if(cigar_op == BAM_CINS) {
+            read_inc = 1;
+        } else if(cigar_op == BAM_CSOFT_CLIP) {
+            read_inc = 1; // special case, do not use read_stride
+        } else if(cigar_op == BAM_CHARD_CLIP) {
+            read_inc = 0;
+        } else {
+            printf("Cigar: %d\n", cigar_op);
+            assert(false && "Unhandled cigar operation");
+        }
+
+        // Iterate over the pairs of aligned bases
+        for(int j = 0; j < cigar_len; ++j) {
+            if(is_aligned) {
+                ref2read.map[ref2read.size] = read_pos;
+            }
+
+            // increment
+            read_pos += read_inc;
+            ref_pos += ref_inc;
+            ref2read.size += ref_inc;
+
+            if(ref2read.size >= ref2read.capacity){
+                ref2read.capacity = ref2read.capacity*+1000;
+                ref2read.map = (int32_t*)realloc(ref2read.map, sizeof(int32_t)*ref2read.capacity);
+                MALLOC_CHK(ref2read.map);
+                for(int i = ref2read.size; i < ref2read.capacity; i++){
+                    ref2read.map[i] = -1;
+                }
+            }
+        }
+    }
+
+    // for(int i = 0; i < ref2read.size; i++){
+    //     fprintf(stderr, "ref2read_map[%d]=%d\n",i+ref2read.ref_pos_start, ref2read.map[i]);
+    // }
+
+    return ref2read;
+}
+
+static inline void sprintf_read_kmer(kstring_t *sp, ref2read_t ref2read, uint64_t ref_pos, char *read, uint32_t kmer_size){
+
+    int32_t map_idx = ref_pos - ref2read.ref_pos_start;
+    assert(map_idx >= 0 && map_idx < ref2read.size);
+    int32_t read_pos = ref2read.map[map_idx];
+
+    // if(read_pos == -1){
+    //     sprintf_append(sp, "\t.\t.");
+    //     return;
+    // } else {
+    char kmer[kmer_size+1];
+    for(uint32_t i = 0; i < kmer_size; i++){
+        int ref_base_idx = map_idx + i;
+        assert(ref_base_idx >= 0 && ref_base_idx < ref2read.size);
+        int read_base_idx=ref2read.map[ref_base_idx];
+        if(read_base_idx == -1){
+            kmer[i] = '.';
+        } else {
+            kmer[i] = read[read_base_idx];
+        }
+    }
+    kmer[kmer_size] = '\0';
+    if(read_pos == -1){
+        sprintf_append(sp, "\t.\t%s", kmer);
+    } else {
+        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)
+                              int8_t print_read_names, int8_t scale_events, int8_t write_samples, int8_t write_signal_index, int8_t collapse, int8_t write_read_kmer,
+                              int64_t read_index, char* read_name, char *ref_name,float sample_rate, float *rawptr, int32_t read_len, char *read, bam1_t* bam_record)
 {
 
     kstring_t str;
     kstring_t *sp = &str;
     str_init(sp, sizeof(char)*alignments.size()*120);
 
+    ref2read_t ref2read = {0};
+    if(write_read_kmer){
+        ref2read = get_ref2read_map(bam_record, read_len);
+    }
+
     size_t n_collapse = 1;
     for(size_t i = 0; i < alignments.size(); i+=n_collapse) {
 
@@ -2167,9 +2296,15 @@ 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());
         }
+        if(write_read_kmer){
+            sprintf_read_kmer(sp, ref2read, ea.ref_position, read, kmer_size);
+        }
         sprintf_append(sp, "\n");
     }
 
+    if(write_read_kmer){
+        free(ref2read.map);
+    }
 
     //str_free(sp); //freeing is later done in free_db_tmp()
     return sp->s;

src/f5c.c

@@ -865,6 +865,7 @@ void eventalign_single(core_t* core, db_t* db, int32_t i){
     int8_t sam_output = (core->opt.flag & F5C_SAM) ? 1 : 0;
     int8_t paf_output = (core->opt.flag & F5C_PAF) ? 1 : 0;
     int8_t m6anet_output = (core->opt.flag & F5C_M6ANET) ? 1 : 0;
+    int8_t write_read_kmer = (core->opt.flag & F5C_READ_KMER) ? 1 : 0;
     int8_t rna = (core->opt.flag & F5C_RNA) ? 1 : 0;
 
     if(paf_output){
@@ -878,8 +879,8 @@ void eventalign_single(core_t* core, db_t* db, int32_t i){
         db->event_alignment_result_str[i] = emit_event_alignment_tsv_m6anet(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);
     } 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->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, write_read_kmer,
+                   db->read_idx[i], qname, contig, db->sig[i]->sample_rate, db->sig[i]->rawptr, db->read_len[i], db->read[i], db->bam_rec[i]);
     }
 }
 

src/f5c.h

@@ -58,6 +58,7 @@
 #define F5C_R10 0x40000 //r10
 #define F5C_PAF 0x80000 //paf (eventalign only)
 #define F5C_M6ANET 0x100000 //m6anet (eventalign only)
+#define F5C_READ_KMER 0x200000 //read-kmer (eventalign only)
 
 /*************************************************************
  * flags for a read status (related to db_t->read_stat_flag) *

src/f5cmisc.h

@@ -78,14 +78,14 @@ void calculate_methylation_for_read(std::map<int, ScoredSite>* site_score_map, c
 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);
+                              int8_t print_read_names, int8_t scale_events, int8_t write_samples, int8_t write_signal_index, int8_t collapse, int8_t write_read_kmer,
+                              int64_t read_index, char* read_name, char *ref_name, float sample_rate, float *samples, int32_t read_len, char *read, bam1_t* bam_record);
 char *emit_event_alignment_tsv_m6anet(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);
-void emit_event_alignment_tsv_header(FILE* fp, int8_t print_read_names, int8_t write_samples, int8_t write_signal_index);
+void emit_event_alignment_tsv_header(FILE* fp, int8_t print_read_names, int8_t write_samples, int8_t write_signal_index, int8_t write_read_kmer);
 void emit_event_alignment_tsv_m6anet_header(FILE* fp, int8_t print_read_names, 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,

src/meth_main.c

@@ -106,6 +106,7 @@ static struct option long_options[] = {
     {"paf",no_argument,0,'c'},                     //47 if print in paf format (only for eventalign)
     {"sam-out-version",required_argument,0,0},     //48 specify the version of the sam output for eventalign (eventalign only)
     {"m6anet",no_argument,0,0},                    //49 m6anet output (eventalign only)
+    {"read-kmer",no_argument,0,0},                 //50 read kmer (eventalign only)
     {0, 0, 0, 0}};
 
 
@@ -454,6 +455,13 @@ int meth_main(int argc, char* argv[], int8_t mode) {
                 exit(EXIT_FAILURE);
             }
             yes_or_no(&opt, F5C_M6ANET, longindex, "yes", 1);
+        } else if (c == 0 && longindex == 50){ //read kmer
+            if(mode!=1){
+                ERROR("%s","Option --read-kmer is available only in eventalign");
+                exit(EXIT_FAILURE);
+            }
+            WARNING("%s", "Option --read-kmer is experimental. Exercise caution.");
+            yes_or_no(&opt, F5C_READ_KMER, longindex, "yes", 1);
         }
     }
 
@@ -528,6 +536,7 @@ int meth_main(int argc, char* argv[], int8_t mode) {
         fprintf(fp_help,"   --signal-index             write the raw signal start and end index values for the event to the tsv output\n");
         fprintf(fp_help,"   --rna                      the dataset is direct RNA\n");
         fprintf(fp_help,"   --collapse-events          collapse events that stays on the same reference k-mer\n");
+        //fprintf(fp_help,"   --read-kmer                print the read k-mer\n");
     }
         fprintf(fp_help,"   --min-recalib-events INT   minimum number of events to recalbrate (decrease if your reads are very short and could not calibrate) [%d]\n",opt.min_num_events_to_rescale);
 
@@ -582,6 +591,7 @@ int meth_main(int argc, char* argv[], int8_t mode) {
         int8_t sam_output =  (core->opt.flag & F5C_SAM) ? 1 : 0 ;
         int8_t paf_output =  (core->opt.flag & F5C_PAF) ? 1 : 0 ;
         int8_t m6anet_output =  (core->opt.flag & F5C_M6ANET) ? 1 : 0 ;
+        int8_t write_read_kmer = (core->opt.flag & F5C_READ_KMER) ? 1 : 0 ;
 
         if(sam_output && paf_output){
             ERROR("%s","-c and --sam cannot be used together");
@@ -603,7 +613,7 @@ int meth_main(int argc, char* argv[], int8_t mode) {
         } else if (m6anet_output){
             emit_event_alignment_tsv_m6anet_header(stdout, print_read_names, write_signal_index);
         } else{
-            emit_event_alignment_tsv_header(stdout, print_read_names, write_samples, write_signal_index);
+            emit_event_alignment_tsv_header(stdout, print_read_names, write_samples, write_signal_index, write_read_kmer);
         }
 
     }