XA: and SA: tags liftover

[zoldrax]

LevioSAM2 lifts over paired alignments accurately, but unfortunately it does not lift over coordinates in XA: and SA: tags. Some structural variant tools rely on these tags. I’m unaware of any tool that can lift over XA: and SA: tags. Is it possible to lift over XA and SA tags?

Example:
SA:Z:chr7:66983788-66996925,515,+,99M52S,60,0;
to
SA:Z:chr7,66984302,+,99M52S,60,0;

[milkschen]

Thank you for the comment. I can see this being a helpful feature. As you said, I'm not aware of other tools supporting lifting over those tags.

I added this to my new feature queue. However, the development of leviosam2 has been very slow due to personal bandwidth limitations. Forking this repo and adding the feature there would probably be a better approach if this feature is urgent to you. I'm happy to collaborate on PRs if there's any.

[zoldrax]

Thank you for the answer. For my urgent tasks, I already use additional simple own python script. It would be great if this feature was added in a planned future release.

[milkschen]

Would you mind sharing your Python script? That would be super helpful 😀

[zoldrax]

It works only for my simple case, is not universal:

import pysam
from pyliftover import LiftOver

def lift_position(lo, chrom, pos, strand):
    try:
        lifted = lo.convert_coordinate(chrom, pos, strand)
        if lifted and len(lifted) > 0:
            return lifted[0][0], lifted[0][1], lifted[0][2]
        else:
            return None, None, None
    except Exception as e:
        return None, None, None

def lift_cigar(cigar, strand):
    # For simplicity, return original CIGAR unless strand-specific adjustments are needed
    return cigar

def lift_sa_tag(lo, sa_tag):
    try:
        sa_parts = sa_tag.split(";")
        lifted_sa = []
        for sa in sa_parts:
            if not sa:
                continue
            fields = sa.split(",")
            if len(fields) != 6:
                continue
            chrom, pos, strand, cigar, mapq, nm = fields[:6]
            pos = int(pos) - 1  # Convert to 0-based
            new_chrom, new_pos, new_strand = lift_position(lo, chrom, pos, strand)
            if new_chrom:
                lifted_sa.append(f"{new_chrom},{new_pos + 1},{new_strand},{cigar},{mapq},{nm}")
        return ";".join(lifted_sa) if lifted_sa else None
    except Exception as e:
        return None

def lift_xa_tag(lo, xa_tag):
    try:
        xa_parts = xa_tag.split(";")
        lifted_xa = []
        for xa in xa_parts:
            if not xa:
                continue
            fields = xa.split(",")
            if len(fields) != 4:
                continue
            chrom, pos, cigar, nm = fields
            pos_val = int(pos)
            strand = '+' if pos_val >= 0 else '-'
            pos = abs(pos_val) - 1  # Convert to 0-based
            new_chrom, new_pos, new_strand = lift_position(lo, chrom, pos, strand)
            if new_chrom:
                new_pos_signed = new_pos + 1 if new_strand == '+' else -(new_pos + 1)
                lifted_xa.append(f"{new_chrom},{new_pos_signed},{cigar},{nm}")
        return ";".join(lifted_xa) if lifted_xa else None
    except Exception as e:
        return None

...
In main function:

               # Lift over SA tag
                if read.has_tag('SA'):
                    sa_tag = read.get_tag('SA')
                    new_sa = lift_sa_tag(lo, sa_tag)
                    if new_sa:
                        read.set_tag('SA', new_sa)

                # Lift over XA tag
                if read.has_tag('XA'):
                    xa_tag = read.get_tag('XA')
                    new_xa = lift_xa_tag(lo, xa_tag)
                    if new_xa:
                        read.set_tag('XA', new_xa)