diff --git a/src/eventalign.c b/src/eventalign.c index 83f5929c..6c835418 100644 --- a/src/eventalign.c +++ b/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 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& 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; diff --git a/src/f5c.c b/src/f5c.c index dbc6f5c6..e8b2c3bf 100644 --- a/src/f5c.c +++ b/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]); } } diff --git a/src/f5c.h b/src/f5c.h index 972fb285..762ed146 100644 --- a/src/f5c.h +++ b/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) * diff --git a/src/f5cmisc.h b/src/f5cmisc.h index 694867cb..2fae0712 100644 --- a/src/f5cmisc.h +++ b/src/f5cmisc.h @@ -78,14 +78,14 @@ void calculate_methylation_for_read(std::map* 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& 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& 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, diff --git a/src/meth_main.c b/src/meth_main.c index 126905f1..53dec18b 100644 --- a/src/meth_main.c +++ b/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); } }