Using partis to trim NGS reads

So as it turns out, many of the NGS reads I have worked with so far are not consistently trimmed causing various problems. An example of what I frequently see is a bunch of sequences that vary in length but all have the same middle part i.e. they are the same sequence but just with different, inconsistent trimming. Currently these badly trimmed reads will be interpreted as different bcr's which is obviously wrong and adding these as extra sequences to partitioning just makes it slower. Other downstream problems also arrises e.g. when read frequency is calculated.

Since partis is doing all the alignment of germline genes it should also enable to trim reads based on anchor points shared by all bcr sequences e.g. position in IGV and conserved cysteine. Having those anchors enables consistent trimming as illustrated below: screen shot 2016-11-18 at 20 20 51

And as Duncan has illustrated it, partis could fix these badly trimmed reads by N-padding them to the same length. Here illustrated with the middle two sequences as badly trimmed and N-padded to equal length: img_20161118_111131_992

If the N-padding functions would be moved to utils it might be a lot easier to access, and then I could try to cook the trimming together.

I started to move them to utils, but realized it probably makes more sense to leave them in waterer.py, but add an option to trim and pad before writing the sw cache file (rather than after, which it has to be for technical reasons, in general to be a valid sw cache file that waterer.py can read later). But you can just get what you need from this faux-sw cache file. So if you run this:

./bin/partis cache-parameters --infname simu.csv --only-smith-waterman --sw-cachefname sw-cache.csv --write-trimmed-and-padded-seqs-to-sw-cachefname --debug 1

with this in simu.csv:

unique_id,reco_id,v_gene,d_gene,j_gene,v_5p_del,v_3p_del,d_5p_del,d_3p_del,j_5p_del,j_3p_del,fv_insertion,vd_insertion,dj_insertion,jf_insertion,cdr3_length,seq,indelfo
-2305064032213891805,-446504016578225677,IGHV3-48*03,IGHD6-6*01,IGHJ3*02,0,0,7,0,0,0,,TCACAGATGCACT,ATGGAT,,60,TTTTTTTTGAGGTGCAGCTGCTGGAGTCTGGGGGAGGCTTGGTACAGCCTGGGCGGTCCCTGAGACTCTCCTGTGCAGCCTCTAGATTCGCCTTCAGTAGTTATGAAATGAACTGGGTCGGCCAGGCTCCAGGGAAGGGGCTGGAGTGGGTTTCATACATTAATAATAGTGGTAGTCCCATATACTACCCAGACTCTGTGAAGGGCCGATTCACCATCTCCAGAGACAACGCCAAGAACTCACTGTATCTGCAAATGAACAGCCTGAGAGCCGAGGACACGGCTGTTTATTACTGTGCGAGAGATCACAGATGCACTGCAGCTCGTCCATAGATTGACGCTTTTGATATCTGGGGCCAAGGGACAATGGTCACCGTCTCTTCAG,"{'reversed_seq': '', 'indels': []}"
-xx2305064032213891805,-446504016578225677,IGHV3-48*03,IGHD6-6*01,IGHJ3*02,0,0,7,0,0,0,,TCACAGATGCACT,ATGGAT,,60,TTTTTTGAGGTGCAGCTGCTGGAGTCTGGGGGAGGCTTGGTACAGCCTGGGCGGTCCCTGAGACTCTCCTGTGCAGCCTCTAGATTCGCCTTCAGTAGTTATGAAATGAACTGGGTCGGCCAGGCTCCAGGGAAGGGGCTGGAGTGGGTTTCATACATTAATAATAGTGGTAGTCCCATATACTACCCAGACTCTGTGAAGGGCCGATTCACCATCTCCAGAGACAACGCCAAGAACTCACTGTATCTGCAAATGAACAGCCTGAGAGCCGAGGACACGGCTGTTTATTACTGTGCGAGAGATCACAGATGCACTGCAGCTCGTCCATAGATTGACGCTTTTGATATCTGGGGCCAAGGGACAATGGTCACCGTCTCTTCAG,"{'reversed_seq': '', 'indels': []}"
3469894191945479084,-446504016578225677,IGHV3-48*03,IGHD6-6*01,IGHJ3*02,0,0,7,0,0,0,,TCACAGATGCACT,ATGGAT,,60,GCTGGTGGAGTCTGGGGGAGGCTTGGTATAGCCTGGGGGGTCCCTGAGACTCTCCTGTGCAGCCTCTGGATTCACCTTCAGTAGTTATGAAATGAACTGGGTCCGCCAGGCTCCAGGGAAGGGGCTGGAGTGGCTTTCATACATTAGTAATAGTAGTAGTCCCTTATGCTACGCAGACTCTGTGAAGGGCCGATTCACCATCTCCAAAGACAACGCCAAGAACTCACTGTATCTGCAAATGAACAGCCTGAGAGCCGAGGACACGGCTGTTTATTACTGTGCGAGAGATCTCACATGCACTGCAGCGCGTCCATGGATTGTCGCTTTTGATATCTGGGGCCAAGGGACAATGGTCACCGTCTCTTCAG,"{'reversed_seq': '', 'indels': []}"
-1682815703800499414,-446504016578225677,IGHV3-48*03,IGHD6-6*01,IGHJ3*02,0,0,7,0,0,0,,TCACAGATGCACT,ATGGAT,,60,GAGGTGCAGCTGGTAGAGTCTGGGGGAGGCTTGGTATAGCCTGGAGGGTCCCTGAGACTCTCCTGTGCAGCCTCTGGATTCACCTTCAGTAGTTATGAAATGAACTGGGTCCGCCAAGCTCCAGGGAAGGGGCTGGAGTGGGTTTCATACATTAGTAGTAGTGGTAGTACCATATACTACGCAGACTCTGTGAACGGCCGATTCACCATCTCCAGAGACAACACCAAGAACTCACTGTATCTGCAAATGAACACCCTGAGAGCCGAGGACACGGCTGTTTATTACTGTGCGAGAGATCACAGATGCACTGCAGCTCGTCCATCGATTGATGCTTTTGATACCTGGGGCCAAGGGACAATGGTCACCGTCTCTTCAG,"{'reversed_seq': '', 'indels': []}"

