[lankage]

Hi Brian,
Came across an instance today where I was attempting to parse bbmap SAM output for cigar string information to observe the mapping of some HIV read sequences against a reference. I was expecting for a cigar string to be present for each record when using the outm=output.sam parameter. In one of my mapping records, I observed an asterisk instead. Am I wrong to assume outm= is intended to include only those mapped reads? I can filter out these reads, but I'd like to make sure my mapping parameters make sense.

From the SAM output:
@HD VN:1.4 SO:unsorted
@SQ SN:NC_001802DRannotations_(modified) LN:9181
@PG ID:BBMap PN:BBMap VN:34.92 CL:java -Djava.library.path=/home/dnanexus/bbmap/jni/ -ea -Xmx10g align2.BBMap build=1 overwrite=true fastareadlen=500 in=reads_file ref=ref_file outm=sam_output minid=.8 strictmaxindel=10 k=8 subfilter=15 -Xmx10g

The offending mapping record:
M01472:214:000000000-AG0YC:1:2108:10755:20410 1:N:0:78 0 NC_001802DRannotations_(modified) 2154 4 * * 0 0 AAGTTAAACAATGGCCATTGACAGAAGAAAAAATAAAAGCATTAGTCATAGTAATATGGGGAAAGACTCCTAAATTTAAATTACCCATACAAAAGGAAACATGGGAAGCATGGTGGACAGAGTATTGGC CCCCCGGGGGGGGGGGGGGGGGGGGGGFDCFGECGGGF<AFEGGFEFGGGGFGGGGGGGGGGGGGGGGGFGFFGGFGEFGGAGFGEAFF<,FGGGGGGGGGGFGGFDGGGFECFGGGGGGGGGGCCCCC AM:i:4

This only occurred after parsing ~8 million mapping records out of a total 50 million.

[Brian Bushnell]

That should not happen... outm should only print mapped reads, except when processing paired data, in which case it will print unmapped reads if the mate is mapped. But according to the flag you are processing single-ended data and that read was mapped, so it should have a cigar string.

It's possible that this is a bug that I've fixed since v34.92, but I have no record of that in my changelog. So I'd like to try to replicate it to see if it's a current bug. I already have the read, but can you send me or direct me to the specific reference you are using?

Thanks,
Brian

[lankage]

>NC_001802DRannotations_(modified)
GGUCUCUCUGGUUAGACCAGAUCUGAGCCUGGGAGCUCUCUGGCUAACUAGGGAACCCACUGCUUAAGCCUCAAUAAAGCUUGCCUUGAGUGCUUCAAGUAGUGUGUGCCCGUCUGUUGUGUGACUCUGGUAACUAGAGAUCCCUCAGACCCUUUUAGUCAGUGUGGAAAAUCUCUAGCAGUGGCGCCCGAACAGGGACCUGAAAGCGAAAGGGAAACCAGAGGAGCUCUCUCGACGCAGGACUCGGCUUGCUGAAGCGCGCACGGCAAGAGGCGAGGGGCGGCGACUGGUGAGUACGCCAAAAAUUUUGACUAGCGGAGGCUAGAAGGAGAGAGAUGGGUGCGAGAGCGUCAGUAUUAAGCGGGGGAGAAUUAGAUCGAUGGGAAAAAAUUCGGUUAAGGCCAGGGGGAAAGAAAAAAUAUAAAUUAAAACAUAUAGUAUGGGCAAGCAGGGAGCUAGAACGAUUCGCAGUUAAUCCUGGCCUGUUAGAAACAUCAGAAGGCUGUAGACAAAUACUGGGACAGCUACAACCAUCCCUUCAGACAGGAUCAGAAGAACUUAGAUCAUUAUAUAAUACAGUAGCAACCCUCUAUUGUGUGCAUCAAAGGAUAGAGAUAAAAGACACCAAGGAAGCUUUAGACAAGAUAGAGGAAGAGCAAAACAAAAGUAAGAAAAAAGCACAGCAAGCAGCAGCUGACACAGGACACAGCAAUCAGGUCAGCCAAAAUUACCCUAUAGUGCAGAACAUCCAGGGGCAAAUGGUACAUCAGGCCAUAUCACCUAGAACUUUAAAUGCAUGGGUAAAAGUAGUAGAAGAGAAGGCUUUCAGCCCAGAAGUGAUACCCAUGUUUUCAGCAUUAUCAGAAGGAGCCACCCCACAAGAUUUAAACACCAUGCUAAACACAGUGGGGGGACAUCAAGCAGCCAUGCAAAUGUUAAAAGAGACCAUCAAUGAGGAAGCUGCAGAAUGGGAUAGAGUGCAUCCAGUGCAUGCAGGGCCUAUUGCACCAGGCCAGAUGAGAGAACCAAGGGGAAGUGACAUAGCAGGAACUACUAGUACCCUUCAGGAACAAAUAGGAUGGAUGACAAAUAAUCCACCUAUCCCAGUAGGAGAAAUUUAUAAAAGAUGGAUAAUCCUGGGAUUAAAUAAAAUAGUAAGAAUGUAUAGCCCUACCAGCAUUCUGGACAUAAGACAAGGACCAAAGGAACCCUUUAGAGACUAUGUAGACCGGUUCUAUAAAACUCUAAGAGCCGAGCAAGCUUCACAGGAGGUAAAAAAUUGGAUGACAGAAACCUUGUUGGUCCAAAAUGCGAACCCAGAUUGUAAGACUAUUUUAAAAGCAUUGGGACCAGCGGCUACACUAGAAGAAAUGAUGACAGCAUGUCAGGGAGUAGGAGGACCCGGCCAUAAGGCAAGAGUUUUGGCUGAAGCAAUGAGCCAAGUAACAAAUUCAGCUACCAUAAUGAUGCAGAGAGGCAAUUUUAGGAACCAAAGAAAGAUUGUUAAGUGUUUCAAUUGUGGCAAAGAAGGGCACACAGCCAGAAAUUGCAGGGCCCCUAGGAAAAAGGGCUGUUGGAAAUGUGGAAAGGAAGGACACCAAAUGAAAGAUUGUACUGAGAGACAGGCUAAUUUUUUAGGGAAGAUCUGGCCUUCCUACAAGGGAAGGCCAGGGAAUUUUCUUCAGAGCAGACCAGAGCCAACAGCCCCACCAGAAGAGAGCUUCAGGUCUGGGGUAGAGACAACAACUCCCCCUCAGAAGCAGGAGCCGAUAGACAAGGAACUGUAUCCUUUAACUUCCCUCAGGUCACUCUUUGGCAACGACCCCUCGUCACAAUAAAGAUAGGGGGGCAACUAAAGGAAGCUCUAUUAGAUACAGGAGCAGAUGAUACAGUAUUAGAAGAAAUGAGUUUGCCAGGAAGAUGGAAACCAAAAAUGAUAGGGGGAAUUGGAGGUUUUAUCAAAGUAAGACAGUAUGAUCAGAUACUCAUAGAAAUCUGUGGACAUAAAGCUAUAGGUACAGUAUUAGUAGGACCUACACCUGUCAACAUAAUUGGAAGAAAUCUGUUGACUCAGAUUGGUUGCACUUUAAAUUUUCCCAUUAGCCCUAUUGAGACUGUACCAGUAAAAUUAAAGCCAGGAAUGGAUGGCCCAAAAGUUAAACAAUGGCCAUUGACAGAAGAAAAAAUAAAAGCAUUAGUAGAAAUUUGUACAGAGAUGGAAAAGGAAGGGAAAAUUUCAAAAAUUGGGCCUGAAAAUCCAUACAAUACUCCAGUAUUUGCCAUAAAGAAAAAAGACAGUACUAAAUGGAGAAAAUUAGUAGAUUUCAGAGAACUUAAUAAGAGAACUCAAGACUUCUGGGAAGUUCAAUUAGGAAUACCACAUCCCGCAGGGUUAAAAAAGAAAAAAUCAGUAACAGUACUGGAUGUGGGUGAUGCAUAUUUUUCAGUUCCCUUAGAUGAAGACUUCAGGAAGUAUACUGCAUUUACCAUACCUAGUAUAAACAAUGAGACACCAGGGAUUAGAUAUCAGUACAAUGUGCUUCCACAGGGAUGGAAAGGAUCACCAGCAAUAUUCCAAAGUAGCAUGACAAAAAUCUUAGAGCCUUUUAGAAAACAAAAUCCAGACAUAGUUAUCUAUCAAUACAUGGAUGAUUUGUAUGUAGGAUCUGACUUAGAAAUAGGGCAGCAUAGAACAAAAAUAGAGGAGCUGAGACAACAUCUGUUGAGGUGGGGACUUACCACACCAGACAAAAAACAUCAGAAAGAACCUCCAUUCCUUUGGAUGGGUUAUGAACUCCAUCCUGAUAAAUGGACAGUACAGCCUAUAGUGCUGCCAGAAAAAGACAGCUGGACUGUCAAUGACAUACAGAAGUUAGUGGGGAAAUUGAAUUGGGCAAGUCAGAUUUACCCAGGGAUUAAAGUAAGGCAAUUAUGUAAACUCCUUAGAGGAACCAAAGCACUAACAGAAGUAAUACCACUAACAGAAGAAGCAGAGCUAGAACUGGCAGAAAACAGAGAGAUUCUAAAAGAACCAGUACAUGGAGUGUAUUAUGACCCAUCAAAAGACUUAAUAGCAGAAAUACAGAAGCAGGGGCAAGGCCAAUGGACAUAUCAAAUUUAUCAAGAGCCAUUUAAAAAUCUGAAAACAGGAAAAUAUGCAAGAAUGAGGGGUGCCCACACUAAUGAUGUAAAACAAUUAACAGAGGCAGUGCAAAAAAUAACCACAGAAAGCAUAGUAAUAUGGGGAAAGACUCCUAAAUUUAAACUGCCCAUACAAAAGGAAACAUGGGAAACAUGGUGGACAGAGUAUUGGCAAGCCACCUGGAUUCCUGAGUGGGAGUUUGUUAAUACCCCUCCCUUAGUGAAAUUAUGGUACCAGUUAGAGAAAGAACCCAUAGUAGGAGCAGAAACCUUCUAUGUAGAUGGGGCAGCUAACAGGGAGACUAAAUUAGGAAAAGCAGGAUAUGUUACUAAUAGAGGAAGACAAAAAGUUGUCACCCUAACUGACACAACAAAUCAGAAGACUGAGUUACAAGCAAUUUAUCUAGCUUUGCAGGAUUCGGGAUUAGAAGUAAACAUAGUAACAGACUCACAAUAUGCAUUAGGAAUCAUUCAAGCACAACCAGAUCAAAGUGAAUCAGAGUUAGUCAAUCAAAUAAUAGAGCAGUUAAUAAAAAAGGAAAAGGUCUAUCUGGCAUGGGUACCAGCACACAAAGGAAUUGGAGGAAAUGAACAAGUAGAUAAAUUAGUCAGUGCUGGAAUCAGGAAAGUACUAUUUUUAGAUGGAAUAGAUAAGGCCCAAGAUGAACAUGAGAAAUAUCACAGUAAUUGGAGAGCAAUGGCUAGUGAUUUUAACCUGCCACCUGUAGUAGCAAAAGAAAUAGUAGCCAGCUGUGAUAAAUGUCAGCUAAAAGGAGAAGCCAUGCAUGGACAAGUAGACUGUAGUCCAGGAAUAUGGCAACUAGAUUGUACACAUUUAGAAGGAAAAGUUAUCCUGGUAGCAGUUCAUGUAGCCAGUGGAUAUAUAGAAGCAGAAGUUAUUCCAGCAGAAACAGGGCAGGAAACAGCAUAUUUUCUUUUAAAAUUAGCAGGAAGAUGGCCAGUAAAAACAAUACAUACUGACAAUGGCAGCAAUUUCACCGGUGCUACGGUUAGGGCCGCCUGUUGGUGGGCGGGAAUCAAGCAGGAAUUUGGAAUUCCCUACAAUCCCCAAAGUCAAGGAGUAGUAGAAUCUAUGAAUAAAGAAUUAAAGAAAAUUAUAGGACAGGUAAGAGAUCAGGCUGAACAUCUUAAGACAGCAGUACAAAUGGCAGUAUUCAUCCACAAUUUUAAAAGAAAAGGGGGGAUUGGGGGGUACAGUGCAGGGGAAAGAAUAGUAGACAUAAUAGCAACAGACAUACAAACUAAAGAAUUACAAAAACAAAUUACAAAAAUUCAAAAUUUUCGGGUUUAUUACAGGGACAGCAGAAAUCCACUUUGGAAAGGACCAGCAAAGCUCCUCUGGAAAGGUGAAGGGGCAGUAGUAAUACAAGAUAAUAGUGACAUAAAAGUAGUGCCAAGAAGAAAAGCAAAGAUCAUUAGGGAUUAUGGAAAACAGAUGGCAGGUGAUGAUUGUGUGGCAAGUAGACAGGAUGAGGAUUAGAACAUGGAAAAGUUUAGUAAAACACCAUAUGUAUGUUUCAGGGAAAGCUAGGGGAUGGUUUUAUAGACAUCACUAUGAAAGCCCUCAUCCAAGAAUAAGUUCAGAAGUACACAUCCCACUAGGGGAUGCUAGAUUGGUAAUAACAACAUAUUGGGGUCUGCAUACAGGAGAAAGAGACUGGCAUUUGGGUCAGGGAGUCUCCAUAGAAUGGAGGAAAAAGAGAUAUAGCACACAAGUAGACCCUGAACUAGCAGACCAACUAAUUCAUCUGUAUUACUUUGACUGUUUUUCAGACUCUGCUAUAAGAAAGGCCUUAUUAGGACACAUAGUUAGCCCUAGGUGUGAAUAUCAAGCAGGACAUAACAAGGUAGGAUCUCUACAAUACUUGGCACUAGCAGCAUUAAUAACACCAAAAAAGAUAAAGCCACCUUUGCCUAGUGUUACGAAACUGACAGAGGAUAGAUGGAACAAGCCCCAGAAGACCAAGGGCCACAGAGGGAGCCACACAAUGAAUGGACACUAGAGCUUUUAGAGGAGCUUAAGAAUGAAGCUGUUAGACAUUUUCCUAGGAUUUGGCUCCAUGGCUUAGGGCAACAUAUCUAUGAAACUUAUGGGGAUACUUGGGCAGGAGUGGAAGCCAUAAUAAGAAUUCUGCAACAACUGCUGUUUAUCCAUUUUCAGAAUUGGGUGUCGACAUAGCAGAAUAGGCGUUACUCGACAGAGGAGAGCAAGAAAUGGAGCCAGUAGAUCCUAGACUAGAGCCCUGGAAGCAUCCAGGAAGUCAGCCUAAAACUGCUUGUACCAAUUGCUAUUGUAAAAAGUGUUGCUUUCAUUGCCAAGUUUGUUUCAUAACAAAAGCCUUAGGCAUCUCCUAUGGCAGGAAGAAGCGGAGACAGCGACGAAGAGCUCAUCAGAACAGUCAGACUCAUCAAGCUUCUCUAUCAAAGCAGUAAGUAGUACAUGUAAUGCAACCUAUACCAAUAGUAGCAAUAGUAGCAUUAGUAGUAGCAAUAAUAAUAGCAAUAGUUGUGUGGUCCAUAGUAAUCAUAGAAUAUAGGAAAAUAUUAAGACAAAGAAAAAUAGACAGGUUAAUUGAUAGACUAAUAGAAAGAGCAGAAGACAGUGGCAAUGAGAGUGAAGGAGAAAUAUCAGCACUUGUGGAGAUGGGGGUGGAGAUGGGGCACCAUGCUCCUUGGGAUGUUGAUGAUCUGUAGUGCUACAGAAAAAUUGUGGGUCACAGUCUAUUAUGGGGUACCUGUGUGGAAGGAAGCAACCACCACUCUAUUUUGUGCAUCAGAUGCUAAAGCAUAUGAUACAGAGGUACAUAAUGUUUGGGCCACACAUGCCUGUGUACCCACAGACCCCAACCCACAAGAAGUAGUAUUGGUAAAUGUGACAGAAAAUUUUAACAUGUGGAAAAAUGACAUGGUAGAACAGAUGCAUGAGGAUAUAAUCAGUUUAUGGGAUCAAAGCCUAAAGCCAUGUGUAAAAUUAACCCCACUCUGUGUUAGUUUAAAGUGCACUGAUUUGAAGAAUGAUACUAAUACCAAUAGUAGUAGCGGGAGAAUGAUAAUGGAGAAAGGAGAGAUAAAAAACUGCUCUUUCAAUAUCAGCACAAGCAUAAGAGGUAAGGUGCAGAAAGAAUAUGCAUUUUUUUAUAAACUUGAUAUAAUACCAAUAGAUAAUGAUACUACCAGCUAUAAGUUGACAAGUUGUAACACCUCAGUCAUUACACAGGCCUGUCCAAAGGUAUCCUUUGAGCCAAUUCCCAUACAUUAUUGUGCCCCGGCUGGUUUUGCGAUUCUAAAAUGUAAUAAUAAGACGUUCAAUGGAACAGGACCAUGUACAAAUGUCAGCACAGUACAAUGUACACAUGGAAUUAGGCCAGUAGUAUCAACUCAACUGCUGUUAAAUGGCAGUCUAGCAGAAGAAGAGGUAGUAAUUAGAUCUGUCAAUUUCACGGACAAUGCUAAAACCAUAAUAGUACAGCUGAACACAUCUGUAGAAAUUAAUUGUACAAGACCCAACAACAAUACAAGAAAAAGAAUCCGUAUCCAGAGAGGACCAGGGAGAGCAUUUGUUACAAUAGGAAAAAUAGGAAAUAUGAGACAAGCACAUUGUAACAUUAGUAGAGCAAAAUGGAAUAACACUUUAAAACAGAUAGCUAGCAAAUUAAGAGAACAAUUUGGAAAUAAUAAAACAAUAAUCUUUAAGCAAUCCUCAGGAGGGGACCCAGAAAUUGUAACGCACAGUUUUAAUUGUGGAGGGGAAUUUUUCUACUGUAAUUCAACACAACUGUUUAAUAGUACUUGGUUUAAUAGUACUUGGAGUACUGAAGGGUCAAAUAACACUGAAGGAAGUGACACAAUCACCCUCCCAUGCAGAAUAAAACAAAUUAUAAACAUGUGGCAGAAAGUAGGAAAAGCAAUGUAUGCCCCUCCCAUCAGUGGACAAAUUAGAUGUUCAUCAAAUAUUACAGGGCUGCUAUUAACAAGAGAUGGUGGUAAUAGCAACAAUGAGUCCGAGAUCUUCAGACCUGGAGGAGGAGAUAUGAGGGACAAUUGGAGAAGUGAAUUAUAUAAAUAUAAAGUAGUAAAAAUUGAACCAUUAGGAGUAGCACCCACCAAGGCAAAGAGAAGAGUGGUGCAGAGAGAAAAAAGAGCAGUGGGAAUAGGAGCUUUGUUCCUUGGGUUCUUGGGAGCAGCAGGAAGCACUAUGGGCGCAGCCUCAAUGACGCUGACGGUACAGGCCAGACAAUUAUUGUCUGGUAUAGUGCAGCAGCAGAACAAUUUGCUGAGGGCUAUUGAGGCGCAACAGCAUCUGUUGCAACUCACAGUCUGGGGCAUCAAGCAGCUCCAGGCAAGAAUCCUGGCUGUGGAAAGAUACCUAAAGGAUCAACAGCUCCUGGGGAUUUGGGGUUGCUCUGGAAAACUCAUUUGCACCACUGCUGUGCCUUGGAAUGCUAGUUGGAGUAAUAAAUCUCUGGAACAGAUUUGGAAUCACACGACCUGGAUGGAGUGGGACAGAGAAAUUAACAAUUACACAAGCUUAAUACACUCCUUAAUUGAAGAAUCGCAAAACCAGCAAGAAAAGAAUGAACAAGAAUUAUUGGAAUUAGAUAAAUGGGCAAGUUUGUGGAAUUGGUUUAACAUAACAAAUUGGCUGUGGUAUAUAAAAUUAUUCAUAAUGAUAGUAGGAGGCUUGGUAGGUUUAAGAAUAGUUUUUGCUGUACUUUCUAUAGUGAAUAGAGUUAGGCAGGGAUAUUCACCAUUAUCGUUUCAGACCCACCUCCCAACCCCGAGGGGACCCGACAGGCCCGAAGGAAUAGAAGAAGAAGGUGGAGAGAGAGACAGAGACAGAUCCAUUCGAUUAGUGAACGGAUCCUUGGCACUUAUCUGGGACGAUCUGCGGAGCCUGUGCCUCUUCAGCUACCACCGCUUGAGAGACUUACUCUUGAUUGUAACGAGGAUUGUGGAACUUCUGGGACGCAGGGGGUGGGAAGCCCUCAAAUAUUGGUGGAAUCUCCUACAGUAUUGGAGUCAGGAACUAAAGAAUAGUGCUGUUAGCUUGCUCAAUGCCACAGCCAUAGCAGUAGCUGAGGGGACAGAUAGGGUUAUAGAAGUAGUACAAGGAGCUUGUAGAGCUAUUCGCCACAUACCUAGAAGAAUAAGACAGGGCUUGGAAAGGAUUUUGCUAUAAGAUGGGUGGCAAGUGGUCAAAAAGUAGUGUGAUUGGAUGGCCUACUGUAAGGGAAAGAAUGAGACGAGCUGAGCCAGCAGCAGAUAGGGUGGGAGCAGCAUCUCGAGACCUGGAAAAACAUGGAGCAAUCACAAGUAGCAAUACAGCAGCUACCAAUGCUGCUUGUGCCUGGCUAGAAGCACAAGAGGAGGAGGAGGUGGGUUUUCCAGUCACACCUCAGGUACCUUUAAGACCAAUGACUUACAAGGCAGCUGUAGAUCUUAGCCACUUUUUAAAAGAAAAGGGGGGACUGGAAGGGCUAAUUCACUCCCAAAGAAGACAAGAUAUCCUUGAUCUGUGGAUCUACCACACACAAGGCUACUUCCCUGAUUAGCAGAACUACACACCAGGGCCAGGGGUCAGAUAUCCACUGACCUUUGGAUGGUGCUACAAGCUAGUACCAGUUGAGCCAGAUAAGAUAGAAGAGGCCAAUAAAGGAGAGAACACCAGCUUGUUACACCCUGUGAGCCUGCAUGGGAUGGAUGACCCGGAGAGAGAAGUGUUAGAGUGGAGGUUUGACAGCCGCCUAGCAUUUCAUCACGUGGCCCGAGAGCUGCAUCCGGAGUACUUCAAGAACUGCUGACAUCGAGCUUGCUACAAGGGACUUUCCGCUGGGGACUUUCCAGGGAGGCGUGGCCUGGGCGGGACUGGGGAGUGGCGAGCCCUCAGAUCCUGCAUAUAAGCAGCUGCUUUUUGCCUGUACUGGGUCUCUCUGGUUAGACCAGAUCUGAGCCUGGGAGCUCUCUGGCUAACUAGGGAACCCACUGCUUAAGCCUCAAUAAAGCUUGCCUUGAGUGCUUC

[GenoMax]

Looks like this (T's changed to U's): http://www.ncbi.nlm.nih.gov/nuccore/9629357/

[Brian Bushnell]

OK, that looks like a bug that is related to the "subfilter" setting. Without subfilter, it gets mapped here:

M01472:214:000000000-AG0YC:1:2108:10755:20410 1:N:0:78  0       NC_001802DRannotations_(modified)   3186        15      1=1X2=1X7=3X2=2X1=1X3=3X9=3X2=2X1=2X34=1X1=1X24=1X21=

With subfilter, that site gets killed because it has too many substitutions, so it chooses a different site that presumably has a worse alignment (more indels) but fewer than 15 substitutions, so it's still valid. However, for some reason the cigar string does not get generated for the new site. I'll find out why and fix that very soon. There's nothing wrong with your command.

-Brian

[Brian Bushnell]

This bug has now been fixed as of 35.14, which I just uploaded. Filtering was essentially not being applied when the best site was filtered and inferior sites passed the minid yet did not have cigar strings (a very rare scenario). But, it works fine now.

-Brian

P.S. Without filtering, just with default parameters, you get this alignment:

M01472:214:000000000-AG0YC:1:2108:10755:20410 1:N:0:78  0       NC_001802DRannotations_(modified)   2154        35      46=1032D34=1X1=1X24=1X21=

