SAM format as input/output

[kclem]

It looks like the core implementation of process_bam should easily be applied to sam files too. The CIGAR string and alignment position need to be modified (implemented properly).

Helpful links:
Explain sam flags:
https://broadinstitute.github.io/picard/explain-flags.html

IVG viewer:
https://software.broadinstitute.org/software/igv/

[kclem]

OK, I pushed some demo code to 2f7d321. If you change to that directory ('CRISPResso2/sam_test') you should be able to run python sam_test.py.

Note that sam_test.py includes a function process_bam that mirrors the code in CRISPRessoCORE. The goal is that once this is working, we can plop process_bam from sam_test.py back into CRISPRessoCORE, so be aware that any changes you make to the functionality or process_bam signature, etc. mean we have to make changes to other parts in CRISPRessoCORE.

if your run python sam_test.py you'll see demo output from a bam file with 50 reads.

I just printed out a couple of things including

• the old cigar string for when this read was aligned to the genome,

old cigar was 164M34D76M4I4M

• The sam line elements (separated by tabs so you can index into them as in the sam specification. Note that CRISPResso adds on the "c2:Z" tags at the end. The rest are from the input sam file.

line els: ['M05634:17:000000000-BKJP5:1:1101:23018:4446', '16', 'chr11', '992', '255', '164M34D76M4I4M', '*', '0', '0', 'CTTACCTCCAGTGTTTGTGCAGCCGCCGCTCCAGAGCCGTGCGAATGGGGCCATGCCGACCAAAGCGCCGATGGATGTGGCGCAGGTAGCGCGCCCACTGCAAGGCCCGGCGCACGGTGGCGGGGTCCCAGGTGCTGACGTAGGTAGTGCTTGAGACCGCCAGAATTCCATGAGGTGCGCGAAGGCCCTACTTCCGCTTTCACCTTGGAGACGGCGACTCTCTGCGTACTGATTGGAACATCCGTCAC', ';9/..9;.9/.A9.FFFF.;-999=>-F99/FBD=A;CBDFFBB9.-:/A9;;B=AFFFF9.D;ABA?:/BFFBF9AA.;@A/FB9B@@DCB.A9F00BC9/GCGD@?DGDC.GGHGCFCGGFDGG/HG0?0EHGGHHGGHFHGBHHGHG2GG?E<GGEGHGFHHHHGHGHHFFCGGGGGFEFEBAGFFBGCGCGHHGE5GEGHHFDFGGGGGFCHHDG2EGGGGHGFDGGGGGGGFFBBBDABAAAB', 'AS:i:-143', 'XN:i:0', 'XM:i:5', 'XO:i:2', 'XG:i:38', 'NM:i:43', 'MD:Z:3G6C2G150^AGCTCGGAAAAGCGATCCAGGTGCTGCAGAAGGG78C0G0', 'YT:Z:UU', 'c2:Z:ALN=Reference CLASS=Reference_MODIFIED MODS=D34;I0;S0 DEL=80(34) INS= SUB= ALN_REF=----CGGATGTTCCAATCAGTACGCAGAGAGTCGCCGTCTCCAAGGTGAAAGCGGAAGTAGGGCCTTCGCGCACCTCATGGAATCCCTTCTGCAGCACCTGGATCGCTTTTCCGAGCTTCTGGCGGTCTCAAGCACTACCTACGTCAGCACCTGGGACCCCGCCACCGTGCGCCGGGCCTTGCAGTGGGCGCGCTACCTGCGCCACATCCATCGGCGCTTTGGTCGG------------------------------------------------------- ALN_SEQ=GTGACGGATGTTCCAATCAGTACGCAGAGAGTCGCCGTCTCCAAGGTGAAAGCGGAAGTAGGGCCTTCGCGCACCTCATGGAAT----------------------------------TCTGGCGGTCTCAAGCACTACCTACGTCAGCACCTGGGACCCCGCCACCGTGCGCCGGGCCTTGCAGTGGGCGCGCTACCTGCGCCACATCCATCGGCGCTTTGGTCGGCATGGCCCCATTCGCACGGCTCTGGAGCGGCGGCTGCACAAACACTGGAGGTAAG']

• the new cigar string when the read is aligned to this amplicon sequence

potentially setting cigar to : 4I80M34D109M55I from alignment

• The alignment between the read (from the input bam) and the reference amplicon (provided by the user):

aln: GTGACGGATGTTCCAATCAGTACGCAGAGAGTCGCCGTCTCCAAGGTGAAAGCGGAAGTAGGGCCTTCGCGCACCTCATGGAAT----------------------------------TCTGGCGGTCTCAAGCACTACCTACGTCAGCACCTGGGACCCCGCCACCGTGCGCCGGGCCTTGCAGTGGGCGCGCTACCTGCGCCACATCCATCGGCGCTTTGGTCGGCATGGCCCCATTCGCACGGCTCTGGAGCGGCGGCTGCACAAACACTGGAGGTAAG
ref: ----CGGATGTTCCAATCAGTACGCAGAGAGTCGCCGTCTCCAAGGTGAAAGCGGAAGTAGGGCCTTCGCGCACCTCATGGAATCCCTTCTGCAGCACCTGGATCGCTTTTCCGAGCTTCTGGCGGTCTCAAGCACTACCTACGTCAGCACCTGGGACCCCGCCACCGTGCGCCGGGCCTTGCAGTGGGCGCGCTACCTGCGCCACATCCATCGGCGCTTTGGTCGG-------------------------------------------------------

• and the alignment position reported by the alignment. These positions were produced by the short read aligner bowtie2 and are all in reference to a very short genome at smallGenome/smallGenome.fa.

alignment position from aligner is 255

[kclem]

You should be able to view the input and output bam files in IGV.

The goal is to have the CRISPResso alignment reflected in the output bam and viewable in the IGV. So for example, many of these reads have a 4bp insertion at the beginning (4I at the beginning of the new CIGAR) and we'd like this to show up in the IGV as well, while modifying the start position to be accurate.

I've created code to create the CIGAR string, but you will have to change the "MD:Z" tag - which the IGV may use for display. You can read about this in the SAM spec.

[Colelyman]

Something I stumbled on that we should be aware of, apparently there is a design flaw in BAM format that CIGAR strings will break when they have >65535 operations (source: https://github.com/lh3/minimap2#working-with-65535-cigar-operations). This will obviously only affect us when we start using long reads, but it wouldn't be hard to do what minimap2 does and have a clipped CIGAR and add the full one to the CG tag.