Skip to content

Add query to ref coord mapping#120

Open
chilampoon wants to merge 4 commits into
nanoporetech:masterfrom
chilampoon:master
Open

Add query to ref coord mapping#120
chilampoon wants to merge 4 commits into
nanoporetech:masterfrom
chilampoon:master

Conversation

@chilampoon

Copy link
Copy Markdown

targeting spliced read issue #118 (comment)

Comment thread src/remora/data_chunks.py Outdated
return knots


def get_ref_coords(ops, lens):

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The dumb version of this function is:

def get_ref_coords(ops, lens):
    ref_coords = []
    s = 1

    for op, l in zip(ops, lens):
        if op in [3]:
            s += l
        elif op in [1,4,5]:
            continue
        elif op in [0,2,7,8]:
            # mapped
            pos_list = list(range(s, s + l))
            ref_coords.extend(pos_list)
            s += l
        
    return ref_coords

Comment thread src/remora/data_chunks.py Outdated
Comment thread src/remora/data_chunks.py Outdated
for op, start, l in zip(ops, starts, lens)
if REFCOORD_OPS[op]
]
ref_coords = np.concatenate(ranges) + 1

@marcus1487 marcus1487 Oct 11, 2023

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I still feel like there is extra logic here. I think the fix here might be as simple as completely dropping the skip cigar ops. Do you get the same results if you switch the REF_OPS skip code value to False without any other code changes?

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was thinking about the same thing I think you're right - just changing the intron operation to False in REF_OPS can work as well... got the same results for genomic reads and didn't raise errors for transcriptomic reads

@chilampoon

Copy link
Copy Markdown
Author

I put the same reference coordinate function to io.py now because it's a way to get a ref coordinate to query mapping with the exact same length as read.ref_to_signal. The get_reference_positions() from pysam doesn't include deletion bases, whereas its get_reference_sequence() includes deletions...
I tested the gapless reads without introns with updated code and got the same outputs, seems like there're many other issues for intronic reads, I'll stick to gapless read analysis for a while

@marcus1487

Copy link
Copy Markdown
Collaborator

I'm not sure I understand the issue. The get_reference_positions function is never used in Remora. Why is support for this function required in Remora?

I'll stick to gapless read analysis for a while

Are you saying that I should close this PR? Or would you like for this PR to be considered still?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants