lh3 / bwa

Burrow-Wheeler Aligner for short-read alignment (see minimap2 for long-read alignment)
GNU General Public License v3.0
1.55k stars 556 forks source link

Increasing the 5' clipping penalty drops alignments with no 5' clipping #296

Open nh13 opened 4 years ago

nh13 commented 4 years ago

I am confused as to why after increasing the 5' clipping penalty, I do not see an alignment with no 5' clipping promoted to be the primary, and in fact, it disappears.

$ cat tmp.fasta
>read
AGCTCACATCCTGTATGCCCTCCCTCCATGCCAAGCAACTTTGGGGAGGATCTTCGTGGCTGTCTTGGTAGCAGCTTTTTAGGTGATCTCAGGGCTGTCTGGAGGCTTCTTTTGGTTTAGCTGCGACTACTGCGAGTACGTGGTATTGGAC

$ bwa mem -Y hs38DH.fa tmp.fasta | grep ^read
read    16  chr1    26317934    60  27S78M46S   *   0   0   GTCCAATACCACGTACTCGCAGTAGTCGCAGCTAAACCAAAAGAAGCCTCCAGACAGCCCTGAGATCACCTAAAAAGCTGCTACCAAGACAGCCACGAAGATCCTCCCCAAAGTTGCTTGGCATGGAGGGAGGGCATACAGGATGTGAGCT *   NM:i:0  MD:Z:78 AS:i:78 XS:i:0  SA:Z:chr1,26318103,-,103S48M,60,0;
read    2064    chr1    26318103    60  103S48M *   0   0   GTCCAATACCACGTACTCGCAGTAGTCGCAGCTAAACCAAAAGAAGCCTCCAGACAGCCCTGAGATCACCTAAAAAGCTGCTACCAAGACAGCCACGAAGATCCTCCCCAAAGTTGCTTGGCATGGAGGGAGGGCATACAGGATGTGAGCT *   NM:i:0  MD:Z:48 AS:i:48 XS:i:21 SA:Z:chr1,26317934,-,27S78M46S,60,0;

$ bwa mem -L1000,5 -Y hs38DH.fa tmp.fasta | grep ^read
read    16  chr1    26317934    51  27S85M33I3M3D3M *   0   0   GTCCAATACCACGTACTCGCAGTAGTCGCAGCTAAACCAAAAGAAGCCTCCAGACAGCCCTGAGATCACCTAAAAAGCTGCTACCAAGACAGCCACGAAGATCCTCCCCAAAGTTGCTTGGCATGGAGGGAGGGCATACAGGATGTGAGCT *   NM:i:38 MD:Z:78A2A6^AGC3    AS:i:78 XS:i:48

So the 5' clipping is turned up to 1000, so I would expect that we get the 103S48M alignment since it has no 5' soft-clipping (sequencing order), but instead we get a strange 27S85M33I3M3D3M alignment.
The latter says it has score 78, which makes no sense:

85M = 85
33I = -6 + (-1 * 33) = -39
3M = 3
3D = -6 + (-1 * 3) = -9
3M = 3
so
85 + (-39) + 3 + (-9) + 3 = 43

unless you look back on the original alignment 27S78M46S where the soft-clips are not penalized (78M score).