alignmentSieve seems to not be working

[gdolsten]

Hi everyone, since I was having a challenge shifting reads manually and using deeptools, I decided to shift reads with alignmentSieve. In this, I truncated a bam file to the first 10 reads to test it out, into a file "sort.test.bam" which is indexed. When I run the following command:

alignmentSieve -b sort.test.bam -o sieved.sort.test.bam --filterMetrics metrics.txt --samFlagInclude 16

I get the sort.test.bam file back. However, when I add a shift:

alignmentSieve -b sort.test.bam -o sieved.sort.test.bam --filterMetrics metrics.txt --samFlagInclude 16 --shift 10 10

I get the header back but nothing else. Do you know why this might be happenign?

[dpryan79]

I'm looking into this now.

[dpryan79]

This is a misunderstanding that may be cleared up by looking at the underlying code:

    start = b.pos
    end = start + b.query_alignment_end
    if b.is_reverse and not b.is_read2:
        end -= args.shift[2]
    elif b.is_reverse and b.is_read2:
        end += args.shift[1]
    elif not b.is_reverse and not b.is_read2:
        start += args.shift[0]
    else:
        start -= args.shift[3]

The various args.shift values are 10 10 10 10 in this case. You're selecting reads that align after reverse complementation, therefore read 1 will have its end decremented by 10 (through changing the CIGAR) and read 2 will have its end increased by 10 (again, by changing the CIGAR).

I do see a bug when it comes to handling soft-clipping though. Reads with 5' soft-clipping can be shifted more due to soft-clipping being ignored. I'll fix that.