...which in my view is better than any alternative, despite the long deletion.

[WhatsOEver]

Hi Brian,
I found a small bug in the v1.4 cigar strings: Whenever there is a "N" in the read sequence, one gets a "M" for this position instead of a "=".

Read:
GATATCGGTTTCATCCTCGCCTTAGCATGATTTATCCTACACTCCAACTCATGAGACCCCANAACAAATAGCCCTTCTAAACGCTAATCCAAGCCTCACCCCCACTACTAGGCCTCCTCCTAGCAGCAGCAGGCAAATCAGCCCAATTGGCTCCACCCCTGACTCCCCTCAG
aligns with the following cigar:
56=1I4=1M36=1I49=1X1=2D22=

[Brian Bushnell]

That's actually intentional, and happens when there's an N in the read or the reference. The SAM spec does not specify what to do in these situations. So, the way I see it, A aligning to A is a match (=), a aligning to C is a substitution (X), and A aligning to N - or even N aligning to N - might be a match, and might be a substitution. Since it's not an insertion or deletion, it's length-neutral, so I call it an alignment match.

Are these causing problems for some downstream application? If so, I can add a flag that changes the behavior.

[sdriscoll]

Hello Brian-

I'm curious if there is an output setting that I'm missing for PE alignments that would maintain a strict grouping of the mates in the output. I'd like to try using bbmap.sh as an aligner to feed RSEM (and maybe eXpress) which was designed around bowtie alignments. RSEM is VERY picky and there is almost no way to get alignments to work with their program without running them through some sort of reformatting script. One thing, however, that would be useful is to have the alignment output be able to keep mates in two-row pairs no matter what. So no matter what if line X is the left mate (mapped or unmapped) then line X+1 is the right mate mapped or unmapped. This way every pair of lines go together even through this will result in one mate or the other's alignment being reported more than once. RSEM flips out if it doesn't find the alignments like this.