you should get the following trimmed and padded seqs in sw-cache.csv:

seqs                                    
GAGGTGCAGCTGCTGGAGTCTGGGGGAGGCTTGGTACAGCCTGGGCGGTCCCTGAGACTCTCCTGTGCAGCCTCTAGATTCGCCTTCAGTAGTTATGAAATGAACTGGGTCGGCCAGGCTCCAGGGAAGGGGCTGGAGTGGGTTTCATACATTAATAATAGTGGTAGTCCCATATACTACCCAGACTCTGTGAAGGGCCGATTCACCATCTCCAGAGACAACGCCAAGAACTCACTGTATCTGCAAATGAACAGCCTGAGAGCCGAGGACACGGCTGTTTATTACTGTGCGAGAGATCACAGATGCACTGCAGCTCGTCCATAGATTGACGCTTTTGATATCTGGGGCCAAGGGACAATGGTCACCGTCTCTTCAG
GAGGTGCAGCTGGTAGAGTCTGGGGGAGGCTTGGTATAGCCTGGAGGGTCCCTGAGACTCTCCTGTGCAGCCTCTGGATTCACCTTCAGTAGTTATGAAATGAACTGGGTCCGCCAAGCTCCAGGGAAGGGGCTGGAGTGGGTTTCATACATTAGTAGTAGTGGTAGTACCATATACTACGCAGACTCTGTGAACGGCCGATTCACCATCTCCAGAGACAACACCAAGAACTCACTGTATCTGCAAATGAACACCCTGAGAGCCGAGGACACGGCTGTTTATTACTGTGCGAGAGATCACAGATGCACTGCAGCTCGTCCATCGATTGATGCTTTTGATACCTGGGGCCAAGGGACAATGGTCACCGTCTCTTCAG
NNNNNNNNGCTGGTGGAGTCTGGGGGAGGCTTGGTATAGCCTGGGGGGTCCCTGAGACTCTCCTGTGCAGCCTCTGGATTCACCTTCAGTAGTTATGAAATGAACTGGGTCCGCCAGGCTCCAGGGAAGGGGCTGGAGTGGCTTTCATACATTAGTAATAGTAGTAGTCCCTTATGCTACGCAGACTCTGTGAAGGGCCGATTCACCATCTCCAAAGACAACGCCAAGAACTCACTGTATCTGCAAATGAACAGCCTGAGAGCCGAGGACACGGCTGTTTATTACTGTGCGAGAGATCTCACATGCACTGCAGCGCGTCCATGGATTGTCGCTTTTGATATCTGGGGCCAAGGGACAATGGTCACCGTCTCTTCAG

https://github.com/psathyrella/partis/commit/496fb886445e412eaba0c4ac568739950dd79e48

upon further discussion, what we want is basically these examples to all be collapsed to one sequence (well, and the crappy one removed, which is taken care of with --skip-unproductive):

