jkbonfield
commented
7 years ago I just had an idea - to use cram_dump -v to dump out the differences in Java produced CRAM files for ref ending AAA vs NNN. This explains why so many bases, aligned in different locations, change. The base substitution statistics change which in turn generates different substitution matrices. The containers in the two cram files show:
Substitution map:
A: CGTN
C: AGTN
G: ACTN
T: ACGN
N: TGCAvs
Substitution map:
A: TNGC
C: AGTN
G: ACTN
T: ACGN
N: ACGTThe main problem though is it is generating base substitutions that are impossible. CRAM BS data series is used for mapping ACGTN to ACGTN subst, as a number to lookup via the table in the header. (This is because we know we cannot ever generate and A->A substitution, thus the numeric solution has lower entropy and the table permits it to be optimised.) However this cannot work with ambiguous bases in either reference or sequence.
Oddly it does this for matches too. So rec 1/30 (aligned correctly, all matches ref, albeit M vs M, R vs R, etc):
FN = 10 (ret 0, out_sz 1)
0: FC = X (ret 0, out_sz 1)
0: FP = 5+0 = 5 (ret 0, out_sz 1)
0: BS = 0 (ret 0)
1: FC = X (ret 0, out_sz 1)
1: FP = 1+5 = 6 (ret 0, out_sz 1)
1: BS = 0 (ret 0)
2: FC = X (ret 0, out_sz 1)
2: FP = 1+6 = 7 (ret 0, out_sz 1)
2: BS = 0 (ret 0)
3: FC = X (ret 0, out_sz 1)
3: FP = 1+7 = 8 (ret 0, out_sz 1)
3: BS = 0 (ret 0)
4: FC = X (ret 0, out_sz 1)
4: FP = 1+8 = 9 (ret 0, out_sz 1)
4: BS = 0 (ret 0)
5: FC = X (ret 0, out_sz 1)
5: FP = 1+9 = 10 (ret 0, out_sz 1)
5: BS = 0 (ret 0)
6: FC = X (ret 0, out_sz 1)
6: FP = 1+10 = 11 (ret 0, out_sz 1)
6: BS = 0 (ret 0)
7: FC = X (ret 0, out_sz 1)
7: FP = 1+11 = 12 (ret 0, out_sz 1)
7: BS = 0 (ret 0)
8: FC = X (ret 0, out_sz 1)
8: FP = 1+12 = 13 (ret 0, out_sz 1)
8: BS = 0 (ret 0)
9: FC = X (ret 0, out_sz 1)
9: FP = 1+13 = 14 (ret 0, out_sz 1)
9: BS = 0 (ret 0)So the first 4 bases (ACGT) matched and didn't generate a feature code, but all subsequent ones produce FC 'X' (substitution) and BS 0 (nonsensical).
Record 2/30, which is slipped by one and hence 100% mismatch, also uses purely FC X and BS 0 (except for the early ACG vs CGT alignment). Samtools/Scramble produced this output instead for rec 1 and 2:
Rec 1/30 at 0,7
FN = 0 (ret 0, out_sz 1)
[Ie entirely matches]
Rec 2/30 at 0,7
FN = 15 (ret 0, out_sz 1)
0: FC = X (ret 0, out_sz 1)
0: FP = 1+0 = 1 (ret 0, out_sz 1)
0: BS = 0 (ret 0)
1: FC = X (ret 0, out_sz 1)
1: FP = 1+1 = 2 (ret 0, out_sz 1)
1: BS = 1 (ret 0)
2: FC = X (ret 0, out_sz 1)
2: FP = 1+2 = 3 (ret 0, out_sz 1)
2: BS = 2 (ret 0)
3: FC = X (ret 0, out_sz 1)
3: FP = 1+3 = 4 (ret 0, out_sz 1)
3: BS = 3 (ret 0)
4: FC = B (ret 0, out_sz 1)
4: FP = 1+4 = 5 (ret 0, out_sz 1)
4: BA/QS(B) = R/-1 (ret 0)
5: FC = B (ret 0, out_sz 1)
5: FP = 1+5 = 6 (ret 0, out_sz 1)
5: BA/QS(B) = Y/-1 (ret 0)
...Ie it uses feature 'X' for the initial ACGT base mismatches and feature 'B' for the rest to store verbatim the corrected sequence. I'm wondering why I chose a series of B instead of a single 'b' string, but maybe that's just historic as 'b' didn't arrive until much later.
This explains why Cramtools decode of the samtools generated CRAM gets the correct sequence for the matching reads (they have no features at all), and why the ACGT bases are correct for the mismatching ones (via valid 'X' features). However it doesn't really explain why all the verbatim 'B' features were replaced with Ns instead.
Admittedly ambiguous bases in the records is pretty rare.
cmnbroad
vadimzalunin
Verify
I'm using cramtools to drive the htsjdk cram code. I see no warnings present from cramtools.
Subject of the issue
Ambiguous bases in sequence and/or reference decode incorrectly when going via CRAM round-trips. Oddly it does this in a baffling way.
My basic example is this as the fasta file:
My SAM file is a sequence of ACGTRYMWSKDHVBN aligned to coordinates 1 to 30 inclusive. Note that this means some alignments are invalid, going beyond the end of the reference, but this itself is an interesting test case which ought to pass or to emit a warning if it cannot be stored. Hence the SAM file is:
Your environment
version of htsjdk The latest checkout that works with cramtools (which is itself the git master). For htsjdk this means 88b6719 as the following commit removes functions used by cramtools. (I assume a fix for cramtools is incoming.)
version of java jdk1.8.0_74
which OS Ubuntu 12.04
Steps to reproduce
The complete round-trip commands are:
Expected behaviour
Expected behaviour is obviously to reproduce the SAM input listed above.
Actual behaviour
Ok so this is where it gets weird...
With ref seq x of ACGTRYMWSKDHVBNacgtrymwskdhvbnAAA as listed above, my SAM->CRAM->SAM round trip ends up as:
Most of the ambiguous bases turned to A apart from a couple which became T. In the latter half, where we know we run off the end of the reference, we see a stripe of "TTT" occuring.
If I edit the last 3 reference bases from AAA to TTT instead the round-trip generates this:
It's obvious that the TTTs in the output became GGGs, but also those isolated Ts earlier became either Gs or Cs too. Think about that for a moment. By changing ref seq bases 30 to 32 the round-trip changed the bases for sequences aligned to base 3..17 (the last T turned into a C).
More "exciting" still is to turn the last 3 reference bases to NNN instead:
We're now really singing in the baffling changes! This is repeatable behaviour, so it's not due to ininitialised data. Samtools decoding the CRAM file shows the same affect (but not when samtools encodeds), which implies this is an issue in the Java CRAM encoder. Indeed Java can mostly decode the samtools generated files, but demonstrates a different bug:
Here Cramtools decode turned all the ambiguous bases to Ns. The bases are explicitly recorded in the CRAM file. Eg read 2 has (from io_lib's "cram_dump -v"):
So the data is there, but is being modified by the decoder somewhere. I think this is perhaps more than one bug.