TL;DR; An option so that every pair of lines in the output SAM represent the first and second mates of a fragment even if this causes redundant reporting of individual mate alignments.

Another question - I wanted to exclude singletons and improper pairs (most specifically pairs with mates mapped to different transcripts) so I enabled pairedonly AND killbadpairs. Each option taken separately works however killbadpairs seems to happen AFTER pairedonly, internally, because killbadpairs creates singletons in the output. Is it possible to have pairedonly apply to the result of killbadpairs as well?

thanks!

[sdriscoll]

Sort of continuing on the subject I posted above. Again I'm looking at paired-end alignments to a transcriptome where I'm allowing multiple alignments per pair. I'm using both ambiguous=all and secondary=t to allow some additional sub-optimal alignments to be reported. I need pairs to be reported as neighbors in the output so when there are 2 alignments for the first mate but only 1 for the second I need that to either be two lines or four lines depending on how those alignments were paired. So in parsing the output of bbmap today I found something that's throwing me off that appears to be a bug.

I've been parsing out reads in cases where there are more than one alignment of either the left or right mate so that I can see how bbmap reports those alignments. The following is one such set (first mate separated from second mate alignments):

#== left ==#
SRR1016945.2666	83	NM_146193	1295	42	100M	=	1199	-196	GAACCTATAGAGATCCTGCCCAATGTCTGTTACACAGCGTGCGCAACACTCAAGGGCCCAGATTCTCACTATGGCACAAAAGGACTGAAGAAAGTAGTGT	DCCDEEDEECCCDDDBDDDDDDEEDEEFFFEFFFHHIIIIIIIJIGGFJIJJJJJIJJJJIJJJJJJIJJJIJJJHJJJJJJJIJJJHHHHHFFFFFCCC	NM:i:1	AM:i:42	NH:i:1
SRR1016945.2666	353	ENSMUST00000121441	1072	40	100M	NM_146193	1295	0	CAGATCATTGAATATGAGAAAAAACAAACCTTGGGACAGAATGATACTGGATCTAGTTGTGATGGGACTGCTAACACATTCAGGGTCATGTTCAAAGAAC	CCCFFFFFHHHGHJJJJJIJJJJJJJJJJJJJJJJJJIIGIJJJJJIJJJJIIGIJIJJHIJJJJJIHHHHHFFFFFDEEDEEDDCCDDDEEEDCCDDDD	NM:i:2	AM:i:40	NH:i:3
SRR1016945.2666	353	ENSMUST00000178360	1054	40	100M	NM_146193	1295	0	CAGATCATTGAATATGAGAAAAAACAAACCTTGGGACAGAATGATACTGGATCTAGTTGTGATGGGACTGCTAACACATTCAGGGTCATGTTCAAAGAAC	CCCFFFFFHHHGHJJJJJIJJJJJJJJJJJJJJJJJJIIGIJJJJJIJJJJIIGIJIJJHIJJJJJIHHHHHFFFFFDEEDEEDDCCDDDEEEDCCDDDD	NM:i:2	AM:i:40	NH:i:3

