Add support for SQ-AN in SAMSequenceDictionary by magicDGS · Pull Request #956 · samtools/htsjdk

magicDGS · 2017-07-28T10:46:25Z

Description

Adding the support for the new @SQ AN tag for alternative sequences (samtools/hts-specs#220).

Checklist

Code compiles correctly
New tests covering changes and new functionality
All tests passing
Extended the README / documentation, if necessary
Is not backward compatible (breaks binary or source compatibility)

codecov-io · 2017-07-28T12:31:37Z

Codecov Report

Merging #956 into master will increase coverage by 0.013%.
The diff coverage is 78.125%.

@@               Coverage Diff               @@
##              master      #956       +/-   ##
===============================================
+ Coverage     65.963%   65.976%   +0.013%     
- Complexity      7663      7675       +12     
===============================================
  Files            537       537               
  Lines          32553     32583       +30     
  Branches        5528      5533        +5     
===============================================
+ Hits           21473     21497       +24     
- Misses          8930      8932        +2     
- Partials        2150      2154        +4

Impacted Files	Coverage Δ	Complexity Δ
...in/java/htsjdk/samtools/SAMSequenceDictionary.java	`73.649% <100%> (+0.921%)`	`46 <1> (+1)`	⬆️
...c/main/java/htsjdk/samtools/SAMSequenceRecord.java	`74.766% <74.074%> (+0.376%)`	`47 <10> (+11)`	⬆️
...sjdk/samtools/util/Md5CalculatingOutputStream.java	`69.231% <0%> (-7.692%)`	`8% <0%> (-1%)`
src/main/java/htsjdk/samtools/BAMFileReader.java	`64.286% <0%> (+0.857%)`	`39% <0%> (+1%)`	⬆️

pshapiro4broad · 2017-07-28T17:49:13Z

+     * @return the contig associated to the 'originalName/altName', with the AN tag including the altName
+     */
+    public SAMSequenceRecord addAlternativeSequenceName(final String originalName,
+            final String altName) {


According to codecov this new code isn't covered by a test.

lbergelson

@magicDGS Thanks for this work. I made a few comments.

I think it's better to store the AN values in the attributes map with the other attributes so it doesn't need special handling all over the place. It makes the add a single alias method gross because it has to do string manipulation now, but I think it's worth it to avoid adding additional complexity everywhere else.

lbergelson · 2017-08-09T15:11:08Z


+    /** Returns unmodifiable set with alternative sequence names. */
+    @XmlTransient //we use the field instead of getter/setter
+    public Set<String> getAlternativeSequeneNames() {


typo: getAlternativeSequeneNames()

lbergelson · 2017-08-09T15:23:23Z

            }
        }
        else {
+            // TODO: should this also check in the alternative sequences names?


It seems like this should match on alternative sequence names if present

Done (although maybe not in the way you were expecting).

lbergelson · 2017-08-09T15:24:44Z

        if (mSequenceLength != that.mSequenceLength) return false;
        if (!attributesEqual(that)) return false;
        if (mSequenceName != that.mSequenceName) return false; // Compare using == since we intern() the name
+        // TODO: should this also compare the alternative names?


it seems strange to leave alternative names out of equals especially given the weaker isSameSequence() also exists

Done by comparing the sets.

lbergelson · 2017-08-09T15:27:04Z

+     * <code>1,chr1,chr01,01,CM000663,NC_000001.10</code> or
+     * <code>MT,chrM</code>.
+     *
+     * @param originalName


there are funny line breaks here for some reason

oh, it matches the weird line breaks above, lets fix those too...

lbergelson · 2017-08-09T15:28:36Z

 {
    public static final long serialVersionUID = 1L; // AbstractSAMHeaderRecord implements Serializable
    private String mSequenceName = null; // Value must be interned() if it's ever set/modified
+    @XmlAttribute(name = "alternative_sequece_names")


typo: alternative_sequece_names

lbergelson · 2017-08-09T15:35:59Z

        fields[1] = SAMSequenceRecord.SEQUENCE_NAME_TAG + TAG_KEY_VALUE_SEPARATOR + sequenceRecord.getSequenceName();
        fields[2] = SAMSequenceRecord.SEQUENCE_LENGTH_TAG + TAG_KEY_VALUE_SEPARATOR + Integer.toString(sequenceRecord.getSequenceLength());
        encodeTags(sequenceRecord, fields, 3);
+        if (sequenceRecord.hasAlternativeSequenceNames()) fields[fields.length - 1] = SAMSequenceRecord.ALTERNATIVE_SEQUENCE_NAME_TAG + TAG_KEY_VALUE_SEPARATOR + String.join(",", sequenceRecord.getAlternativeSequeneNames());


Why does this need special handling? it seems like the alternate contig names should be added to the attributes map and handled the same way as the other attributes.

lbergelson · 2017-08-09T15:40:22Z

@@ -220,6 +223,9 @@ public boolean equals(Object o) {
     * <code>1,chr1,chr01,01,CM000663,NC_000001.10</code> e.g:


I think the comment on the equals method above needs updating to say that alias's are not considered but Alternative Sequences are.

lbergelson · 2017-08-09T15:44:30Z

    public static final long serialVersionUID = 1L; // AbstractSAMHeaderRecord implements Serializable
    private String mSequenceName = null; // Value must be interned() if it's ever set/modified
+    @XmlAttribute(name = "alternative_sequece_names")
+    private Set<String> mAlternativeSequenceName = new LinkedHashSet<>();


This should be handled the same as the other attributes, stored in the attributes map. It means doing annoying string processing when you add a single alias to an existing set of aliases, but I think making it consistent is worthwhile

lbergelson · 2017-08-09T15:45:45Z

+
+    /** Adds an alternative sequence name if it is not the same as the sequence name or it is not present already. */
+    public void addAlternativeSequenceName(final String name) {
+        if (!mSequenceName.equals(name) && ! mAlternativeSequenceName.contains(name) ) {


These should change to store the names in the attributes map so that you don't need special handling for all the places that handle the tags.

lbergelson · 2017-08-09T15:46:38Z

+
+    /** Adds an alternative sequence name if it is not the same as the sequence name or it is not present already. */
+    public void addAlternativeSequenceName(final String name) {
+        if (!mSequenceName.equals(name) && ! mAlternativeSequenceName.contains(name) ) {


this and the set method should validate that the alternative names are legal according to the regex provided in the sam spec.

magicDGS · 2017-08-21T14:19:45Z

I implement your comments in two different PRs, @lbergelson. I am not sure if it is a good idea for efficiency to store the AN tag in the attributes map, but I did in the second commit as you requested, just in case we should revert that changes because you didn't mean that.

Back to you @lbergelson!

magicDGS · 2017-09-15T13:19:29Z

Friendly reminder here, @lbergelson!

lbergelson · 2017-09-18T17:28:41Z

@magicDGS Hey, the changes you made here look good. I've been debating with myself what the right solution is though. I hadn't thought through what a large potential performance impact the changes I asked for would be. I'm torn because it's much cleaner than the alternative, but introducing a gross performance regression for minor code cleanlyness seems wrong.

magicDGS · 2017-09-19T09:21:54Z

Thanks for the reply @lbergelson - should I revert the last commit?

magicDGS · 2017-10-09T10:12:58Z

Friendly ping here, @lbergelson!

jmarshall · 2018-05-28T10:06:21Z

ALTERNATIVE_SEQUENCE_NAME_REGEXP is missing _ and | compared to the definition of AN in the SAM spec. Is there a reason for this?

magicDGS · 2018-05-28T10:18:17Z

I can't remember - maybe it was still under discussion if they should be allowed or I just forgot to add them. Changed in next commit matching the current definition.

jmarshall · 2018-06-06T14:50:23Z

Apologies and FYI: After further contemplation of the issue (see samtools/hts-specs#124 (comment)), I suspect we may end up allowing colons in AN after all.

magicDGS · 2018-06-06T14:59:43Z

@jmarshall - also in the AN? I think that not allowing it in the alternative-name is a good solution for unambiguously define contigs with colons: although the main contig name has colons, if the AN does not then it is safe to use to define intervals...

jmarshall · 2018-09-05T17:49:17Z

Yes, also in AN. See the draft PR samtools/hts-specs#333.

I used to think as you do, but… For better or worse these things exist in the wild, so IMHO we have to give up on thinking AN might be a transition mechanism away from allowing them. Reread that comment, which contains pseudocode showing how tools should parse their interval notation even in the presence of stupid colons in names (and I plan to turn that into an appendix for the spec).

jmarshall · 2019-02-04T12:59:49Z

FYI RNAME changes have now landed in the SAM spec: in particular the allowed characters (which include colons!) for SN and AN are now the same, so a separate ALTERNATIVE_SEQUENCE_NAME_REGEXP is unnecessary.

nh13 · 2019-03-09T23:21:25Z

    private static char[] WHITESPACE_CHARS = {' ', '\t', '\n', '\013', '\f', '\r'}; // \013 is vertical tab

+    // alternative sequence name regexp
+    private static final Pattern ALTERNATIVE_SEQUENCE_NAME_REGEXP = Pattern.compile("[0-9A-Za-z][0-9A-Za-z*+.@_|-]*");


no longer needed, use SAMSequenceRecord.LEGAL_RNAME_PATTERN

yfarjoun · 2019-06-22T03:54:24Z

let's resolve this and move it past the finish line....

nh13 · 2019-06-22T05:12:45Z

@yfarjoun agreed. I am already using Picard and this PR to store alternative coffee tog names, and I actually use them! fulcrumgenomics/fgbio#467

pshapiro4broad · 2019-07-08T18:36:05Z

+    public void addAlternativeSequenceName(final String name) {
+        validateAltRegExp(name);
+        final Set<String> altSequences = getAlternativeSequenceNames();
+        if (!mSequenceName.equals(name) && ! altSequences.contains(name) ) {


inconsistent whitespace

Suggested change

if (!mSequenceName.equals(name) && ! altSequences.contains(name) ) {

if (!mSequenceName.equals(name) && !altSequences.contains(name)) {

yfarjoun · 2019-07-08T18:35:09Z

+                // comma-separated
+                {"chr1,alt"},
+                // coordinate-like
+                {"chr1:1000"}, {"chr1:100-200"}


add a testcase with "{" and "}"

pshapiro4broad · 2019-07-08T18:37:48Z

+    @XmlAttribute(name = "alternative_sequence_names")
+    public Set<String> getAlternativeSequenceNames() {
+        final String anTag = getAttribute(ALTERNATIVE_SEQUENCE_NAME_TAG);
+        return (anTag == null) ? new LinkedHashSet<>() : new LinkedHashSet<>(Arrays.asList(anTag.split(ALTERNATIVE_SEQUENCE_NAME_SEPARATOR)));


This shouldn't return a modifiable set, should it? E.g. I would use Collections.emptySet() for the null case.

yfarjoun · 2020-04-23T21:08:31Z

I'm closing this PR as I don't know how to push to daniel's fork....I'll open it with a new branch...sorry for the inconvenience.

* Add support for SQ-AN in SAMSequenceDictionary Co-authored-by: Daniel Gómez-Sánchez <daniel.gomez.sanchez@hotmail.es>

pshapiro4broad reviewed Jul 28, 2017

View reviewed changes

lbergelson requested changes Aug 9, 2017

View reviewed changes

magicDGS force-pushed the dgs_sq_an_tag branch from 9e34333 to ccad70a Compare October 9, 2017 10:12

magicDGS mentioned this pull request Mar 1, 2018

CreateSequenceDictionary support for alternative names (@SQ-AN) broadinstitute/picard#1127

Merged

5 tasks

magicDGS added 6 commits May 29, 2018 10:11

Add support for SQ-AN in SAMSequenceDictionary

f3f5996

Fix XML seralization

f4d4bfa

Add test to uncovered new code

e437a3f

Address minor comments

a540558

Handle AN tag in the attributes map

a8685ee

Fix regular expression (and add tests failing previously)

bfe2206

magicDGS force-pushed the dgs_sq_an_tag branch from dd431e5 to bfe2206 Compare May 29, 2018 08:11

magicDGS mentioned this pull request Oct 28, 2018

add SAM 1.6 to supported versions #1186

Open

nh13 reviewed Mar 9, 2019

View reviewed changes

pshapiro4broad reviewed Jul 8, 2019

View reviewed changes

yfarjoun reviewed Jul 8, 2019

View reviewed changes

yfarjoun self-assigned this Jul 8, 2019

pshapiro4broad reviewed Jul 8, 2019

View reviewed changes

yfarjoun added the Waiting for revisions This PR has received comments from reviewers and is waiting for the Author to respond label Sep 30, 2019

yfarjoun closed this Apr 23, 2020

yfarjoun mentioned this pull request Apr 23, 2020

Add support for SQ-AN in SAMSequenceDictionary #956 #1474

Merged

yfarjoun pushed a commit that referenced this pull request Sep 2, 2020

Add support for SQ-AN in SAMSequenceDictionary #956 (#1474)

5240c40

* Add support for SQ-AN in SAMSequenceDictionary Co-authored-by: Daniel Gómez-Sánchez <daniel.gomez.sanchez@hotmail.es>

		@@ -220,6 +223,9 @@ public boolean equals(Object o) {
		* <code>1,chr1,chr01,01,CM000663,NC_000001.10</code> e.g:

	if (!mSequenceName.equals(name) && ! altSequences.contains(name) ) {
	if (!mSequenceName.equals(name) && !altSequences.contains(name)) {

Conversation

magicDGS commented Jul 28, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Checklist

Uh oh!

codecov-io commented Jul 28, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lbergelson left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

magicDGS commented Aug 21, 2017

Uh oh!

magicDGS commented Sep 15, 2017

Uh oh!

lbergelson commented Sep 18, 2017

Uh oh!

magicDGS commented Sep 19, 2017

Uh oh!

magicDGS commented Oct 9, 2017

Uh oh!

jmarshall commented May 28, 2018

Uh oh!

magicDGS commented May 28, 2018

Uh oh!

jmarshall commented Jun 6, 2018

Uh oh!

magicDGS commented Jun 6, 2018

Uh oh!

jmarshall commented Sep 5, 2018

Uh oh!

jmarshall commented Feb 4, 2019

Uh oh!

Choose a reason for hiding this comment

magicDGS commented Jul 28, 2017 •

edited

Loading

codecov-io commented Jul 28, 2017 •

edited

Loading