Low complexity read filtering step in FindJunctions.cpp is broken for long reads

[mbeavitt]

I'm seeing a lot of valid reads get filtered out in this step in long read .bams, and I was wondering if I could have some help understanding the logic of this section of the code:

PsiCLASS/FindJunction.cpp

Lines 536 to 584 in 04e2a11

	int count[5] = { 0, 0, 0, 0, 0 } ;
	int pos = 0 ;
	for ( i = 0 ; i < ccnt ; ++i )
	{
	switch ( cigarSeg[i].type )
	{
	case 'S':
	pos += cigarSeg[i].len ;
	case 'M':
	case 'I':
	{
	for ( j = 0 ; j < cigarSeg[i].len ; ++j )
	++count[ (unsigned char) nucToNum[ col[9][pos + j] - 'A' ] ] ;
	pos += j ;
	} break ;
	case 'N':
	{
	int max = 0 ;
	int sum = 0 ;
	for ( j = 0 ; j < 5 ; ++j )
	{
	if ( count[j] > max )
	max = count[j] ;
	sum += count[j] ;
	}
	if ( max > 0.8 * sum ) {
	validRead = false ;
	}
	count[0] = count[1] = count[2] = count[3] = count[4] = 0 ;
	} break ;
	case 'H':
	case 'P':
	case 'D':
	default: break ;
	}
	}
	int max = 0 ;
	int sum = 0 ;
	for ( j = 0 ; j < 5 ; ++j )
	{
	if ( count[j] > max ) {
	max = count[j] ;
	}
	sum += count[j] ;
	}
	if ( max > 0.8 * sum )
	validRead = false ;

First off, there is an error - no break statement in case 'S' means there is a fallthrough for soft-clipped reads to case 'I', which offsets the pos variable. This later causes reads into incorrect data left over in the string from another alignment, or simply reads into garbage memory. Somehow nucToNum handles junk input silently, i.e. nucToNum[-65] == 0, so this does not cause a crash/error and simply churns out incorrect data. Below is an example, where the program just starts reading zeroes after a certain point. In combination with the reset of "count" when 'N' is reached in the cigar string, this read gets incorrectly filtered. Below is an example of the issue with the program loaded into a debugger.

The numbers are the numerical values of the characters being read from the sequence string, with a -'A' offset. (i.e. col[9][pos + j] - 'A' ). You can see that after a certain point, it starts reading zeroes (which appear as -65 due to -'A')

Apart from some soft clipped poly-T at the start of the read, the content seems decently complex, but because the program is reading junk from memory, the 126 matches register as the same nucleotide.

❯ samtools view lr_test.bam | grep 465S1446M534N59M1206N77M435N126M
SRR14298432.320776      16      Contig10        10022   0       465S1446M534N59M1206N77M435N126M        *   00       TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTAGCTAAAATAAGTATTACTTAGAAAACATCTGAAGTTAACATATTTTGGAAGTCGACGAAGATCAAGATGTACATTAATTTTTATAGAAATAACTGCTTTCTTCTCCACACATCTTATATCTCACAATACTGCTGCCGCACAACACCTAAGGCCCTTCAGACGCTCAAAAACTTGTTTCATTTCTGGTCGTTTCTTTTGTTCATCATCCATACACTGTCTTGCAAGTTCTGTAAGTTTACGTCCATCCCGACGATCAAAGTAGGTATCATCCTCAAGACTTTCGTCGACATATGAACGCCTGTGCTTGTACTCATTCGCGGCCCACTTATCAATTTCACACCCTGTATATTTTTCCTTGTCAAGAGCTCTCTTGCTTAGAAGATTGTAGAGCGAAATACCGAATGCATAAACATCAGACTTGACATTGTATGCAGCTGTGATCACATCGTCCGAGGTCTCCGCATCTGCTCTAGCTGGAGCTGCATAAAAGTGGGCCTGGACTTCCCTACCCTGTGCATCGGCTCTGCCCCCTAATCTACCTCCTCGGGGACCACCCCGTATACCCTGGTGACCACCCCTACGGTCCTGTCCCTGGCCTCTGCCTCTAGCTGGAGGTGTTGGCGCATGCCTAACTGCCTGTGGAGCGGGTAAGGCTAACCTGCGATGGGGACACTCACGGGCCCAATGATCTAACTCCCCACACTCATAACACCCCTGGGACTGTCGACTGAGAGCTGGAGGTCGGCTACCCCGGCCAGAAGAACCCATCTGTGAACTATCTCGGCCTCGACTAGGACCGGCACTAAAACCAATATAGCGGTGCTGATCACCCTCTGATGCCTGGAGAGAAGCCTGGAAGGGGCTACTAGATTGACCTGACTGAAACCTCTGCTGAGATCTATCATAAGACTCTCGGGGTCTGAAACGGCGCCTATCATGGCTATCCTCCTGTCTAGCTCTCTTGCCGCTGCCCCCTTGGGCCTCGCGACGGAGCTGCTCGATAGTGCGGGCAAGGTCCACTACATCCAGAAACGAGCGGCCTGCTGAGACTAAAGACTAAGTCTCAATCCGCAACTGTAGTCTCAAACCCCGAACAAAACAACGAACCCTCTCCTCCTCGGTAGGTAATATCATGGCGGCATGCCTGGAGAGCTCGTGGAACCGGGCCTCATATTCTGATACTGTCATGGAGCCCTGCTCTAACCAAGCAAACTGATCTCTGAGCTGATCTCTGATGCTCCTGGGAATAAATCGGGCCAAGAAGGCCTCTGAAAACTCAGTCCATGTGACCGGTGGAGACCCAGCTGGCCTGGTCTCTAAGTATGTACGCCACCACTGCCTGGCGGGACCGTCTAACTGGTAAGCGGTGAAGTCTGCCCCTCTCGACTCCACTAAACCCAAAGTCTGTAGCCGCTCGCGGCAAGTAACTAGAAACTCGTGGGCATCCTCCCCCACTGCCCCAGAAAACCTCGGCGGTGATAGTCTAAGGAACACCCCGAGCATCTTCTGCTCCCGTAATGGCATGACGCCCTCCAAATCGTCAGGCTGAACAACCGGTGCCGCTAGCTGCTGAACAACCGGTACCGCTGGTGTAGGCTGGCCCTGCGGTACTGGAACGACTGGAGCTGGGGCCTGCTGCATACCCCCAATGGCAGCTAGTATCTGTACGAGCATCGCCTGTATATCTGGAGTTGCCACAGGCTCGGGCGGTGGCTGGGCTGGTCCCGGGCACATCGGCTGCTGCTGTGCTGGAGCTGCTGGTGGCTCTATCACCGGCTCCAGAGTTTGCTCTCGGTCCTGGGCTGGTGCTGTAACTCTGGCACGACCTCGAGGTCGGCCTCTACCCCTAGCTGGGGACCTACCACGCGTGCGAGCCATCCCCCCGTCTCCGGATAAGTCTACTCCAAGCTAGATCAACGCTGGCTACTAGTGCTCGGACCTGAAAGCCGATCACGCGTGCTCCAACTCACGGTACCGGATCTTTTCCTCTCCGAACGGTCTCCAAATGCTTCCCCACTCAAATGAAGCCTCCCCTCGCGTTTCCGAGATCACCGATAGATTCCCCTCGAAATTCTCTTGAAGTTTGCCTTTTGAATCACTTAATTTGGACTCAAAATGAAGGATTAATCTCAAGTTGAAGTT       ,~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~w~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~u~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~w~~~~~~~~~~~~~~~~~~}~~~~~~~n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~L~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~N~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~       NM:i:39 ms:i:1513       AS:i:1441       nn:i:0  ts:A:+  tp:A:P  cm:i:374        s1:i:1488       s2:i:1487    de:f:0.0228     SA:Z:chr01,87709743,+,1705S439M2D29S,1,8;       rl:i:73

I'll write a fix for this, but even after this is corrected, some valid reads are still being filtered - this is due to the count resetting when an 'N' is reached in the CIGAR string, so a read that is complex but is terminated by a low complexity region will often get filtered out in the case that there is an N near the end of the alignment. This is an example of a read that was incorrectly filtered:

CIGAR:1S93M1D184M979N143M354N185M236N91M143N334M151N58M71N402M13319N17M16S
READ:GATTTGCTGTTGTTAGAAGCATTTAGTATTATAAACTATTTCTCTTTCATTCCTTTCTCACTTTTCTTTGCTTTTGGGTTTCTGGGGTTTGTTGTTTTTTTTTCTAGTTTTGTTTGAGTTTGATCTGCTTCTTGGGGTTCTTTAAATTCTCTTGAAATCTTGAATTTTGGGTGTTTCTTGAAAATTCATAGAGATGGAGAGGGACTTTATGGGATTGAATATTAAAGATTCTTTAGTTGTAGTCAAGGAAGAACATGTTGAAAGCTCCAAAGACTCTGGGTTTCGCTGGCCAATGTCGAGTAAGGTTGGTGTACCTCATTTCATGTCCTTGAACTCTGCTCAAGATGAAAACACCTTCAAAGCTCTATCTGCCACAGATGGAGTCGATGCTGGTCTCAAACGTCAGTCTGGTGAACTCCAGAATGTGCATGCAATGCATCTTCCGTATGATGTTAAGATGCTTCCATTTAACATGAACAATCCCTCCTACAAGACTCATTTTGGTGGTTCTGGCCAGATGAAGCAAGTTCTTGGTGGAATTCCTGTTACAGCTCCTCATTCAATGCTTCCATCGTGTGGCTCTGTGGCTGGGACAACCGAACCTTGGTTTAATTCCAAAAGTTCTGCAGCACCTGTTCAACTGACCATCTTCTATGGTGGGACGGTCAATGTCTTCGAGGATATCTCCCCTGAGAAGGCACAGGCTATTATGTTTTTGGCTGGACATGGCTGTGCTCCACCTAATGTGGTGCAACCAAGGTTTCAACTTCAGGCATCTGCATCAAAACCTGCTGCTGCGGATGGTGTTTGTGTGAACCAAACCCCAAACATGCTGCCTGCCTCGGGTCTCTCTAGCCCTATGTCTGTTTCTTCCCATCCCATTGGTCAATCTGATGGCAGTTCTGGAAACAAAGATGACATGAAGATGTCTAAAACTGCAAACATTTCAGTGACTCCCCTTGTCAAATTGGACACTTCAAAGATTGTGACATCACTAGGACCTGTTGGGGCGACTACCATAATGACAGCAGCGGTTCCACAAGCTAGGAAAGCATCTCTGGCTCGTTTTCTGGAGAAGCGTAAGGAAAGGGTGATGAACTTAGCACCATATGGCCTCAGCAAGAAATCGCCTGAGTGCTTCACCCCCGAGTCTAATGGAGTTGGTTTCTCTGCGACTTCCAGTCCTCTGTTAGCCGGTAAGGAGACCTAGCTATGATCTCTGAAGACTGCTTAACGCTGTAAATCAGTGAAGCATAATCTATTTTATTCAGCAAGTTAAAAATCTGTGTAGCTATCAAGTTTTGTTTGTTGTAATTTGTAATGTCAGCTTTGGTTTAGCTCTAGCTAGACATCCTAGTATCCAATTGGTTTGGTCTTTTTTTTTCTGTAACCTGATGGAATTTATCCCCATCCAAATAACTTCTTCAATTTGTAACTATTTATTTTGTTACATTGTGTTTCACAAGTATATATATATATTAATGTTATTTAGTTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

The CIGAR ends with 13319N17M16S. 13319N stops the count, resetting it, 17M reads into the start of the poly A: GTTAAAAAAAAAAAAAA and finally the rest of the poly-A is softclipped. As a result, the whole read is filtered based on this tail section.

What's the best way of dealing with this? I guess the easiest option is just to not reset when N is reached, but then this introduces the assumption that the whole read has a similar complexity. Alternatively, you could argue that the aligner didn't clip accurately, I guess, and then this is not a PsiCLASS issue.

I guess the ideal solution is to take some combined measure of the complexities of each chunk, somehow? I'd love some input, but in any case I'll be headscratching on this one for a while 😅.

Here's the test data I'm using (minimap2 bam):

https://drive.google.com/file/d/1nT7th6WKVPs1Jb5mhbm9CRAaHENP9MRd/view?usp=drive_link

[mbeavitt]

Ok, so the solution I've come up with is to take the complexity of each skip-separated chunk (the proportion taken up by one type of nucleotide, so a chunk AAAACGT would be ~= 0.57), where skip is 'N' in the CIGAR string, and add it to a weighted total complexity score. The weight is simply the size of the chunk.

We can then get a weighted average complexity by dividing the total complexity by the sum of the chunk sizes.

Formally:

• $n$: Total number of chunks.
• $s_i$: Size (number of bases) of the $i$-th chunk.
• $f_i$: Frequency of the most common base in the $i$-th chunk.
• $c_i$: Complexity of the $i$-th chunk, given by $c_i = \frac{f_i}{s_i}$.
• $T_c$: Total complexity score.
• $\bar{C}$: Weighted average complexity.

1. Chunk Complexity:

$c_i = \frac{f_i}{s_i}$

2. Total Complexity Score:

$T_c = \sum_{i=1}^{n} c_i \cdot s_i$

3. Weighted Average Complexity:

$\bar{C} = \frac{T_c}{\sum_{i=1}^{n} s_i}$

This works quite well in practice with slightly lower thresholds than the previous threshold of 0.8. A threshold of 0.5 filters out repetitive sequences without losing too many good alignments, with the biological meaning of 0.5 being along the lines of "on average, the chunks are at least 50% composed of one type of nucleotide across this alignment". Please see the log attached, which is the result of the new filter applied to this bam file: https://drive.google.com/file/d/1STrUGw6xASnWPaBM8huQ7qgp8jNkEZkN/view?usp=drive_link.

This is the log of the new filter:
filter_log.txt

And for comparison, this is the same log of alignments filtered in the old code:
old_code_filtering_log.txt

I've applied this to the code already on my fork, but I will wait for feedback to make a PR. I have actually not observed any impact on performance after testing, which is nice. This is the modified section:

https://github.com/mbeavitt/PsiCLASS/blob/99073d89e000832771e8f05d9e0f38acf06d2a0a/FindJunction.cpp#L543-L647

[mbeavitt]

Actually I've just realised that it's a bit confusing that the "complexity score" is higher when the read is less complex. Maybe I should think of a different term to use.

Also, I am aware that there are almost certainly other ways of determining DNA sequence complexity, but I'm not sure if there is a pre-existing algorithm for segmented RNA complexity like this. I would love to know if one exists already, and I'd happily change my approach.

[mbeavitt]

Sorry for the wait on this - I'm happy with this implementation now, please have a look over when you have time.

I renamed complexity score -> repetitiveness score for clarity