#== right ==#
SRR1016945.2666	163	NM_146193	1199	42	100M	=	1295	196	CAGATCATTGAATATGAGAAAAAACAAACCTTGGGACAGAATGATACTGGATCTAGTTGTGATGGGACTGCTAACACATTCAGGGTCATGTTCAAAGAAC	CCCFFFFFHHHGHJJJJJIJJJJJJJJJJJJJJJJJJIIGIJJJJJIJJJJIIGIJIJJHIJJJJJIHHHHHFFFFFDEEDEEDDCCDDDEEEDCCDDDD	NM:i:1	AM:i:42	NH:i:3

The first three alignments all have 0x40 set indicating they are the first mate. The last one has 0x80 set indicating it is the second mate. The first+second pair with flags 83 and 163 are properly paired and confirmed with bowtie2. So what are the other two first mate alignments? Those aren't actually first mates...and the so-called RNEXT and PNEXT fields indicate they are mated to the first first-mate alignment. Also the sequences and qualities match up with those of the second-mate read. So those should actually have the 0x80 flag set since they are NOT first-mate reads at all.

This appears to be happening randomly because sometimes I'll see multiple first-mate mappings and a single second-mate mapping but all of the 0x40 and 0x80 flags are set correctly.

Another interesting case:

#== left ==#
SRR1016945.1906	83	NM_011653	367	44	100M	=	266	-201	AGCAGCTCATCACAGGCAAGGAGGATGCTGCCAATAACTATGCTCGTGGCCACTACACCATTGGCAAGGAGATCATTGACCTTGTCCTGGACAGGATTCG	DDDDDCEDDDC@DDDDDCDEEDEEEFFFFHGHHGHIIJIIIJJJJJIIJIJJIHHGDJIGGGJJJJJIJJJJIIGHFIGJJJJJJJJHHHHHFFFFFCCC	NM:i:0	AM:i:44	NH:i:3
SRR1016945.1906	337	NM_009448	800	42	100M	NM_011653	266	0	AGCAGCTCATCACAGGCAAGGAGGATGCTGCCAATAACTATGCTCGTGGCCACTACACCATTGGCAAGGAGATCATTGACCTTGTCCTGGACAGGATTCG	DDDDDCEDDDC@DDDDDCDEEDEEEFFFFHGHHGHIIJIIIJJJJJIIJIJJIHHGDJIGGGJJJJJIJJJJIIGHFIGJJJJJJJJHHHHHFFFFFCCC	NM:i:1	AM:i:42	NH:i:3
SRR1016945.1906	337	ENSMUST00000122404	269	42	100M	NM_011653	266	0	AGCAGCTCATCACAGGCAAGGAGGATGCTGCCAATAACTATGCTCGTGGCCACTACACCATTGGCAAGGAGATCATTGACCTTGTCCTGGACAGGATTCG	DDDDDCEDDDC@DDDDDCDEEDEEEFFFFHGHHGHIIJIIIJJJJJIIJIJJIHHGDJIGGGJJJJJIJJJJIIGHFIGJJJJJJJJHHHHHFFFFFCCC	NM:i:1	AM:i:42	NH:i:3
SRR1016945.1906	353	NM_009448	699	42	100M	NM_011653	367	0	AGGAGCTGGCAAGCATGTGCCCCGGGCAGTGTTCGTAGACCTGGAACCCACGGTCATCGATGAAGTTCGCACCGGCACCTACCGCCAGCTCTTCCACCCT	B@CDDDFDHHHHHIJIJI>AHGHIGIBGIGHDHGIBFHGJJJJJIJFI6D;BHDFFFCEEDEDDDACDDBBB=@BBBDBDDCB@BBBDBDCCDDDDD<<A	NM:i:1	AM:i:42	NH:i:2
#== right ==#
SRR1016945.1906	163	NM_011653	266	44	100M	=	367	201	AGGAGCTGGCAAGCATGTGCCCCGGGCAGTGTTCGTAGACCTGGAACCCACGGTCATCGATGAAGTTCGCACCGGCACCTACCGCCAGCTCTTCCACCCT	B@CDDDFDHHHHHIJIJI>AHGHIGIBGIGHDHGIBFHGJJJJJIJFI6D;BHDFFFCEEDEDDDACDDBBB=@BBBDBDDCB@BBBDBDCCDDDDD<<A	NM:i:0	AM:i:44	NH:i:2