>too_long_left_side
GGCCCAGGACTGGTGAAGCCTTCACAGACCCTGTCCCTCACCTGCACTGTCTCTGGTGGCTCCATCAGCAGTGGTGGTTACTACTGGAGCTGGATCCGCCAGCACCCAGGGAAGGGCCTGGAGTGGATTGGGTACATCTATTACAGTGGGAGCACCTACTACAACCCGTCCCTCAAGAGTCGAGTTACCATATCAGTAGACAAGTCCAAGAACCAGTTCTCCCTGAAGCTGAGCTCTGTGACCGCCGCGGACACGGCCGTGTATTACTGTGCGAGAGATGAGGAGGGGATCGGATTTGACTACTGGGGC
>too_long_right_side
CCAGGACTGGTGAAGCCTTCACAGACCCTGTCCCTCACCTGCACTGTCTCTGGTGGCTCCATCAGCAGTGGTGGTTACTACTGGAGCTGGATCCGCCAGCACCCAGGGAAGGGCCTGGAGTGGATTGGGTACATCTATTACAGTGGGAGCACCTACTACAACCCGTCCCTCAAGAGTCGAGTTACCATATCAGTAGACAAGTCCAAGAACCAGTTCTCCCTGAAGCTGAGCTCTGTGACCGCCGCGGACACGGCCGTGTATTACTGTGCGAGAGATGAGGAGGGGATCGGATTTGACTACTGGGGCCAA
>too_short_right_side_discard
CCAGGACTGGTGAAGCCTTCACAGACCCTGTCCCTCACCTGCACTGTCTCTGGTGGCTCCATCAGCAGTGGTGGTTACTACTGGAGCTGGATCCGCCAGCACCCAGGGAAGGGCCTGGAGTGGATTGGGTACATCTATTACAGTGGGAGCACCTACTACAACCCGTCCCTCAAGAGTCGAGTTACCATATCAGTAGACAAGTCCAAGAACCAGTTCTCCCTGAAGCTGAGCTCTGTGACCGCCGCGGACACGGCCGTGTATTACTGTGCGAGAGATGAGGAGGGGATCGGATTTGA
>good_length1
CCAGGACTGGTGAAGCCTTCACAGACCCTGTCCCTCACCTGCACTGTCTCTGGTGGCTCCATCAGCAGTGGTGGTTACTACTGGAGCTGGATCCGCCAGCACCCAGGGAAGGGCCTGGAGTGGATTGGGTACATCTATTACAGTGGGAGCACCTACTACAACCCGTCCCTCAAGAGTCGAGTTACCATATCAGTAGACAAGTCCAAGAACCAGTTCTCCCTGAAGCTGAGCTCTGTGACCGCCGCGGACACGGCCGTGTATTACTGTGCGAGAGATGAGGAGGGGATCGGATTTGACTACTGGGGC
>good_length2
CCAGGACTGGTGAAGCCTTCACAGACCCTGTCCCTCACCTGCACTGTCTCTGGTGGCTCCATCAGCAGTGGTGGTTACTACTGGAGCTGGATCCGCCAGCACCCAGGGAAGGGCCTGGAGTGGATTGGGTACATCTATTACAGTGGGAGCACCTACTACAACCCGTCCCTCAAGAGTCGAGTTACCATATCAGTAGACAAGTCCAAGAACCAGTTCTCCCTGAAGCTGAGCTCTGTGACCGCCGCGGACACGGCCGTGTATTACTGTGCGAGAGATGAGGAGGGGATCGGATTTGACTACTGGGGC
>good_length3
CCAGGACTGGTGAAGCCTTCACAGACCCTGTCCCTCACCTGCACTGTCTCTGGTGGCTCCATCAGCAGTGGTGGTTACTACTGGAGCTGGATCCGCCAGCACCCAGGGAAGGGCCTGGAGTGGATTGGGTACATCTATTACAGTGGGAGCACCTACTACAACCCGTCCCTCAAGAGTCGAGTTACCATATCAGTAGACAAGTCCAAGAACCAGTTCTCCCTGAAGCTGAGCTCTGTGACCGCCGCGGACACGGCCGTGTATTACTGTGCGAGAGATGAGGAGGGGATCGGATTTGACTACTGGGGC

I think this is now accomplished (https://github.com/psathyrella/partis/commit/6c152e0ca322765b66f154fbf0493b0f4b5671c4):

at least, running this: ./bin/partis cache-parameters --infname trimming_example.fa --debug 1 --skip-unproductive --also-remove-duplicate-sequences-with-different-lengths

results in I just getting good_length2, which I think is what we want?

I would've turned on --also-remove-duplicate-sequences-with-different-lengths but I think it'll be too slow, and in most data sets I don't think there's enough sequences that are substrings of each other to make it worth while.

Thank you Duncan, this is totally what I wanted, with the only addition that I would like to have a count on how many that have been collapsed so the read frequencies are not compromised.

ah, that's great.

The number of duplicates should be showing up in the 'duplicates' column, which is fairly recent. Let me know if that's not working, though.

psathyrella / partis

Using partis to trim NGS reads #224