Added base modification interfaces.

[jkbonfield]

First draft of the base modification (methylation etc) API for comment on interface.

I'll work on samtools to add an option to get this into the mpileup text output too, once I've figured out exactly what we want it to look like. (Probably like one of the test/pileup_mod choices.)

---

The general principle is

1. Allocate a base modification state
2. Parse mods
3. Iterate over mods or query positions
4. Destroy state.

An example usage is as follows (minus error checking for clarity):

hts_base_mod_state *m = hts_base_mod_state_alloc();
sam_read1(fp, hdr, b);
bam_parse_basemod(b, m);

hts_base_mod mod;
while ((n = bam_next_basemod(b, m, &mod, 1, &pos) > 0))
    printf("Modified %c to %c at %d\n",
        mod.canonical_base, mod.modified_base, pos);
}
hts_base_mod_state_free(m);

If being called from an iterator that marches along the query sequence
we may wish to have a loop calling bam_mods_at_next_pos instead (pos
0, 1, 2, 3 etc).

A simple example of this is in test/test_mod.c

It may also be used from within a pileup iterator which marches along
query position instead of sequence position. In this case call
bam_mods_at_qpos instead with pileup->qpos; which is normally +1 each
time but may be more or less with indels.

A simple example of this is in test/pileup_mod.c although it doesn't
look for modifications within indels. (That would need to mix in
bam_mods_at_next_pos calls for each base in the indel.)

Note currently bam_mods_at_qpos is a brain-dead loop around
bam_mods_at_next_pos which means it may be excessive when huge
soft-clips are being used, however this optimisation may happen later.

Right now the overhead of pileup with modifications isn't vast. On a
~1Kb long seqs it's around 70% slower than the same code without these
calls. On ~50Kb seqs it drops to aroudn 20% slower. Hence the
significant cost in the initial bam_parse_basemod call.

[jmarshall]

rmd?

[jkbonfield]

Lol cut, paste ;)

[whitwham]

Removed a couple of my comments as I think the problem is elsewhere. Though as far as I can tell, bam_next_basemod can never return -1 for an error.

[cjw85]

@jkbonfield Thanks for this.

I've been using the branch in anger recently with a little program to summarize human 5mC data to bedMethyl format. I've not observed any issue. The performance is fine in this use case particularly writing out all contexts; writing the records to file becomes an issue. When filtering to CpG sites I can process 10s of fold coverage of chr11 in less than half a minute --- easily sufficient for practical purposes.

[jkbonfield]

Thanks. It's always useful to get comments "from the wild" rather than just internal to the group.

It's long overdue merging, and I think we're pretty much there.

[jmarshall]

PR #1320 adds some minor source file fixes that were not noticed before this was merged.