In the set of 0x40 flagged alignments only the fourth one should have the 0x80 flag. I didn't catch this until I looked at alignments for this pair from bowtie2 but the second first-mate and the fourth first-mate (the two mapped to NM_009448) should actually be mated to one another (i.e. bowtie reports that exact pair as a properly paired alignment with insert size 201 which is the same as the primary here to NM_011653). The alignment at NM_011653 is perfect (no mismatches) however the alignment at NM_009448 has 1 mismatch per mate (from bowtie2). So that should fall under the category of secondary alignments in bbmap. Also, curiously, bowtie2 has reported 6 different pairings for this fragment with between 0 and 2 mismatches. I allowed secondary mappings to be within 0.95 of the primary score so maybe I need to lower that in order to see alignments with > 1 mismatch when the best has 0 mismatches however there are two additional sequences that bowtie2 aligned this fragment to with only a single mismatch and they aren't reported in bbmap's output so it seems it may be missing some valid mappings.

[Brian Bushnell]

Woah, I replied to your first post and it never registered. To summarize:

I was unable to reproduce an interaction of the "pairedonly" and "killbadpairs" flag generating singletons, so if you can send me a read pair and reference for which that occurs, it would be helpful. In my testing, when pairedonly=true, I get identical results with killbadpairs true or false. But, I don't use killbadpairs anymore - when I want to produce only properly paired mapped reads, I set pairedonly and use the outm stream, e.g.

bbmap.sh in=reads.fq outm=mapped_paired.fq pairedonly

For RSEM - the default behaviour of BBMap is to produce R1, R2, R1, R2, etc. in sam files, so it should be fine. It will do this with the flags "ambig=best" (default), "ambig=random", and "ambig=toss". The only time it will not do that is with "ambig=all", in which it will print all R1 alignments, then all R2 alignments. Enforcing strict interleaving when printing all alignments would probably result in output that is either incorrect or violates the sam spec, which makes me wonder where it ever occurs... so, I'm a bit reluctant to add the option, but I already have other flags like "keepnames" that (as noted) violate the sam spec but are still useful, so I'm willing to add it if useful. Does RSEM require all mapping locations, rather than a random location, for multimapping reads?

As for your second post - that seems much more concerning (to me). I'll examine it in detail tomorrow, and thanks as always for your thorough bug reports. I'm still working on the Seal qhdist issue, by the way; I've just been preoccupied with upgrading some programs to support unlimited max kmer lengths, which is now mostly finished.

[sdriscoll]

I'll have to doublecheck that pairedonly/killbadpairs thing I thought I was seeing. I like the pairedonly strategy and I can deal with mates mapped to different references in other ways.

Yes, tools like RSEM and eXpress need to have as many alignments as possible per read within reason. In fact the default alignment for RSEM is with bowtie (1) allowing 2 mismatches in the seed and then up to the entire rest of the read to be mismatches and all alignments that fall under that threshold. Multiple best alignments is also super useful for the program because they it can see which reads are totes ambiguous (i.e. the ones that map exactly the same to mulitple isoforms - those that hit an exon that's shared among multiple isoforms). So that information is necessary but it's also necessary to have what most aligners consider suboptimal alignments because those programs don't necessarily look at each alignment's quality but the overall solution to the abundance problem. The need to see which target references a read can map to within some reason. I don't like to allow as much wiggle room as the default RSEM bowtie setting - I typically allow up to 10% error. In some cases it would be entirely reasonable to only accept the best alignment...but I would still need ambig=all to satisfy the requirements for RSEM/eXpress. Obviously I could just use bowtie2 which works great but it is slooowwwwwww. STAR works pretty good and is a great substitute for bowtie2. I want to see if I can find the right settings for BBMap too and see how it ranks. The variable between programs tends to be the number of alignments found. bowtie2 finds many more PE alignments than bowtie1, for example. STAR comes close to bowtie2 while bwa falls short of them both. Basically whichever aligner combines finding the 'most' mappings with speed is my favorite for this job. I haven't specifically looked to see how well BBMap stacks up but I'll do that. The kind of test I'm talking about is to generate transcriptome based simulated reads and then map them with different aligners. I just check, per read/pe-fragment, if the aligner successfully reported the original correct alignment along with, if any, additional alignments. Again bowtie2 was the winner between bowtie1, STAR, bwa (aln and mem), and I think subread. Unfortunately it is slow in the alignment mode where it reports more than a single alignment (-k or -a mode).

FYI the default behavior for STAR, bowtie and bowtie2 when allowing multiple mappings per pair is to report them in R1, R2, R1, R2 strict formatting. So if bowtie2 finds 12 locations for a single fragment it will print out 24 lines in R1,R2 pairs. All of the secondary locations get the 0x100 flag set. Basically that just makes parsing and interpreting them really easy since, as you know, if you can count on a certain formatting it means you can write less code. Obviously the RSEM program could do this but they chose to put the burden of formatting on the user - probably to keep the number of decisions going on behind the scenes to a minimum. I was parsing the BBMap alignments by watching the read name and piling up the 0x40 mates in one array with 0x80 mates in a second array and then afterwards I ran all possible combinations of them through a function that validates whether or not they are a valid pair. It would make their program a little more powerful, in my opinion, if they could just incorporate something like that but alas...I guess that's what keeps Perl in business.

[sdriscoll]

As for your second post - that seems much more concerning (to me). I'll examine it in detail tomorrow, and thanks as always for your thorough bug reports. I'm still working on the Seal qhdist issue, by the way; I've just been preoccupied with upgrading some programs to support unlimited max kmer lengths, which is now mostly finished.

By the way this was when I had ambig=all and aligning PE reads. The odd behavior didn't seem to care whether or not I had 'secondary=t' set.

[rkanwar]

Aligning reads after bbmerge

Hello,

I have PE reads, many of which overlap. I used bbmerge to merge the overlapping reads. In the end I have 3 fastq files:

1. Merged single end reads
2. End 1 of no-merged reads
3. End 2 of no-merged reads

How do I specify these 3 fastq files to bbmap [what should I give for in and in2] ? Thanks.

bye,
Rahul

[Brian Bushnell]

Hi Rahul,

For a paired end dataset and a single-end dataset, you need to run mapping twice. However, you can do it in one command with BBWrap like this:

bbwrap.sh ref=reference.fa in=read1.fq,merged.fq in2=read2.fq,null out=mapped.sam append=t

This would be equivalent to:

bbmap.sh ref=reference.fa
bbmap.sh in=read1.fq in2=read2.fq out=mapped_pairs.sam
bbmap.sh in=merged.fq out=mapped_merged.sam

...except that BBWrap would only produce a single sam output file.