Skip to content

Commit 2178567

Browse files
imsunethimsuneth
committed
fix varmod: tag CpG entries with strand
1 parent 4c693f2 commit 2178567

26 files changed

Lines changed: 108 additions & 527 deletions

‎src/minimod.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,7 @@ typedef struct {
9595
int var_idx; // index into vars->vars[]
9696
int8_t is_insertion_only; // (cg_offsets[o] > ref_len && cg_offsets[o] <= alt_len)
9797
int8_t is_compound; // 1 if CG spans >=2 variants on the same hap (phase-only)
98+
char strand; // '+' for the C cytosine, '-' for the G (reverse-strand) cytosine
9899
int seq; // append order for sorting
99100
} cg_entry_t;
100101

‎src/varmod.c‎

Lines changed: 39 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -432,6 +432,9 @@ void load_var_map(const char* vcf_file, const char* sample_name, khash_t(varm)*
432432
vars->cg_entries[vars->cg_entries_len].var_idx = var_idx;
433433
vars->cg_entries[vars->cg_entries_len].is_insertion_only = is_ins;
434434
vars->cg_entries[vars->cg_entries_len].is_compound = 0;
435+
// o_val points at the C ('+' cytosine) or the G ('-' cytosine) of the CpG
436+
vars->cg_entries[vars->cg_entries_len].strand =
437+
(after_site[o_val] == 'C' || after_site[o_val] == 'c') ? '+' : '-';
435438
vars->cg_entries[vars->cg_entries_len].seq = vars->cg_entries_len;
436439
vars->cg_entries_len++;
437440
}
@@ -570,37 +573,48 @@ void load_var_map(const char* vcf_file, const char* sample_name, khash_t(varm)*
570573
if (is_reference_cpg(ref_seq, ref->ref_seq_length,
571574
hap_ref_pos[p], hap_ref_pos[p+1],
572575
hap_is_ins[p], hap_is_ins[p+1])) continue;
573-
int ref_cg_pos = hap_ref_pos[p];
574-
int8_t is_ins = hap_is_ins[p];
575576
// attribute to the owning variant of the C; if C is a ref base, fall back to the
576577
// owner of the G, else the nearest preceding variant base.
577578
int attr = hap_owner[p];
578579
if (attr < 0) attr = hap_owner[p+1];
579580
if (attr < 0) { for (long q = p; q >= 0; q--) { if (hap_owner[q] >= 0) { attr = hap_owner[q]; break; } } }
580581
if (attr < 0) attr = hap_idx[0];
581582

582-
// dedup against existing entries at same (ref_cg_pos, is_ins) on this hap;
583-
// for insertion entries also require matching var.pos (matches scan filter)
584-
int dup = 0;
585-
for (int e = 0; e < vars->cg_entries_len; e++) {
586-
if (vars->cg_entries[e].ref_cg_pos != ref_cg_pos) continue;
587-
if (vars->cg_entries[e].is_insertion_only != is_ins) continue;
588-
if (is_ins && vars->vars[vars->cg_entries[e].var_idx].pos != vars->vars[attr].pos) continue;
589-
if (vars->vars[vars->cg_entries[e].var_idx].hap == hap) { dup = 1; break; }
590-
}
591-
if (dup) continue;
583+
// emit both cytosines of the CpG: C at p ('+' strand) and G at p+1 ('-' strand),
584+
// each at its own reference coordinate so pileup matches the correct strand.
585+
int cg_pos[2] = { hap_ref_pos[p], hap_ref_pos[p+1] };
586+
int8_t cg_ins[2] = { hap_is_ins[p], hap_is_ins[p+1] };
587+
char cg_strand[2] = { '+', '-' };
588+
for (int s = 0; s < 2; s++) {
589+
int ref_cg_pos = cg_pos[s];
590+
int8_t is_ins = cg_ins[s];
591+
char strand = cg_strand[s];
592+
593+
// dedup against existing entries at same (ref_cg_pos, is_ins, strand) on this hap;
594+
// for insertion entries also require matching var.pos (matches scan filter)
595+
int dup = 0;
596+
for (int e = 0; e < vars->cg_entries_len; e++) {
597+
if (vars->cg_entries[e].ref_cg_pos != ref_cg_pos) continue;
598+
if (vars->cg_entries[e].is_insertion_only != is_ins) continue;
599+
if (vars->cg_entries[e].strand != strand) continue;
600+
if (is_ins && vars->vars[vars->cg_entries[e].var_idx].pos != vars->vars[attr].pos) continue;
601+
if (vars->vars[vars->cg_entries[e].var_idx].hap == hap) { dup = 1; break; }
602+
}
603+
if (dup) continue;
592604

593-
if (vars->cg_entries_len >= vars->cg_entries_cap) {
594-
vars->cg_entries_cap *= 2;
595-
vars->cg_entries = (cg_entry_t*)realloc(vars->cg_entries, sizeof(cg_entry_t) * vars->cg_entries_cap);
596-
MALLOC_CHK(vars->cg_entries);
605+
if (vars->cg_entries_len >= vars->cg_entries_cap) {
606+
vars->cg_entries_cap *= 2;
607+
vars->cg_entries = (cg_entry_t*)realloc(vars->cg_entries, sizeof(cg_entry_t) * vars->cg_entries_cap);
608+
MALLOC_CHK(vars->cg_entries);
609+
}
610+
vars->cg_entries[vars->cg_entries_len].ref_cg_pos = ref_cg_pos;
611+
vars->cg_entries[vars->cg_entries_len].var_idx = attr;
612+
vars->cg_entries[vars->cg_entries_len].is_insertion_only = is_ins;
613+
vars->cg_entries[vars->cg_entries_len].is_compound = 1;
614+
vars->cg_entries[vars->cg_entries_len].strand = strand;
615+
vars->cg_entries[vars->cg_entries_len].seq = vars->cg_entries_len;
616+
vars->cg_entries_len++;
597617
}
598-
vars->cg_entries[vars->cg_entries_len].ref_cg_pos = ref_cg_pos;
599-
vars->cg_entries[vars->cg_entries_len].var_idx = attr;
600-
vars->cg_entries[vars->cg_entries_len].is_insertion_only = is_ins;
601-
vars->cg_entries[vars->cg_entries_len].is_compound = 1;
602-
vars->cg_entries[vars->cg_entries_len].seq = vars->cg_entries_len;
603-
vars->cg_entries_len++;
604618
}
605619

606620
free(hap_seq); free(hap_ref_pos); free(hap_is_ins); free(hap_owner);
@@ -1000,6 +1014,7 @@ void varviewfreq_single(core_t * core, db_t *db, int32_t bam_i) {
10001014
for (; ei < vars->cg_entries_len && vars->cg_entries[ei].ref_cg_pos == lookup_pos; ei++) {
10011015
if (vars->cg_entries[ei].is_insertion_only != want_ins) continue;
10021016
if (vars->cg_entries[ei].is_compound && !core->opt.haplotypes) continue;
1017+
if (vars->cg_entries[ei].strand != strand) continue;
10031018
var_t var = vars->vars[vars->cg_entries[ei].var_idx];
10041019
if (want_ins && var.pos != ins_start) continue;
10051020
// phase-aware filter: only when --haplotypes is on; phased ALT (var.hap > 0) only counts reads with matching HP tag
@@ -1060,6 +1075,7 @@ void varviewfreq_single(core_t * core, db_t *db, int32_t bam_i) {
10601075
for (; ei < vars->cg_entries_len && vars->cg_entries[ei].ref_cg_pos == lookup_pos; ei++) {
10611076
if (vars->cg_entries[ei].is_insertion_only != want_ins) continue;
10621077
if (vars->cg_entries[ei].is_compound && !core->opt.haplotypes) continue;
1078+
if (vars->cg_entries[ei].strand != strand) continue;
10631079
var_t var = vars->vars[vars->cg_entries[ei].var_idx];
10641080
if (want_ins && var.pos != skip_ins_start) continue;
10651081
if (core->opt.haplotypes && var.hap > 0 && (int)read_hp != var.hap) continue;
@@ -1107,6 +1123,7 @@ void varviewfreq_single(core_t * core, db_t *db, int32_t bam_i) {
11071123
for (; ei < vars->cg_entries_len && vars->cg_entries[ei].ref_cg_pos == lookup_pos; ei++) {
11081124
if (vars->cg_entries[ei].is_insertion_only != want_ins) continue;
11091125
if (vars->cg_entries[ei].is_compound && !core->opt.haplotypes) continue;
1126+
if (vars->cg_entries[ei].strand != strand) continue;
11101127
var_t var = vars->vars[vars->cg_entries[ei].var_idx];
11111128
if (want_ins && var.pos != skip_ins_start) continue;
11121129
uint16_t offset = want_ins ? (uint16_t)skip_ins_offset : REF_OFFSET(out_pos, var.pos);

‎test/expected/varmod/dna_4mC_5mC_mm_chr22.mm.varfreq.bedmethyl‎

Lines changed: 0 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,33 +1,22 @@
11
chr22 11211064 11211065 m 7 + 11211064 11211065 255,0,0 7 0.000000 11211064 1/1 T C 0
2-
chr22 11211065 11211066 m 4 + 11211065 11211066 255,0,0 4 0.000000 11211064 1/1 T C 1
32
chr22 11211065 11211066 m 3 - 11211065 11211066 255,0,0 3 100.000000 11211064 1/1 T C 1
43
chr22 11211066 11211067 m 4 + 11211066 11211067 255,0,0 4 0.000000 11211066 1/1 T C 0
5-
chr22 11211067 11211068 m 6 + 11211067 11211068 255,0,0 6 0.000000 11211066 1/1 T C 1
64
chr22 11211105 11211106 m 4 + 11211105 11211106 255,0,0 4 0.000000 11211105 1/1 T C 0
7-
chr22 11211224 11211225 m 2 + 11211224 11211225 255,0,0 2 0.000000 11211224 1/1 C G 0
85
chr22 11211224 11211225 m 9 - 11211224 11211225 255,0,0 9 0.000000 11211224 1/1 C G 0
96
chr22 11211273 11211274 m 7 + 11211273 11211274 255,0,0 7 0.000000 11211273 0/1 T C 0
10-
chr22 11211274 11211275 m 2 + 11211274 11211275 255,0,0 2 0.000000 11211273 0/1 T C 1
117
chr22 11211522 11211523 m 7 + 11211522 11211523 255,0,0 7 0.000000 11211523 0/1 A G -1
12-
chr22 11211523 11211524 m 6 + 11211523 11211524 255,0,0 6 0.000000 11211523 0/1 A G 0
138
chr22 11211523 11211524 m 2 - 11211523 11211524 255,0,0 2 100.000000 11211523 0/1 A G 0
149
chr22 11211648 11211649 m 2 + 11211648 11211649 255,0,0 2 50.000000 11211648 1/1 T C 0
1510
chr22 11212145 11212146 m 5 + 11212145 11212146 255,0,0 5 0.000000 11212145 1/1 T C 0
1611
chr22 11212146 11212147 m 3 - 11212146 11212147 255,0,0 3 0.000000 11212145 1/1 T C 1
1712
chr22 11212977 11212978 m 2 + 11212977 11212978 255,0,0 2 0.000000 11212977 1/1 G C 0
18-
chr22 11212977 11212978 m 1 - 11212977 11212978 255,0,0 1 0.000000 11212977 1/1 G C 0
19-
chr22 11212978 11212979 m 2 + 11212978 11212979 255,0,0 2 0.000000 11212977 1/1 G C 1
2013
chr22 11212978 11212979 m 1 - 11212978 11212979 255,0,0 1 0.000000 11212977 1/1 G C 1
2114
chr22 11213000 11213001 m 6 + 11213000 11213001 255,0,0 6 0.000000 11213000 1/1 T C 0
22-
chr22 11213001 11213002 m 1 + 11213001 11213002 255,0,0 1 0.000000 11213000 1/1 T C 1
2315
chr22 11213183 11213184 m 7 + 11213183 11213184 255,0,0 7 0.000000 11213184 0/1 A G -1
24-
chr22 11213184 11213185 m 2 + 11213184 11213185 255,0,0 2 0.000000 11213184 0/1 A G 0
2516
chr22 11213184 11213185 m 1 - 11213184 11213185 255,0,0 1 100.000000 11213184 0/1 A G 0
2617
chr22 11213334 11213335 m 6 + 11213334 11213335 255,0,0 6 0.000000 11213334 1/1 A C 0
2718
chr22 11214662 11214663 m 5 + 11214662 11214663 255,0,0 5 0.000000 11214662 0/1 T C 0
28-
chr22 11214663 11214664 m 4 + 11214663 11214664 255,0,0 4 0.000000 11214662 0/1 T C 1
2919
chr22 11214868 11214869 m 4 + 11214868 11214869 255,0,0 4 0.000000 11214869 1/1 C G -1
30-
chr22 11214869 11214870 m 4 + 11214869 11214870 255,0,0 4 0.000000 11214869 1/1 C G 0
3120
chr22 11214869 11214870 m 7 - 11214869 11214870 255,0,0 7 0.000000 11214869 1/1 C G 0
3221
chr22 11215132 11215133 m 5 + 11215132 11215133 255,0,0 5 0.000000 11215133 1/1 A G -1
3322
chr22 11215133 11215134 m 7 - 11215133 11215134 255,0,0 7 0.000000 11215133 1/1 A G 0

0 commit comments

Comments
 (0)