Naive sequence reconstruction is introducing a stop codon

[matsen]

Naive sequence reconstruction is introducing a stop codon in a clonal family that originally had no stops.
I'm guessing that the reconstruction doesn't know about frame.
No big deal if you don't want to address this, but it was a surprise for me.

Issue Summary

File: /fh/fast/matsen_e/shared/bcr-mut-sel/working/data/loris/jaffe-paired/jaffe_donor2_memory.csv.gz
Family: 100150-IGK-100150 (donor 2, single-member family)
Problem: Naive reconstruction introduces TAG stop codon at nucleotide positions 303-305

Details

The original memory sequence translates cleanly with no stop codons:

Memory sequence CDR3: ...CARDSVQWDRLIGYFQHW...

But the reconstructed naive sequence contains a premature stop:

Naive sequence CDR3:  ...CARDSV*WELLLGYFQHW...  (stop codon at position 101)

The stops_heavy: False field in the output correctly reflects the memory sequence, but the naive sequence itself contains the TAG stop codon.

Full CSV row

d2,100150-IGK-100150,AGATTGCAGCCCAGCT-1-1279054_contig_h,CAGGTGCACCTGGTGGAGTCTGGGGGAGGCGTGGTCCAGCCTGGGAGGTCCCTGAGACTCTCTTGTGCAGCCTCTGGATTCACCATAAGTACTTATGCTATGCACTGGGTCCGCCAGGCTCCAGGCAAGGGGCTGGAGTGGGTGGCAGTTACATTATATGATGGAAGCAATACATACTATGGAGATTCCGTGAAGGGCCGATTCACCATCTCCAGAGACAATTCCAAGAACACCCTGTATCTGCAAATGAACGGCCTGAGAGCTGAGGACACGGCTCTATATTATTGTGCGAGAGATAGCGTACAGTGGGACCGACTAATAGGCTACTTCCAGCACTGGGGCCAGGGCACCCTGGTCACCGTCTCCTCA,CAGGTGCACCTGGTGGAGTCTGGGGGAGGCGTGGTCCAGCCTGGGAGGTCCCTGAGACTCTCTTGTGCAGCCTCTGGATTCACCATAAGTACTTATGCTATGCACTGGGTCCGCCAGGCTCCAGGCAAGGGGCTGGAGTGGGTGGCAGTTACATTATATGATGGAAGCAATACATACTATGGAGATTCCGTGAAGGGCCGATTCACCATCTCCAGAGACAATTCCAAGAACACCCTGTATCTGCAAATGAACGGCCTGAGAGCTGAGGACACGGCTCTATATTATTGTGCGAGAGATAGCGTACAGTGGGACCGACTAATAGGCTACTTCCAGCACTGGGGCCAGGGCACCCTGGTCACCGTCTCCTCA,CAGGTGCAGCTGGTGGAGTCTGGGGGAGGCGTGGTCCAGCCTGGGAGGTCCCTGAGACTCTCCTGTGCAGCCTCTGGATTCACCTTCAGTAGCTATGGCATGCACTGGGTCCGCCAGGCTCCAGGCAAGGGGCTGGAGTGGGTGGCAGTTATATCATATGATGGAAGTAATAAATACTATGCAGACTCCGTGAAGGGCCGATTCACCATCTCCAGAGACAATTCCAAGAACACGCTGTATCTGCAAATGAACAGCCTGAGAGCTGAGGACACGGCTGTGTATTACTGTGCGAGAGATAGCGTATAGTGGGAGCTACTACTAGGCTACTTCCAGCACTGGGGCCAGGGCACCCTGGTCACCGTCTCCTCA,IGH,IGHV3-30*03,IGHJ1*01,CARDSVQWDRLIGYFQHW,True,False,False,0,75,96,150,171,288,333,75,87,156,165,297,330,23,AGATTGCAGCCCAGCT-1-1279054_contig_l,AAAATTGTGTTGACACAGTCTCCAGCCACCCTGTCTTTGTCTCCAGGGGAAAGAGCCACCCTCTCCTGCAGGGCCAGTCAGAGTGTTAGCAGCTACTTAGCCTGGTACCAACAGAAGCCTGGCCAGCCTCCCAGGCTCCTCATCTATGATGCATCCTACAGGGCCACTGGCATCCCAGCCAGGTTCAGTGGCAGTGGGTCTGGGACAGACTTCACTCTCACCATCAGCAGCCTAGAGCCTGAAGATTTTGCAGTTTATTACTGTCAGCAGCGTAGCAGCTGGCCTCTCACTTTCGGCGGAGGGACCAAGGTGGAGATCAAA,AAAATTGTGTTGACACAGTCTCCAGCCACCCTGTCTTTGTCTCCAGGGGAAAGAGCCACCCTGTCAGGGCCAGTCAGAGTGTTAGCAGCTACTTAGCCTGGTACCAACAGAAGCCTGGCCAGCCTCCCAGGCTCCTCATCTATGATGCATCCTACAGGGCCACTGGCATCCCAGCCAGGTTCAGTGGCAGTGGGTCTGGGACAGACTTCACTCTCACCATCAGCAGCCTAGAGCCTGAAGATTTTGCAGTTTATTACTGTCAGCAGCGTAGCAGCTGGCCTCTCACTTTCGGCGGAGGGACCAAGGTGGAGATCAAA,GAAATTGTGTTGACACAGTCTCCAGCCACCCTGTCTTTGTCTCCAGGGGAAAGAGCCACCCTCTCCTGCAGGGCCAGTCAGAGTGTTAGCAGCTACTTAGCCTGGTACCAACAGAAACCTGGCCAGGCTCCCAGGCTCCTCATCTATGATGCATCCAACAGGGCCACTGGCATCCCAGCCAGGTTCAGTGGCAGTGGGTCTGGGACAGACTTCACTCTCACCATCAGCAGCCTAGAGCCTGAAGATTTTGCAGTTTATTACTGTCAGCAGCGTAGCAACTGGCCTCTCACTTTCGGCGGAGGGACCAAGGTGGAGATCAAA,IGK,IGKV3-11*01,IGKJ4*01,CQQRSSWPLTF,True,False,False,0,78,93,147,153,264,288,75,93,147,132,270,285,5.0

[matsen]

Here's how often this happened in the jaffe data:

• Donor 1: 674 families filtered out (from 69,922 → 69,248 families)
• Donor 2: 303 families filtered out (from 25,506 → 25,203 families)
• Donor 3: 338 families filtered out (from 33,730 → 33,392 families)
• Donor 4: 464 families filtered out (from 38,401 → 37,937 families)

[psathyrella]

Yes, this is (unfortunately) expected behavior. In the past (Overbaugh and Goo data I think?) we've just ignored those alternative naive sequences with stop codons, or have arbitrarily changed a base so it wasn't a stop.

I'm not sure how we could make the HMM aware of stop codons -- it's calculating probabilities for V/D/J gene sequences and insertion and deletion lengths, and only arrives at the full sequence (which you need to know if there's a stop codon) as a final output of the viterbi calculation. So the bit of the HMM that's calculating which V gene + V deletion is most likely only knows the length it's working with (recall the HMM "factorization" scheme).

I'm not sure that we can use the "alternative annotations" as a general method, since they're only meaningful for larger families, and the per-gene support framework also, to me, doesn't seem useful here.

Which, I think, means we'd have to come up with a post-hoc scheme for modifying naive sequences with stop codons. Perhaps the easiest would be reporting all naive sequences separate by single nucleotide differences from the "most likely" one that has a stop codon. Or perhaps we could rank these alternatives by their hamming distance to all sequences in the family.