Hi, I think that option was used initially in oxford Nanopore pipelines to remove internal priming in the cDNA-PCR protocols, before the addition (in the new protocols) of the ligation step. However it seems that when using this command, seqkit does not apply correctly the search at the 3' of the alignment. Theoretically it should look for the stretch of As, depending if the bam flag is 0 or 16 (forward or reverse), at the end or at the beginning of the alignment, after performing the reverse complement. However when I tested it, on the reverse reads (the ones with 16 flag) it was looking at the biological 5' end of the read instead of the 3' (the log also shows the sequence searched was inverted from start and end).

Looking into the code I hypothesized the bug might be here:

func GetSamReverse(r *sam.Record) bool {
    return (r.Flags&sam.Reverse == 0)  
}

seems to return FALSE when flag 16 is set (shouldn't it be the contrary?)
But then:

if GetSamReverse(r) {
    strand = -1
    if stranded {
        startSeq, endSeq = RevCompDNA(endSeq), RevCompDNA(startSeq)
    }
}

if GetSamReverse is TRUE (so the flag is not 16), the code assigns strand=-1.

Could you help me?

Best,
Camilla

seqkit bam #200

Description

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions