Mismatched read bases rendered incorrectly #523

Closed
09140930 opened this Issue Oct 17, 2014 · 7 comments

Comments

Projects
None yet
2 participants
@09140930

Hi - apologies in advance if this stems from a misunderstanding.

JBrowse 1.11.5.
I'm looking at 4 reads overlapping a position where a variation exists. The reference is "A", and the four observed bases are G/C/C/C/C. They are highlighted as being mismatches (which I expect), but are rendered as T/A/A/A/A, which makes no sense to me. Maybe I'm misunderstanding something?

The attached screenshot demonstrates what I mean. The popup refers to the bottom-most of the 4 visible reads.

I've checked the underlying BAM and have confirmed that the bases are indeed G/C/C/C/C.

Perhaps this is a bug?

screen shot 2014-10-17 at 19 42 27

@cmdcolin

This comment has been minimized.

Show comment
Hide comment
@cmdcolin

cmdcolin Oct 20, 2014

Contributor

Hi there,
I see on your displayed read that the MD tag displays "88A1" which means 88 matches, an A mismatch, and then a match. Overall the CIGAR string is combined with the MD tag to produce the visualizaation (the CIGAR string contains the extra information about the gapped alignment for example)

What kind of aligner is being used?

Contributor

cmdcolin commented Oct 20, 2014

Hi there,
I see on your displayed read that the MD tag displays "88A1" which means 88 matches, an A mismatch, and then a match. Overall the CIGAR string is combined with the MD tag to produce the visualizaation (the CIGAR string contains the extra information about the gapped alignment for example)

What kind of aligner is being used?

@09140930

This comment has been minimized.

Show comment
Hide comment
@09140930

09140930 Oct 21, 2014

Hi,
This BAM was created with TopHat 2.0.10.

I think the MD is OK - the SAM specification (v1) implies the the MD should show the reference base, not the observed base. In this case, the reference is A, so the MD is fine. The actual observed base on the read is C.

So I guess something else must be going on here?

Thanks.

From spec:
"The MD field aims to achieve SNP/indel calling without looking at the reference. For example, a string '10A5^AC6'
means from the leftmost reference base in the alignment, there are 10 matches followed by an A on the reference which
is di�fferent from the aligned read base"

Hi,
This BAM was created with TopHat 2.0.10.

I think the MD is OK - the SAM specification (v1) implies the the MD should show the reference base, not the observed base. In this case, the reference is A, so the MD is fine. The actual observed base on the read is C.

So I guess something else must be going on here?

Thanks.

From spec:
"The MD field aims to achieve SNP/indel calling without looking at the reference. For example, a string '10A5^AC6'
means from the leftmost reference base in the alignment, there are 10 matches followed by an A on the reference which
is di�fferent from the aligned read base"

@cmdcolin

This comment has been minimized.

Show comment
Hide comment
@cmdcolin

cmdcolin Oct 21, 2014

Contributor

Oh wow! thanks for the correction. I guess in that case I am still confused where jbrowse is getting the 'A' that it is displaying as the SNP. Would you happen to have a small piece of your BAM file that we could check out?

Thanks

Contributor

cmdcolin commented Oct 21, 2014

Oh wow! thanks for the correction. I guess in that case I am still confused where jbrowse is getting the 'A' that it is displaying as the SNP. Would you happen to have a small piece of your BAM file that we could check out?

Thanks

@09140930

This comment has been minimized.

Show comment
Hide comment
@09140930

09140930 Oct 23, 2014

Can't attach BAM here; sent to Colin's gmail.

Can't attach BAM here; sent to Colin's gmail.

@cmdcolin

This comment has been minimized.

Show comment
Hide comment
@cmdcolin

cmdcolin Oct 23, 2014

Contributor

Thank you very much for providing the data. I think this is actually similar to another recently reported issue

#516

The following additional patch seems to restore functionality for this case

diff --git a/src/JBrowse/Store/SeqFeature/_MismatchesMixin.js b/src/JBrowse/Store/SeqFeature/_MismatchesMixin.js
index b59e54e..6628cb5 100644
--- a/src/JBrowse/Store/SeqFeature/_MismatchesMixin.js
+++ b/src/JBrowse/Store/SeqFeature/_MismatchesMixin.js
@@ -129,7 +129,7 @@ return declare( null, {
                 if( op == 'S' || op == 'I' ) { 
                     templateOffset += len;
                 }   
-                else if( op == 'D' || op == 'N' || op == 'P' ) { 
+                else if( op == 'D' || op == 'P' ) { 
                     refOffset += len;
                 }   
                 else {

I think it may be wise to go over this code more closely to verify the logic in all cases now!

screenshot-localhost 2014-10-22 21-35-51

Contributor

cmdcolin commented Oct 23, 2014

Thank you very much for providing the data. I think this is actually similar to another recently reported issue

#516

The following additional patch seems to restore functionality for this case

diff --git a/src/JBrowse/Store/SeqFeature/_MismatchesMixin.js b/src/JBrowse/Store/SeqFeature/_MismatchesMixin.js
index b59e54e..6628cb5 100644
--- a/src/JBrowse/Store/SeqFeature/_MismatchesMixin.js
+++ b/src/JBrowse/Store/SeqFeature/_MismatchesMixin.js
@@ -129,7 +129,7 @@ return declare( null, {
                 if( op == 'S' || op == 'I' ) { 
                     templateOffset += len;
                 }   
-                else if( op == 'D' || op == 'N' || op == 'P' ) { 
+                else if( op == 'D' || op == 'P' ) { 
                     refOffset += len;
                 }   
                 else {

I think it may be wise to go over this code more closely to verify the logic in all cases now!

screenshot-localhost 2014-10-22 21-35-51

@09140930

This comment has been minimized.

Show comment
Hide comment
@09140930

09140930 Oct 23, 2014

Excellent, thanks for looking into this Colin!

Looking forward to the next release :)

Excellent, thanks for looking into this Colin!

Looking forward to the next release :)

@cmdcolin cmdcolin added this to the 1.11.6 milestone Jan 23, 2015

@cmdcolin cmdcolin closed this in 2734ae1 Jan 26, 2015

@cmdcolin

This comment has been minimized.

Show comment
Hide comment
@cmdcolin

cmdcolin Jan 26, 2015

Contributor

I would like to add unit testing for this eventually

Contributor

cmdcolin commented Jan 26, 2015

I would like to add unit testing for this eventually

cmdcolin added a commit that referenced this issue Feb 9, 2015

cmdcolin added a commit that referenced this issue Feb 10, 2015

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment